text:
"""
Clamped cylinder, transient vibration.

The cylinder is banged at the initial time (i.e. it is given an initial velocity).
"""
module clamp_cyl_forc_damp_expl_examples

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.AssemblyModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using SymRCM
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
using Gnuplot; #@gp "clear"
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection
using ThreadedSparseCSR
using UnicodePlots
using InteractiveUtils
using BenchmarkTools

const E = 200e9;
const nu = 0.3;
const rho = 7800.0
const R = 0.1
const L = 0.8
const ksi = 0.0
const omegad = 2*pi*1000
const color = "red"
const omegaf = 2*pi*10^5
const tend = 25 * (2*pi/omegaf)
const visualize = true

function loop!(M, K, ksi, U0, V0, tend, dt, force!, peek)
    U = deepcopy(U0)
    V = deepcopy(U0)
    A = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invMC .* force!(F, t);
    peek(0, U, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U += dt*V + ((dt^2)/2)*A;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-=  mul!(E, K, U)
        @inbounds @simd for i in eachindex(U)
            _Fi = F[i]; _Ai = A[i]; _Vi = V[i]
            # Add damping forces
            _Fi -= C[i] * (_Vi + (dt/2) * _Ai)
            # Compute the new acceleration.
            _A1i = invMC[i] * _Fi
            A[i] = _A1i
            # Update the velocity
            V[i] += (dt/2)* (_Ai + _A1i);
        end
        t = t + dt
        peek(step, U, t)
    end
end

struct ThreadBuffer{KVT}
    kcolumns::KVT
    uv::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}
    result::Vector{Float64}
end

struct ThreadBufferVOM{RT}
    rang::RT
    result::Vector{Float64}
end

# function parmul!(R, tbs)
#     tasks = [];
#     for th in 1:length(tbs)
# tb = tbs[th]
# push!(tasks, Threads.@spawn begin
#     mul!(tb.result, tb.kcolumns, tb.uv)
# end);
#     end
#     Threads.wait(tasks[1]);
#     R .= tbs[1].result
#     for th in 2:length(tbs)
# Threads.wait(tasks[th]);
# R .+= tbs[th].result
#     end
#     R
# end

function parmul!(R, tbs)
    Threads.@threads for i in 1:length(tbs)
        if i == 1
            mul!(R, tbs[i].kcolumns, tbs[i].uv)
        else
            mul!(tbs[i].result, tbs[i].kcolumns, tbs[i].uv)
        end
    end
    for i in 2:length(tbs)
        R .+= tbs[i].result
    end
    R
end

function parloop!(M, K, ksi, U0, V0, tend, dt, force!, peek, nthr)
    U = deepcopy(U0)
    V = deepcopy(U0)
    A = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"
    chunk = Int(floor(length(U) / nth))
    threadbuffs = ThreadBuffer[];
    for th in 1:nth
        colrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:length(U))
        push!(threadbuffs, ThreadBuffer(K[:, colrange], view(U, colrange), deepcopy(U)));
    end

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invMC .* force!(F, t);
    peek(0, U, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U += dt*V + ((dt^2)/2)*A;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-=  parmul!(E, threadbuffs)
        @inbounds @simd for i in eachindex(U)
            _Fi = F[i]; _Ai = A[i]; _Vi = V[i]
            # Add damping forces
            _Fi -= C[i] * (_Vi + (dt/2) * _Ai)
            # Compute the new acceleration.
            _A1i = invMC[i] * _Fi
            A[i] = _A1i
            # Update the velocity
            V[i] += (dt/2)* (_Ai + _A1i);
        end
        t = t + dt
        peek(step, U, t)
    end
end

function parloop_csr!(M, K, ksi, U0, V0, tend, dt, force!, peek, nthr)
    U = deepcopy(U0)
    V = deepcopy(U0)
    A = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invMC .* force!(F, t);
    peek(0, U, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U += dt*V + ((dt^2)/2)*A;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-= ThreadedSparseCSR.bmul!(E, K, U)
        @inbounds @simd for i in eachindex(U)
            _Fi = F[i]; _Ai = A[i]; _Vi = V[i]
            # Add damping forces
            _Fi -= C[i] * (_Vi + (dt/2) * _Ai)
            # Compute the new acceleration.
            _A1i = invMC[i] * _Fi
            A[i] = _A1i
            # Update the velocity
            V[i] += (dt/2)* (_Ai + _A1i);
        end
        t = t + dt
        peek(step, U, t)
    end
end

function parvommul!(R, vom, P, threadbuffs, U)
    Threads.@threads for i in 1:length(threadbuffs)
        if i == 1
            vommul!(R, vom, P, threadbuffs[i].rang, U)
        else
            vommul!(threadbuffs[i].result, vom, P, threadbuffs[i].rang, U)
        end
    end
    for i in 2:length(threadbuffs)
        R .+= threadbuffs[i].result
    end
    R
end

function vommul!(R, vom, P, rang, U)
    z = zero(eltype(R))
    R .= z
    # v = Kasm.buffer
    # P = Kasm.elem_mat_dim
    # r = fill(z, P)
    # u = fill(z, P)
    for i in rang
        d = vom[i][1]
        m = vom[i][2]
        @inbounds for i in 1:P
            acc = z
            @inbounds for j in 1:P
                acc += m[i, j] * U[d[j]]
            end
            R[d[i]] += acc
        end
    end
    R
end

function parloop_vom!(M, Kasm, threadbuffs, ksi, U0, V0, tend, dt, force!, peek, nthr)
    U = deepcopy(U0)
    V = deepcopy(U0)
    A = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"

    vom = Kasm.buffer
    P = Kasm.elem_mat_dim

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invMC .* force!(F, t);
    peek(0, U, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U += dt*V + ((dt^2)/2)*A;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-= parvommul!(E, vom, P, threadbuffs, U)
        @inbounds @simd for i in eachindex(U)
            _Fi = F[i]; _Ai = A[i]; _Vi = V[i]
            # Add damping forces
            _Fi -= C[i] * (_Vi + (dt/2) * _Ai)
            # Compute the new acceleration.
            _A1i = invMC[i] * _Fi
            A[i] = _A1i
            # Update the velocity
            V[i] += (dt/2)* (_Ai + _A1i);
        end
        t = t + dt
        peek(step, U, t)
    end
end

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ); r[2] = 0.0;
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function _execute_parallel(n = 64, thickness = 0.01, nthr = 0)
    color = "red"
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    K = FEMMShellT3FFModule.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    # Solve
    function pwr(K, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            mul!(w, K, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        sqrt((v' * K * v) / (v' * M * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[0 -L/4 R])[1]
    cpoint = selectnode(fens; nearestto=[0 0 R])[1]

    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]
    cpointdof6 = dchi.dofnums[cpoint, 6]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[qpoint, 1]
        dy = XYZ[2] - fens.xyz[qpoint, 2]
        dz = XYZ[3] - fens.xyz[qpoint, 3]
        mag = 10^5*exp(-20/R*sqrt(dx^2+dy^2+dz^2))
        r = vec(XYZ); r[2] = 0.0
        r .= vec(r)/norm(vec(r))
        theta = atan(r[3], r[1])
        forceout[1:3] .= r*mag
        forceout[4:6] .= 0.0
        return forceout
    end

    # Sinusoidal loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        if t < 2 * (2*pi/omegaf)
            F .= sin(omegaf*t) .* Fmag
        else
            F .= 0.0
        end
        # @show norm(F)
        return F
    end

    # Suddenly applied constant force at a node
    # Fmag = fill(0.0, dchi.nfreedofs)
    # Fmag[qpointdof] = -1.0
    # function force!(F, t)
    #     F .= Fmag
    #     return F
    # end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        cdeflections[step+1] = U[cpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    parloop!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek, nthr)

    if visualize
        # @gp  "set terminal windows 0 "  :-
        # @gp "clear"
        @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-

        @gp  :- "set xlabel 'Time'" :-
        @gp  :- "set ylabel 'Deflection'" :-
        @gp  :- "set title 'Free-floating plate'"

        # Visualization
        @info "Dumping visualization"
        times = Float64[]
        vectors = []
        for i in 1:length(displacements)
            scattersysvec!(dchi, displacements[i])
            push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
            push!(times, i*dt*nbtw)
        end
        vtkwritecollection("clamp_cyl_forc_damp_expl_$n", fens, fes, times; vectors = vectors)
    end
end

function _execute_parallel_csr(n = 64, thickness = 0.01, nthr = 0)
    color = "green"
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    K = FEMMShellT3FFModule.stiffness(femm, SM.SysmatAssemblerSparseCSRSymm(0.0), geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    # Solve
    function pwr(K, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            ThreadedSparseCSR.bmul!(w, K, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        sqrt((v' * (K * v)) / (v' * M * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[0 -L/4 R])[1]
    cpoint = selectnode(fens; nearestto=[0 0 R])[1]

    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]
    cpointdof6 = dchi.dofnums[cpoint, 6]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[qpoint, 1]
        dy = XYZ[2] - fens.xyz[qpoint, 2]
        dz = XYZ[3] - fens.xyz[qpoint, 3]
        mag = 10^5*exp(-20/R*sqrt(dx^2+dy^2+dz^2))
        r = vec(XYZ); r[2] = 0.0
        r .= vec(r)/norm(vec(r))
        theta = atan(r[3], r[1])
        forceout[1:3] .= r*mag
        forceout[4:6] .= 0.0
        return forceout
    end

    # Sinusoidal loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        if t < 2 * (2*pi/omegaf)
            F .= sin(omegaf*t) .* Fmag
        else
            F .= 0.0
        end
        # @show norm(F)
        return F
    end

    # Suddenly applied constant force at a node
    # Fmag = fill(0.0, dchi.nfreedofs)
    # Fmag[qpointdof] = -1.0
    # function force!(F, t)
    #     F .= Fmag
    #     return F
    # end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        cdeflections[step+1] = U[cpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    parloop_csr!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek, nthr)

    if visualize
        # @gp  "set terminal windows 0 "  :-
        # @gp "clear"
        @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-

        @gp  :- "set xlabel 'Time'" :-
        @gp  :- "set ylabel 'Deflection'" :-
        @gp  :- "set title 'Free-floating plate'"

        # Visualization
        @info "Dumping visualization"
        times = Float64[]
        vectors = []
        for i in 1:length(displacements)
            scattersysvec!(dchi, displacements[i])
            push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
            push!(times, i*dt*nbtw)
        end
        vtkwritecollection("clamp_cyl_forc_damp_expl_$n", fens, fes, times; vectors = vectors)
    end
end

function _execute_serial(n = 64, thickness = 0.01)
    color = "blue"

    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    K = FEMMShellT3FFModule.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    # Solve
    function pwr(K, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            mul!(w, K, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        sqrt((v' * K * v) / (v' * M * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[0 -L/4 R])[1]
    cpoint = selectnode(fens; nearestto=[0 0 R])[1]

    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]
    cpointdof6 = dchi.dofnums[cpoint, 6]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[qpoint, 1]
        dy = XYZ[2] - fens.xyz[qpoint, 2]
        dz = XYZ[3] - fens.xyz[qpoint, 3]
        mag = 10^5*exp(-20/R*sqrt(dx^2+dy^2+dz^2))
        r = vec(XYZ); r[2] = 0.0
        r .= vec(r)/norm(vec(r))
        theta = atan(r[3], r[1])
        forceout[1:3] .= r*mag
        forceout[4:6] .= 0.0
        return forceout
    end

    # Sinusoidal loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        if t < 2 * (2*pi/omegaf)
            F .= sin(omegaf*t) .* Fmag
        else
            F .= 0.0
        end
        # @show norm(F)
        return F
    end

    # Suddenly applied constant force at a node
    # Fmag = fill(0.0, dchi.nfreedofs)
    # Fmag[qpointdof] = -1.0
    # function force!(F, t)
    #     F .= Fmag
    #     return F
    # end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        cdeflections[step+1] = U[cpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    loop!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek)

    if visualize
        # @gp  "set terminal windows 0 "  :-
        # @gp "clear"
        @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-

        @gp  :- "set xlabel 'Time'" :-
        @gp  :- "set ylabel 'Deflection'" :-
        @gp  :- "set title 'Free-floating plate'"

        # Visualization
        @info "Dumping visualization"
        times = Float64[]
        vectors = []
        for i in 1:length(displacements)
            scattersysvec!(dchi, displacements[i])
            push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
            push!(times, i*dt*nbtw)
        end
        vtkwritecollection("clamp_cyl_forc_damp_expl_$n", fens, fes, times; vectors = vectors)
    end
end

function _execute_parallel_vom(n = 64, thickness = 0.01, nthr = 0)
    color = "magenta"
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    Kasm = SM.SysmatAssemblerVecOfMatSymm(0.0)
    FEMMShellT3FFModule.stiffness(femm, Kasm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"
    N = Kasm.buffer_length
    chunk = Int(floor(N / nth))
    threadbuffs = ThreadBufferVOM[];
    for th in 1:nth
        rang = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:N)
        push!(threadbuffs, ThreadBufferVOM(rang, fill(0.0, dchi.nfreedofs)));
    end

    # Solve
    function pwr(vom, P, threadbuffs, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            parvommul!(w, vom, P, threadbuffs, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        parvommul!(w, vom, P, threadbuffs, v)
        sqrt((v' * w) / (v' * M * v))
    end
    vom = Kasm.buffer
    P = Kasm.elem_mat_dim
    @time omega_max = pwr(vom, P, threadbuffs, M)
    @show omega_max

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[0 -L/4 R])[1]
    cpoint = selectnode(fens; nearestto=[0 0 R])[1]

    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]
    cpointdof6 = dchi.dofnums[cpoint, 6]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[qpoint, 1]
        dy = XYZ[2] - fens.xyz[qpoint, 2]
        dz = XYZ[3] - fens.xyz[qpoint, 3]
        mag = 10^5*exp(-20/R*sqrt(dx^2+dy^2+dz^2))
        r = vec(XYZ); r[2] = 0.0
        r .= vec(r)/norm(vec(r))
        theta = atan(r[3], r[1])
        forceout[1:3] .= r*mag
        forceout[4:6] .= 0.0
        return forceout
    end

    # Sinusoidal loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        if t < 2 * (2*pi/omegaf)
            F .= sin(omegaf*t) .* Fmag
        else
            F .= 0.0
        end
        # @show norm(F)
        return F
    end

    # Suddenly applied constant force at a node
    # Fmag = fill(0.0, dchi.nfreedofs)
    # Fmag[qpointdof] = -1.0
    # function force!(F, t)
    #     F .= Fmag
    #     return F
    # end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        cdeflections[step+1] = U[cpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    parloop_vom!(M, Kasm, threadbuffs, ksi, U0, V0, nsteps*dt, dt, force!, peek, nthr)

    if visualize
        # @gp  "set terminal windows 0 "  :-
        # @gp "clear"
        @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-

        @gp  :- "set xlabel 'Time'" :-
        @gp  :- "set ylabel 'Deflection'" :-
        @gp  :- "set title 'Free-floating plate'"

        # Visualization
        @info "Dumping visualization"
        times = Float64[]
        vectors = []
        for i in 1:length(displacements)
            scattersysvec!(dchi, displacements[i])
            push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
            push!(times, i*dt*nbtw)
        end
        vtkwritecollection("clamp_cyl_forc_damp_expl_$n", fens, fes, times; vectors = vectors)
    end
end

function test_serial(ns = [4*64])
    @info "Clamped cylinder: serial"
    for n in ns
        _execute_serial(n, 0.01)
    end
    return true
end

function test_parallel(ns = [4*64], nthr = 0)
    @info "Clamped cylinder: parallel"
    for n in ns
        _execute_parallel(n, 0.01, nthr)
    end
    return true
end

function test_parallel_csr(ns = [4*64], nthr = 0)
    @info "Clamped cylinder: parallel CSR"
    for n in ns
        _execute_parallel_csr(n, 0.01, nthr)
    end
    return true
end

function test_parallel_vom(ns = [4*64], nthr = 0)
    @info "Clamped cylinder: parallel VOM"
    for n in ns
        _execute_parallel_vom(n, 0.01, nthr)
    end
    return true
end

function _explore_csr(n = 64, thickness = 0.01, nthr = 0)
    color = "green"
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);
    # @show dchi.dofnums

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    K = FEMMShellT3FFModule.stiffness(femm, SM.SysmatAssemblerSparseCSRSymm(0.0), geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    I, J, V = findnz(K)
    starts = fill(Inf, size(K, 1))
    finishes = fill(0, size(K, 1))
    nzs = fill(0, size(K, 1))

    for i in 1:length(I)
        nzs[I[i]] += 1
        starts[I[i]] = min(starts[I[i]], J[i])
        finishes[I[i]] = max(finishes[I[i]], J[i])
    end
    bw = finishes .- starts .+ 1
    display(lineplot(1:length(bw), bw))
        @show median(bw)
        @show minimum(bw)
        @show maximum(bw)
        @show median(nzs)
        @show minimum(nzs)
        @show maximum(nzs)
    display(spy(K))

    nothing
end

function allrun(ns = [4*64])
    println("#####################################################")
    println("# test_serial ")
    test_serial(ns)
    println("#####################################################")
    println("# test_parallel ")
    test_parallel(ns)
    println("#####################################################")
    println("# test_parallel_csr ")
    test_parallel_csr(ns)
    println("#####################################################")
    println("# test_parallel_vom ")
    test_parallel_vom(ns)
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing

# #####################################################
# # test_serial
# [ Info: Clamped cylinder: serial
# [ Info: Mesh 262400 nodes, 524288 elements
#   2.341943 seconds (40 allocations: 97.752 MiB)
# omega_max = 1.41323e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27368e-07
# [ Info: 1963 steps
# 146.272046 seconds (4.32 k allocations: 1.147 GiB, 0.05% gc time)
# [ Info: Dumping visualization
# #####################################################
# # test_parallel
# [ Info: Clamped cylinder: parallel
# [ Info: Mesh 262400 nodes, 524288 elements
#   2.376553 seconds (40 allocations: 97.752 MiB)
# omega_max = 1.41333e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27359e-07
# [ Info: 1963 steps
# [ Info: 2 threads used
#  84.598355 seconds (176.37 k allocations: 1.158 GiB, 0.16% gc time, 0.10% compilation time)
# [ Info: Dumping visualization
# #####################################################
# # test_parallel_csr
# [ Info: Clamped cylinder: parallel CSR
# [ Info: Mesh 262400 nodes, 524288 elements
#   2.114689 seconds (118.45 k allocations: 103.899 MiB, 29.08% compilation time)
# omega_max = 1.41333e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27359e-07
# [ Info: 1963 steps
# [ Info: 2 threads used
#  82.601808 seconds (20.03 k allocations: 1.148 GiB, 0.25% gc time)
# [ Info: Dumping visualization
# #####################################################
# # test_parallel_vom
# [ Info: Clamped cylinder: parallel VOM
# [ Info: Mesh 262400 nodes, 524288 elements
# [ Info: 2 threads used
#   4.372176 seconds (172.21 k allocations: 95.321 MiB, 1.67% compilation time)
# omega_max = 1.41354e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27340e-07
# [ Info: 1963 steps
# [ Info: 2 threads used
# 268.287995 seconds (30.07 k allocations: 1.150 GiB, 0.04% gc time)
# [ Info: Dumping visualization
# true

# #####################################################
# # test_serial
# [ Info: Clamped cylinder: serial
# [ Info: Mesh 262400 nodes, 524288 elements
#   2.360393 seconds (40 allocations: 97.752 MiB)
# omega_max = 1.41322e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27369e-07
# [ Info: 1963 steps
# 147.455232 seconds (4.32 k allocations: 1.147 GiB, 0.06% gc time)
# [ Info: Dumping visualization
# #####################################################
# # test_parallel
# [ Info: Clamped cylinder: parallel
# [ Info: Mesh 262400 nodes, 524288 elements
#   2.412289 seconds (40 allocations: 97.752 MiB)
# omega_max = 1.41318e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27373e-07
# [ Info: 1963 steps
# [ Info: 4 threads used
#  80.991000 seconds (205.23 k allocations: 1.160 GiB, 0.18% gc time, 0.13% compilation time)
# [ Info: Dumping visualization
# #####################################################
# # test_parallel_csr
# [ Info: Clamped cylinder: parallel CSR
# [ Info: Mesh 262400 nodes, 524288 elements
#   1.988305 seconds (158.39 k allocations: 106.171 MiB, 34.71% compilation time)
# omega_max = 1.41331e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27361e-07
# [ Info: 1963 steps
# [ Info: 4 threads used
#  72.222234 seconds (20.03 k allocations: 1.148 GiB, 0.35% gc time)
# [ Info: Dumping visualization
# #####################################################
# # test_parallel_vom
# [ Info: Clamped cylinder: parallel VOM
# [ Info: Mesh 262400 nodes, 524288 elements
# [ Info: 4 threads used
#   3.393447 seconds (172.64 k allocations: 95.362 MiB, 2.69% compilation time)
# omega_max = 1.41302e+07
# dt = Float64((9.00000e-01 * 2) / omega_max) * (sqrt(1 + ksi ^ 2) - ksi) = 1.27387e-07
# [ Info: 1963 steps
# [ Info: 4 threads used
# 220.470659 seconds (54.74 k allocations: 1.152 GiB, 0.05% gc time)
# [ Info: Dumping visualization
# true


text:
# Exceptions
@testset "Exceptions" begin
    data = dataset("Ecdat", "Crime")
    data.AvgSen2 = 2 * data.AvgSen
    # Rank-deficient model
    @test isa(
        fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + AvgSen2),
            data,
        ),
        EconometricModel,
    )
    @test isa(
        fit(
            RandomEffectsEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
            data,
            panel = :County,
            time = :Year,
        ),
        EconometricModel,
    )
    # Model specifics
    @test_throws(
        ArgumentError("The between estimator requires the panel identifier."),
        fit(
            BetweenEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
            data,
        )
    )
    @test_throws(
        ArgumentError("Absorbing features is only implemented for least squares."),
        fit(
            BetweenEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
            data,
            panel = :County,
        )
    )
    @test_throws(
        ArgumentError("The random effects estimator requires the panel and temporal identifiers."),
        fit(
            RandomEffectsEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
            data,
            panel = :County,
        )
    )
    @test_throws(
        ArgumentError("Absorbing features is only implemented for least squares."),
        fit(
            RandomEffectsEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
            data,
            panel = :County,
            time = :Year,
        )
    )
    data = dataset("datasets", "iris")
    @test_throws(
        ConvergenceException,
        fit(EconometricModel, @formula(Species ~ SepalLength + SepalWidth), data)
    )
    # Robust variance covariance estimators
    @test_throws(
        ArgumentError("Robust variance covariance estimators only available for continous response models."),
        fit(
            EconometricModel,
            @formula(Species ~ SepalLength + SepalWidth),
            data,
            vce = HC1,
        )
    )
    @test_throws(
        ArgumentError("When absorbing features, only HC2-HC4 are unavailable."),
        fit(
            EconometricModel,
            @formula(SepalLength ~ SepalLength + absorb(Species)),
            data,
            vce = HC2,
        )
    )
end

ext: jl


text:
using sbl
using Test

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


text:
module Compat

import LinearAlgebra
using LinearAlgebra: Adjoint, Diagonal, Transpose, UniformScaling, RealHermSymComplexHerm

include("compatmacro.jl")

# https://github.com/JuliaLang/julia/pull/29440
if VERSION < v"1.1.0-DEV.389"
    Base.:(:)(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
        CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end

# https://github.com/JuliaLang/julia/pull/29442
if VERSION < v"1.1.0-DEV.403"
    Base.oneunit(::CartesianIndex{N}) where {N} = oneunit(CartesianIndex{N})
    Base.oneunit(::Type{CartesianIndex{N}}) where {N} = CartesianIndex(ntuple(x -> 1, Val(N)))
end

# https://github.com/JuliaLang/julia/pull/29679
if VERSION < v"1.1.0-DEV.472"
    export isnothing
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

# https://github.com/JuliaLang/julia/pull/29749
if VERSION < v"1.1.0-DEV.792"
    export eachrow, eachcol, eachslice
    eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1))
    eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))
    @inline function eachslice(A::AbstractArray; dims)
        length(dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
        dim = first(dims)
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
        return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
    end
end

function rangeargcheck(;step=nothing, length=nothing, kwargs...)
    if step===nothing && length===nothing
        throw(ArgumentError("At least one of `length` or `step` must be specified"))
    end
end

if VERSION < v"1.1.0-DEV.506"
    function Base.range(start, stop; kwargs...)
        rangeargcheck(;kwargs...)
        range(start; stop=stop, kwargs...)
    end
end

# https://github.com/JuliaLang/julia/pull/30496
if VERSION < v"1.2.0-DEV.272"
    Base.@pure hasfield(::Type{T}, name::Symbol) where T =
        Base.fieldindex(T, name, false) > 0
    export hasfield
    hasproperty(x, s::Symbol) = s in propertynames(x)
    export hasproperty
end

# https://github.com/JuliaLang/julia/pull/29259
if VERSION < v"1.1.0-DEV.594"
    Base.merge(a::NamedTuple, b::NamedTuple, cs::NamedTuple...) = merge(merge(a, b), cs...)
    Base.merge(a::NamedTuple) = a
end

# https://github.com/JuliaLang/julia/pull/33129
if VERSION < v"1.4.0-DEV.142"
    export only

    Base.@propagate_inbounds function only(x)
        i = iterate(x)
        @boundscheck if i === nothing
            throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
        end
        (ret, state) = i
        @boundscheck if iterate(x, state) !== nothing
            throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
        end
        return ret
    end

    # Collections of known size
    only(x::Ref) = x[]
    only(x::Number) = x
    only(x::Char) = x
    only(x::Tuple{Any}) = x[1]
    only(x::Tuple) = throw(
        ArgumentError("Tuple contains $(length(x)) elements, must contain exactly 1 element")
    )
    only(a::AbstractArray{<:Any, 0}) = @inbounds return a[]
    only(x::NamedTuple{<:Any, <:Tuple{Any}}) = first(x)
    only(x::NamedTuple) = throw(
        ArgumentError("NamedTuple contains $(length(x)) elements, must contain exactly 1 element")
    )
end

# https://github.com/JuliaLang/julia/pull/32628
if VERSION < v"1.3.0-alpha.8"
    Base.mod(i::Integer, r::Base.OneTo) = mod1(i, last(r))
    Base.mod(i::Integer, r::AbstractUnitRange{<:Integer}) = mod(i-first(r), length(r)) + first(r)
end

# https://github.com/JuliaLang/julia/pull/32739
# This omits special methods for more exotic matrix types, Triangular and worse.
if VERSION < v"1.4.0-DEV.92" # 2425ae760fb5151c5c7dd0554e87c5fc9e24de73

    # stdlib/LinearAlgebra/src/generic.jl
    LinearAlgebra.dot(x, A, y) = LinearAlgebra.dot(x, A*y) # generic fallback

    function LinearAlgebra.dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
        (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
        T = typeof(LinearAlgebra.dot(first(x), first(A), first(y)))
        s = zero(T)
        i₁ = first(eachindex(x))
        x₁ = first(x)
        @inbounds for j in eachindex(y)
            yj = y[j]
            if !iszero(yj)
                temp = zero(adjoint(A[i₁,j]) * x₁)
                @simd for i in eachindex(x)
                    temp += adjoint(A[i,j]) * x[i]
                end
                s += LinearAlgebra.dot(temp, yj)
            end
        end
        return s
    end
    LinearAlgebra.dot(x::AbstractVector, adjA::Adjoint, y::AbstractVector) =
         adjoint(LinearAlgebra.dot(y, adjA.parent, x))
    LinearAlgebra.dot(x::AbstractVector, transA::Transpose{<:Real}, y::AbstractVector) =
        adjoint(LinearAlgebra.dot(y, transA.parent, x))

    # stdlib/LinearAlgebra/src/diagonal.jl
    function LinearAlgebra.dot(x::AbstractVector, D::Diagonal, y::AbstractVector)
        mapreduce(t -> LinearAlgebra.dot(t[1], t[2], t[3]), +, zip(x, D.diag, y))
    end

    # stdlib/LinearAlgebra/src/symmetric.jl
    function LinearAlgebra.dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
        require_one_based_indexing(x, y)
        (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
        data = A.data
        r = zero(eltype(x)) * zero(eltype(A)) * zero(eltype(y))
        if A.uplo == 'U'
            @inbounds for j = 1:length(y)
                r += LinearAlgebra.dot(x[j], real(data[j,j]), y[j])
                @simd for i = 1:j-1
                    Aij = data[i,j]
                    r += LinearAlgebra.dot(x[i], Aij, y[j]) +
                        LinearAlgebra.dot(x[j], adjoint(Aij), y[i])
                end
            end
        else # A.uplo == 'L'
            @inbounds for j = 1:length(y)
                r += LinearAlgebra.dot(x[j], real(data[j,j]), y[j])
                @simd for i = j+1:length(y)
                    Aij = data[i,j]
                    r += LinearAlgebra.dot(x[i], Aij, y[j]) +
                        LinearAlgebra.dot(x[j], adjoint(Aij), y[i])
                end
            end
        end
        return r
    end

    # stdlib/LinearAlgebra/src/uniformscaling.jl
    LinearAlgebra.dot(x::AbstractVector, J::UniformScaling, y::AbstractVector) =
        LinearAlgebra.dot(x, J.λ, y)
    LinearAlgebra.dot(x::AbstractVector, a::Number, y::AbstractVector) =
        sum(t -> LinearAlgebra.dot(t[1], a, t[2]), zip(x, y))
    LinearAlgebra.dot(x::AbstractVector, a::Union{Real,Complex}, y::AbstractVector) =
        a*LinearAlgebra.dot(x, y)
end

# https://github.com/JuliaLang/julia/pull/30630
if VERSION < v"1.2.0-DEV.125" # 1da48c2e4028c1514ed45688be727efbef1db884
    require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError(
        "offset arrays are not supported but got an array with index other than 1"))
# At present this is only used in Compat inside the above dot(x,A,y) functions, #32739
elseif VERSION < v"1.4.0-DEV.92"
    using Base: require_one_based_indexing
end

# https://github.com/JuliaLang/julia/pull/33568
if VERSION < v"1.4.0-DEV.329"
    Base.:∘(f, g, h...) = ∘(f ∘ g, h...)
end

# https://github.com/JuliaLang/julia/pull/33128
if VERSION < v"1.4.0-DEV.397"
    export pkgdir
    function pkgdir(m::Module)
        rootmodule = Base.moduleroot(m)
        path = pathof(rootmodule)
        path === nothing && return nothing
        return dirname(dirname(path))
    end
end

# https://github.com/JuliaLang/julia/pull/33736/
if VERSION < v"1.4.0-DEV.493"
    Base.Order.ReverseOrdering() = Base.Order.ReverseOrdering(Base.Order.Forward)
end

# https://github.com/JuliaLang/julia/pull/32968
if VERSION < v"1.4.0-DEV.551"
    Base.filter(f, xs::Tuple) = Base.afoldl((ys, x) -> f(x) ? (ys..., x) : ys, (), xs...)
    Base.filter(f, t::Base.Any16) = Tuple(filter(f, collect(t)))
end

# https://github.com/JuliaLang/julia/pull/34652
if VERSION < v"1.5.0-DEV.247"
    export ismutable
    ismutable(@nospecialize(x)) = (Base.@_pure_meta; typeof(x).mutable)
end

# https://github.com/JuliaLang/julia/pull/28761
export uuid5
if VERSION < v"1.1.0-DEV.326"
    import SHA
    import UUIDs: UUID
    function uuid5(ns::UUID, name::String)
        nsbytes = zeros(UInt8, 16)
        nsv = ns.value
        for idx in Base.OneTo(16)
            nsbytes[idx] = nsv >> 120
            nsv = nsv << 8
        end
        hash_result = SHA.sha1(append!(nsbytes, convert(Vector{UInt8}, codeunits(unescape_string(name)))))
        # set version number to 5
        hash_result[7] = (hash_result[7] & 0x0F) | (0x50)
        hash_result[9] = (hash_result[9] & 0x3F) | (0x80)
        v = zero(UInt128)
        #use only the first 16 bytes of the SHA1 hash
        for idx in Base.OneTo(16)
            v = (v << 0x08) | hash_result[idx]
        end
        return UUID(v)
    end
else
    using UUIDs: uuid5
end

# https://github.com/JuliaLang/julia/pull/34773
if VERSION < v"1.5.0-DEV.301"
    Base.zero(::AbstractIrrational) = false
    Base.zero(::Type{<:AbstractIrrational}) = false

    Base.one(::AbstractIrrational) = true
    Base.one(::Type{<:AbstractIrrational}) = true
end

# https://github.com/JuliaLang/julia/pull/32753
if VERSION < v"1.4.0-DEV.513"
    function evalpoly(x, p::Tuple)
        if @generated
            N = length(p.parameters)
            ex = :(p[end])
            for i in N-1:-1:1
                ex = :(muladd(x, $ex, p[$i]))
            end
            ex
        else
            _evalpoly(x, p)
        end
    end

    evalpoly(x, p::AbstractVector) = _evalpoly(x, p)

    function _evalpoly(x, p)
        N = length(p)
        ex = p[end]
        for i in N-1:-1:1
            ex = muladd(x, ex, p[i])
        end
        ex
    end

    function evalpoly(z::Complex, p::Tuple)
        if @generated
            N = length(p.parameters)
            a = :(p[end])
            b = :(p[end-1])
            as = []
            for i in N-2:-1:1
                ai = Symbol("a", i)
                push!(as, :($ai = $a))
                a = :(muladd(r, $ai, $b))
                b = :(muladd(-s, $ai, p[$i]))
            end
            ai = :a0
            push!(as, :($ai = $a))
            C = Expr(:block,
            :(x = real(z)),
            :(y = imag(z)),
            :(r = x + x),
            :(s = muladd(x, x, y*y)),
            as...,
            :(muladd($ai, z, $b)))
        else
            _evalpoly(z, p)
        end
    end
    evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]

    evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)

    function _evalpoly(z::Complex, p)
        length(p) == 1 && return p[1]
        N = length(p)
        a = p[end]
        b = p[end-1]

        x = real(z)
        y = imag(z)
        r = 2x
        s = muladd(x, x, y*y)
        for i in N-2:-1:1
            ai = a
            a = muladd(r, ai, b)
            b = muladd(-s, ai, p[i])
        end
        ai = a
        muladd(ai, z, b)
    end
    export evalpoly
end

# https://github.com/JuliaLang/julia/pull/35304
if VERSION < v"1.5.0-DEV.574"
    Base.similar(A::PermutedDimsArray, T::Type, dims::Base.Dims) = similar(parent(A), T, dims)
end

# https://github.com/JuliaLang/julia/pull/34548
if VERSION < v"1.5.0-DEV.314"
    macro NamedTuple(ex)
        Meta.isexpr(ex, :braces) || Meta.isexpr(ex, :block) ||
            throw(ArgumentError("@NamedTuple expects {...} or begin...end"))
        decls = filter(e -> !(e isa LineNumberNode), ex.args)
        all(e -> e isa Symbol || Meta.isexpr(e, :(::)), decls) ||
            throw(ArgumentError("@NamedTuple must contain a sequence of name or name::type expressions"))
        vars = [QuoteNode(e isa Symbol ? e : e.args[1]) for e in decls]
        types = [esc(e isa Symbol ? :Any : e.args[2]) for e in decls]
        return :(NamedTuple{($(vars...),), Tuple{$(types...)}})
    end

    export @NamedTuple
end

# https://github.com/JuliaLang/julia/pull/34296
if VERSION < v"1.5.0-DEV.182"
    export mergewith, mergewith!
    _asfunction(f::Function) = f
    _asfunction(f) = (args...) -> f(args...)
    mergewith(f, dicts...) = merge(_asfunction(f), dicts...)
    mergewith!(f, dicts...) = merge!(_asfunction(f), dicts...)
    mergewith(f) = (dicts...) -> mergewith(f, dicts...)
    mergewith!(f) = (dicts...) -> mergewith!(f, dicts...)
end

# https://github.com/JuliaLang/julia/pull/32003
if VERSION < v"1.4.0-DEV.29"
    hasfastin(::Type) = false
    hasfastin(::Union{Type{<:AbstractSet},Type{<:AbstractDict},Type{<:AbstractRange}}) = true
    hasfastin(x) = hasfastin(typeof(x))
else
    const hasfastin = Base.hasfastin
end

# https://github.com/JuliaLang/julia/pull/34427
if VERSION < v"1.5.0-DEV.124"
    const FASTIN_SET_THRESHOLD = 70

    function isdisjoint(l, r)
        function _isdisjoint(l, r)
            hasfastin(r) && return !any(in(r), l)
            hasfastin(l) && return !any(in(l), r)
            Base.haslength(r) && length(r) < FASTIN_SET_THRESHOLD &&
                return !any(in(r), l)
            return !any(in(Set(r)), l)
        end
        if Base.haslength(l) && Base.haslength(r) && length(r) < length(l)
            return _isdisjoint(r, l)
        end
        _isdisjoint(l, r)
    end

    export isdisjoint
end

include("deprecated.jl")

end # module Compat


text:
using Plots

function makeCosPlots()
    x = 0:0.01:1.00
    yn2, yn4, yn8, yn16 = [], [], [], []
    for i in x
        push!(yn2, cosx(i, 2))
    end
    for i in x
        push!(yn4, cosx(i, 4))
    end
    for i in x
        push!(yn8, cosx(i, 8))
    end
    for i in x
        push!(yn16, cosx(i, 16))
    end
    p = plot(x, yn2)
    plot!(p, x, yn4)
    plot!(p, x, yn8)
    plot!(p, x, yn16)
    display(p)
end


text:
# coding:utf-8

doc=
"""
  Dialogue Act Classification for Infomation-Seeking Dialogue System
  修正予定
  2017/02/05 ~
"""

include("switch_depwarn.jl")
switch_depwarn!(false)

using DataFrames
using DataStructures
using JSON
using PyCall
using MLBase
using ScikitLearn

using ScikitLearn: fit!, predict, transform
@sk_import ensemble: RandomForestClassifier
@sk_import ensemble: GradientBoostingClassifier
@sk_import ensemble: RandomTreesEmbedding
@sk_import metrics: classification_report
@sk_import neighbors: KNeighborsClassifier
@sk_import neural_network: MLPClassifier
@sk_import gaussian_process: GaussianProcessClassifier
@sk_import svm: SVC

# python module
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport prep
@pyimport numpy as np

# 追加素性作成
# normarized sequence number & message length
function add_feature(d, session_ids)
  for sid in session_ids
    normarized_sequence, message_length = [], []
    for (i,cont) in enumerate(d[sid][:unigram])
      push!(normarized_sequence, d[sid][:sequence][i]/length(cont))
      push!(message_length, length(cont))
    end
    d[sid][:norm_sequence] = normarized_sequence
    d[sid][:message_length] = message_length
  end
  return d
end

# 出現回数でフィルタリング
function convert_sequence(word_sequence, word_d, n_count)
  # rank_dict
  rank =  Dict(w[1] => i for (i, w) in enumerate(sort(collect(word_d), by=x->x[2], rev = true)))
  for (i, seq) in enumerate(word_sequence)
      for (j,w) in enumerate(seq)
        if word_d[w] > n_count
          word_sequence[i][j] = rank[w]
        else
          word_sequence[i][j] = 0
        end
      end
    end
  return word_sequence
end

# 素性作成
function make_feature(d, word_d, session_ids, bowlen, boblen, poplen)
  unigram, bigram, sequence, norm_sequence, message_length, pos, speaker = [], [], [], [], [], [], []
  da_function, da_domain = [], []
  da_function_label, da_domain_label = [], []
  @printf "push ...\n"
  for sid in session_ids
    for i in collect(1:length(d[sid][:utterance]))
      push!(unigram, d[sid][:unigram][i]), push!(bigram, d[sid][:bigram][i]), push!(sequence, d[sid][:sequence][i])
      push!(norm_sequence, d[sid][:norm_sequence][i]), push!(message_length, d[sid][:message_length][i])
      push!(da_function, d[sid][:function_da][i]), push!(da_domain, d[sid][:domain_da][i]), push!(pos, d[sid][:pos][i])
      push!(speaker, d[sid][:speaker][i] - 1)
      push!(da_function_label, d[sid][:function_da_label][i][1]), push!(da_domain_label, d[sid][:domain_da_label][i][1])
    end
  end
  @printf "to sequence ...\n"
  unigram = convert_sequence(unigram, word_d[:unigram], 0)
  bigram = convert_sequence(bigram, word_d[:bigram], 0)
  pos = convert_sequence(pos, word_d[:pos], 0)
  @printf "to bag-of-words ...\n"
  bow  = prep.sequences_to_matrix(unigram, bowlen, "count")
  @printf "to bag-of-bigram ...\n"
  bob = prep.sequences_to_matrix(bigram, boblen, "count")
  @printf "to bag-of-pos ...\n"
  bop = prep.sequences_to_matrix(pos, poplen, "count")
  return Dict(:unigram => unigram, :bigram =>  bigram, :sequence => sequence, :norm_sequence => norm_sequence, :da_domain_label => da_domain_label, :da_function_label => da_function_label,
   :message_length => message_length, :bow => bow, :bob => bob, :da_domain => da_domain, :da_function => da_function, :pos => pos, :bop => bop, :speaker => speaker)
end

# 単語出現頻度
function word_dict(d, session_ids)
  word_d = Dict(:unigram => [], :bigram => [], :pos => [])
  for ids in session_ids
     for i in  collect(length(d[ids][:unigram]))
       for x in d[ids][:unigram][i]; push!(word_d[:unigram], x) ; end
       for x in d[ids][:bigram][i]; push!(word_d[:bigram], x) ; end
       for x in d[ids][:pos][i]; push!(word_d[:pos], x) ; end
    end
  end
  return Dict(:unigram => counter(word_d[:unigram]),
    :bigram => counter(word_d[:bigram]), :pos => counter(word_d[:pos]))
end

# 論文記載以外のDAが定義された発話をスキップ
function skip_label(feature)
  skipindex = []
  skipindex! = []
  for i in collect(1:length(feature[:da_function]))
    in(feature[:da_function][i], [201, 192]) || feature[:da_domain_label][i] == "NULL"  && push!(skipindex, i)
    !(in(feature[:da_function][i], [201, 192]) || feature[:da_domain_label][i] == "NULL") && push!(skipindex!, i)
  end
  for k in keys(feature)
    feature[k] = feature[k][skipindex!, :]
  end
  return feature, skipindex!
end

# 素性作成
function get_features()
  @printf "read corpus ...\n"
  d = open(deserialize, "./data/annotated_corpus.dic")
  @printf "add feature ...\n"
  session_ids = filter(x -> !(in(x, [:unigram, :bigram])), collect(keys(d)))
  d = add_feature(d, session_ids)
  @printf "make feature ...\n"
  word_d = word_dict(d, session_ids)
  feature = make_feature(d, word_d, session_ids, int(length(keys(word_d[:unigram]))), int(length(keys(word_d[:bigram]))), int(length(keys(word_d[:pos]))))
  return d, word_d, feature
end

# baseline
function baseline_randomforest(tgt)
  @printf "preprocessing ..\n"
  corpus_d, word_d, feature = get_features()
  @printf "skip label ...\n"
  feature, skipindex! = skip_label(feature)
  @printf "make kfold validation dataset ...\n"
  data = prep.cross_validation(collect(1:length(feature[tgt])), 10)
  @printf "start 10 fold validation ...\n"
  for line_id in zip(data[:,1], data[:, 2])
    train_id, test_id = line_id[1], line_id[2]
    @printf "split feature ...\n"
    train_feature, test_feature = Dict(), Dict()
    for k in keys(feature)
      train_feature[k]  = feature[k][train_id, :]
      test_feature[k] =  feature[k][test_id, :]
    end
    @printf "train model ...\n"
    train = hcat(train_feature[:bow], train_feature[:bob], train_feature[:speaker], train_feature[:norm_sequence], train_feature[:message_length])
    test = hcat(test_feature[:bow], test_feature[:bob], test_feature[:speaker], test_feature[:norm_sequence], test_feature[:message_length])
    model = RandomForestClassifier(n_estimators=200, max_depth=100, max_features="auto", criterion="gini") # n_estimator = 200 max_deps=200 auto
    model = fit!(model, train, vec(train_feature[tgt]))
    @printf "predict ...\n"
    output = predict(model, test)
    println(classification_report(vec(test_feature[tgt]), output, digits=4))
    #return
  end
end

# proposal
function proposal_neuralnet(tgt)
  maxlen = 200
  @printf "preprocessing ..\n"
  corpus_d, word_d, feature = get_features()
  @printf "skip label ...\n"
  feature, skipindex! = skip_label(feature)
  @printf "pad sequence ...\n"
  feature[:unigram] = prep.pad_sequence(feature[:unigram], maxlen)
  feature[:bigram] = prep.pad_sequence(feature[:bigram], maxlen)
  feature[:pos] = prep.pad_sequence(feature[:pos], maxlen)
  @printf "dummy variable...\n"
  feature[tgt] = vec(feature[tgt])
  feature[tgt], d = prep.get_dummies(feature[tgt])
  @printf "make kfold validation dataset ...\n"
  data = prep.cross_validation(collect(1:length(feature[tgt][:, 1])), 10)
  @printf "start 10 fold validation ...\n"
  for line_id in zip(data[:,1], data[:, 2])
    train_id, test_id = line_id[1], line_id[2]
    @printf "split feature ...\n"
    train_feature, test_feature = Dict(), Dict()
    for k in keys(feature)
      train_feature[k]  = feature[k][train_id, :]
      test_feature[k] =  feature[k][test_id, :]
    end
    @printf "read model ...\n"
    # dense_dim, embedding_dim, sequence_maxlen, max_features_uni, max_features_pos, max_features_bi, bowlen, boblen, boplen, glove_embeddings
    bowlen, boblen, boplen = length(feature[:bow][1, :]), length(feature[:bob][1, :]), length(feature[:bop][1, :])
    glove_embeddings = ""
    #glove_embeddings = Dict(i => lowercase(w[1]) for (i, w) in enumerate(sort(collect(word_d[:unigram]), by=x->x[2], rev = true)))
    model = prep.nn_classifier(length(feature[tgt][1, :]), 200, maxlen, bowlen, boblen, boplen, bowlen, boblen, boplen, glove_embeddings)
    @printf "model fit ...\n"
    batch_size, nb_epoch = 128, 30 # 216 25 dim: 200
    feature_list =  ["bow", "bob", "speaker", "norm_sequence", "message_length"]
    model = prep.best_fit(model, 1, train_feature, test_feature, string(tgt), feature_list, batch_size, nb_epoch)
    @printf "predict ...\n"
    prep.predict(model, test_feature, tgt, feature_list)
    #println(classification_report(vec(test_feature[:da_function]), output, digits=4))
  end
end

# baseline
# baseline_randomforest(:da_domain)

# proposal
proposal_neuralnet(:da_domain_label)

# baseline
# acc .8574
# acc


text:
######### StanSample Bernoulli example  ###########

using StanSample, DataFrames, Test

ProjDir = @__DIR__

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir across multiple runs to prevent re-compilation
tmpdir = joinpath(ProjDir, "tmp")
isdir(tmpdir) &&  rm(tmpdir; recursive=true)
sm = SampleModel("bernoulli", bernoulli_model, tmpdir);
rc1 = stan_sample(sm; data);

if success(rc1)
  st = read_samples(sm)
  #display(DataFrame(st))
end
@test size(DataFrame(st), 1) == 4000

isdir(tmpdir) && rm(tmpdir; recursive=true)
sm = SampleModel("bernoulli", bernoulli_model, tmpdir);
rc2 = stan_sample(sm; use_cpp_chains=true, data);

if success(rc2)
  st = read_samples(sm)
  #display(DataFrame(st))
end
@test size(DataFrame(st), 1) == 4000

isdir(tmpdir) && rm(tmpdir; recursive=true)
sm = SampleModel("bernoulli", bernoulli_model, tmpdir);
rc3 = stan_sample(sm; use_cpp_chains=true, check_num_chains=false,
    num_cpp_chains=2, num_julia_chains=2, data);

if success(rc3)
  st = read_samples(sm)
  #display(DataFrame(st))
end
@test size(DataFrame(st), 1) == 4000

isdir(tmpdir) && rm(tmpdir; recursive=true)
sm = SampleModel("bernoulli", bernoulli_model, tmpdir);
rc4 = stan_sample(sm; use_cpp_chains=true, check_num_chains=false,
  num_cpp_chains=4, num_julia_chains=4, data);

if success(rc4)
  st = read_samples(sm)
  #display(DataFrame(st))
end
@test size(DataFrame(st), 1) == 16000

isdir(tmpdir) && rm(tmpdir; recursive=true)
sm = SampleModel("bernoulli", bernoulli_model, tmpdir);
rc4 = stan_sample(sm; use_cpp_chains=true, check_num_chains=false,
  num_cpp_chains=1, num_julia_chains=4, data);

if success(rc4)
  st = read_samples(sm)
  #display(DataFrame(st))
end
@test size(DataFrame(st), 1) == 4000

isdir(tmpdir) && rm(tmpdir; recursive=true)
sm = SampleModel("bernoulli", bernoulli_model, tmpdir);
rc4 = stan_sample(sm; use_cpp_chains=true, check_num_chains=false,
  num_cpp_chains=4, num_julia_chains=1, data);

if success(rc4)
  st = read_samples(sm)
  #display(DataFrame(st))
end
@test size(DataFrame(st), 1) == 4000

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
###############################################################################
#
# Initialize drought legacy
#
###############################################################################
"""
    inititialize_legacy!(
                hs::Union{LeafHydraulics,RootHydraulics,StemHydraulics})
    inititialize_legacy!(tree::GrassLikeOrganism)
    inititialize_legacy!(tree::PalmLikeOrganism)
    inititialize_legacy!(tree::TreeLikeOrganism)
    inititialize_legacy!(tree::TreeSimple)

Initialize the drought legacy effects in the xylem, given
- `hs` [`AbstractHydraulicOrgan`] type struct
- `tree` [`TreeSimple`](@ref) type struct
"""
function inititialize_legacy!(
            hs::Union{LeafHydraulics,RootHydraulics,StemHydraulics}
)
    hs.k_history .= 1;
    hs.p_history .= 0;

    return nothing
end

function inititialize_legacy!(tree::GrassLikeOrganism)
    for root in tree.roots
        inititialize_legacy!(root);
    end
    for leaf in tree.leaves
        inititialize_legacy!(leaf);
    end

    return nothing
end

function inititialize_legacy!(tree::PalmLikeOrganism)
    for root in tree.roots
        inititialize_legacy!(root);
    end
    inititialize_legacy!(tree.trunk);
    for leaf in tree.leaves
        inititialize_legacy!(leaf);
    end

    return nothing
end

function inititialize_legacy!(tree::TreeLikeOrganism)
    for root in tree.roots
        inititialize_legacy!(root);
    end
    inititialize_legacy!(tree.trunk);
    for stem in tree.branch
        inititialize_legacy!(stem);
    end
    for leaf in tree.leaves
        inititialize_legacy!(leaf);
    end

    return nothing
end

function inititialize_legacy!(tree::TreeSimple)
    inititialize_legacy!(tree.root);
    inititialize_legacy!(tree.stem);
    inititialize_legacy!(tree.leaf);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
cd("/Users/adrienpajon/Documents/all.nosync/git_repo/pulse_analysis/2020-11-12 xp pulse all seuil/Julia_analysis")

using CSV
using DataFrames
using JLD

# include("functions/prony.jl")
include("functions/pronyAlgo.jl")

#using Pyplot
pyplot()

Ts_MTiXsens=0.01
Ts_IBS = 0.02
Ts_Xsensor=1/55.2; #correctif dans la pratique on est plus vers 55.2Hz
convertlbpo2Nm = 0.11298482933333
#MTiXsens
#filenameMTiXsens : namefile of each experiment
#filenameMTiXsensID : name ID of each sensor with keys that refers to the position on the body
#sensorXsensNames : keys of each sensor equivalent to their position on the body
# include("imports/Import_filename_MTiXsens.jl")
filenameMTiXsens=load("imports/import_filename_MTiXsens.jld")["filenameMTiXsens"]
sensorXsensNames=load("imports/import_filename_MTiXsens.jld")["sensorXsensNames"]
filenameMTiXsensID=load("imports/import_filename_MTiXsens.jld")["filenameMTiXsensID"]

#Xsensor
#filenameXsensor : namefile of each experiment
# include("imports/Import_filename_Xsensor.jl")
filenameXsensor=load("imports/Import_filename_Xsensor.jld")["filenameXsensor"]

#Instrumented_backseat
#filenameIBS : namefile of each experiment
# include("imports/Import_filename_IBS.jl")
filenameIBS=load("imports/Import_filename_IBS.jld")["filenameIBS"]

println("FileName import finished")

####################################################################
#script to read CSV and store in
# MTiXsens
# Xsensor
# IBS
encoderFolderPath="../measures/encoder_file/20201026_test_1p/"
IBSFolderPath="../measures/IBS/"
XsensorFolderPath="../measures/Xsensor/"
MTiXsensFolderPath="../measures/MTiXsens/"
include("scripts/script_read_CSV.jl")

####################################################################
#script to extract data from MTiXsens and Xsensor in
# AccFree_MTiXsens
# Xsensor_reduce
include("scripts/script_data_extract.jl")
include("scripts/script_reduce_IBS.jl")

# include("imports/Import_synchro.jl")
# sync_OLD=sync

sync=Dict()

#MTiXsens est en Packet counter

sync[1]=(MTiXsens=0,Xsensor=0,Instrumented_backseat=0)
sync[2]=(MTiXsens=0,Xsensor=0,Instrumented_backseat=0)
#MTiXsens est en Packet counter

########
index_test_MTiXsens=searchAllSynchroMTiXsens(AccFree_MTiXsens,filenameMTiXsens,20,50)
println(index_test_MTiXsens)
# plot(pTotal[1])

######
if size(Xsensor_reduce[filenameXsensor[1,1]],2)>6
    m=10
else
    m=4
end
index_test_Xsensor=searchAllSynchro(Xsensor_reduce,m,filenameXsensor,20,1*10^4)
println(index_test_Xsensor)

######
index_test_IBS=searchAllSynchro(IBS,2,filenameIBS,20,1*10^5)
println(index_test_IBS)

# sync = Dict{String, Int64}()
sync=Dict()
sync_=Dict()
for k=1:length(index_test_MTiXsens)
    sync[k]=(MTiXsens=index_test_MTiXsens[k]-100,
                Xsensor=index_test_Xsensor[k]-55,
                IBS=index_test_IBS[k]-50)
    sync_[k]=(MTiXsens=index_test_MTiXsens[k]-100,Xsensor=1,IBS=1)
end
step_max_perturb=400

save("imports/import_synchro.jld", "sync", sync)
save("imports/Import_synchro_perturb.jld", "sync_", sync_,
                                            "step_max_perturb",step_max_perturb)


text:
# This script solves two capital adjustment models in two different ways. Adding
# no structure to the problem versus adding as much structure as possible.
# 	- no structure: state-action formulation
# 	- structure: add as much information as possible
# todo: benchmark compared to POMDPs.jl instead of my own package
# https://stackoverflow.com/questions/52885318/porting-an-example-from-quantecon-jl-to-pomdps-jl

using StructuralDDP
using DataFrames, Gadfly
using InteractiveUtils, BenchmarkTools

include("../test/models/CapitalAdjustModel.jl")
include("../test/models/CapitalAdjustModel2.jl")
createmodel(model::Symbol; kwargs...) = eval(model)(; kwargs...)

par_1 = Dict(:nK => 150, :nz => 5, :γ => 2., :β=>0.9)
p_1_sa = createmodel(:CapitalAdjustModel; par_1..., intdim=:All)
p_1_str = createmodel(:CapitalAdjustModel; par_1..., intdim=:Separable_ExogStates)

# precompile
solve(p_1_sa; concavity = false, monotonicity = false, rewardcall=:pre)
solve(p_1_str, concavity = true, monotonicity = true,
	rewardcall=:pre_partial)

# measure
t_1_sa = @timed solve(p_1_sa)
t_1_str = @timed solve(p_1_str, concavity = true, monotonicity = true,
	rewardcall=:pre_partial)

par_2 = Dict(:nK => 15, :nN=>15, :nz=>3, :β=>0.9)
p_2_sa = createmodel(:CapitalAdjustModel2; par_2..., intdim=:All)
p_2_str = createmodel(:CapitalAdjustModel2; par_2..., intdim=:Separable_ExogStates)

solve(p_2_sa, concavity = false, monotonicity = false, rewardcall=:pre)
solve(p_2_str, concavity = true, monotonicity = true, rewardcall=:pre_partial)

t_2_sa = @timed solve(p_2_sa,
	concavity = false, monotonicity = false, rewardcall=:pre)
t_2_str = @timed solve(p_2_str,
	concavity = true, monotonicity = true, rewardcall=:pre_partial)

test1 = "one choice variable"
test2 = "two choice variables"
method1 = "standard"
method2 = "structured"
df = DataFrame(test = test1, method = method1, time = t_1_sa[2])
	push!(df, [test1, method2, t_1_str[2]])
	push!(df, [test2, method1, t_2_sa[2]])
	push!(df, [test2, method2, t_2_str[2]])

plt = Gadfly.plot(df,
    x = :test,
    y = :time,
	color = :method,
    Scale.y_log10,
    Guide.ylabel("seconds (orders of magnitude)"),
    Guide.xlabel(nothing),
    # Coord.Cartesian(xmin=1,xmax=13.3,ymin=-0.5,ymax=4.2),
    Theme(
        guide_title_position = :left,
        colorkey_swatch_shape = :circle,
        minor_label_font = "Georgia",
        major_label_font = "Georgia"),
	)

draw(SVGJS("compare.svg"), plt)


text:
using Adapt
using CUDA: CuArray
using Oceananigans.Fields: fill_halo_regions!
using Oceananigans.Architectures: arch_array

import Oceananigans.TurbulenceClosures: ivd_upper_diagonal,
                                        ivd_lower_diagonal

abstract type AbstractGridFittedBoundary <: AbstractImmersedBoundary end

const GFIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:AbstractGridFittedBoundary}

#####
##### GridFittedBoundary
#####

struct GridFittedBoundary{S} <: AbstractGridFittedBoundary
    mask :: S
end

@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary) = ib.mask(node(c, c, c, i, j, k, underlying_grid)...)

#####
##### GridFittedBottom
#####

"""
    GridFittedBottom(bottom)

Return an immersed boundary.
"""
struct GridFittedBottom{B} <: AbstractGridFittedBoundary
    bottom :: B
end

@inline function is_immersed(i, j, k, underlying_grid, ib::GridFittedBottom)
    x, y, z = node(c, c, c, i, j, k, underlying_grid)
    return z < ib.bottom(x, y)
end

@inline function is_immersed(i, j, k, underlying_grid, ib::GridFittedBottom{<:AbstractArray})
    x, y, z = node(c, c, c, i, j, k, underlying_grid)
    return @inbounds z < ib.bottom[i, j]
end

const ArrayGridFittedBottom = GridFittedBottom{<:Array}
const CuArrayGridFittedBottom = GridFittedBottom{<:CuArray}

function ImmersedBoundaryGrid(grid, ib::Union{ArrayGridFittedBottom, CuArrayGridFittedBottom})
    # Wrap bathymetry in an OffsetArray with halos
    arch = grid.architecture
    bottom_field = Field{Center, Center, Nothing}(grid)
    bottom_data = arch_array(arch, ib.bottom)
    bottom_field .= bottom_data
    fill_halo_regions!(bottom_field, arch)
    offset_bottom_array = dropdims(bottom_field.data, dims=3)
    new_ib = GridFittedBottom(offset_bottom_array)
    return ImmersedBoundaryGrid(grid, new_ib)
end

Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom))

#####
##### Implicit vertical diffusion
#####

####
#### For a center solver we have to check the interface "solidity" at faces k+1 in both the Upper diagonal and the Lower diagonal
#### (because of tridiagonal convention where lower_diagonal on row k is found at k-1)
#### Same goes for the face solver, where we check at centers k in both Upper and lower diagonal
####

@inline ivd_immersed_solid_interface(LX, LY, ::Center, i, j, k, ibg) = solid_interface(LX, LY, Face(), i, j, k+1, ibg)
@inline ivd_immersed_solid_interface(LX, LY, ::Face, i, j, k, ibg)   = solid_interface(LX, LY, Center(), i, j, k, ibg)

# extending the upper and lower diagonal functions of the batched tridiagonal solver

for location in (:upper_, :lower_)
    alt_func = Symbol(:_ivd_, location, :diagonal)
    func     = Symbol(:ivd_ , location, :diagonal)
    @eval begin
        @inline function $alt_func(i, j, k, ibg::GFIBG, LX, LY, LZ, clock, Δt, interp_κ, κ)
            return ifelse(ivd_immersed_solid_interface(LX, LY, LZ, i, j, k, ibg),
                          zero(eltype(ibg.grid)),
                          $func(i, j, k, ibg.grid, LX, LY, LZ, clock, Δt, interp_κ, κ))
        end
        @inline $func(i, j, k, ibg::GFIBG, LX, LY, LZ::Face, clock, Δt, interp_κ, κ) =
                $alt_func(i, j, k, ibg::GFIBG, LX, LY, LZ, clock, Δt, interp_κ, κ)
        @inline $func(i, j, k, ibg::GFIBG, LX, LY, LZ::Center, clock, Δt, interp_κ, κ) =
                $alt_func(i, j, k, ibg::GFIBG, LX, LY, LZ, clock, Δt, interp_κ, κ)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "libass"
version = v"0.14.0"
# Collection of sources required to build LibassBuilder
sources = [
    "https://github.com/libass/libass/releases/download/0.14.0/libass-0.14.0.tar.xz" =>
    "881f2382af48aead75b7a0e02e65d88c5ebd369fe46bc77d9270a94aa8fd38a2",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
apk add nasm
cd libass-0.14.0/
sed -i 's/9.10.3/2.9.1/' configure.ac
autoreconf
# Grumble-grumble apple grumble-grumble broken linkers...
if [[ ${target} == *-apple-* ]]; then
    export AR=/opt/${target}/bin/ar
fi
./configure --prefix=$prefix --host=$target --disable-require-system-font-provider
make -j${ncore}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libass", :libass),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-4/build_FreeType2.v2.10.0.jl",
    "https://github.com/SimonDanisch/FribidiBuilder/releases/download/0.14.0/build_fribidi.v0.14.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Bzip2-v1.0.6-0/build_Bzip2.v1.0.6.jl",
    "https://github.com/ianshmean/ZlibBuilder/releases/download/v1.2.11/build_Zlib.v1.2.11.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
type ProbabilityDensityDistribution{V} <: PDistribution
	states::Array{V,1}
	pdfs::Array{Function,1}

	function ProbabilityDensityDistribution{V}(_states::Array{V,1}, _pdfs::Array{Function,1})
		if length(_states) == length(_pdfs) && valid_pdfs(_pdfs)
			new(_states,_pdfs)
		else
			throw("Mismatch between number of states and pdfs when constructing ProbabilityDensityDistribution")
		end
	end
end

ProbabilityDensityDistribution{V}(_states::Array{V,1}, _pdfs::Array{Function,1}) = ProbabilityDensityDistribution{V}(_states,_pdfs)

function valid_pdfs(pdfs::Array{Function,1})
	for f in pdfs
		if !verify_real_to_real(f)
			return false
		end
	end
	true
end

states(pdd::ProbabilityDensityDistribution) = pdd.states
pdfs(pdd::ProbabilityDensityDistribution) = pdd.pdfs


text:
m = @model y ~ x1 fe = pn  vcov = cluster(pid1)
x = reg(df, m, drop_singletons = false)


text:
export graph2pdf,
	graph2ps

################################################################################

graph2ps(file_in, file_out="out.ps") = graph2format(file_in, "ps", file_out)
graph2pdf(file_in, file_out="out.pdf") = graph2format(file_in, "pdf", file_out)

function graph2format(file_in, format, file_out)
	g = load_graphs_from_file(file_in)[1] # Only take the first start point of the graph
	s = graph2dot(g)

	open(`dot -T$(format) -o $(file_out)`, write = true) do io
		print(io,s)
	end
end

function graph2dot(g::QuellerGraph)
	body = mapreduce(node2dot, *, g.nodes)

	return """
		digraph {
			rankdir=TB

		$(body)
		}
		"""
end

################################################################################

node2dot(n::Start) =
	"""
		$(getid(n)) [shape=ellipse, style=filled, fillcolor=green, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::End) =
	"""
		$(getid(n)) [shape=diamond, style=filled, fillcolor=red, label="$(escape_string(getmsg(n)))"];

	"""

node2dot(n::Dummy) =
	"""
		$(getid(n)) [shape=ellipse, style=filled, fillcolor=lightblue];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::JumpToGraph) =
	"""
		$(getid(n)) [shape=octagon, style=filled, fillcolor=grey, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::ReturnFromGraph) =
	"""
		$(getid(n)) [shape=octagon, style=filled, fillcolor=grey, label="$(escape_string(getmsg(n)))"];

	"""

################################################################################

node2dot(n::PerformAction) =
	"""
		$(getid(n)) [shape=box, style=filled, fillcolor=purple, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::BinaryCondition) =
	"""
		$(getid(n)) [shape=box, style=filled, fillcolor=yellow, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.n_true)) [label = "True"];
		$(getid(n)) -> $(getid(n.n_false)) [label = "False"];

	"""

function node2dot(n::MultipleChoice)
	dot_node = """
		$(getid(n)) [shape=box, style=filled, fillcolor=yellow, label="$(escape_string(getmsg(n)))"];
	"""

	for (i, nx) in enumerate(n.nexts)
		dot_node *= """
			$(getid(n)) -> $(getid(nx)) [label = "$(i)"];
		"""
	end
	return dot_node*'\n'
end

################################################################################

node2dot(n::CheckStrategy) =
	"""
		$(getid(n)) [shape=ellipse, style=filled, fillcolor=orange, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.n_true)) [label="True"];
		$(getid(n)) -> $(getid(n.n_false)) [label="False"];

	"""

node2dot(n::SetStrategy) =
	"""
		$(getid(n)) [shape=box, style=filled, fillcolor=orange, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::CheckActiveDie) =
	"""
		$(getid(n)) [shape=ellipse, style=filled, fillcolor=pink, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.n_true)) [label="True"];
		$(getid(n)) -> $(getid(n.n_false)) [label="False"];

	"""

node2dot(n::SetActiveDie) =
	"""
		$(getid(n)) [shape=box, style=filled, fillcolor=pink, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next)) [label="Die Available"];
		$(getid(n)) -> $(getid(n.no_die)) [label="No Matching Die"];

	"""

node2dot(n::UseActiveDie) =
	"""
		$(getid(n)) [shape=hexagon, style=filled, fillcolor=pink, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::GetAvailableDice) =
	"""
		$(getid(n)) [shape=ellipse, style=filled, fillcolor=pink, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::SetRingAvailable) =
	"""
		$(getid(n)) [shape=box, style=filled, fillcolor=pink, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

node2dot(n::SetMoDTAvailable) =
	"""
		$(getid(n)) [shape=box, style=filled, fillcolor=pink, label="$(escape_string(getmsg(n)))"];
		$(getid(n)) -> $(getid(n.next));

	"""

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
function calc_mean_Δv(λ::AbstractVector{T1}, λ_min::Real, λ_max::Real ) where { T1<:Real }
    idx_min = searchsortedfirst(λ,λ_min)
    idx_max = searchsortedlast(λ,λ_max)
    @assert 1 <= idx_min <= length(λ)
    @assert 1 <= idx_max <= length(λ)
    Δv = log(λ[idx_max]/λ[idx_min])/(idx_max-idx_min) * RvSpectMLBase.speed_of_light_mps
end

""" `find_ranges_with_tellurics( spectra )`
Return Dataframe of non-overlapping, sorted wavelength ranges that were marked as having tellurics in any of provided spectra.
Calls instrument-specific `find_ranges_with_tellurics` on each spectrum.
"""
function find_ranges_with_tellurics(spectra::AS; min_order::Integer = 1, max_order::Integer = size(first(spectra).λ,2), min_Δv_clean::Real = default_min_Δv_clean, telluric_threshold::Real = 1 ) where { ST<:AbstractSpectra, AS<:AbstractArray{ST} }
    inst_mod = get_inst_module(first(spectra).inst)
    list_of_ranges = mapreduce(s->inst_mod.find_ranges_with_tellurics(s, min_order=min_order, max_order=max_order, min_Δv_clean=min_Δv_clean, telluric_threshold=telluric_threshold), append!, spectra)
    #return list_of_ranges
    if size(list_of_ranges,1)<=1   return list_of_ranges   end
    #sort!(list_of_ranges, by=x->x[1] )
    sort!(list_of_ranges, :lambda_lo )
    non_overlapping_ranges = merge_sorted_wavelength_ranges(list_of_ranges, min_Δv_clean=min_Δv_clean)
    return non_overlapping_ranges
end

""" `merge_sorted_wavelength_ranges( df; min_Δv_clean )`
Return Dataframe of non-overlapping, sorted wavelength ranges (keys lambda_lo and lambda_hi) that is (approximately) the union of all the wavelength ranges in the input dataframe.
min_Δv_clean limits results included in output.
Input dataframe are assumed to be sorted by :lambda_lo
"""
function merge_sorted_wavelength_ranges(df::DataFrame; min_Δv_clean::Real = default_min_Δv_clean) where { T<:Real }
    @assert hasproperty(df, :lambda_lo)
    @assert hasproperty(df, :lambda_hi)
    @assert issorted(df[!,:lambda_lo])
    if size(df,1) <= 1    return df    end
    c_mps = RvSpectMLBase.speed_of_light_mps
    non_overlapping_ranges = DataFrame(:lambda_lo=>[], :lambda_hi=>[])
    last_range_lo = df[1,:lambda_lo]
    last_range_hi = df[1,:lambda_hi]
    for i in 2:size(df,1)
        this_range = df[i,:]
        if this_range.lambda_lo - last_range_hi > min_Δv_clean/c_mps*last_range_hi
            push!(non_overlapping_ranges, Dict(:lambda_lo=>last_range_lo, :lambda_hi=>last_range_hi))
            last_range_lo = this_range.lambda_lo
            last_range_hi = this_range.lambda_hi
        else
            last_range_hi = this_range.lambda_hi
        end
    end
    push!(non_overlapping_ranges, Dict(:lambda_lo=>last_range_lo, :lambda_hi=>last_range_hi))
    return non_overlapping_ranges
end

""" `calc_complement_wavelength_ranges( df_in ; lambda_start, lambda_stop )`
Return DataFrame with the compliment of wavelength ranges in input DataFrame.
Inputs:
- `df_in`: DataFrame containing `lambda_lo` and `lambda_hi`
Optional Inputs:
- `λ_start`:   Output range should start at `λ_start` rather than first entry in df
- `λ_stop`:   Output range should end at `λ_stop` rather than first entry in df
Returns:
- `df_out`: DataFrame containing `lambda_lo` and `lambda_hi`
"""
function calc_complement_wavelength_ranges(df_in::DataFrame; λ_start::Real = df_in[1,:lambda_lo], λ_stop::Real=df_in[end,:lambda_hi] )
    @assert hasproperty(df_in,:lambda_lo)
    @assert hasproperty(df_in,:lambda_hi)
    df_out = λ_start < df_in[1,:lambda_lo] ? DataFrame( :lambda_lo=>[λ_start], :lambda_hi=>[df_in[1,:lambda_lo]] ) :
                                                     DataFrame( :lambda_lo=>[], :lambda_hi=>[] )
    append!(df_out, DataFrame( :lambda_lo=>df_in[1:end-1,:lambda_hi],  :lambda_hi=>df_in[2:end,:lambda_lo] ) )
    if  df_in[end,:lambda_hi] < λ_stop
        append!(df_out, Dict( :lambda_lo=>df_in[end,:lambda_hi], :lambda_hi=>λ_stop ) )
    end
    df_out
end

#global already_printed_stuff = 0

""" `make_ranges_without_tellurics( telluric_list ; lambda_start, lambda_stop, min_Δv, max_Δv )`
Return DataFrame with the compliment of wavelength ranges in input DataFrame.
Inputs:
- `telluric_list`: DataFrame containing `lambda_lo` and `lambda_hi` for each wavelength range to be excluded due to tellurics
Optional Inputs:
- `min_Δv`:  Don't include chunks that are smaller than `min_Δv` (2*max barycentric correction)
- `max_Δv`:  Split up any chunks larger than `max_Δv` (Inf)
- `λ_start`:   Output range should start at `λ_start` rather than first entry in telluric_list
- `λ_stop`:   Output range should end at `λ_stop` rather than first entry in telluric_list
"""
function make_ranges_without_tellurics(telluric_list::DataFrame; min_Δv::Real = 2*RvSpectMLBase.max_bc, max_Δv::Real =Inf, λ_start::Real = telluric_list[1,:lambda_lo], λ_stop::Real=telluric_list[end,:lambda_hi] )
    @assert issorted(telluric_list.lambda_lo)
    c = RvSpectMLBase.speed_of_light_mps
    df_full = calc_complement_wavelength_ranges(telluric_list,λ_start=λ_start,λ_stop=λ_stop)
    df_large_chunks = df_full |> @filter((_.lambda_hi-_.lambda_lo)*2/(_.lambda_hi+_.lambda_lo) >= min_Δv/c) |> DataFrame
    df_nice_sized_chunks =  df_large_chunks |> @filter((_.lambda_hi-_.lambda_lo)*2/(_.lambda_hi+_.lambda_lo) <= max_Δv/c) |> DataFrame
    if max_Δv < Inf
        df_too_large_chunks = df_large_chunks |> @filter((_.lambda_hi-_.lambda_lo)*2/(_.lambda_hi+_.lambda_lo) > max_Δv/c) |> DataFrame
        for row in eachrow(df_too_large_chunks)
            #global already_printed_stuff
            #if already_printed_stuff > 100   continue   end
            #  TODO: Pay attention to where spectrum has natural breaks (order boundaries, no lines, etc.)
            Δv_to_be_divided = (row.lambda_hi-row.lambda_lo)*2/(row.lambda_hi+row.lambda_lo) * c
            println("# row: λ= ", row.lambda_lo, " - ", row.lambda_hi, " Δv_to_be_divided= ",Δv_to_be_divided)
            Δv_to_be_divided = log(row.lambda_hi/row.lambda_lo)  * c
            num_chunks = ceil(Int, Δv_to_be_divided/max_Δv )
            Δv_new_chunks = Δv_to_be_divided/num_chunks
            println("# Δv_to_be_divided = ", Δv_to_be_divided, " num_chunks = ", num_chunks, " Δv_new_chunks= ",Δv_new_chunks)
            new_λ_lo = exp.(log(row.lambda_lo) .+ Δv_new_chunks./c .* (0:(num_chunks-1)))
            new_λ_hi = exp.(log(row.lambda_lo) .+ Δv_new_chunks./c .* (1:num_chunks))
            append!(df_nice_sized_chunks,Dict( :lambda_lo=>new_λ_lo, :lambda_hi=>new_λ_hi ) )
            #already_printed_stuff += 1
        end
        sort!(df_nice_sized_chunks, :lambda_lo )
    end
    return df_nice_sized_chunks
    #return df_large_chunks
end

function calc_complement_index_ranges(idx_all::AR1, idx_bad::AA1 ) where { AR1 <: AbstractRange, AR2 <: AbstractRange, AA1 <: AbstractArray{AR2} }
    function is_in_any_range(x::eltype(AR1), r::eltype(AA1) )
        any(map(b->b.start<=x<=b.stop, idx_bad))
    end
    is_not_in_any_range(x::eltype(AR1), r::eltype(AA1) ) = !is_in_any_range(x,r)
    idx_out = Array{typeof(idx_all),1}()
    start = 0
    stop = 0
    while stop < last(idx_all)
        start = findnext(x->is_not_in_any_range(x,idx_bad), idx_all,stop+1)
        if start == nothing
            break
        end
        stop = findnext(x->is_in_any_range(x,idx_bad), idx_all, start)
        if stop == nothing
            stop = length(idx_all)
        else
            stop -= 1
        end
        push!(idx_out,idx_all[start]:idx_all[stop])
    end
    return idx_out
end

function make_bad_pixel_list_into_ranges(x::V) where { T<:Integer, V<:AbstractVector{T} }
    if length(x) == 0
        return []
    elseif length(x) == 1
        return [x[1]:x[1]]
    else
        tmp = findall(x[2:end].-x[1:end-1].!=1)
        pix_start = vcat( x[1], x[tmp.+1] )
        pix_stop = vcat( x[tmp], x[length(x)] )
        return map(p->p[1]:p[2], zip(pix_start,pix_stop) )
    end
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
# To test
# compile with ragel:
#   ~/Apps/ragel-julia/src/ragel -U ~/Dropbox/Workspace/julia/julia-ragel/test-ragel-moabs-2.rl
# the run the test-ragel-moabs-3 which calles the parse moabs function:
#   ~/Apps/julia-0.3/julia ~/Dropbox/Workspace/julia/julia-ragel/test-ragel-moabs-3.jl
#
module MOABSParserImpl
using Lumberjack
using Switch

type MOABSCoverage
    coverdict::Dict{ASCIIString,Int32}
    id::ASCIIString
    start_pos::ASCIIString
    stop_pos::ASCIIString
    score::ASCIIString
    tot_c::Int32
    bis_c::Int32
end

function parse_moabs(data::String, dict::Dict{ASCIIString,Int32}; gcdisable=false)
    #gc_disable()
    if gcdisable
         println("DISABLE GARBAGE COLLECTION")
        gc_disable()
    end
    moabs_dict= MOABSCoverage( dict, "","","","",0,0)
    mark=0
    p=0
    pe=length(data)

_Test2_actions = Int8  [ 0, 1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1, 10, 1, 11, 1, 12, 0  ]
_Test2_key_offsets = Int8  [ 0, 0, 7, 9, 12, 14, 17, 22, 28, 30, 33, 35, 38, 46, 55, 0  ]
_Test2_trans_keys = Int8  [ 9, 48, 57, 65, 90, 97, 122, 48, 57, 9, 48, 57, 48, 57, 9, 48, 57, 46, 97, 110, 48, 57, 9, 46, 97, 110, 48, 57, 48, 57, 9, 48, 57, 48, 57, 9, 48, 57, 9, 43, 45, 46, 48, 57, 65, 90, 9, 10, 43, 45, 46, 48, 57, 65, 90, 48, 57, 65, 90, 97, 122, 0  ]
_Test2_single_lengths = Int8  [ 0, 1, 0, 1, 0, 1, 3, 4, 0, 1, 0, 1, 2, 3, 0, 0  ]
_Test2_range_lengths = Int8  [ 0, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 0  ]
_Test2_index_offsets = Int8  [ 0, 0, 5, 7, 10, 12, 15, 20, 26, 28, 31, 33, 36, 42, 49, 0  ]
_Test2_trans_cond_spaces = Int8  [ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0  ]
_Test2_trans_offsets = Int8  [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 0  ]
_Test2_trans_lengths = Int8  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0  ]
_Test2_cond_keys = Int8  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0  ]
_Test2_cond_targs = Int8  [ 2, 1, 1, 1, 0, 3, 0, 4, 3, 0, 5, 0, 6, 5, 0, 7, 7, 7, 7, 0, 8, 7, 7, 7, 7, 0, 9, 0, 10, 9, 0, 11, 0, 12, 11, 0, 13, 13, 13, 13, 13, 0, 13, 14, 13, 13, 13, 13, 0, 1, 1, 1, 0, 0  ]
_Test2_cond_actions = Int8  [ 3, 0, 0, 0, 0, 5, 0, 7, 0, 0, 9, 0, 11, 0, 0, 13, 13, 13, 13, 0, 15, 0, 0, 0, 0, 0, 17, 0, 19, 0, 0, 21, 0, 23, 0, 0, 0, 0, 0, 0, 0, 0, 0, 25, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0  ]
Test2_start  :: Int  = 14
Test2_first_final  :: Int  = 14
Test2_error  :: Int  = 0
Test2_en_main  :: Int  = 14

begin
cs = Test2_start;
	end

begin
_acts = 0
_keys = 0
_ckeys = 0
if p == pe
    @goto _test_eof

end
if cs == 0
    @goto _out

end
@label _resume
_keys = int(_Test2_key_offsets[1 + cs] );
	_trans = convert(Uint ,_Test2_index_offsets[1 + cs]);

	_klen = convert(Int ,_Test2_single_lengths[1 + cs]);
	if _klen > 0
    begin
_lower = 0
_mid = 0
_upper = 0
_lower = _keys;
		_upper = _keys + _klen - 1;
		while ( true  )
    begin
if _upper < _lower
    break;

end
_mid = _lower + ((_upper-_lower) >> 1);
			if ( data[1 + p ]) < _Test2_trans_keys[1 + _mid ]
    _upper = _mid - 1;

elseif ( ( data[1 + p ]) > _Test2_trans_keys[1 + _mid ] )
    _lower = _mid + 1;

else
    begin
_trans += convert(Uint ,(_mid - _keys));
				@goto _match
end

end
end

end
_keys += _klen;
		_trans += convert(Uint ,_klen);
	end

end
_klen = convert(Int ,_Test2_range_lengths[1 + cs]);
	if _klen > 0
    begin
_lower = 0
_mid = 0
_upper = 0
_lower = _keys;
		_upper = _keys + (_klen<<1) - 2;
		while ( true  )
    begin
if _upper < _lower
    break;

end
_mid = _lower + (((_upper-_lower) >> 1) & ~1);
			if ( data[1 + p ]) < _Test2_trans_keys[1 + _mid ]
    _upper = _mid - 2;

elseif ( ( data[1 + p ]) > _Test2_trans_keys[1 + _mid + 1 ] )
    _lower = _mid + 2;

else
    begin
_trans += convert(Uint ,((_mid - _keys)>>1));
				@goto _match
end

end
end

end
_trans += convert(Uint ,_klen);
	end

end
@label _match
_ckeys = int(_Test2_trans_offsets[1 + _trans] );
	_klen = convert(Int ,_Test2_trans_lengths[1 + _trans]);
	_cond = convert(Uint ,_Test2_trans_offsets[1 + _trans]);

	_cpc = 0;
	@switch _Test2_trans_cond_spaces[1 + _trans]  begin
end
begin
_lower = 0
_mid = 0
_upper = 0
_lower = _ckeys;
		_upper = _ckeys + _klen - 1;
		while ( true  )
    begin
if _upper < _lower
    break;

end
_mid = _lower + ((_upper-_lower) >> 1);
			if _cpc < convert(Int ,_Test2_cond_keys[1 + _mid ])
    _upper = _mid - 1;

elseif ( _cpc > convert(Int ,_Test2_cond_keys[1 + _mid ]) )
    _lower = _mid + 1;

else
    begin
_cond += convert(Uint ,(_mid - _ckeys));
				@goto _match_cond
end

end
end

end
cs = 0;
		@goto _again
end
@label _match_cond
cs = convert(Int ,_Test2_cond_targs[1 + _cond]);

	if _Test2_cond_actions[1 + _cond] == 0
    @goto _again

end
_acts = int(_Test2_cond_actions[1 + _cond] );
	_nacts = convert(Uint ,_Test2_actions[1 + _acts ]); _acts += 1
;
	while ( _nacts > 0  )
    begin
@switch _Test2_actions[1 + _acts ] begin
@case 0
begin
mark=p
    end

	break;
	@case 1
begin
moabs_dict.id=data.data[mark+1:p]
    end

	break;
	@case 2
begin
mark=p
    end

	break;
	@case 3
begin
moabs_dict.start_pos= ASCIIString( data.data[mark+1:p])
    end

	break;
	@case 4
begin
mark=p
    end

	break;
	@case 5
begin
moabs_dict.stop_pos= ASCIIString( data.data[mark+1:p])
    end

	break;
	@case 6
begin
mark=p
    end

	break;
	@case 7
begin
moabs_dict.score=  ASCIIString( data.data[mark+1:p])
    end

	break;
	@case 8
begin
mark=p
    end

	break;
	@case 9
begin
moabs_dict.tot_c=parseint( Int32, ASCIIString( data.data[mark+1:p]) )
    end

	break;
	@case 10
begin
mark=p
    end

	break;
	@case 11
begin
moabs_dict.bis_c=parseint( Int32, ASCIIString( data.data[mark+1:p]) )
    end

	break;
	@case 12
begin
str=string( moabs_dict.id,".", moabs_dict.start_pos,".", moabs_dict.stop_pos )
       if !haskey( moabs_dict.coverdict, str)
           moabs_dict.coverdict[ str ] = moabs_dict.tot_c
       else
           moabs_dict.coverdict[ str ] += moabs_dict.tot_c
       end
    end

	break;
		end
_nacts -= 1
;
		_acts += 1
;
	end

end
@label _again
if cs == 0
    @goto _out

end
p += 1
;
	if p != pe
    @goto _resume

end
@label _test_eof
begin
end
@label _out
begin
end
end
  Lumberjack.info("garbage clean")
  #if gcdisable
  #    gc_enable()
  #    gc()
  #end
  gc_enable()
  gc()
  Lumberjack.info("done garbage clean")
return moabs_dict
end

end


text:
abstract type CategoricalWeightLearner <: WeightLearner end
abstract type CategoricalWeightLearnerFit <: WeightLearner end

struct GBHLearner <: CategoricalWeightLearner
    pi0estimator::Pi0Estimator
end

GBHLearner() = GBHLearner(Storey(0.5))

function learn_weights(
    gbh::GBHLearner,
    Ps_train,
    Xs_train::AbstractVector{<:CategoricalValue},
    Xs_test::AbstractVector{<:CategoricalValue},
    α,
    multiple_testing_method,
)
    #TODO: check levels are the same

    function gbh_wt(Ps)
        pi0 = MultipleTesting.estimate(Ps, gbh.pi0estimator)
        (1 - pi0) / pi0
    end

    groupwise_wt = map(gbh_wt, group(Xs_train, Ps_train))
    ws = ones(size(Xs_test, 1))
    @inbounds for i = 1:size(Xs_test, 1)
        ws[i] = get(groupwise_wt, Xs_test[i], 0.0)
    end
    ws, groupwise_wt
end

struct GrenanderLearner <: CategoricalWeightLearner end

function _Fdr(ts, Fs, ms)
    num = sum(ts .* ms)
    denom = sum(Fs .* ms)
    iszero(num) ? 0.0 : num / denom
end

function _Fdr_linearized(ts, Fs, ms, α)
    num = sum(ts .* ms)
    denom = sum(Fs .* ms)
    num - α * denom
end

function lagrange_balance(λ, grs, ms, α, ::BenjaminiHochberg; linearized = false)
    ts = IndependentHypothesisWeighting.invert_subgradient.(grs, λ)
    ts_L = first.(ts)
    ts_R = last.(ts)
    Fs_L = cdf.(grs, ts_L)
    Fs_R = cdf.(grs, ts_R)

    if linearized
        Fdr_L = _Fdr_linearized(ts_L, Fs_L, ms, α)
        Fdr_R = _Fdr_linearized(ts_R, Fs_R, ms, α)
    else
        Fdr_L = _Fdr(ts_L, Fs_L, ms)
        Fdr_R = _Fdr(ts_R, Fs_R, ms)
    end
    # wt_sum = Interval( ((sum( ms .* first.(ts)),sum( ms .* last.(ts))) .* m ./ 0.1 ./m)...) #want this to be 1
    Interval(Fdr_L, Fdr_R)
end

function learn_weights(
    ::GrenanderLearner,
    Ps_train,
    Xs_train::AbstractVector{<:CategoricalValue},
    Xs_test::AbstractVector{<:CategoricalValue},
    α,
    ::BenjaminiHochberg,
)
    # check levels are the same

    train_data = group(Xs_train, Ps_train)
    test_ms = groupcount(Xs_test)
    train_data_view = view(train_data, keys(test_ms))
    train_grenander = fit.(Ref(Grenander()), train_data_view)

    all_λs = sort(vcat([gr.fs for gr in train_grenander]...), rev = true)
    α = 0.1
    λ_idx = findfirst(
        λ -> α ∈ lagrange_balance(λ, train_grenander, test_ms, α, BenjaminiHochberg()),
        all_λs,
    )
    λ_opt = all_λs[λ_idx]

    ts = invert_subgradient.(train_grenander, λ_opt)
    ts_L = first.(ts)
    ts_R = last.(ts)
    Fs_L = cdf.(train_grenander, ts_L)
    Fs_R = cdf.(train_grenander, ts_R)

    lin_Fdr = lagrange_balance(
        λ_opt,
        train_grenander,
        test_ms,
        α,
        BenjaminiHochberg();
        linearized = true,
    )
    tmp_comb = last(lin_Fdr) ./ Intervals.span(lin_Fdr)

    ts_mixed = tmp_comb .* ts_L .+ (1 - tmp_comb) .* ts_R
    ts_mixed = sum(test_ms) .* ts_mixed ./ sum(ts_mixed .* test_ms)

    ws = ones(size(Xs_test, 1))
    @inbounds for i = 1:size(Xs_test, 1)
        ws[i] = get(ts_mixed, Xs_test[i], 0.0)
    end
    ws, ts_mixed
end


text:
struct NodeItem
    min_x::Float64
    min_y::Float64
    max_x::Float64
    max_y::Float64
    offset::UInt64
end

mutable struct FlatGeobuffer
    header::Header
    rtree::Vector{NodeItem}
    io::IO
    offset::Int
    offsets::Vector{Int}
    filtered::Bool
end

Base.show(io::IO, fgb::FlatGeobuffer) = print(io, "$(fgb.filtered ? "Filtered " : "")FlatGeobuffer with $(length(fgb)) $(fgb.header.geometry_type) features.")