text:
function statsheader(s::AbstractSimulation)
    simulationheader = "iteration,time,u1,u2,u3,u1p2,u2p2,u3p2,kh,kv,kaniso,k,dissh,dissv,diss,nlh,nlv,pressh,pressv,tnlh,tnlv,tdiss"
    header = join(Iterators.filter(x->x !== "",
      (simulationheader,statsheader.(getfield.(Ref(s),sim_fields))...,"\n")),
      ",","")
    return header
end

function writeheader(s::AbstractSimulation)
    open("Stats.txt","w") do f
        write(f,statsheader(s))
    end

    if (SPEC_FIL_D2 != 0.0)
        open("Stats_fil.txt","w") do f
            write(f,statsheader(s))
        end
    end

    open("Scales.txt","w") do f
        write(f,"iteration,time,Lk,Lt,lamg,laml,lamt,Lhh,Lvv,Lhv,Lvh,lam12,lam21,L11,L12,L13,L21,L22,L23,L31,L32,L33")
        if hasdensity(s)
            write(f,",Lb,Lo,Lom,Tk,Tt,Thh,Tb,Tbm\n")
        else
            write(f,",Tk,Tt,Thh\n")
        end
    end

    open("Numbers.txt","w") do f
        write(f,"iteration,time,Reh,Reh1,Relamg,Relaml,Relamt,Relamg1,Relaml1,Relamt1")
        if hasdensity(s)
            write(f,",Reb,Rebm,Frh,Frh1\n")
        else
            write(f,"\n")
        end
    end
end

function writestats(s::AbstractSimulation)
    init = s.iteration[]
    time = s.time[]
    results, lens, anum, results_fil = stats(s)

    open("Stats.txt","a+") do file
        join(file,(init, time, results..., "\n"), ",","")
    end

    if (STAT_FIL)
        open("Stats_fil.txt","a+") do file
            join(file,(init, time, results_fil..., "\n"), ",","")
        end
    end

    open("Scales.txt","a+") do file
        join(file,(init, time, lens..., "\n"), ",","")
    end

    open("Numbers.txt","a+") do file
        join(file,(init, time, anum..., "\n"), ",","")
    end

end

function stats(s::AbstractSimulation)
    vstats, vstats_fil = velocity_stats(s)

    u1p2 = vstats[4]
    u2p2 = vstats[5]
    u3p2 = vstats[6]
    uh = sqrt((u1p2+u2p2)/2)
    urms = sqrt((u1p2 + u2p2 + u3p2)/3)
    k = vstats[10]

    L11 = integral_lenght_x(s.u.c.x,u1p2)
    L12 = integral_lenght_y(s.u.c.x,u1p2)
    L13 = integral_lenght_z(s.u.c.x,u1p2)

    #λ12 = taylor_lenght_y(s.yspec,s.hspec,s.u.c.x,u1p2)
    du1dyp2 = dy_squared_mean(s.reductionh,s.u.c.x)
    λ12 = sqrt(2*u1p2/du1dyp2)
    #λ11 = taylor_lenght_x(s.xspec,s.u.c.x,u1p2)
    du1dxp2 = dx_squared_mean(s.reductionh,s.u.c.x)
    λ11 = sqrt(2*u1p2/du1dxp2)

    L21 = integral_lenght_x(s.u.c.y,u2p2)
    L22 = integral_lenght_y(s.u.c.y,u2p2)
    L23 = integral_lenght_z(s.u.c.y,u2p2)

    #λ21 = taylor_lenght_x(s.xspec,s.u.c.y,u2p2)
    du2dyp2 = dy_squared_mean(s.reductionh,s.u.c.y)
    λ22 = sqrt(2*u2p2/du2dyp2)
    #λ22 = taylor_lenght_y(s.yspec,s.hspec,s.u.c.y,u2p2)
    du2dxp2 = dx_squared_mean(s.reductionh,s.u.c.y)
    λ21 = sqrt(2*u2p2/du2dxp2)

    L31 = integral_lenght_x(s.u.c.z,u3p2)
    L32 = integral_lenght_y(s.u.c.z,u3p2)
    L33 = integral_lenght_z(s.u.c.z,u3p2)

    Lhh = (L11 + L22)/2
    Thh = Lhh/urms
    #λl = (λ11 + λ22)/2
    λl = sqrt(2*(u1p2+u2p2)/(du1dxp2+du2dyp2))
    #λt = (λ12 + λ21)/2
    λt = sqrt(2*(u1p2+u2p2)/(du1dyp2+du2dxp2))
    Lhv = (L13 + L23)/2
    Lvv = L33
    Lvh = (L31 + L32)/2

    _otherstats = stats.(getfield.(Ref(s),sim_fields),Ref(s))
    otherstats = getindex.(_otherstats,1)
    otherstats_fil = getindex.(_otherstats,2)

    tdiss = vstats[13]
    tdiss_fil = vstats_fil[13]

    hashyperviscosity(s) && (tdiss += otherstats[5][3])
    hashyperviscosity(s) && (tdiss_fil += otherstats_fil[5][3])
    hasles(s) && (tdiss += otherstats[3][4])
    hasles(s) && (tdiss_fil += otherstats_fil[3][4])

    Lt = sqrt(k*k*k)/tdiss # outer scale
    Tt = k/tdiss # outer time scale
    Lk = sqrt(sqrt(ν*ν*ν/tdiss)) # Kolmogorov length scale
    Tk = sqrt(ν/tdiss) # Kolmogorov timescale
    λg = sqrt(15*ν/tdiss)*urms

    lensv = (Lk,Lt,λg, λl,λt,Lhh,Lvv,Lhv,Lvh,λ12,λ21,L11,L12,L13,L21,L22,L23,L31,L32,L33)

    Reh = uh*Lhh/ν
    Reh1 = (Lhh/Lk)^(4/3)

    Reλg = urms*λg/ν
    Reλg1 = (λg/Lk)^(4/3)

    Reλl = urms*λl/ν
    Reλl1 = (λl/Lk)^(4/3)

    Reλt = urms*λt/ν
    Reλt1 = (λt/Lk)^(4/3)

    if hasdensity(s)
        Nb = calculate_N(s.densitystratification)
        Lb = 2π*urms/Nb # Buoyancy lenght scale
        Lo = sqrt(tdiss/(Nb*Nb*Nb)) # Ozmidov lenght scale
        Lom = sqrt(tdiss/(Nb*Nb*Nb/(8*pi*pi*pi))) # Ozmidov lenght scale
        Tb = 2π/Nb # Buoyancy time scale
        Tbm = 1/Nb # Buoyancy time scale
        Reb = (Lo/Lk)^(4/3)
        Rebm = (Lom/Lk)^(4/3)
        Frh = 2π*uh/(Nb*Lhh)
        Frh1 = Lhv/Lhh
        lens = (lensv...,Lb,Lo,Lom,Tk,Tt,Thh,Tb,Tbm)
        anums = (Reh,Reh1, Reλg, Reλl, Reλt, Reλg1, Reλl1, Reλt1, Reb, Rebm, Frh,Frh1)
    else
        lens = (lensv...,Tk,Tt,Thh)
        anums = (Reh, Reh1, Reλl, Reλt, Reλl1, Reλt1)
    end

    results = (vstats..., tdiss, flatten(otherstats)...)
    results_fil = (vstats_fil..., tdiss_fil, flatten(otherstats_fil)...)
    return results, lens, anums, results_fil
end
#  (velocity_stats(s)..., stats(s.passivescalar,s)..., stats(s.densitystratification,s)..., stats(s.lesmodel,s)..., stats(s.forcing,s)...)

@par function velocity_stats(s::@par(AbstractSimulation))
    if !STAT_FIL
        u1 = real(s.u.c.x[1,1,1])
        u2 = real(s.u.c.y[1,1,1])
        u3 = real(s.u.c.z[1,1,1])

        u12 = squared_mean(s.reductionh, s.u.c.x)
        u22 = squared_mean(s.reductionh, s.u.c.y)
        u32 = squared_mean(s.reductionh, s.u.c.z)
        kh = (u12+u22)/2
        kv = u32/2
        k = (kh+kv)
        kaniso = 2*kv/kh

        εh,εv,ε = viscosity_stats(s.reductionh,s.reductionv,s.u)

        nlh, nlv = s.nlstats[]

        pressh, pressv = s.pressstats[]

        return (u1, u2, u3, u12, u22, u32, kh, kv, kaniso, k, εh, εv, ε, nlh, nlv, pressh, pressv, -(nlh+pressh),(nlv+pressv)), ntuple(x->0.0,19)
    else
        u1 = real(s.u.c.x[1,1,1])
        u2 = real(s.u.c.y[1,1,1])
        u3 = real(s.u.c.z[1,1,1])

        u1_fil,u2_fil,u3_fil = (u1,u2,u3)

        u12,u12_fil = squared_mean(s.reductionh,s.reductionv, s.u.c.x)
        u22,u22_fil = squared_mean(s.reductionh,s.reductionv, s.u.c.y)
        u32,u32_fil = squared_mean(s.reductionh,s.reductionv, s.u.c.z)
        kh = (u12+u22)/2
        kh_fil = (u12_fil+u22_fil)/2
        kv = u32/2
        kv_fil = u32_fil/2
        k = (kh+kv)
        k_fil = (kh_fil+kv_fil)
        kaniso = 2*kv/kh
        kaniso_fil = 2*kv_fil/kh_fil

        (εh,εv,ε),(εh_fil,εv_fil,ε_fil) = viscosity_stats(s.reductionh,s.reductionv,s.reductionh_fil,s.reductionv_fil,s.u)

        #nlh, nlv = non_linear_stats(s.reductionh,s.reductionv,s.u,s.rhs)
        nlh, nlv,nlh_fil, nlv_fil = s.nlstats[]

        #pressh, pressv = pressure_stats(s.reductionh,s.reductionv,s)
        pressh, pressv,pressh_fil, pressv_fil = s.pressstats[]

        return ((u1, u2, u3, u12, u22, u32, kh, kv, kaniso, k, εh, εv, ε, nlh, nlv, pressh, pressv, -(nlh+pressh),(nlv+pressv)),
               (u1_fil, u2_fil, u3_fil, u12_fil, u22_fil, u32_fil, kh_fil, kv_fil, kaniso_fil, k_fil, εh_fil, εv_fil, ε_fil, nlh_fil, nlv_fil, pressh_fil, pressv_fil, -(nlh_fil+pressh_fil),(nlv_fil+pressv_fil)))
    end
end

@par function scalar_stats(ρ,s1,s::@par(AbstractSimulation))
    if !STAT_FIL
        rho = real(ρ[1,1,1])
        rho2 = squared_mean(s1.reduction,ρ)

        drd1 = dx_squared_mean(s1.reduction, ρ)
        drd2 = dy_squared_mean(s1.reduction, ρ)
        drd3 = dz_squared_mean(s1.reduction, ρ)

        rhodiss = diffusivity(s1)*(drd1 + drd2 + drd3)

        return (rho, rho2, drd1, drd2, drd3, rhodiss), ntuple(x->0.0,6)
    else
        rho = real(ρ[1,1,1])
        rho_fil = rho

        rho2, rho2_fil = squared_mean(s1.reduction,s.reductionh,ρ)

        drd1, drd1_fil = dx_squared_mean(s1.reduction,s.reductionh, ρ)
        drd2, drd2_fil = dy_squared_mean(s1.reduction,s.reductionh, ρ)
        drd3, drd3_fil = dz_squared_mean(s1.reduction,s.reductionh, ρ)

        rhodiss = diffusivity(s1)*(drd1 + drd2 + drd3)
        rhodiss_fil = diffusivity(s1)*(drd1_fil + drd2_fil + drd3_fil)

        return (rho, rho2, drd1, drd2, drd3, rhodiss), (rho_fil, rho2_fil, drd1_fil, drd2_fil, drd3_fil, rhodiss_fil)
    end
end

#ape(s::AbstractSimulation) = tmean(x->x^2,parent(real(s.ρ)),s)

#@par function tmean(f::F,x::AbstractArray{T,3},s::@par(AbstractSimulation)) where {F<:Function,T<:Number}
#
#    result = fill!(s.reduction,0.0)
#    @mthreads for k in ZRANGE
#        for j in YRANGE
#            @inbounds @msimd for i in RXRANGE
#                result[Threads.threadid()] += f(x[i,j,k])
#            end
#        end
#    end
#
#    return sum(result)/(NRX*NY*NZ)
#end

#tmean(x::AbstractArray,s::AbstractSimulation) = tmean(identity,x,s)

#@generated function getind(f::F,x::NT,i::Int,j::Int,k::Int) where {F<:Function,N,NT<:NTuple{N,AbstractArray}}
#    args = Array{Any,1}(undef,N)
#    for l=1:N
#       args[l] = :(x[$l][i,j,k])
#    end
#    ex = Expr(:call,:f,args...)
#    return Expr(:block,Expr(:meta,:inline),Expr(:meta,:propagate_inbounds),ex)
#end
#
#@par function tmean(f::F,x::NT,reduction::L) where {F<:Function,N,NT<:NTuple{N,AbstractArray},L}
#    result = fill!(reduction,0.0)
#    @mthreads for k in ZRANGE
#        for j in YRANGE
#            @inbounds @msimd for i in RXRANGE
#                result[Threads.threadid()] += getind(f,x,i,j,k)
#            end
#        end
#    end
#    return sum(result)/(NRX*NY*NZ)
#end

function squared_mean(reduction,u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = abs2(u[i,j,l])
                    ee += (1 + (i>1))*magsq
                end
            result[ii] += ee
        end
    end
    return sum(result)
end

function squared_mean(reduction,reduction_fil, u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = abs2(u[i,j,l])
                    ee += (1 + (i>1))*magsq
                    ee_fil += (1 + (i>1))*magsq*Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    return sum(result), sum(result_fil)
end

function dx_squared_mean(reduction,u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = abs2(im*KX[i]*u[i,j,l])
                    ee += (1 + (i>1))*magsq
                end
            result[ii] += ee
        end
    end
    return sum(result)
end

function dx_squared_mean(reduction, reduction_fil, u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = abs2(im*KX[i]*u[i,j,l])
                    lee = (1 + (i>1))*magsq
                    ee += lee
                    ee_fil = lee*Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    return sum(result), sum(result_fil)
end

function dy_squared_mean(reduction,u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ii = Threads.threadid()
                ky = KY[j]
                @simd for i in XRANGE
                    magsq = abs2(im*ky*u[i,j,l])
                    ee += (1 + (i>1))*magsq
                end
            result[ii] += ee
        end
    end
    return sum(result)
end

function dy_squared_mean(reduction,reduction_fil, u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
                ky = KY[j]
                @simd for i in XRANGE
                    magsq = abs2(im*ky*u[i,j,l])
                    lee = (1 + (i>1))*magsq
                    ee += lee
                    ee_fil += lee*Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    return sum(result), sum(result_fil)
end

function dz_squared_mean(reduction,u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            kz = KZ[l]
            ee = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = abs2(im*kz*u[i,j,l])
                    ee += (1 + (i>1))*magsq
                end
            result[ii] += ee
        end
    end
    return sum(result)
end

function dz_squared_mean(reduction, reduction_fil, u::ScalarField{T}) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            kz = KZ[l]
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = abs2(im*kz*u[i,j,l])
                    lee = (1 + (i>1))*magsq
                    ee += lee
                    ee_fil += lee*Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    return sum(result), sum(result_fil)
end

@inline proj(u::Vec{<:Complex},v::Vec{<:Complex}) = proj(u.x,v.x) + proj(u.y,v.y) + proj(u.z,v.z)
@inline proj(a::SymTen{<:Complex},b::SymTen{<:Complex}) =
    proj(a.xx,b.xx) + proj(a.yy,b.yy) + proj(a.zz,b.zz) + 2*(proj(a.xy,b.xy) + proj(a.xz,b.xz) + proj(a.yz,b.yz))

function proj_mean(reduction,u::AbstractField{T},v::AbstractField{T}) where {T}
    isrealspace(u) && fourier!(u)
    isrealspace(v) && fourier!(v)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = proj(u[i,j,l],v[i,j,l])
                    ee += (1 + (i>1))*magsq
                end
            result[ii] += ee
        end
    end
    return sum(result)
end

function proj_mean(reduction,reduction_fil,u::AbstractField{T},v::AbstractField{T}) where {T}
    isrealspace(u) && fourier!(u)
    isrealspace(v) && fourier!(v)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    magsq = (1 + (i>1))*proj(u[i,j,l],v[i,j,l])
                    ee += magsq
                    ee_fil += magsq*Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    return sum(result), sum(result_fil)
end

@inline mag2(u::Vec{T}) where T<:Complex = muladd(u.x.re, u.x.re, muladd(u.x.im, u.x.im,
                                   muladd(u.y.re, u.y.re, muladd(u.y.im, u.y.im,
                                   muladd(u.z.re, u.z.re, u.z.im*u.z.im)))))

@inline mag2(u::Vec{T}) where T<:Real = u⋅u

function hyperviscosity_stats(reductionh,reductionv,u::VectorField{T},s) where {T}
    isrealspace(u) && fourier!(u)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            eeh = zero(T)
            eev = zero(T)
            ii = Threads.threadid()
            M::Int = get_hyperviscosity_exponent(s)
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    uh = u[i,j,l]
                    eeh += (1 + (i>1)) * (k2^M)*(proj(uh.x,uh.x)+proj(uh.y,uh.y))
                    eev += (1 + (i>1)) * (k2^M)*proj(uh.z,uh.z)
                end
            resulth[ii] += eeh
            resultv[ii] += eev
        end
    end
    mν = nuh(s)
    eh = mν*sum(resulth)
    ev = mν*sum(resultv)
    return eh,ev,eh+ev
end

function hyperviscosity_stats(reductionh,reductionv,reductionh_fil,reductionv_fil,u::VectorField{T},s) where {T}
    isrealspace(u) && fourier!(u)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    resulth_fil = fill!(reductionh_fil,zero(T))
    resultv_fil = fill!(reductionv_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            eeh = zero(T)
            eev = zero(T)
            eeh_fil = zero(T)
            eev_fil = zero(T)
            ii = Threads.threadid()
            M::Int = get_hyperviscosity_exponent(s)
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    uh = u[i,j,l]
                    leeh = (1 + (i>1)) * (k2^M)*(proj(uh.x,uh.x)+proj(uh.y,uh.y))
                    leev = (1 + (i>1)) * (k2^M)*proj(uh.z,uh.z)
                    eeh += leeh
                    eev += leev

                    G2 = Gaussfilter(SPEC_FIL_D2,k2)^2
                    eeh_fil += leeh*G2
                    eev_fil += leev*G2
                end
            resulth[ii] += eeh
            resultv[ii] += eev
            resulth_fil[ii] += eeh_fil
            resultv_fil[ii] += eev_fil
        end
    end
    mν = nuh(s)
    eh = mν*sum(resulth)
    ev = mν*sum(resultv)
    eh_fil = mν*sum(resulth_fil)
    ev_fil = mν*sum(resultv_fil)
    return (eh,ev,eh+ev),(eh_fil,ev_fil,eh_fil+ev_fil)
end

function viscosity_stats(reductionh,reductionv,u::VectorField{T}) where {T}
    isrealspace(u) && fourier!(u)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            eeh = zero(T)
            eev = zero(T)
            ii = Threads.threadid()
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    uh = u[i,j,l]
                    eeh += (1 + (i>1)) * k2*(proj(uh.x,uh.x)+proj(uh.y,uh.y))
                    eev += (1 + (i>1)) * k2*proj(uh.z,uh.z)
                end
            resulth[ii] += eeh
            resultv[ii] += eev
        end
    end
    eh = ν*sum(resulth)
    ev = ν*sum(resultv)
    return eh,ev,eh+ev
end

function viscosity_stats(reductionh,reductionv,reductionh_fil,reductionv_fil,u::VectorField{T}) where {T}
    isrealspace(u) && fourier!(u)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    resulth_fil = fill!(reductionh_fil,zero(T))
    resultv_fil = fill!(reductionv_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            eeh = zero(T)
            eev = zero(T)
            eeh_fil = zero(T)
            eev_fil = zero(T)
            ii = Threads.threadid()
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    uh = u[i,j,l]
                    leeh = (1 + (i>1)) * k2*(proj(uh.x,uh.x)+proj(uh.y,uh.y))
                    leev = (1 + (i>1)) * k2*proj(uh.z,uh.z)
                    eeh += leeh
                    eev += leev
                    G2 = Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                    eeh_fil += leeh*G2
                    eev_fil += leev*G2
                end
            resulth[ii] += eeh
            resultv[ii] += eev
            resulth_fil[ii] += eeh_fil
            resultv_fil[ii] += eev_fil
        end
    end
    eh = ν*sum(resulth)
    ev = ν*sum(resultv)
    eh_fil = ν*sum(resulth_fil)
    ev_fil = ν*sum(resultv_fil)
    return (eh,ev,eh+ev),(eh_fil,ev_fil,eh_fil+ev_fil)
end

function les_stats(reductionh,reductionv,τ::AbstractField{T},u::VectorField{T}) where {T}
    isrealspace(u) && fourier!(u)
    isrealspace(τ) && fourier!(τ)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            eeh = zero(T)
            eev = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    kh = K[i,j,l]
                    out = vecouterproj((im*kh)⋅τ[i,j,l],u[i,j,l])
                    eeh += (1 + (i>1))*(out.x + out.y)
                    eev += (1 + (i>1))*out.z
                end
            resulth[ii] += eeh
            resultv[ii] += eev
        end
    end
    prh = -sum(resulth)
    prv = -sum(resultv)
    return prh,prv,(2prv)/prh,prh+prv
end

function les_stats(reductionh,reductionv,reductionh_fil,reductionv_fil,τ::AbstractField{T},u::VectorField{T}) where {T}
    isrealspace(u) && fourier!(u)
    isrealspace(τ) && fourier!(τ)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    resulth_fil = fill!(reductionh_fil,zero(T))
    resultv_fil = fill!(reductionv_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            eeh = zero(T)
            eev = zero(T)
            eeh_fil = zero(T)
            eev_fil = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    kh = K[i,j,l]
                    out = vecouterproj((im*kh)⋅τ[i,j,l],u[i,j,l])
                    leeh = (1 + (i>1))*(out.x + out.y)
                    leev = (1 + (i>1))*out.z
                    eeh += leeh
                    eev += leev

                    G2 = Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                    eeh_fil += leeh*G2
                    eev_fil += leev*G2
                end
            resulth[ii] += eeh
            resultv[ii] += eev
            resulth_fil[ii] += eeh_fil
            resultv_fil[ii] += eev_fil
        end
    end
    prh = -sum(resulth)
    prv = -sum(resultv)
    prh_fil = -sum(resulth_fil)
    prv_fil = -sum(resultv_fil)
    return (prh,prv,(2prv)/prh,prh+prv), (prh_fil,prv_fil,(2prv_fil)/prh_fil,prh_fil+prv_fil)
end

@par function pressure_stats(resulth,resultv,s::@par(Simulation))
    resulth = fill!(resulth,0.0)
    resultv = fill!(resultv,0.0)

    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        pressure_stats(j,k,resulth,resultv,s)
    end

    return sum(resulth), sum(resultv), 0.0, 0.0
end

@par function pressure_stats(j::Int,k::Int,hr,vr,s::A) where {A<:@par(Simulation)}
    eeh = 0.0
    eev = 0.0
    ii = Threads.threadid()

    u = s.u.c
    rhsv = s.rhs.c
    ke = s.kspec.c

    if hasles(A)
        τ = s.lesmodel.tau.c
    end

    if hasdensity(A)
        ρ = s.densitystratification.ρ.field
        g = gravity(s.densitystratification)
    end

        @inbounds @msimd for i in XRANGE

            kh = K[i,j,k]
            K2 = kh⋅kh
            v = u[i,j,k]
            rhs = rhsv[i,j,k] - (im*ke[i,j,k])*kh

            if hasles(A)
                if !(!is_SandP(A) && is_FakeSmagorinsky(A))
                    rhs += (im*kh)⋅τ[i,j,k]
                end
            end

            if hasdensity(A)
                rhs += ρ[i,j,k]*g
            end

            p1 = ifelse(k==j==i==1,zero(ComplexF64),-(kh⋅rhs)/K2)
            pressure = p1*kh
            out = vecouterproj(v,pressure)

            outxy = out.x + out.y
            outz = out.z

            eeh += (1 + (i>1))*outxy
            eev += (1 + (i>1))*outz
        end

    hr[ii] += eeh
    vr[ii] += eev

    return nothing
end

@par function pressure_stats(resulth,resultv,resulth_fil,resultv_fil,s::@par(Simulation))
    resulth = fill!(resulth,0.0)
    resultv = fill!(resultv,0.0)
    resulth_fil = fill!(resulth_fil,0.0)
    resultv_fil = fill!(resultv_fil,0.0)

    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        pressure_stats(j,k,resulth,resultv,resulth_fil,resultv_fil,s)
    end

    return sum(resulth), sum(resultv), sum(resulth_fil), sum(resultv_fil)
end

@par function pressure_stats(j::Int,k::Int,hr,vr,hr_fil,vr_fil,s::A) where {A<:@par(Simulation)}
    eeh = 0.0
    eev = 0.0
    eeh_fil = 0.0
    eev_fil = 0.0
    ii = Threads.threadid()

    u = s.u.c
    rhsv = s.rhs.c
    ke = s.kspec.c

    if hasles(A)
        τ = s.lesmodel.tau.c
    end

    if hasdensity(A)
        ρ = s.densitystratification.ρ.field
        g = gravity(s.densitystratification)
    end

        @inbounds @msimd for i in XRANGE

            kh = K[i,j,k]
            K2 = kh⋅kh
            v = u[i,j,k]
            rhs = rhsv[i,j,k] - (im*ke[i,j,k])*kh

            if hasles(A)
                if !(!is_SandP(A) && is_FakeSmagorinsky(A))
                    rhs += (im*kh)⋅τ[i,j,k]
                end
            end

            if hasdensity(A)
                rhs += ρ[i,j,k]*g
            end

            p1 = ifelse(k==j==i==1,zero(ComplexF64),-(kh⋅rhs)/K2)
            pressure = p1*kh
            out = vecouterproj(v,pressure)

            outxy = out.x + out.y
            outz = out.z

            leeh = (1 + (i>1))*outxy
            leev = (1 + (i>1))*outz

            eeh += leeh
            eev += leev

            G2 = Gaussfilter(SPEC_FIL_D2,K2)^2

            eeh += leeh*G2
            eev += leev*G2
        end

    hr[ii] += eeh
    vr[ii] += eev
    hr_fil[ii] += eeh_fil
    vr_fil[ii] += eev_fil

    return nothing
end

function non_linear_stats(hout,vout,u,nl,ke)
    fill!(hout,0.0)
    fill!(vout,0.0)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        eeh = 0.0
        eev = 0.0
        ii = Threads.threadid()
            @simd for i in XRANGE
                ∇ = im*K[i,j,k]
                out = vecouterproj(nl[i,j,k] - ∇*ke[i,j,k],u[i,j,k])
                eeh += (1 + (i>1))*(out.x + out.y)
                eev += (1 + (i>1))*out.z
            end
        hout[ii]+=eeh
        vout[ii]+=eev
    end
    return sum(hout), sum(vout), 0.0, 0.0
end

function non_linear_stats(hout,vout,hout_fil,vout_fil,u,nl,ke)
    fill!(hout,0.0)
    fill!(vout,0.0)
    fill!(hout_fil,0.0)
    fill!(vout_fil,0.0)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        eeh = 0.0
        eev = 0.0
        eeh_fil = 0.0
        eev_fil = 0.0
        ii = Threads.threadid()
            @simd for i in XRANGE
                ∇ = im*K[i,j,k]
                out = vecouterproj(nl[i,j,k] - ∇*ke[i,j,k],u[i,j,k])
                leeh = (1 + (i>1))*(out.x + out.y)
                leev = (1 + (i>1))*out.z
                eeh += leeh
                eev += leev
                G2 = Gaussfilter(SPEC_FIL_D2,i,j,k)^2
                eeh_fil += leeh*G2
                eev_fil += leev*G2
            end
        hout[ii]+=eeh
        vout[ii]+=eev
        hout_fil[ii]+=eeh_fil
        vout_fil[ii]+=eev_fil
    end
    return sum(hout), sum(vout), sum(hout_fil), sum(vout_fil)
end

function integral_lenght_x(u,u1p2)
    r = zero(typeof(proj(u[1],u[1])))
    @inbounds for k in ZRANGE
        for j in YRANGE
            r += proj(u[1,j,k],u[1,j,k])
        end
    end
    return LX*r*pi/u1p2
end

function integral_lenght_y(u,u1p2)
    r = zero(typeof(proj(u[1],u[1])))
    @inbounds for k in ZRANGE
        @msimd for i in XRANGE
            r += (1+(i>1))*proj(u[i,1,k],u[i,1,k])
        end
    end
    return LY*r*pi/u1p2
end

function integral_lenght_z(u,u1p2)
    r = zero(typeof(proj(u[1],u[1])))
    @inbounds for j in YRANGE
        @msimd for i in XRANGE
            r += (1+(i>1))*proj(u[i,j,1],u[i,j,1])
        end
    end
    return LZ*r*pi/u1p2
end

function taylor_lenght_x(re,u,u1p2)

    fill!(re.data,0.0)

    @inbounds for k in ZRANGE
        for j in YRANGE
            for i in XRANGE
                re[i] += proj(u[i,j,k],u[i,j,k])/u1p2
            end
        end
    end

    @inbounds for i in XRANGE
        re[i] *= -KX[i]*KX[i]
    end

    brfft!(re)

    return sqrt(-2/re.data[1])
end

function taylor_lenght_y(re,hpart,u,u1p2)

    fill!(re.data,0.0)
    @inbounds for j in YRANGE
        hpart[1,j,1] = 0.0
    end

    @inbounds for k in ZRANGE
        for j in YRANGE
            for i in XRANGE
                hpart[1,j,1] += proj(u[i,j,k],u[i,j,k])/u1p2
            end
        end
    end

    re[1] = 0.0 + 0.0im

    ep = div(NY,2) + 1

    re[ep] = KY[ep]^2*hpart[1,ep,1]
    @inbounds for i in 2:(ep-1)
        re[i] = (-KY[NY-i]*KY[NY-i])*(hpart[1,i,1] + hpart[1,NY-i,1])/2
    end

    brfft!(re)

    return sqrt(-2/re.data[1])
end

function scalar_les_stats(reduction,f::AbstractField{T},ρ::ScalarField{T}) where {T}
    isrealspace(ρ) && fourier!(ρ)
    isrealspace(f) && fourier!(f)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    kh = K[i,j,l]
                    out = proj((im*kh)⋅f[i,j,l],ρ[i,j,l])
                    ee += (1 + (i>1))*out
                end
            result[ii] += ee
        end
    end
    pr = -sum(result)
    return pr
end

function scalar_les_stats(reduction,reduction_fil, f::AbstractField{T},ρ::ScalarField{T}) where {T}
    isrealspace(ρ) && fourier!(ρ)
    isrealspace(f) && fourier!(f)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
                @simd for i in XRANGE
                    kh = K[i,j,l]
                    out = proj((im*kh)⋅f[i,j,l],ρ[i,j,l])
                    lee = (1 + (i>1))*out
                    ee += lee
                    G2 = Gaussfilter(SPEC_FIL_D2,kh⋅kh)^2
                    ee_fil += lee*G2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    pr = -sum(result)
    pr_fil = -sum(result_fil)
    return pr, pr_fil
end

function scalar_hvis_stats(reduction,u::ScalarField{T},s::HyperViscosity) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ii = Threads.threadid()
            M::Int = get_hyperviscosity_exponent(s)
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    uh = u[i,j,l]
                    ee += (1 + (i>1)) * (k2^M)*proj(uh,uh)
                end
            result[ii] += ee
        end
    end
    mν = nuh(s)
    eh = mν*sum(result)
    return eh
end

function scalar_hvis_stats(reduction,reduction_fil, u::ScalarField{T},s::HyperViscosity) where {T}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
            M::Int = get_hyperviscosity_exponent(s)
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    uh = u[i,j,l]
                    lee = (1 + (i>1)) * (k2^M)*proj(uh,uh)
                    ee += lee
                    ee_fil += lee*Gaussfilter(SPEC_FIL_D2,i,j,l)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    mν = nuh(s)
    eh = mν*sum(result)
    eh_fil = mν*sum(result_fil)
    return eh, eh_fil
end

function scalar_hvis_stats(reduction,u::ScalarField{T},hv::SpectralBarrier{ini,cut,F}) where {T,F,ini,cut}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            f = hv.func
            ee = zero(T)
            ii = Threads.threadid()
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    k = sqrt(k2)
                    uh = u[i,j,l]
                    fk = f(k)
                    fk = ifelse(fk == -Inf, 0.0,fk)
                    ee += (1 + (i>1)) * fk*proj(uh,uh)
                end
            result[ii] += ee
        end
    end
    diss = -sum(result)
    return diss
end

function scalar_hvis_stats(reduction, reduction_fil, u::ScalarField{T},hv::SpectralBarrier{ini,cut,F}) where {T,F,ini,cut}
    isrealspace(u) && fourier!(u)
    result = fill!(reduction,zero(T))
    result_fil = fill!(reduction_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            f = hv.func
            ee = zero(T)
            ee_fil = zero(T)
            ii = Threads.threadid()
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    k = sqrt(k2)
                    uh = u[i,j,l]
                    fk = f(k)
                    fk = ifelse(fk == -Inf, 0.0,fk)
                    lee = (1 + (i>1)) * fk*proj(uh,uh)
                    ee += lee
                    ee_fil += lee*Gaussfilter(SPEC_FIL_D2,k2)^2
                end
            result[ii] += ee
            result_fil[ii] += ee_fil
        end
    end
    diss = -sum(result)
    diss_fil = -sum(result_fil)
    return diss, diss_fil
end


text:
# # Multinomial logistic regression

using TransformVariables, LogDensityProblems, DynamicHMC, DynamicHMC.Diagnostics
using MCMCDiagnostics
using Parameters, Statistics, Random, Distributions
using NNlib: softmax
import ForwardDiff # use for automatic differentiation (AD)

"""
Multinomial logistic regression.

For each draw, ``Pr(yᵢ) ∼ softmax(Xᵢ β)``. Uses a `β ∼ MultivariateNormal(0, σI)` prior.

`X` is supposed to include the `1`s for the intercept.
"""
struct MultinomialLogisticRegression{Ty, TX, Tσ}
    y::Ty
    X::TX
    σ::Tσ
end

function (problem::MultinomialLogisticRegression)(θ)
    @unpack y, X, σ = problem
    @unpack β = θ
    num_rows, num_covariates = size(X)
    num_classes = size(β, 2) + 1
    η = X * hcat(zeros(num_covariates), β) # the first column of all zeros corresponds to the base class
    μ = softmax(η; dims=2)
    loglik = sum([logpdf(Multinomial(1, μ[i, :]), y[i, :]) for i = 1:num_rows])
    logpri = sum([logpdf(MultivariateNormal(num_classes - 1, σ), β[i, :]) for i = 1:num_covariates])
    return loglik + logpri
end

# Make up parameters, generate data using random draws.
N = 10_000
# There are two covariates. The first one (the column of all ones) is the intercept.
X = hcat(ones(N), randn(N));
# If we have C classes, then for each covariate we need (C - 1) coefficients.
# we consider the first class to be the "base class"
# and then for each of the other classes, we have a coefficient comparing that class to the base class
# In this example, we have four classes, so we need three coefficients for each covariate.
# There are two covariates, so we will have six coefficients in total.
# the rows of β correspond to the covariates
# e.g. the first row of β corresponds to the first covariate (the intercept)
# e.g. the second row of β corresponds to the second covariate
# the columns of β correspond to classes
# recall that we set the first class as our "base class"
# then the jth column of β contains the coefficients comparing the (j+1) class against the first class
# e.g. the first column of β contains coefficients comparing the second class against the first class
# e.g. the second column of β contains coefficients comparing the third class against the first class
# e.g. the third column of β contains coefficients comparing the fourth class against the first class
β_true = [1.0 2.0 3.0; 4.0 5.0 6.0]
η = X * hcat(zeros(2), β_true);
μ = softmax(η; dims=2);
# y has N rows and C columns
# the rows of y correspond to observations
# the columns of y correspond to classes
y = vcat([rand(Multinomial(1, μ[i,:]))' for i in 1:N]...);

# Create a problem, apply a transformation, then use automatic differentiation.
p = MultinomialLogisticRegression(y, X, 100.0) # data and (vague) priors
t = as((β = as(Array, size(β_true)), )) # identity transformation, just to get the dimension
P = TransformedLogDensity(t, p) # transformed
∇P = ADgradient(:ForwardDiff, P) # use ForwardDiff for automatic differentiation (AD)

# Sample using NUTS, random starting point.
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1_000);

# Extract the posterior. (Here the transformation was not really necessary).
β_posterior = first.(transform.(t, results.chain));

# Check that we recover the parameters.
mean(β_posterior)

function _median(x)
	n = length(x)
	result = similar(first(x))
	for i in eachindex(result)
		result[i] = median([x[j][i] for j = 1:n])
	end
	return result
end

_median(β_posterior)

# Quantiles
qs = [0.05, 0.25, 0.5, 0.75, 0.95]
quantile([β_posterior[i][1, 1] for i in 1:length(β_posterior)], qs)
quantile([β_posterior[i][1, 2] for i in 1:length(β_posterior)], qs)
quantile([β_posterior[i][1, 3] for i in 1:length(β_posterior)], qs)
quantile([β_posterior[i][2, 1] for i in 1:length(β_posterior)], qs)
quantile([β_posterior[i][2, 2] for i in 1:length(β_posterior)], qs)
quantile([β_posterior[i][2, 3] for i in 1:length(β_posterior)], qs)

# Check that mixing is good.
ess = vec(mapslices(effective_sample_size, reduce(hcat, [vec(a) for a in β_posterior]); dims = 2))

ext: jl


text:
@testset "smarts.smarts" begin

@testset "chain" begin
    nullmol = parse(SMARTS, "")
    @test nodecount(nullmol) == 0
    @test edgecount(nullmol) == 0

    aliphatic = parse(SMARTS, "C")
    @test nodeattr(aliphatic, 1).query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :C),
        QueryFormula(:isaromatic, false)
    ]))
    @test edgecount(aliphatic) == 0

    carbonyl = parse(SMARTS, "[CX3]=[OX1]")
    @test nodeattr(carbonyl, 1).query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :C),
        QueryFormula(:isaromatic, false),
        QueryFormula(:connectivity, 3)
    ]))
    @test nodeattr(carbonyl, 2).query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :O),
        QueryFormula(:isaromatic, false),
        QueryFormula(:connectivity, 1)
    ]))

    ether = parse(SMARTS, "[#6][OD2][#6]")
    @test nodeattr(ether, 1).query == QueryFormula(:atomsymbol, :C)
    @test nodeattr(ether, 2).query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :O),
        QueryFormula(:isaromatic, false),
        QueryFormula(:nodedegree, 2)
    ]))

    notH = parse(SMARTS, "[!#1]")
    @test nodeattr(notH, 1).query == QueryFormula(:not, QueryFormula(:atomsymbol, :H))

    imine = parse(SMARTS, raw"[$([CX3]([#6])[#6]),$([CX3H][#6])]=[$([NX2][#6]),$([NX2H])]")
    @test nodeattr(imine, 1).query == QueryFormula(:or, Set([
        QueryFormula(:recursive, "[CX3]([#6])[#6]"),
        QueryFormula(:recursive, "[CX3H][#6]")
    ]))
    @test nodeattr(imine, 2).query == QueryFormula(:or, Set([
        QueryFormula(:recursive, "[NX2][#6]"),
        QueryFormula(:recursive, "[NX2H]")
    ]))
end

end # smarts.SMARTS


text:
Adapt.adapt_structure(to, prob::ODEProblem{<:Any, <:Any, iip}) where {iip} =
    ODEProblem{iip,true}(
        adapt(to, prob.f),
        adapt(to, prob.u0),
        adapt(to, prob.tspan),
        adapt(to, prob.p);
        adapt(to, prob.kwargs)...
    )

## Fixed TimeStep Integrator
mutable struct GPUTsit5Integrator{IIP, S, T, P, F} <: DiffEqBase.AbstractODEIntegrator{GPUTsit5, IIP, S, T}
    f::F                  # eom
    uprev::S              # previous state
    u::S                  # current state
    tmp::S                # dummy, same as state
    tprev::T              # previous time
    t::T                  # current time
    t0::T                 # initial time, only for reinit
    dt::T                 # step size
    tdir::T
    p::P                  # parameter container
    u_modified::Bool
    k1::S                 #intepolants
    k2::S
    k3::S
    k4::S
    k5::S
    k6::S
    k7::S
    cs::SVector{6, T}     # ci factors cache: time coefficients
    as::SVector{21, T}    # aij factors cache: solution coefficients
    rs::SVector{22, T}    # rij factors cache: interpolation coefficients
end
const GPUT5I = GPUTsit5Integrator

DiffEqBase.isinplace(::GPUT5I{IIP}) where {IIP} = IIP

## Adaptive TimeStep Integrator

mutable struct GPUATsit5Integrator{IIP, S, T, P, F, N, TOL, Q} <: DiffEqBase.AbstractODEIntegrator{GPUTsit5, IIP, S, T}
    f::F                  # eom
    uprev::S              # previous state
    u::S                  # current state
    tmp::S                # dummy, same as state
    tprev::T              # previous time
    t::T                  # current time
    t0::T                 # initial time, only for reinit
    tf::T
    dt::T                 # step size
    dtnew::T
    tdir::T
    p::P                  # parameter container
    u_modified::Bool
    k1::S         # interpolants of the algorithm
    k2::S
    k3::S
    k4::S
    k5::S
    k6::S
    k7::S
    cs::SVector{6, T}     # ci factors cache: time coefficients
    as::SVector{21, T}    # aij factors cache: solution coefficients
    btildes::SVector{7,T}
    rs::SVector{22, T}    # rij factors cache: interpolation coefficients
    qold::Q
    abstol::TOL
    reltol::TOL
    internalnorm::N       # function that computes the error EEst based on state
end

const GPUAT5I = GPUATsit5Integrator

#######################################################################################
# Initialization of Integrators
#######################################################################################
@inline function gputsit5_init(f::F, IIP::Bool, u0::S, t0::T, dt::T, p::P) where {F, P, T, S<:AbstractArray{T}}

    cs, as, rs = SimpleDiffEq._build_tsit5_caches(T)

    !IIP && @assert S <: SArray

    integ = GPUT5I{IIP, S, T, P, F}(f, copy(u0), copy(u0), copy(u0), t0, t0, t0, dt, sign(dt), p, true,
                                    copy(u0), copy(u0), copy(u0), copy(u0), copy(u0), copy(u0), copy(u0), cs, as, rs)
end

@inline function gpuatsit5_init(f::F, IIP::Bool, u0::S, t0::T, tf::T, dt::T, p::P, abstol::TOL, reltol::TOL, internalnorm::N) where {F, P, S, T, N, TOL}

    cs, as, btildes, rs = SimpleDiffEq._build_atsit5_caches(T)

    !IIP && @assert S <: SArray

    qoldinit = eltype(S)(1e-4)

    integ = GPUAT5I{IIP, S, T, P, F, N, TOL, typeof(qoldinit)}(f, copy(u0), copy(u0), copy(u0), t0, t0, t0, tf, dt, dt, sign(tf-t0), p, true, copy(u0), copy(u0), copy(u0), copy(u0), copy(u0), copy(u0), copy(u0), cs, as, btildes, rs, qoldinit, abstol, reltol, internalnorm)
end
## GPU solver

function vectorized_solve(prob::ODEProblem, ps::CuVector, alg::GPUSimpleTsit5;
                          dt, saveat = nothing,
                          save_everystep = true,
                          debug = false, kwargs...)
    # if saveat is specified, we'll use a vector of timestamps.
    # otherwise it's a matrix that may be different for each ODE.
    if saveat === nothing
        if save_everystep
            len = length(prob.tspan[1]:dt:prob.tspan[2])
        else
            len = 2
        end
        ts = CuMatrix{typeof(dt)}(undef, (len, length(ps)))
        us = CuMatrix{typeof(prob.u0)}(undef, (len, length(ps)))
    else
        error("Not fully implemented yet")  # see the TODO in the kernel
        ts = saveat
        us = CuMatrix{typeof(prob.u0)}(undef, (length(ts), length(ps)))
    end

    kernel = @cuda launch=false tsit5_kernel(prob, ps, us, ts, dt,
                                             Val(saveat !== nothing), Val(save_everystep))
    if debug
        @show CUDA.registers(kernel)
        @show CUDA.memory(kernel)
    end

    config = launch_configuration(kernel.fun)
    threads = min(length(ps), config.threads)
    # XXX: this kernel performs much better with all blocks active
    blocks = max(cld(length(ps), threads), config.blocks)
    threads = cld(length(ps), blocks)
    kernel(prob, ps, us, ts, dt; threads, blocks)

    # we build the actual solution object on the CPU because the GPU would create one
    # containig CuDeviceArrays, which we cannot use on the host (not GC tracked,
    # no useful operations, etc). That's unfortunate though, since this loop is
    # generally slower than the entire GPU execution, and necessitates synchronization
    #EDIT: Done when using with DiffEqGPU
    ts,us
end

@inline function step!(integ::GPUT5I{false, S, T}) where {T, S}

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dt; t = integ.t; p = integ.p
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as

    tmp = integ.tmp; f = integ.f
    integ.uprev = integ.u; uprev = integ.u

    if integ.u_modified
      k1 = f(uprev, p, t)
      integ.u_modified=false
    else
      @inbounds k1 = integ.k7;
    end

    tmp = uprev+dt*a21*k1
    k2 = f(tmp, p, t+c1*dt)
    tmp = uprev+dt*(a31*k1+a32*k2)
    k3 = f(tmp, p, t+c2*dt)
    tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
    k4 = f(tmp, p, t+c3*dt)
    tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
    k5 = f(tmp, p, t+c4*dt)
    tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
    k6 = f(tmp, p, t+dt)

    integ.u = uprev+dt*((a71*k1+a72*k2+a73*k3+a74*k4)+a75*k5+a76*k6)
    k7 = f(integ.u, p, t+dt)

    @inbounds begin # Necessary for interpolation
        integ.k1 = k7; integ.k2 = k2; integ.k3 = k3
        integ.k4 = k4; integ.k5 = k5; integ.k6 = k6
        integ.k7 = k7
    end

    integ.tprev = t
    integ.t += dt

    return  nothing
end

# saveat is just a bool here:
#  true: ts is a vector of timestamps to read from
#  false: each ODE has its own timestamps, so ts is a vector to write to
function tsit5_kernel(_prob, ps, _us, _ts, dt,
                      ::Val{saveat}, ::Val{save_everystep}) where {saveat, save_everystep}
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    i >= length(ps) && return

    # get the actual problem for this thread
    p = @inbounds ps[i]
    prob = remake(_prob; p)

    # get the input/output arrays for this thread
    ts = if saveat
        _ts
    else
        @inbounds view(_ts, :, i)
    end
    us = @inbounds view(_us, :, i)

    integ = gputsit5_init(prob.f,false,prob.u0, prob.tspan[1], dt, prob.p)
    # TODO: optimize contiguous view to return a CuDeviceArray

    # FSAL
    for i in 2:length(ts)
        step!(integ)
        if saveat
            # TODO
        elseif save_everystep
            @inbounds us[i] = integ.u
            @inbounds ts[i] = integ.t
        end
    end

    if !saveat && !save_everystep
        @inbounds us[2] = integ.u
        @inbounds ts[2] = integ.t
    end

    return nothing
end

#############################Adaptive Version#####################################

function vectorized_asolve(prob::ODEProblem, ps::CuVector, alg::GPUSimpleATsit5;
    dt=0.1f0, saveat=nothing,
    save_everystep=false,
    abstol=1.0f-6, reltol=1.0f-3,
    debug=false, kwargs...)
    # if saveat is specified, we'll use a vector of timestamps.
    # otherwise it's a matrix that may be different for each ODE.
    if saveat === nothing
        if save_everystep
            error("Don't use adaptive version with saveat == nothing and save_everystep = true")
        else
            len = 2
        end
        ts = CuMatrix{typeof(dt)}(undef, (len, length(ps)))
        us = CuMatrix{typeof(prob.u0)}(undef, (len, length(ps)))
    else
        ts = saveat
        us = CuMatrix{typeof(prob.u0)}(undef, (length(ts), length(ps)))
    end

    kernel = @cuda launch = false atsit5_kernel(prob, ps, us, ts, dt, abstol, reltol,
        Val(saveat !== nothing), Val(save_everystep))
    if debug
        @show CUDA.registers(kernel)
        @show CUDA.memory(kernel)
    end

    config = launch_configuration(kernel.fun)
    threads = min(length(ps), config.threads)
    # XXX: this kernel performs much better with all blocks active
    blocks = max(cld(length(ps), threads), config.blocks)
    threads = cld(length(ps), blocks)
    kernel(prob, ps, us, ts, dt, abstol, reltol; threads, blocks)

    # we build the actual solution object on the CPU because the GPU would create one
    # containig CuDeviceArrays, which we cannot use on the host (not GC tracked,
    # no useful operations, etc). That's unfortunate though, since this loop is
    # generally slower than the entire GPU execution, and necessitates synchronization
    #EDIT: Done when using with DiffEqGPU
    ts, us
end

function build_adaptive_tsit5_controller_cache(::Type{T}) where {T}

    beta1 = T(7 / 50)
    beta2 = T(2 / 25)
    qmax = T(10.0)
    qmin = T(1 / 5)
    gamma = T(9 / 10)
    qoldinit = T(1e-4)
    qold = qoldinit

    return beta1, beta2, qmax, qmin, gamma, qoldinit, qold
end

@inline function step!(integ::GPUAT5I{false, S, T}) where {S, T}

    beta1, beta2, qmax, qmin, gamma, qoldinit , _ = build_adaptive_tsit5_controller_cache(eltype(integ.u))
    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dtnew; t = integ.t; p = integ.p; tf = integ.tf
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
    btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes

    tmp = integ.tmp; f = integ.f
    integ.uprev = integ.u; uprev = integ.u

    qold = integ.qold
    abstol = integ.abstol
    reltol = integ.reltol

    if integ.u_modified
      k1 = f(uprev, p, t)
      integ.u_modified=false
    else
      @inbounds k1 = integ.k7;
    end

    EEst = Inf

    while EEst > 1.0
        dt < 1e-14 && error("dt<dtmin")

        tmp = uprev + dt * a21 * k1
        k2 = f(tmp, p, t + c1 * dt)
        tmp = uprev + dt * (a31 * k1 + a32 * k2)
        k3 = f(tmp, p, t + c2 * dt)
        tmp = uprev + dt * (a41 * k1 + a42 * k2 + a43 * k3)
        k4 = f(tmp, p, t + c3 * dt)
        tmp = uprev + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4)
        k5 = f(tmp, p, t + c4 * dt)
        tmp = uprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)
        k6 = f(tmp, p, t + dt)
        u = uprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)
        k7 = f(u, p, t + dt)

        tmp = dt * (btilde1 * k1 + btilde2 * k2 + btilde3 * k3 + btilde4 * k4 +
                    btilde5 * k5 + btilde6 * k6 + btilde7 * k7)
        tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
        EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)

        if iszero(EEst)
            q = inv(qmax)
        else
            q11 = EEst^beta1
            q = q11 / (qold^beta2)
        end

        if EEst > 1
            dt = dt / min(inv(qmin), q11 / gamma)
        else # EEst <= 1
            q = max(inv(qmax), min(inv(qmin), q / gamma))
            qold = max(EEst, qoldinit)
            dtnew = dt/q #dtnew
            dtnew = min(abs(dtnew),abs(tf-t-dt))

            @inbounds begin # Necessary for interpolation
                integ.k1 = k1; integ.k2 = k2; integ.k3 = k3
                integ.k4 = k4; integ.k5 = k5; integ.k6 = k6
                integ.k7 = k7;
            end

            integ.dt = dt
            integ.dtnew = dtnew
            integ.qold = qold
            integ.tprev = t
            integ.u = u

            if (tf - t - dt) < 1e-14
              integ.t = tf
            else
              integ.t += dt
            end
        end
    end
    return nothing
end

function atsit5_kernel(_prob, ps, _us, _ts, dt, abstol, reltol,
    ::Val{saveat}, ::Val{save_everystep}) where {saveat,save_everystep}
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i >= length(ps) && return

    # get the actual problem for this thread
    p = @inbounds ps[i]
    prob = remake(_prob; p)

    # get the input/output arrays for this thread
    ts = if saveat
        _ts
    else
        @inbounds view(_ts, :, i)
    end
    us = @inbounds view(_us, :, i)
    # TODO: optimize contiguous view to return a CuDeviceArray

    u0 = prob.u0
    tspan = prob.tspan
    f = prob.f
    p = prob.p

    t = tspan[1]
    tf = prob.tspan[2]

    cur_t = 0
    if saveat !== nothing
        cur_t = 1
        if tspan[1] == ts[1]
            cur_t += 1
            @inbounds us[1] = u0
        end
    else
        @inbounds ts[1] = tspan[1]
        @inbounds us[1] = u0
    end

    integ = gpuatsit5_init(prob.f, false, prob.u0, prob.tspan[1], prob.tspan[2], dt, prob.p, abstol, reltol, DiffEqBase.ODE_DEFAULT_NORM)

    while integ.t < tspan[2]
        step!(integ)
        if saveat === nothing && save_everystep
            error("Do not use saveat == nothing & save_everystep = true in adaptive version")
        else saveat !== nothing
            while cur_t <= length(ts) && ts[cur_t] <= integ.t
                savet = ts[cur_t]
                θ = (savet - integ.tprev) / integ.dt
                b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = SimpleDiffEq.bθs(integ.rs, θ)
                us[cur_t] = integ.uprev + integ.dt * (
                    b1θ * integ.k1 + b2θ * integ.k2 + b3θ * integ.k3 + b4θ * integ.k4 + b5θ * integ.k5 + b6θ * integ.k6 + b7θ * integ.k7)
                cur_t += 1
            end
        end
    end

    if !saveat && !save_everystep
        @inbounds us[2] = integ.u
        @inbounds ts[2] = integ.t
    end

    return nothing
end


text:
# doc.jl - infrastructure

"""
    process(entries::Array{JournalEntry, 1}, path_journal="./journal.txt", path_ledger="./ledger.txt", )

It adds a `JournalEntry` to the file `path_journal,` converts it into `Record`'s, and adds them to the file `path_ledger.`

@see also [`read_from_file`](@ref)

# Example

```
julia> include("./src/infrastructure/infrastructure.jl")
process (generic function with 1 method)

julia> entry1 = create_journal_entry("20200700", "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20200700", 3, 2020-03-03T13:31:07.953, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> entry2 = create_journal_entry("20201500", "PRG", "A1002", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20201500", 3, 2020-03-03T13:31:07.955, "PRG", "A1002", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> process([entry1, entry2])
```
"""
function process end

"""
    function read_from_file(path::String)::Array{Any, 1}

It reads (structured) data back from a file.

# Example

```
julia> include("./src/infrastructure/infrastructure.jl")
process (generic function with 1 method)

julia> entry1 = create_journal_entry("20200700", "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20200700", 3, 2020-03-03T13:32:46.7, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> entry2 = create_journal_entry("20201500", "PRG", "A1002", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20201500", 3, 2020-03-03T13:32:46.702, "PRG", "A1002", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> process([entry1, entry2])

julia> r = read_from_file("./test_journal.txt")
2-element Array{Any,1}:
 JournalEntry("20200700", 3, 2020-03-03T13:32:46.7, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
 JournalEntry("20201500", 3, 2020-03-03T13:32:46.702, "PRG", "A1002", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
```
"""
function read_from_file end

"""
    add_to_file(path::String, data::Array{T, 1} where T <: Any)

Adds an object to a file.

# Example

```
julia> struct Person
           name:: String
       end

julia> donald = Person("Donald Duck")
Person("Donald Duck")

julia> add_to_file("./persons.txt", [donald])

julia> read_from_file("./persons.txt")
1-element Array{Any,1}:
 Person("Donald Duck")
```
"""
function add_to_file end


text:
struct ZeroPolynomialSet{DT <: AbstractSemialgebraicSet,
                         BT <: MB.AbstractPolynomialBasis
                        } <: MOI.AbstractVectorSet
    domain::DT
    basis::BT
end

# `x`-in-`PlusMinusSet(set)` iff `x`-in-`set` and `-x`-in-`set`.
struct PlusMinusSet{ST <: MOI.AbstractVectorSet} <: MOI.AbstractVectorSet
    set::ST
end


text:
using Test, bad

pnull = .3
pmcr  = .4
p     = Beta(6, 5)
α, β  = .05, .2
problem = Problem(
    minimise( SampleSize(p) ),
    Power(p  | pnull) <= α,
    Power(p >= pmcr)  >= 1 - β
)
design = optimise(problem; verbosity = 0)

mle      = MaximumLikelihoodEstimator()
ordering = EstimatorOrdering(mle)
ci = ClopperPearsonInterval(ordering, design, α)
ci([0], [0])
coverage_probability(ci, [.5])
mean_width(ci, .3)
@test !compatible(ci, design, pnull)["compatible"]

pci = PosteriorCredibleInterval(Beta(3, 7), design, α)
pci(0, 0)
coverage_probability(pci, .5)
mean_width(pci, .3)
@test !compatible(pci, design, pnull)["compatible"]

jprior = JeffreysPrior(design)
jci = PosteriorCredibleInterval(jprior, design, α)
jci(0, 0)
coverage_probability(jci, .5)
mean_width(jci, .3)
@test compatible(jci, design, pnull)["compatible"]

cmle  = CompatibleMLE(design)
cci = ClopperPearsonInterval(EstimatorOrdering(cmle), design, α)
cci([0], [0])
coverage_probability(cci, [.5])
mean_width(cci, .3)
@test compatible(cci, design, pnull)["compatible"]


text:
# ------------------------------------------------------------------------------------------
# ## Get Julia running locally
#
#
# ## Local package installations
#
# If you'd like to run these tutorial notebooks locally, you'll want to install all the
# packages used in them. Since installation can take some time, you may want to run this
# notebook before getting started with the tutorial, rather than trying to install
# everything as you go.
#
# #### Installations
# ------------------------------------------------------------------------------------------

using Pkg
Pkg.add(["Example", "Colors", "Primes", "Plots", "BenchmarkTools"])

# ------------------------------------------------------------------------------------------
# #### Loading all packages
# ------------------------------------------------------------------------------------------

using Example, Colors, Plots, BenchmarkTools, Primes

# ------------------------------------------------------------------------------------------
# #### Tests
#
# `plot` should generate a plot,
# ------------------------------------------------------------------------------------------

plot(x -> x^2, -10:10)

# ------------------------------------------------------------------------------------------
# `RGB(0, 0, 0)` should return a black square,
# ------------------------------------------------------------------------------------------

RGB(0, 0, 0)

# ------------------------------------------------------------------------------------------
# and `@btime primes(1000000);` should report an execution time in ms and memory used. For
# example, on one computer, this yielded "2.654 ms (5 allocations: 876.14 KiB)".
# ------------------------------------------------------------------------------------------

@btime primes(1000000);

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using IntervalLinearAlgebra, LazySets, Plots

A = [2..4 -1..1;-1..1 2..4]
b = [-2..2, -1..1]

Xenclose = solve(A, b)
polytopes = solve(A, b, LinearOettliPrager())

plot(UnionSetArray(polytopes), ratio=1, label="solution set", legend=:top)
plot!(IntervalBox(Xenclose), label="enclosure")


text:
# This file holds all the functions necessary to save a system configuration to
#   the SQLite database.

export saveSimConfigToDatabase,
       readParFileToDatabase

"""
```
saveSimConfigToDatabase( mpSim::ManpowerSimulation )
```
This function stores the configuration parameters of the manpower simulation
`mpSim` in table `config` of the database defined in the simulation.

This function returns `nothing`.
"""
function saveSimConfigToDatabase( mpSim::ManpowerSimulation )::Void

    saveSimConfigToDatabase( mpSim, mpSim.dbName )

end  # saveSimConfigToDatabase( mpSim::ManpowerSimulation )

"""
```
saveSimConfigToDatabase( mpSim::ManpowerSimulation,
                         dbName::String,
                         configName::String = "config" )
```
This function stores the configuration parameters of the manpower simulation
`mpSim` in the database with filename `dbName`, in the table with name
`configName`. If the database filename doesn't have the proper extension,
`.sqlite`, it will be appended to the name.

This function returns `nothing`.
"""
function saveSimConfigToDatabase( mpSim::ManpowerSimulation, dbName::String,
    configName::String = "config" )::Void

    tmpDBname = ( dbName == "" ) || endswith( dbName, ".sqlite" ) ?
        dbName : dbName * ".sqlite"
    configDB = SQLite.DB( tmpDBname )

    createConfigTable( configName, configDB )

    command = "BEGIN TRANSACTION"
    SQLite.execute!( configDB, command )

    # try
        readGeneralParsFromSim( mpSim, configDB, configName )
        readAttributesFromSim( mpSim, configDB, configName )
        readAttritionFromSim( mpSim, configDB, configName )
        readStatesFromSim( mpSim, configDB, configName )
        readCompoundStatesFromSim( mpSim, configDB, configName )
        readTransTypesFromSim( mpSim, configDB, configName )
        readTransitionsFromSim( mpSim, configDB, configName )
        readRecruitmentFromSim( mpSim, configDB, configName )
        readRetirementFromSim( mpSim, configDB, configName )
    # catch errType
    #     command = "COMMIT"
    #     SQLite.execute!( configDB, command )
    #     error( errType )
    # end

    command = "COMMIT"
    SQLite.execute!( configDB, command )

    return

end  # saveSimConfigToDatabase( mpSim, dbName, configName )

"""
```
readParFileToDatabase( fileName::String,
                       dbName::String,
                       configName::String = "config" )
```
This function reads the Excel file with filename `fileName`, processes the
parameters, and stores them in the `config` of the SQLite database with filename
`dbName`. If these filenames do not have the proper extension, `.xlsx` for the
Excel and `.sqlite` for the database, it will be apended to the name.

This function returns `nothing`. If the Excel file doesn't exist, this function
will throw an error.
"""
function readParFileToDatabase( fileName::String, dbName::String,
    configName::String = "config" )::Void

    newMPsim = ManpowerSimulation( fileName )
    saveSimConfigToDatabase( newMPsim, dbName, configName )
    return

end  # readParFileToDatabase( fileName, dbName )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

relationEntries = Dict{Function, String}( Base.:(==) => "IS",
    Base.:(!=) => "IS NOT", Base.:∈ => "IN", Base.:∉ => "NOT IN",
    Base.:> => ">", Base.:(>=) => ">=", Base.:< => "<", Base.:(<=) => "<=" )

function createConfigTable( configName::String, configDB::SQLite.DB )::Void

    SQLite.drop!( configDB, configName, ifexists = true )
    command = "CREATE TABLE $configName(
        parName VARCHAR( 32 ),
        parType VARCHAR( 32 ),
        intPar MEDIUMINT,
        realPar FLOAT,
        boolPar VARCHAR( 5 ),
        strPar TEXT
    )"
    SQLite.execute!( configDB, command )

    return

end  # createConfigTable( configName, configDB )

function readGeneralParsFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    dbCommits = Int( mpSim.simLength / mpSim.commitFrequency )
    command = "INSERT INTO config
        (parName, parType, intPar, realPar, strPar) VALUES
        ('Config file', 'General', NULL, NULL, '$(mpSim.parFileName)'),
        ('DB name', 'General', NULL, NULL, '$(mpSim.dbName)'),
        ('Catalogue name', 'General', NULL, NULL, '$(mpSim.catFileName)'),
        ('Sim name', 'General', NULL, NULL, '$(mpSim.simName)'),
        ('ID key', 'General', NULL, NULL, '$(mpSim.idKey)'),
        ('Personnel target', 'General', '$(mpSim.personnelTarget)', NULL, NULL),
        ('Start date', 'General', NULL, NULL, '$(mpSim.simStartDate)'),
        ('Sim length', 'General', NULL, '$(mpSim.simLength)', NULL),
        ('Sim time', 'General', NULL, $(min( now( mpSim ), mpSim.simLength)), NULL),
        ('DB commits', 'General', '$dbCommits', NULL, NULL)"
    SQLite.execute!( configDB, command )

    return

end  # readGeneralParsFromSim( mpSim, configDB, configName )

function readAttributesFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Don't read anything if there aren't any attributes defined.
    if isempty( mpSim.initAttrList ) && isempty( mpSim.otherAttrList )
        return
    end  # if isempty( mpSim.initAttrList ) && ...

    # Read every single attribute.
    attrInserts = map( attr -> saveAttributeToDatabase( attr ),
        vcat( mpSim.initAttrList, mpSim.otherAttrList ) )
    command = "INSERT INTO $configName
        (parName, parType, boolPar, strPar) VALUES
        $(join( attrInserts, ", " ))"
    SQLite.execute!( configDB, command )
    return

end  # readAttributesFromSim( mpSim, configDB, configName )

function saveAttributeToDatabase( attr::PersonnelAttribute )::String

    valueList = "[" * join( attr.possibleValues, "," ) * "];[" *
        join( map( val -> "$val:$(attr.values[ val ])", keys( attr.values ) ),
        "," ) * "]"
    return "('$(attr.name)', 'Attribute', '$(attr.isFixed)', '$valueList')"

end  # saveAttrToDatabase( attr )

function readAttritionFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    attrInserts = map( attrName -> saveAttritionToDatabase(
        mpSim.attritionSchemes[ attrName ]), keys( mpSim.attritionSchemes ) )
    push!( attrInserts, saveAttritionToDatabase(
        mpSim.defaultAttritionScheme ) )
    command = "INSERT INTO $configName
        (parName, parType, realPar, strPar) VALUES
        $(join( attrInserts, ", " ))"
    SQLite.execute!( configDB, command )
    return

end  # readAttritionFromSim( mpSim, configDB, configName )

function saveAttritionToDatabase( attrScheme::Attrition )::String

    attrCurve = join( map( ii -> "$(attrScheme.attrCurvePoints[ ii ]):$(attrScheme.attrRates[ ii ])",
        eachindex( attrScheme.attrCurvePoints ) ), "," )
    return "('$(attrScheme.name)', 'Attrition', $(attrScheme.attrPeriod), '$attrCurve')"

end  # saveAttritionToDatabase( attrScheme )

function readStatesFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Don't read anything if there aren't any states defined.
    if isempty( mpSim.stateList )
        return
    end  # if isempty( mpSim.stateList )

    # Read every single state.
    stateInserts = map( state -> saveStateToDatabase( mpSim.stateList[ state ],
        mpSim ), keys( mpSim.stateList ) )
    command = "INSERT INTO $configName
        (parName, parType, intPar, realPar, boolPar, strPar) VALUES
        $(join( stateInserts, ", " ))"
    SQLite.execute!( configDB, command )
    return

end  # readStatesFromSim( mpSim, configDB, configName )

function saveStateToDatabase( state::State, mpSim::ManpowerSimulation )::String

    reqList = collect( keys( state.requirements ) )
    map!( attr -> "$attr:$(join( state.requirements[ attr ], "//" ))", reqList,
        reqList )
    stateEntry = "('$(state.name)', 'State', $(state.stateTarget), "
    stateEntry *= "$(state.stateRetAge), '$(state.isInitial)', '["
    stateEntry *= join( reqList, "," ) * "];"

    if state.attrScheme == mpSim.defaultAttritionScheme
        stateEntry *= "default"
    elseif haskey( mpSim.attritionSchemes, state.attrScheme.name )
        stateEntry *= state.attrScheme.name
    else
        attrScheme = state.attrScheme
        stateEntry *= "$(attrScheme.attrPeriod):$(attrScheme.attrRates[ 1 ])"
    end

    stateEntry *= "')"

    return stateEntry

end  # saveStateToDatabase( state, attrName, configDB, configName )

function readCompoundStatesFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Don't read anything if there aren't any compound states defined.
    if isempty( mpSim.compoundStates ) && isempty( mpSim.compoundStateList )
        return
    end  # if isempty( mpSim.stateList )

    # Read every single compound state.
    compStateInserts = vcat( map( compState -> saveCompStateToDatabase(
        mpSim.compoundStates[ compState ], mpSim ),
        keys( mpSim.compoundStates ) ), map(
        compState -> saveCompStateToDatabase(
        mpSim.compoundStateList[ compState ], mpSim ),
        keys( mpSim.compoundStateList ) ) )
    command = "INSERT INTO $configName
        (parName, parType, intPar, boolPar, strPar) VALUES
        $(join( compStateInserts, ", " ))"
    SQLite.execute!( configDB, command )
    return

end  # readCompoundStatesFromSim( mpSim, configDB, configName )

function saveCompStateToDatabase( state::State,
    mpSim::ManpowerSimulation )::String

    reqList = collect( keys( state.requirements ) )
    map!( attr -> "$attr:$(state.requirements[ attr ][ 1 ])", reqList, reqList )
    stateEntry = "('$(state.name)', 'Compound state', $(state.stateTarget), "
    stateEntry *= "'false', '" * join( reqList, ";" ) * "')"
    return stateEntry

end  # saveCompStateToDatabase( state, mpSim )

function saveCompStateToDatabase( compState::CompoundState,
    mpSim::ManpowerSimulation )::String

    stateEntry = "('$(compState.name)', 'Compound state', "
    stateEntry *= "$(compState.stateTarget), 'false', '"
    stateEntry *= join( compState.stateList, ";" ) * "')"
    return stateEntry

end  # saveCompStateToDatabase( state, mpSim )

function readTransTypesFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Don't do anything if there aren't any transition types defined.
    if isempty( mpSim.transList )
        return
    end  # if isempty( transInserts )

    transTypeInserts = map( keys( mpSim.transList ) ) do trans
        return "('$trans', 'Transition type', $(mpSim.transList[ trans ]))"
    end  # map( keys( transList ) ) do

    command = "INSERT INTO $configName
        (parName, parType, intPar) VALUES
        $(join( transTypeInserts, ", " ))"
    SQLite.execute!( configDB, command )

    return

end  # readTransTypesFromSim( mpSim, configDB, configName )

function readTransitionsFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    transInserts = Vector{String}()
    stateList = merge( mpSim.initStateList, mpSim.otherStateList )

    for state in keys( stateList )
        for trans in stateList[ state ]
            entry = saveTransitionToDatabase( trans )
            push!( transInserts, saveTransitionToDatabase( trans ) )
        end  # for trans in stateList[ state ]
    end  # for state in keys( stateList )

    # Don't do anything if there aren't any transitions defined.
    if isempty( transInserts )
        return
    end  # if isempty( transInserts )

    command = "INSERT INTO $configName
        (parName, parType, strPar) VALUES
        $(join( transInserts, ", " ))"
    SQLite.execute!( configDB, command )

    return

end  # readTransitionsFromSim( mpSim, configDB, configName )

function saveTransitionToDatabase( trans::Transition )::String

    transEntry = "('$(trans.name)', 'Transition', '"
    transEntry *= "$(trans.startState.name);$(trans.endState.name);"
    transEntry *= "$(trans.freq);$(trans.offset);["
    condList = Vector{String}( length( trans.extraConditions ) )

    for ii in eachindex( trans.extraConditions )
        cond = trans.extraConditions[ ii ]
        condRel = relationEntries[ cond.rel ]
        condList[ ii ] = "$(cond.attr):$condRel:"

        if isa( cond.val, Vector )
            condList[ ii ] *= join( cond.val, "//" )
        else
            condList[ ii ] *= string( cond.val )
        end  # if isa( cond.val, Vector )
    end  # for cond in trans.extraConditions

    transEntry *= join( condList, "," ) * "];["
    changeList = Vector{String}( length( trans.extraChanges ) )

    for ii in eachindex( trans.extraChanges )
        attr = trans.extraChanges[ ii ]
        changeList[ ii ] = "$(attr.name):" * collect( keys( attr.values ) )[ 1 ]
    end  # for ii in eachindex( trans.extraChanges )

    transEntry *= join( changeList, "," ) * "];["
    transEntry *= join( trans.probabilityList, "," ) * "];"
    transEntry *= "$(trans.maxAttempts);$(trans.isFiredOnFail);"
    transEntry *= "$(trans.maxFlux);$(trans.hasPriority)')"
    return transEntry

end  # saveTransitionToDatabase( trans::Transition )

function readRecruitmentFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Don't do anything if there aren't any recruitment schemes defined.
    if isempty( mpSim.recruitmentSchemes )
        return
    end  # if isempty( mpSim.recruitmentSchemes )

    # Read every single recruitment scheme.
    recInserts = map( recScheme -> saveRecruitmentToDatabase( recScheme ),
        mpSim.recruitmentSchemes )
    command = "INSERT INTO $configName
        (parName, parType, boolPar, strPar) VALUES
        $(join( recInserts, ", " ))"
    SQLite.execute!( configDB, command )
    return

end  # readRecruitmentFromSim( mpSim, configDB, configName )

function saveRecruitmentToDatabase( recScheme::Recruitment )::String

    recEntry = "('$(recScheme.name)', 'Recruitment', '$(recScheme.isAdaptive)',"
    recEntry *= "'$(recScheme.recruitFreq);$(recScheme.recruitOffset);"
    recEntry *= "$(recScheme.recState);"

    if recScheme.isAdaptive
        recEntry *= "$(recScheme.minRecruit);$(recScheme.maxRecruit);"
    else
        recEntry *= "$(recScheme.recDistType);["
        recDistList = Vector{Any}( collect( keys( recScheme.recDistNodes ) ) )

        for ii in eachindex( recDistList )
            node = recDistList[ ii ]
            recDistList[ ii ] = "$node:$(recScheme.recDistNodes[ node ])"
        end  # for ii in eachindex( recDistList )

        recEntry *= join( recDistList, "," ) * "];"
    end  # if recScheme.isAdpative

    recEntry *= "$(recScheme.ageDistType);["
    ageDistList = Vector{Any}( collect( keys( recScheme.ageDistNodes ) ) )

    for ii in eachindex( ageDistList )
        node = ageDistList[ ii ]
        ageDistList[ ii ] = "$(node):$(recScheme.ageDistNodes[ node ])"
    end  # for ii in eachindex( ageDistList )

    recEntry *= join( ageDistList, "," ) * "]')"
    return recEntry

end  # saveAttrToDatabase( attr, configDB, configName )

function readRetirementFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Only process if there is a retirement scheme.
    if mpSim.retirementScheme === nothing
        return
    end  # if mpSim.retirementScheme === nothing

    retScheme = mpSim.retirementScheme

    command = "INSERT INTO $configName
        (parName, parType, boolPar, strPar) VALUES
        ('Retirement', 'Retirement', '$(retScheme.isEither)',
            '$(retScheme.maxCareerLength);$(retScheme.retireAge);$(retScheme.retireFreq);$(retScheme.retireOffset)')"
    SQLite.execute!( configDB, command )

    return

end  # readRetirementFromSim( mpSim, configDB, configName )


text:
abstract type Hello end

struct HelloWorld <: Hello
    name::String
    HelloWorld(s) = new(s)
end

struct HelloTime <: Hello
    name::String
    tnew::DateTime
    HelloTime(s) = new(s, now())
end

sayhello(hlo) = println("Hello to this world, $(hlo.name)!")

sayhello(hlo::HelloTime) = println("It is now $(now()). Hello from back in $(hlo.tnew), $(hlo.name)!")

h1 = HelloWorld("world")
h2 = HelloTime("new world")

sayhello(h1)
sayhello(h2)

fh = open("objects.dat", "w")
serialize(fh, h1)
serialize(fh,h2)
close(fh)

sleep(10)

fh = open("objects.dat", "r")
hh1 = deserialize(fh)
hh2 = deserialize(fh)
close(fh)

sayhello(hh1)
sayhello(hh2)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Covariance Matrix Adaptation Evolution Strategy Implementation: (μ/μ_I,λ)-CMA-ES

The constructor takes following keyword arguments:

- `μ` is the number of parents
- `λ` is the number of offspring
- `τ` is a time constant for a direction vector `s`
- `τ_c` is a time constant for a covariance matrix `C`
- `τ_σ` is a time constant for a global step size `σ`
"""
@kwdef struct CMAES{TT} <: AbstractOptimizer
    μ::Int = 1
    λ::Int = μ+1
    τ::TT    = NaN
    τ_c::TT  = NaN
    τ_σ::TT  = NaN
end
population_size(method::CMAES) = method.μ
default_options(method::CMAES) = (iterations=1500, abstol=1e-10)

mutable struct CMAESState{T, TI, TT} <: AbstractOptimizerState
    N::Int
    τ::TT
    τ_c::TT
    τ_σ::TT
    fitpop::Vector{T}
    C::Matrix{T}
    s::Vector{T}
    s_σ::Vector{T}
    σ::T
    parent::TI
    fittest::TI
    function CMAESState(N::Int, τ::T1, τ_c::T2, τ_σ::T3, fitpop::Vector{T}, C::Matrix{T},
        s::Vector{T}, s_σ::Vector{T}, σ::T, parent::TI, fittest::TI) where {T, TI, T1, T2, T3}
        TP = promote_type(T1,T2,T3)
        new{T,TI,TP}(N, TP(τ), TP(τ_c), TP(τ_σ), fitpop, C,
            s, s_σ, σ, parent, fittest)
    end
end
value(s::CMAESState) = first(s.fitpop)
minimizer(s::CMAESState) = s.fittest

"""Initialization of CMA-ES algorithm state"""
function initial_state(method::CMAES, options, objfun, population)
    @unpack μ,λ,τ,τ_c,τ_σ = method
    @assert μ < λ "Offspring population must be larger then parent population"

    T = typeof(value(objfun))
    individual = first(population)
    N = length(individual)

    # setup time constraints
    τ = isnan(τ) ? sqrt(N) : τ
    τ_c = isnan(τ_c) ? N^2 : τ_c
    τ_σ = isnan(τ_σ) ? sqrt(N) : τ_σ

    # setup initial state
    return CMAESState(N, τ, τ_c, τ_σ,
            fill(convert(T, Inf), μ),
            diagm(0=>ones(T,N)),
            zeros(T, N), zeros(T, N), one(T),
            copy(individual), copy(individual) )
end

function update_state!(objfun, state, population::AbstractVector{IT}, method::CMAES) where {IT}
    @unpack μ,λ,τ,τ_c,τ_σ = method
    N,σ,τ,τ_c,τ_σ = state.N, state.σ, state.τ, state.τ_c, state.τ_σ

    E = zeros(N, λ)
    W = zeros(N, λ)
    offspring = Array{IT}(undef, λ)
    fitoff = fill(Inf, λ)

    SqrtC = (state.C + state.C') / 2.0
    try
        SqrtC = cholesky(SqrtC).U
    catch ex
        @error "Break on Cholesky: $ex: $(state.C)"
        return true
    end

    for i in 1:λ
        # offspring are generated by transforming standard normally distributed random vectors using a transformation matrix
        E[:,i] = randn(N)
        W[:,i] = σ * (SqrtC * E[:,i])
        offspring[i] = state.parent + W[:,i]   # (L1)
        fitoff[i] = value(objfun, offspring[i]) # Evaluate fitness
    end

    # Select new parent population
    idx = sortperm(fitoff)[1:μ]
    for i in 1:μ
        population[i] = offspring[idx[i]]
        state.fitpop[i] = fitoff[idx[i]]
    end

    w = vec(mean(W[:,idx], dims=2))
    ɛ = vec(mean(E[:,idx], dims=2))
    state.parent += w     #  forming recombinant perent for next generation (L2)
    state.s = (1.0 - 1.0/τ)*state.s + (sqrt(μ/τ * (2.0 - 1.0/τ))/σ)*w     # (L3)
    state.C = (1.0 - 1.0/τ_c).*state.C + (state.s./τ_c)*state.s'          # (L4)
    state.s_σ = (1.0 - 1.0/τ_σ)*state.s_σ + sqrt(μ/τ_σ*(2.0 - 1.0/τ_σ))*ɛ # (L5)
    state.σ = σ*exp(((state.s_σ'*state.s_σ)[1] - N)/(2*N*sqrt(N)))

    state.fittest = population[1]

    return false
end


text:
@testset "URIAuthority" begin
    import Spirit.URIAuthority
    tests = Dict((
        "google.com" => [nothing, "google.com", nothing],
        "google.com:8080" => [nothing, "google.com", UInt(8080)],
        "me@google.com:8080" => ["me", "google.com", UInt(8080)],
        "user:pass@danger.somewhere.com" => ["user:pass", "danger.somewhere.com", nothing],
        "255.255.255.255" => [nothing, IPv4("255.255.255.255"), nothing],
        "255.255.255.256" => [nothing, "255.255.255.256", nothing],
        "[::1]" => [nothing, IPv6("::1"), nothing],
        "[::1]:456" => [nothing, IPv6("::1"), UInt(456)]
    ))

    for (str, params) in tests
        authority = URIAuthority(str)
        expected = URIAuthority(params...)
        @test authority.userinfo == expected.userinfo
        @test authority.host == expected.host
        @test authority.port == expected.port

        @test string(authority) == str
    end
end

@testset "URI" begin
    import Spirit.URI
    import Spirit.@uri_str

    uri = uri"https://bro@www.coolsite.com/path/to/awesome?pie=yes#partytime"
    @test uri.scheme == "https"
    @test uri.authority.userinfo == "bro"
    @test uri.authority.host == "www.coolsite.com"
    @test uri.authority.port == nothing
    @test uri.path == "/path/to/awesome"
    @test uri.query == "pie=yes"
    @test uri.fragment == "partytime"
    @test string(uri) == "https://bro@www.coolsite.com/path/to/awesome?pie=yes#partytime"

    uri = uri"urn:starwars:battles:hoth"
    @test uri.scheme == "urn"
    @test uri.authority == nothing
    @test uri.path == "starwars:battles:hoth"
    @test uri.query == nothing
    @test uri.fragment == nothing
    @test string(uri) == "urn:starwars:battles:hoth"

    uri = uri"/jump/to/conclusions?mat=ofcourse?"
    @test uri.scheme == nothing
    @test uri.authority == nothing
    @test uri.path == "/jump/to/conclusions"
    @test uri.query == "mat=ofcourse?"
    @test uri.fragment == nothing
    @test string(uri) == "/jump/to/conclusions?mat=ofcourse?"
end


text:
#common definitions for general usage

const Float = Float64
# const Int = Int64
# const UInt = UInt64

# float and int types for storing all-pairs shortest paths data
# can reduce precision in order to reduce memory used
const FloatSpTime = Float64 # precision for storing shortest path travel times
const IntRNode = Int16 # precision for storing number of nodes in rGraph
const IntFadj = Int8 # precision for storing maximum number of nodes adjacent to any node in rGraph (= length of longest vector in network.rGraph.fadjList)

const sourcePath = @__DIR__

# run modes
const debugMode = false
const checkMode = true # for data checking, e.g. assertions that are checked frequently

# file chars
const delimiter = ','
const newline = "\r\n"

# misc null values
const nullIndex = -1
const nullX = -200 # outside lat/lon range
const nullY = -200
const nullTime = -1.0
const nullDist = -1.0
const nullValue = -1.0
const nullString = ""

nullFunction() = nothing

max_stars_repo_name: jprendo/TaaSSim


text:
# A simple example to test out Convex.jl

using Convex
using GLPK
using Printf

A = [1.0 1.0 1.0]
c = [1.0, 2.0, -1.0]
b = [1.0]

# the variables that get optimized by Convex.jl
y = Variable(1)
s = Variable(3)

# set up the objective (thing to optimize)
objective   = dot(b, y)

# set up the constrains (things that the Variable instances must obey)
constraints = [s >= 0, s == c - transpose(A)*y]

# create the problem and solve it.
problem = Convex.maximize(objective, constraints)
solve!(problem, GLPK.Optimizer)

# get the results
println(problem.optval)
println(y.value)
println(s.value)


text:
module RVFLMod

using DataFrames
using LinearAlgebra
using Random
using ..DataMod

export RVFL, enhanced, radbas

randomrange(rng::AbstractRNG, C::Number, dims::Integer...)::VecOrMat{Float64} =
    2C * rand(rng, dims...) .- C
enames(n::Integer)::Vector{Symbol} = Symbol.('e', 1:n)

# RVFL Input layer
struct RVFL
    m::Integer
    A::Matrix{Float64}
    b::Vector{Float64}
    g::Function

    RVFL(in::Integer, enh::Integer, rng::AbstractRNG = MersenneTwister(); range::Float64 = 1.0, act::Function=radbas) =
        new(enh, randomrange(rng, range, enh, in), randomrange(rng, range, enh), act)
end

enhanced(m::RVFL, x::Union{Vector{Float64}, SubArray{Float64}}) =
    [x; map(m.g, [dot(m.A[i, :], x) for i in 1:m.m] .+ m.b)]

enhanced(m::RVFL, x::Matrix{Float64}) =
    mapslices(row -> enhanced(m, row), Matrix(x), dims = 2)

enhanced(m::RVFL, xdf::DataFrame) =
    DataFrame(enhanced(m, Matrix(xdf)), [propertynames(xdf); enames(m.m)])

enhanced(m::RVFL, d::Data) =
    Data(d, hcat(enhanced(m, d.xdf), d.ydf), [propertynames(d.xdf); enames(m.m)])

# Activation functions
radbas(x::Float64) = exp(-x^2)

end

using .RVFLMod

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
struct CPAParam <: EoSParam
    a::PairParam{Float64}
    b::PairParam{Float64}
    c1::SingleParam{Float64}
    Tc::SingleParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type CPAModel <: EoSModel end
@newmodel CPA CPAModel CPAParam

export CPA
function CPA(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params,sites = getparams(components, ["SAFT/CPA", "properties/molarmass.csv","properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    Mw  = params["Mw"]
    k  = params["k"]
    Tc = params["Tc"]
    c1 = params["c1"]
    params["a"].values .*= 1E-1
    params["b"].values .*= 1E-3
    a  = epsilon_LorentzBerthelot(params["a"], k)
    b  = sigma_LorentzBerthelot(params["b"])
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]
    packagedparams = CPAParam(a, b, c1, Tc, epsilon_assoc, bondvol,Mw)
    references = ["10.1021/ie051305v"]

    model = CPA(packagedparams, sites, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::CPAModel, V, T, z)
    return @f(a_SRK) + @f(a_assoc) + log(V)  # + f(x)
end

function a_SRK(model::CPAModel, V, T, z)
    x = z/∑(z)
    n = ∑(z)
    a = model.params.a.values
    b = model.params.b.values
    Tc = model.params.Tc.values
    c1 = model.params.c1.values

    α = @. (1+c1*(1-√(T/Tc)))^2

    āᾱ = ∑(a .* .√(α * α') .* (x * x'))
    b̄ = ∑(b .* (x * x'))

    return -log(V-n*b̄) - āᾱ/(R̄*T*b̄)*log(1+n*b̄/V)
end
#same as SRK
function ab_consts(model::CPAModel)
    Ωa =  1/(9*(2^(1/3)-1))
    Ωb = (2^(1/3)-1)/3
    return Ωa,Ωb
end

# function cubic_ab(model::CPAModel, T, z)
#     x = z/∑(z)
#     n = ∑(z)
#     a = model.params.a.values
#     b = model.params.b.values
#     Tc = model.params.Tc.values
#     c1 = model.params.c1.values

#     α = @. min((1+c1*(1-√(T/Tc)))^2,one(T))

#     āᾱ = ∑(a .* .√(α * α') .* (x * x'))
#     b̄ = ∑(b .* (x * x'))

#     return āᾱ ,b̄
# end

function a_assoc(model::CPAModel, V, T, z)
    x = z/∑(z)
    n_sites = model.sites.n_flattenedsites
    X_ = @f(X)
    return ∑(x[i]*∑(n_sites[i][a] * (log(X_[i][a])+(1-X_[i][a])/2) for a in @sites(i)) for i in @comps)
end

function X(model::CPAModel, V, T, z)
    _1 = one(V+T+first(z))
    x = z/∑(z)
    ρ = ∑(z)/V
    n_sites = model.sites.n_flattenedsites
    itermax = 100
    dampingfactor = 0.5
    error = 1.
    tol = model.absolutetolerance
    iter = 1
    X_ = [[_1 for a ∈ @sites(i)] for i ∈ @comps]
    X_old = deepcopy(X_)
    while error > tol
        iter > itermax && error("X has failed to converge after $itermax iterations")
        for i ∈ @comps, a ∈ @sites(i)
            rhs = 1/(1+∑(ρ*x[j]*∑(n_sites[j][b]*X_old[j][b]*@f(Δ,i,j,a,b) for b ∈ @sites(j)) for j ∈ @comps))
            X_[i][a] = (1-dampingfactor)*X_old[i][a] + dampingfactor*rhs
        end
        error = sqrt(∑(∑((X_[i][a] - X_old[i][a])^2 for a ∈ @sites(i)) for i ∈ @comps))
        X_old = deepcopy(X_)
        iter += 1
    end
    return X_
end

function Δ(model::CPAModel, V, T, z, i, j, a, b)
    ϵ_associjab = model.params.epsilon_assoc.values[i,j][a,b] * 1e2/R̄
    βijab = model.params.bondvol.values[i,j][a,b] * 1e-3
    x = z/∑(z)
    b = model.params.b.values
    b̄ = ∑(b .* (x * x'))
    η = b̄*∑(z)/(4*V)
    g = (1-η/2)/(1-η)^3
    bij = (b[i,i]+b[j,j])/2
    return g*(exp(ϵ_associjab/T)-1)*βijab*bij
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

l = 1.0       # Total length [m]
N = 5         # Number of nodes
els = N - 1   # Number of finite elements (in x direction)
nod = 2       # Number of nodes per finite elements
nodof = 1     # Degrees of freedom for each node
np_types = 1  # Number of property types
EA = 1.0e5    # Strain stiffness
nip = 1       # Number of integration points

data = Dict(
  # StructuralElement(nxe, np_types, nip, FiniteElement(nod, nodof))
  :struc_el => Rod(els, np_types, nip, Line(nod, nodof)),
  :properties => [EA;],
  # Compute x_coords using length l and number of elements, els
  :x_coords => 0.0:l/els:l,
  # Define a support for node N
  # Fix the single dof (x direction displacement) of node N
  :support => [(N, [0])],
  # External forces are applied to nodes 1 to 5.
  :loaded_nodes => [
      (1, [-0.625]),
      (2, [-1.25]),
      (3, [-1.25]),
      (4, [-1.25]),
      (5, [-0.625])
    ]
);

# Display the data dictionary
data |> display
println()

# Solve the FEM model
@time fem, dis_df, fm_df = p41(data)
println()

display(fem.displacements)
println()
display(fem.actions)
println()

display(dis_df)
println()
display(fm_df)
println()

using Plots
gr(size=(400,500))

x = 0.0:l/els:l
u = convert(Array, dis_df[:, :x_translation])

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
titles = ["PtFEM Ex41.1 u(x)", "PtFEM Ex41.1 N(x)"]

p[1]=plot(ylim=(0.0, 1.0), xlim=(0.0, 5.0),
  yflip=true, xflip=false, xlab="Normal force [N]",
  ylab="x [m]", title=titles[2]
)
vals = convert(Array, fm_df[:, :normal_force_2])
for i in 1:els
    plot!(p[1],
      [vals[i], vals[i]],
      [(i-1)*l/els, i*l/els], color=:blue,
      fill=true, fillalpha=0.1, leg=false
    )
    delta = abs(((i-1)*l/els) - (i*l/els)) / 20.0
    y1 = collect(((i-1)*l/els):delta:(i*l/els))
    for j in 1:length(y1)
      plot!(p[1],
        [vals[i], 0.0],
        [y1[j], y1[j]], color=:blue, alpha=0.5
      )
    end
end

p[2] = plot(u, x, xlim=(-0.00003, 0.0), yflip=true,
  xlab="Displacement [m]", ylab="x [m]",
  fillto=0.0, fillalpha=0.1, leg=false, title=titles[1])

plot(p..., layout=(1, 2))
savefig(ProjDir*"/Ex41.1.png")


text:
# functions related to Poisson distribution

# R implementations
# For pdf and logpdf we use the Julia implementation
using .RFunctions:
    poiscdf,
    poisccdf,
    poislogcdf,
    poislogccdf,
    poisinvcdf,
    poisinvccdf,
    poisinvlogcdf,
    poisinvlogccdf

# Julia implementations
poispdf(λ::Real, x::Real) = exp(poislogpdf(λ, x))

poislogpdf(λ::Real, x::Real) = poislogpdf(promote(λ, x)...)
function poislogpdf(λ::T, x::T) where {T <: Real}
    val = xlogy(x, λ) - λ - loggamma(x + 1)
    return x >= 0 && isinteger(x) ? val : oftype(val, -Inf)
end


text:
using Preferences: load_preference

@testset "Dispatcher" begin
    if load_preference(Vk, "USE_DISPATCH_TABLE", "true") == "true"
        disp = Vk.global_dispatcher[]
        @test haskey(disp.instance_tables, instance)
        @test haskey(disp.instance_tables[instance].pointers, :vkCreateDevice)
        @test disp.instance_tables[instance].pointers[:vkCreateDevice] ≠ :C_NULL
        @test haskey(disp.device_tables, device)
        @test haskey(disp.device_tables[device].pointers, :vkCreateGraphicsPipelines)
        @test disp.device_tables[device].pointers[:vkCreateGraphicsPipelines] ≠ C_NULL
        @test_throws ErrorException("Could not retrieve function pointer for 'vkCreateSwapchainKHR'. This can be caused by an extension not being enabled for a function that needs it; see the help with `?` or the documentation for more information.") function_pointer(disp, device, :vkCreateSwapchainKHR)
    else
        @test !haskey(disp.instance_tables, instance)
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
module TupleTests

using BenchmarkTools

@inline to_tuple1(x::AbstractVector) = Tuple(x)
@inline to_tuple2(x::AbstractVector) = (x...,)
@inline function to_tuple3(x::AbstractVector)
    n = length(x)
    @assert eachindex(x) == 1:n
    ntuple(i->x[i], n)
end

# The cutoff at n = 10 below reflects what is used by `ntuple`.  This value is
# somewhat arbitrary, on the machines where I tested the code, the explicit
# unrolled expression for n = 10 is still about 44 times faster than `(x...,)`.
# Calling `ntuple` for n ≤ 10 is about twice slower; for n > 10, `ntuple` is
# slower than `(x...,)`.

function to_tuple4(x::AbstractVector)
    n = length(x)
    @inbounds begin
        n == 0 ? () :
        n > 10 || firstindex(x) != 1 ? (x...,) :
        n == 1 ? (x[1],) :
        n == 2 ? (x[1], x[2]) :
        n == 3 ? (x[1], x[2], x[3]) :
        n == 4 ? (x[1], x[2], x[3], x[4]) :
        n == 5 ? (x[1], x[2], x[3], x[4], x[5]) :
        n == 6 ? (x[1], x[2], x[3], x[4], x[5], x[6]) :
        n == 7 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7]) :
        n == 8 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8]) :
        n == 9 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9]) :
        (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
    end
end

# As to_tuple4 but inlined.
@inline function to_tuple5(x::AbstractVector)
    n = length(x)
    @inbounds begin
        n == 0 ? () :
        n > 10 || firstindex(x) != 1 ? (x...,) :
        n == 1 ? (x[1],) :
        n == 2 ? (x[1], x[2]) :
        n == 3 ? (x[1], x[2], x[3]) :
        n == 4 ? (x[1], x[2], x[3], x[4]) :
        n == 5 ? (x[1], x[2], x[3], x[4], x[5]) :
        n == 6 ? (x[1], x[2], x[3], x[4], x[5], x[6]) :
        n == 7 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7]) :
        n == 8 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8]) :
        n == 9 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9]) :
        (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
    end
end

for n in ((0:20)...,30,40,50)
    println("\nmaking tuples from a vector of $n elements:")
    let x = rand(n)
        @assert to_tuple1(x) === Tuple(x)
        @assert to_tuple2(x) === Tuple(x)
        @assert to_tuple3(x) === Tuple(x)
        @assert to_tuple4(x) === Tuple(x)
        @assert to_tuple5(x) === Tuple(x)
        print("                    Tuple(x): ")
        @btime to_tuple1($x)
        print("                    (x, ...): ")
        @btime to_tuple2($x)
        print("  ntuple(i->x[i], length(x)): ")
        @btime to_tuple3($x)
        print("                to_tuple4(x): ")
        @btime to_tuple4($x)
        print("                to_tuple5(x): ")
        @btime to_tuple5($x)
    end
end

end # module

nothing


text:
#######################################################################
# An example of a simple Excel chart with XlsxWriter.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart.xlsx")

	ws = add_worksheet!(wb)

	chart = add_chart!(wb, Dict("type"=>"column"))

	data = [
	[1 2 3 4 5];
	[2 4 6 8 10];
	[3 6 9 12 15];
	]

	write_column!(ws, "A1", data[1,:])
	write_column!(ws, "B1", data[2,:])
	write_column!(ws, "C1", data[3,:])

	add_series!(chart, Dict("values"=>"=Sheet1!\$A\$1:\$A\$5"))
	add_series!(chart, Dict("values"=>"=Sheet1!\$B\$1:\$B\$5"))
	add_series!(chart, Dict("values"=>"=Sheet1!\$C\$1:\$C\$5"))

	insert_chart!(ws, "A7", chart)

	close(wb)
	isfile("chart.xlsx")
end

test()