text:
using RobustAndOptimalControl, ControlSystems, MonteCarloMeasurements

d = δr()
@test d.val == 0
@test d.radius == 1

d2 = d+d
@test d2.val == 0
@test d2.radius == 2
@test d2.name === d.name

d0 = d - d
@test d0.val == 0
@test d0.radius == 0
@test d2.name === d.name

d_ = δr()
d0_ = d - d_
@test d0_ isa UncertainSS

sys = uss(d)
@test sys.nx == 0
@test sys.D == [0 1; d.radius d.val] # lft representation of uncertain scalar

D = RobustAndOptimalControl.Δ(2, δc)
@test rand(D, 100) isa Matrix{Complex{Particles{Float64, 100}}}

# Example from Mu tools user guide
d = 2.4 + 0.4δr()
sys = uss(d)
@test sys.D == [0 1; 0.4 2.4]
isys = inv(sys)
@test isys.D ≈ [-1/6 1/2.4; -1/6 1/2.4] # Mu user guide seems wrong here

P = uss([δc(), δc()])
@test iszero(P.D11)
@test P.D21 == I
@test P.D12 == I
@test iszero(P.D22)
# W = tf(1.0, [1, 1])
# W*P

δel = δr(2, 1, :δel)
η = δr(6, 1, :η)
ρ = δr(-1, 1, :ρ)

s1 = uss(δel)
s2 = uss(η)

@test (1/η).D ≈ [
    -0.1667    0.1667
    -0.1667    0.1667
] atol=1e-3

@test (δel*η).D ≈ [
    0     1     6
    0     0     1
    1     2    12
] atol = 1e-3

@test (δel*η + ρ).D == [
    0     1     0     6
    0     0     0     1
    0     0     0     1
    1     2     1    11
]

@test (δel/η).D ≈ [
    0   -0.1667    0.1667
    0   -0.1667    0.1667
1.0000   -0.3333    0.3333
] atol=1e-3

@test svdvals([η+3+δel δel/η].D)[1] ≈ 11.1896 rtol=1e-2

A = [[3+δel+η δel/η]; [7+ρ ρ+δel*η]]
@test svdvals(A.D)[1] ≈ 15.3545 rtol=0.01

# lft(A.D, A.Δ) not yet supported
# A.Δ
# A2 = lft(ss(A.sys), ss(Diagonal(A.Δ)), :u)
# @test A2 ≈ A

δel_mcm = rand(δel, 100)
η_mcm = rand(η, 100)
ρ_mcm = rand(ρ, 100)

A_mcm = [[3+δel_mcm+η_mcm δel_mcm/η_mcm]; [7+ρ_mcm ρ_mcm+δel_mcm*η_mcm]]

As = rand(A, 100)
@test all(As.D .≈ A_mcm)

##
for i = 1:2, j = 1:2, k = 1:2
    d1 = δss(i,j)
    d2 = δss(j,k)
    p = d1*d2
    @test p.ny == d1.ny
    @test p.nu == d2.nu
    @test size(p.M) == (j+k, i+j)
end

##
@test uss(I(2)).M.D == I
@test sum(δss(2, 2).D) == 4
@test sum(δss(4, 4).D) == 8

##
H = δss(2, 3)
W = tf([1, .1],[.1, 1]) .* I(2)
WH = W*H
@test WH.D == [zeros(3,2) I(3); 10I(2) zeros(2,3)]

##

w = 2π .* exp10.(LinRange(-2, 2, 500))
W = makeweight(0.40,15,3)
Pn = tf(1, [1/60, 1]) |> ss
d = δss(1,1)

Gu = rand(d, 100)
@test Gu.nx == d.Δ[1].nx
@test Gu.ny == d.Δ[1].ny
@test Gu.nu == d.Δ[1].nu

@test d.ny == 1
@test d.nu == 1
@test d.nz == 1
@test d.nw == 1

@test d.zinds == 1:1
@test d.winds == 1:1
@test d.uinds == 2:2
@test d.yinds == 2:2

temp = (W*d)
@test temp.nu == temp.ny == 1
@test temp.nz == temp.nw == 1

temp = I(1) + W*d
@test temp.nu == temp.ny == 1
@test temp.nz == temp.nw == 1

@test length(d.Δ) == 1

Pd = Pn*(I(1) + W*d)

usyss = sminreal(system_mapping(Pd))
@test !iszero(usyss.B)
@test !iszero(usyss.C)
@test (usyss.C*usyss.B)[] == 60
@test usyss == Pn

Pp = rand(Pd, 50)
@test Pp.nx == 1+1+2 # == nom, W, 2 sample unc

if isinteractive()
    bodeplot(Pp, w, ylims=(1e-1, 1e1))
    bodeplot!(Pn, w)
end

Pn = ssrand(3,4,5)
@test δss(4,4, bound=0.2).D == [0I(4) sqrt(0.2)*I(4); sqrt(0.2)*I(4) 0I(4)]
mu = I(4) + δss(4,4, bound=0.2)

@test mu.D == [0I(4) sqrt(0.2)*I(4); sqrt(0.2)*I(4) I(4)]
@test size(Pn) == size(convert(UncertainSS, Pn))
@test convert(UncertainSS, Pn).D22 == Pn.D

P = Pn*mu
Pn2 = system_mapping(P)
@test Pn2 == Pn
@test size(P) == (7,8)
@test size(P.M) == (4,4)

## this time mimo real
delta = uss([δr(), δr()])
a = 1
P = ss([0 a; -a -1], I(2), [1 a; 0 1], 0)* (I(2) + delta)
K = ss(I(2))

G = lft(P, -K)
hn = norm(G, Inf)

w = exp10.(LinRange(-5, 1, 100))
M = freqresp(G.M, w)
M = permutedims(M, (2,3,1))
# mu = mussv_DG(M)
# maximum(mu)
# # maximum(structured_singular_value(M))
# @test 1/maximum(mu) ≈ √(2) atol=0.01

## Test block structure
d = δr()
@test RobustAndOptimalControl.makeblock([d]) == [-1, 0]
@test RobustAndOptimalControl.makeblock([d,d]) == [-2, 0]

d = δc()
@test RobustAndOptimalControl.makeblock([d]) == [1,0]
@test RobustAndOptimalControl.makeblock([d,d]) == [2, 0]
@test RobustAndOptimalControl.makeblock(δss(2,3).Δ) == [2, 3]

Δ = [δr(0,1,:qkarne), δr(0,1,:chuck_norris)]
@test RobustAndOptimalControl.block_structure(Δ)[1] == [[-1, 0], [-1, 0]]

Δ = [δr(0,1,:qkarne), δc(0,1,:chuck_norris)]
@test RobustAndOptimalControl.block_structure(Δ)[1] == [[1, 0], [-1, 0]] # The two elements should have been sorted and 1,-1 are therefore switched

Δ = [δr(0,1,:a), δss(2,3).Δ[], δr(0,1,:a), δc(0,1,:chuck_norris), δr(0,1,:a)]

@test RobustAndOptimalControl.block_structure(Δ)[1] == [
    [2, 3], # The gensym name will come first
    [-3, 0],
    [1, 0],
]

perm = RobustAndOptimalControl.block_structure(Δ)[2]
@test perm == [2,1,3,5,4]

P = partition(ssrand(8, 7, 2), 6, 7)
P = UncertainSS(P, Δ)

blocks, M = RobustAndOptimalControl.blocksort(P)
@test blocks == [
    [2, 3], # The gensym name will come first
    [-3, 0],
    [1, 0],
]

@test M[7,6] == P.M[6, 5] # test that the permutation was correctly applied, only test for the complex permutation
@test M[1:3,1:2] == P.M[2:4, 2:3]# test that the permutation was correctly applied, only test for the matrix block

## Test structured_singular_value with different kinds of blocks
w = exp10.(LinRange(-2, 2, 500))
delta = δss(1,1)
P = (tf(1,[1, .2, 1])) * (1+0.2*delta)
s = tf("s")
K = ss(1 + 2/s + 0.9s/(0.1s+1))
Gcl = lft(P, -K)

# These should all be equivalent
# Full block complex uncertainty gives μ = hinfnorm
@test norm(Gcl, Inf) ≈ 1.1861236982687908 rtol=1e-3
mu = structured_singular_value(Gcl, w)
@test 1/norm(mu, Inf) ≈ 1/1.1861236982687908 atol=0.005
@test RobustAndOptimalControl.robstab(Gcl) ≈ 0.84 atol=0.005
@test 1/structured_singular_value(Gcl) ≈ 0.84 atol=0.005

## same as above but with scalar instead of 1×1 system
w = exp10.(LinRange(-2, 2, 500))
delta = δc()
P = (tf(1,[1, .2, 1])) * (1+0.2*delta)
s = tf("s")
K = ss(1 + 2/s + 0.9s/(0.1s+1))
Gcl = lft(P, -K)

# These should all be equivalent
# Full block complex uncertainty gives μ = hinfnorm
@test norm(Gcl, Inf) ≈ 1.1861236982687908 rtol=1e-3
mu = structured_singular_value(Gcl, w)
@test 1/norm(mu, Inf) ≈ 1/1.1861236982687908 atol=0.005
@test RobustAndOptimalControl.robstab(Gcl) ≈ 0.84 atol=0.005
@test 1/structured_singular_value(Gcl) ≈ 0.84 atol=0.005

## this time mimo complex
delta = uss([δc(), δc()])
P = ss([0 1; 0 0], I(2), [1 0], 0) * (I(2) + delta)
# diagonal input uncertainty

K = ss([1;1])
G = lft(P, -K)
hn = norm(G, Inf)
@test robstab(G) ≈ 0.4824 atol=0.001
dm = diskmargin(system_mapping(P), K, 0, w)
dmm = argmin(dm->dm.margin, dm.simultaneous_input)
@test dmm.margin ≈ 0.5335 atol = 0.001


text:
#hline!/vline!
for (fname, axis, axis2, axis3, range) in [(:vline!, :y, :y2, :x, :height), (:hline!, :x, :x2, :y, :width)]

        @eval begin function ($fname)(v::VegaVisualization; value::Real = 0, strokeWidth::Real = 1.5, stroke::AbstractString = "#000", strokeDash::Real = 0)
            line = VegaMark(_type = "rule",
                            properties = VegaMarkProperties(enter = VegaMarkPropertySet($axis = VegaValueRef(field = Dict{Any, Any}("group" => $"$axis")),
                                                                                $axis2 = VegaValueRef(field = Dict{Any, Any}("group" => $"$range")),
                                                                                $axis3 = VegaValueRef(scale = $"$axis3", value = value),
                                                                                strokeWidth = VegaValueRef(value = strokeWidth),
                                                                                stroke = VegaValueRef(value = stroke),
                                                                                strokeDash = VegaValueRef(value = strokeDash != 0? [strokeDash] : [])
                                                                                )
                                        )
                            )

            push!(v.marks, line)

            return v
        end
    end
end

ext: jl

max_stars_repo_name: jpthor/vega


text:
module StrongStabilityPreservingRungeKuttaMethod
export StrongStabilityPreservingRungeKutta, updatedt!
export StrongStabilityPreservingRungeKutta33
export StrongStabilityPreservingRungeKutta34

using Requires

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    using .CuArrays
    using .CuArrays.CUDAnative
    using .CuArrays.CUDAnative.CUDAdrv

    include("StrongStabilityPreservingRungeKuttaMethod_cuda.jl")
end

using ..ODESolvers
ODEs = ODESolvers
using ..SpaceMethods

"""
    StrongStabilityPreservingRungeKutta(f, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

        Q̇ = f(Q)

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses either 2nd or 3rd order Strong-Stability-Preserving (SS) time-integrators.

### References
S. Gottlieb and C.-W. Shu,Total variation diminishing Runge-Kutta schemes, Math .Comp .,67 (1998), pp .73–85
"""

struct StrongStabilityPreservingRungeKutta{T, AT, Nstages} <: ODEs.AbstractODESolver
"time step"
dt::Array{T,1}
"time"
t::Array{T,1}
"rhs function"
rhs!::Function
"Storage for RHS during the StrongStabilityPreservingRungeKutta update"
Rstages::NTuple{Nstages,AT}
"Storage for RHS during the StrongStabilityPreservingRungeKutta update"
Qstages::NTuple{Nstages,AT}
"RK coefficient vector A (rhs scaling)"
RKA::Array{T,2}
"RK coefficient vector B (rhs add in scaling)"
RKB::Array{T,1}
"RK coefficient vector C (time scaling)"
RKC::Array{T,1}

function StrongStabilityPreservingRungeKutta(rhs!::Function, RKA, RKB, RKC, Q::AT; dt=nothing, t0=0) where {AT<:AbstractArray}

    @assert dt != nothing

    T = eltype(Q)
    dt = [T(dt)]
    t0 = [T(t0)]
    new{T, AT, length(RKB)}(dt, t0, rhs!, ntuple(i->similar(Q), length(RKB)), ntuple(i->similar(Q), length(RKB)), RKA, RKB, RKC)
end
end

function StrongStabilityPreservingRungeKutta(spacedisc::AbstractSpaceMethod, RKA, RKB, RKC,Q; dt=nothing, t0=0)
    rhs! = (x...) -> SpaceMethods.odefun!(spacedisc, x...)
    StrongStabilityPreservingRungeKutta(rhs!, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

function StrongStabilityPreservingRungeKutta33(spacedisc, Q; dt=nothing,t0=0)
    T=eltype(Q)
    RKA = [ T(1) T(0); T(3//4) T(1//4); T(1//3) T(2//3) ]
    RKB = [T(1), T(1//4), T(2//3)]
    RKC = [ T(0), T(1//4), T(2//3) ]
    StrongStabilityPreservingRungeKutta(spacedisc, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

function StrongStabilityPreservingRungeKutta34(spacedisc, Q; dt=nothing,t0=0)
    T=eltype(Q)
    RKA = [ T(1) T(0); T(0) T(1); T(2//3) T(1//3); T(0) T(1) ]
    RKB = [ T(1//2); T(1//2); T(1//6); T(1//2) ]
    RKC = [ T(0); T(1//4); T(2//3); T(3//3) ]
    StrongStabilityPreservingRungeKutta(spacedisc, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

"""
    updatedt!(ssp::StrongStabilityPreservingRungeKutta, dt)

Change the time step size to `dt` for `ssp.
"""
updatedt!(ssp::StrongStabilityPreservingRungeKutta, dt) = ssp.dt[1] = dt

function ODEs.dostep!(Q, ssp::StrongStabilityPreservingRungeKutta, timeend, adjustfinalstep)
    time, dt = ssp.t[1], ssp.dt[1]
    if adjustfinalstep && time + dt > timeend
        dt = timeend - time
        @assert dt > 0
    end
    RKA, RKB, RKC = ssp.RKA, ssp.RKB, ssp.RKC
    rhs! = ssp.rhs!
    Rstages, Qstages = ssp.Rstages, ssp.Qstages
    for s = 1:length(RKB)
        Qstages[s] .= Q
        Rstages[s].Q .= 0
        rhs!(Rstages[s], Qstages[s], time)
        update!(Val(size(Q,2)), Val(size(Q,1)), Rstages[s].Q, Qstages[1].Q, Q.Q, Q.realelems, RKA[s,1], RKA[s,2], RKB[s], dt)
        time += RKC[s] * dt
    end
    if dt == ssp.dt[1]
        ssp.t[1] += dt
    else
        ssp.t[1] = timeend
    end

end

# {{{ Update solution (for all dimensions)
function update!(::Val{nstates}, ::Val{Np}, rhs::Array{T,3}, Q0, Q, elems, rka1, rka2, rkb, dt) where {nstates, Np, T}
    @inbounds for e = elems, s = 1:nstates, i = 1:Np
        Q[i, s, e] = rka1*Q0[i, s, e] + rka2*Q[i, s, e] + dt*rkb*rhs[i, s, e]
    end
end
# }}}

end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
### MODEL ###
# initialize all JuAFEM objects
#const CVT = Union{Tuple{CellVectorValues,CellScalarValues},Tuple{CellVectorValues}}
"""
    Model{dim,N,D,V,E,P,S,MS,CVT,FV,DH,CH,NBC,BF,SP,C}

A `Model` instance contains all the informations needed to run a simulation.
Type parameters associated to it are used to dispatch methods with respect to model specificities, especially {dim,N,D,V,E,P,S} (see description below).

# Fields
- `grid::Grid{dim}` : JuAFEM.Grid instance
- `dofhandler::DH` : JuAFEM.DofHandler instance containing information about primitive variables discretization over the grid
- `dirichlet_bc::CH` : JuAFEM.ConstraintHandler instance containing dirichlet boundary conditions
- `neumann_bc::NBC` : dictionnary with `set_names`, `traction(x::Vector)` pairs
- `body_forces::Vector{BC}` : body forces for every cell
- `cellvalues_tuple::CVT` : A tuple with as many elements as primitive variables in the model, with each one containing in cell integrations related informations
- `facevalues::FV` : Same as cellvalues but for tractions integration on faces
- `material_properties::Vector{Rheology{Float64,D,V,E,P}}` : A Vector of `Rheology`s instance associated to each cell
- `material_state::Vector{Vector{MS}}` : `MaterialState`s for each integration point
- `K::SP` : global stiffness matrix
- `RHS::Vector{Float64}` : global right hand side vector (external force vector or residual vector in incremental formulation)
- `clock::C` : `Clock` instance containing all time related informations
- `solver::S` : `Solver` subtype instance containing informations about the solving procedure
- `multithreading::Bool` : `True` to use multithreaded assembly

# Type parameters
- `dim::Int` : spatial dimensions of the model
- `N::Int` : number of primitives variables. Commonly 1 (displacement) or 2 (displacement and pressure)
- `D<:Damage` : micromechanical properties of the rheology
- `V<:Viscosity` : viscous properties of the rheology
- `E<:Elasticity` : elastic properties of the rheology
- `P<:Plasticity` : plastic properties of the rheology
- `S<:AbstractSolver` : linear or non linear solver
"""
struct Model{dim,N,D,V,E,P,S,MS,CVT,FV,DH,CH,NBC,BF,SP,C}
    grid::Grid{dim}
    dofhandler::DH
    dirichlet_bc::CH
    neumann_bc::NBC # TODO : use a struct similar to ConstraintHandler may be relevant
    body_forces::Vector{BF}
    cellvalues_tuple::CVT
    facevalues::FV
    material_properties::Vector{Rheology{Float64,D,V,E,P}}
    material_state::Vector{Vector{MS}}
    K::SP
    RHS::Vector{Float64}
    clock::C
    solver::S
    multithreading::Bool
    function Model(grid::Grid{dim},
                   dh::DH,
                   dbc::CH,
                   nbc::NBC,
                   bf::Vector{BF},
                   cv_tuple::CVT,
                   fv::FV,
                   mp::Vector{Rheology{Float64,D,V,E,P}},
                   ms::Vector{Vector{MS}},
                   K::SP,
                   RHS::Vector{Float64},
                   clock::C,
                   solver::S,
                   multithreading::Bool) where {dim,D,V,E,P,S,MS,CVT,FV,DH,CH,NBC,BF,SP,C}

        N = length(cv_tuple)
        return new{dim,N,D,V,E,P,S,MS,CVT,FV,DH,CH,NBC,BF,SP,C}(grid,dh,dbc,nbc,bf,cv_tuple,fv,mp,ms,K,RHS,clock,solver,multithreading)
    end
end

# Model{dim,N,D,V,E,P,S,MS,CVT,FV,DH,CH,NBC,SP,C}(grid,dh,dbc,nbc,cv_tuple,fv,mp,ms,K,RHS,clock,solver,multithreading)  = Model(grid::Grid{dim},
#                dh::DH,
#                dbc::CH,
#                nbc::NBC,
#                cv_tuple::CVT,
#                fv::FV,
#                mp::Vector{Rheology{Float64,D,V,E,P}},
#                ms::Vector{Vector{MS}},
#                K::SP,
#                RHS::Vector{Float64},
#                clock::C,
#                solver::S,
#                multithreading::Bool) where {dim,D,V,E,P,S,MS,CVT,FV,DH,CH,NBC,SP,C}
#

# Constructor using input file user defined variables
"""
    Model(grid::Grid, variables::PrimitiveVariables, quad_order::Int, quad_type::Symbol, bc::BoundaryConditions, body_forces::BodyForces, rheology::Rheology, clock::Clock, solver::Solver, multithreading::Bool = false)

`Model` constructor. Keyword arguments can also be used but all arguments are mandatory.

# Fields
- `grid::Grid{dim}` : JuAFEM.Grid instance
- `variables::PrimitiveVariables` : objet containing informations about primitive variables and their interpolation
- `quad_order::Int` : quadrature order to evaluate the element integrals
- `quad_type::Symbol` : quadrature type. Can be `:legendre` or `:lobatto`
- `bc::BoundaryConditions` : contains dictionaries mapping boundary condition functions to sets of nodes/elements/faces.
- `body_forces::BodyForces` : contains body forces as an homogeneous vectorial quantity or as a function of a space location vector.
- `rheology::Rheology` : contains material properties
- `clock::Clock` : `Clock` instance containing all time related informations
- `solver::AbstractSolver` : `Solver` subtype instance containing informations about the solving procedure
- `multithreading::Bool` : `True` to use multithreaded assembly

"""
function Model(grid::Grid, variables::PrimitiveVariables,
               quad_order::Int, quad_type::Symbol,
               bc::BoundaryConditions, body_forces::BodyForces, rheology::Rheology,
               initial_state, clock::Clock, solver::AbstractSolver, multithreading)

    dh, dbc, cv_tuple, fv, mp, ms, bf, K, RHS = setup_model(grid, variables, quad_order, quad_type, bc, body_forces, initial_state, rheology)
    nbc = bc.neumann
    return Model(grid,dh,dbc,nbc,bf,cv_tuple,fv,mp,ms,K,RHS,clock,solver,multithreading)
end

Model(; grid, variables, quad_order, quad_type, bc, body_forces, rheology, initial_state, clock, solver, multithreading = false) =
         Model(grid, variables, quad_order, quad_type, bc, body_forces, rheology, initial_state, clock, solver, multithreading)

# Ask for a non linear solver if rheology is damaged and/or plastic
# NLRheology = Union{Rheology{T,Nothing,Nothing,E,P} where {T,E,P<:Plasticity},
#                    Rheology{T,D,Nothing,E,P} where {T,D<:Damage,E,P<:Plasticity},
#                    Rheology{T,D,Nothing,E,Nothing} where {T,E,D<:Plasticity}}
# Model(grid, variables, quad_order, quad_type, bc, rheology::NLRheology, clock, solver::S, multithreading) where {S<:LinearSolver} = @error "Non linear rheology, please provide a non linear solver"

Base.show(io::IO, grid::Grid) = Base.show(io::IO, MIME"text/plain"(), grid::Grid)

function Base.show(io::IO, ::MIME"text/plain", model::Model{dim,D,V,E,P,N,S}) where {dim,D,V,E,P,N,S}
    print(io, "⌗  Model to be solved for $(string(model.dofhandler.field_names)) using $(model.dofhandler.field_interpolations) on a grid $(model.grid) \n")
    #print(io,"\n")

end


text:
using ESDL
using Test
using Dates
import Base.Iterators
using Distributed
using Statistics
addprocs(2)
@everywhere using ESDL, Statistics

@everywhere function sub_and_return_mean(xout1,xout2,xin)
    m=mean(skipmissing(xin))
    for i=1:length(xin)
        xout1[i]=xin[i]-m
    end
    xout2[1]=m
end
function sub_and_return_mean(c::ESDL.AbstractCubeData)
  mapCube(sub_and_return_mean,c,
  indims=InDims("Time"),
  outdims=(OutDims("Time"),OutDims()))
end

function doTests()
  # Test simple Stats first
  c=Cube()

  d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))

  dmem=readcubedata(d)

  @testset "Simple statistics using mapslices" begin
  # Basic statistics
  m=mapslices(mean∘skipmissing,d,dims = "Time")

  @test isapprox(readcubedata(m).data,[282.04703 281.84494 281.65927 282.2035;
    281.90726 281.59995 281.69135 282.38644;
    281.84158 281.35898 281.66083 282.49042;
    281.92462 281.12613 281.50513 282.4352])

  #Test Spatial meann along laitutde axis
  d1=subsetcube(c,variable="gross_primary_productivity",time=(Date("2002-01-01"),Date("2002-01-01")),lon=(30,30))

  dmem=readcubedata(d1)
  mtime=mapslices(mean∘skipmissing,dmem,dims = ("lon","lat"))

  end
  # Test Mean seasonal cycle retrieval
  @testset "Seasonal cycle statistics and anomalies" begin
  cdata=subsetcube(c,variable="soil_moisture",lon=10,lat=50.75)
  d=readcubedata(cdata)

  x2=getMSC(d)

  x3=getMedSC(d)

  a = d[:]
  a = a[3:46:end]
  @test isapprox(mean(skipmissing(a)),x2[3])
  @test isapprox(median(skipmissing(a)),x3[3])

  # Test gap filling
  cube_filled=readcubedata(gapFillMSC(d))
  imiss=findfirst(i->ismissing(i),d.data)
  @test !ismissing(cube_filled.data[imiss])
  its=mod(imiss-1,46)+1
  @test cube_filled.data[imiss]≈readcubedata(x2).data[its]
  @test !any(ismissing(cube_filled.data))

  # Test removal of MSC

  cube_anomalies=readcubedata(removeMSC(cube_filled))
  @test isapprox(cube_anomalies.data[47:92],(cube_filled.data[47:92].-readcubedata(x2).data[1:46]))

  # Test normalization
  anom_normalized=normalizeTS(cube_anomalies)[:]
  #@show cube_anomalies[:,:,:]
  @test mean(anom_normalized)<1e7
  @test 1.0-1e-6 <= std(anom_normalized) <= 1.0+1e-6
  end

  d1=subsetcube(c,variable=["gross_primary_productivity","net_ecosystem_exchange"],lon=(10,10),lat=(50,50))
  d2=subsetcube(c,variable=["gross_primary_productivity","air_temperature_2m"],lon=(10,10),lat=(50,50))

  @testset "Multiple output cubes" begin
  #Test onvolving multiple output cubes
  c1=subsetcube(c,variable="gross_primary_productivity",lon=(10,11),lat=(50,51),time=Date(2001)..Date(2010))

  c2=readcubedata(c1)

  cube_wo_mean,cube_means=sub_and_return_mean(c2)

  @test isapprox(permutedims(c2[:,:,:].-mean(c2[:,:,:],dims=3),(3,1,2)),readcubedata(cube_wo_mean)[:,:,:])
  @test isapprox(mean(c2[:,:,:],dims=3)[:,:,1],cube_means[:,:])
  end
end

@testset "Parallel processing" begin
doTests()
end
rmprocs(workers())

@testset "Single proc processing" begin
doTests()
end


text:
module PopGenSims

using DataFrames, PooledArrays, StaticArrays
using StatsBase: sample
using PopGen:
    PopData,
    Genotype,
    GenoArray,
    get_genotypes,
    read_from,
    write_to,
    sort

include("Cross.jl")
include("Samples.jl")
include("Sibship.jl")
include("Utils.jl")

end


text:
"""
Overview: structures and functions associated to mesh
import and management.
"""

"""
> mutable struct Domain
- Sen : number of supported elements
- Set : type of supported elements
- Senn : number of nodes of the supported elements
- ne : number of elements of each type in the domain
- e2n : volume element connectivity
- mat : material associated to the domain
"""
mutable struct Domain
  #
  Sen::Int64
  Set::Vector{Symbol}
  Senn::Vector{Int64}
  eskip::Vector{Int64}
  ne::Vector{Int64}
  e2n::Vector{Int64}
  mat::material
  #
  Domain() = new()
end

"""
> mutable struct Boundary
- Sen : number of supported elements
- Set : type of supported elements
- Senn : number of nodes of the supported elements
- ne : number of elements of each type in the domain
- e2n : surface element connectivity
- bcdofs : which dofs are constrained
- bcvals : value associated to the boundary condition
"""
mutable struct Boundary
  #
  Sen::Int64
  Set::Vector{Symbol}
  Senn::Vector{Int64}
  eskip::Vector{Int64}
  ne::Vector{Int64}
  e2n::Vector{Int64}
  bcdofs::Vector{Int64}
  bcvals::Vector{Float64}
  #
  Boundary() = new()
  #
end

"""
> mutable struct grid
- netype : number of supperted types of elements
- etype : array that stores the number of nodes for each element
- nn : number of nodes
- nΩ : number of domains
- nΓ : number of boundaries
- n2c : nodes coordinates
- Ω : domains
- Γ : boundaries
- tcn : total cells number
- grid() : empty initialization method
"""
mutable struct Grid
  #
  nn::Int64
  nΩ::Int64
  nΓ::Int64
  #
  n2c::Vector{Float64}
  Ω::Vector{Domain}
  Γ::Vector{Boundary}
  #
  tcn::Int64
  tnn::Int64
  #
  Grid() = new()
  #
end

"""
> read_mesh( mesh_file::String, Ω_list::Vector{Vector{Int64}}, Γ_list::Vector{Vector{Int64}}, mat::Vector{material}, bc_dof::Vector{Vector{Int64}}, bc_vals::Vector{Vector{Float64}})
It reads the mesh file and it inializes the grid data structure
- mesh_file : name of the mesh file to read
- Ω_list : list of domains in which the structures is organised
- Γ_list : list of boundaries
- mat : list of materials associated to each domain
- bc_dof : list of constrained degrees of freedom of each boundary
- bc_vals : boundary condition values
"""
function read_mesh(mesh_file::String, Ω_list,
                   Γ_list, mat::Vector{String},
                   bc_dof::Vector{Vector{Int64}}, bc_vals::Vector{Vector{Float64}})
  mesh = Grid()
  if contains(mesh_file,".med")
    import_mesh!(mesh,mesh_file,Ω_list,Γ_list,mat,bc_dof,bc_vals,Val{:MED})
  elseif contains(mesh_file,".mphtxt")
    import_mesh!(mesh,mesh_file,Ω_list,Γ_list,mat,bc_dof,bc_vals,Val{:MPHTXT})
  else
    println("Mesh format not supported.")
  end
  return mesh
end

"""
> get_coor!(mesh::Grid, conn::Vector{Int64}, X::Vector{Float64}, nn::Int64)
it retrieves element nodes coordinates.
- mesh : Grid
- conn : connectivity of the element
- X : coordinates of the nodes in {x₁, y₁, z₁, ..., xₙ, yₙ, zₙ} format
- nn : number of nodes
"""
function get_coor!(mesh::Grid, conn::Vector{Int64}, X::Vector{Float64}, nn::Int64)
  #
  @inbounds for i = 1:nn
    X[1+(i-1)*dim:i*dim] = mesh.n2c[1+(conn[i]-1)*dim:conn[i]*dim]
  end
  #
  return nothing
end

"""
> get_coor!(mesh::Grid, conn::SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}, X::Vector{Float64}, nn::Int64)
it retrieves element nodes coordinates.
- mesh : Grid
- conn : connectivity of the element
- X : coordinates of the nodes in {x₁, y₁, z₁, ..., xₙ, yₙ, zₙ} format
- nn : number of nodes
"""
function get_coor!(mesh::Grid, conn::SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true},
                   X::Vector{Float64}, nn::Int64)
  #
  @inbounds for i = 1:nn
    X[1+(i-1)*dim:i*dim] = mesh.n2c[1+(conn[i]-1)*dim:conn[i]*dim]
  end
  #
  return nothing
end

function import_mesh!(mesh::Grid, mesh_file::String,
                      Ω_list,
                      Γ_list,
                      mat::Vector{String}, bc_dof::Vector{Vector{Int64}},
                      bc_vals::Vector{Vector{Float64}}, ::Type{Val{:MPHTXT}})
  #
  # it instantiates a mesh_type structure from a COMSOL .mphtxt file
  #
  local neT6,e2nT6,e2gT6
  local neQ9,e2nQ9,e2gQ9
  local neT10,e2nT10,e2gT10
  local neP18,e2nP18,e2gP18
  local neH27,e2nH27,e2gH27
  #
  fname = "./mesh/"*mesh_file
  # extract mesh infos
  open(fname,"r") do fhand
    while !eof(fhand)
      line = readline(fhand)
      if (contains(line,"# number of mesh vertices"))
        line = split(line)
        mesh.nn = Meta.parse(line[1])
        mesh.n2c = Vector{Float64}(undef,mesh.nn*dim)
      elseif (contains(line,"# Mesh vertex coordinates"))
        for i = 1:mesh.nn
          line = split(readline(fhand))
          for j = 1:dim
            mesh.n2c[j+(i-1)*dim] = Meta.parse(line[j])
          end
        end
      elseif (contains(line,"tri2 # type name"))
        line = readline(fhand)
        line = readline(fhand)
        line = readline(fhand)
        line = split(readline(fhand))
        neT6 = Meta.parse(line[1])
        line = readline(fhand)
        e2nT6 = Vector{Int64}(undef,neT6*TR6n)
        e2gT6 = Vector{Int64}(undef,neT6)
        for i = 1:neT6
          line = split(readline(fhand))
          for j = 1:TR6n
            e2nT6[j+(i-1)*TR6n] = Meta.parse(line[j])+1
          end
        end
        while (true)
          line = readline(fhand)
          if (contains(line,"# Geometric entity indices"))
            for i = 1:neT6
              line = readline(fhand)
              e2gT6[i] = Meta.parse(line)+1
            end
            break
          end
        end
      elseif (contains(line,"quad2 # type name"))
        line = readline(fhand)
        line = readline(fhand)
        line = readline(fhand)
        line = split(readline(fhand))
        neQ9 = Meta.parse(line[1])
        line = readline(fhand)
        e2nQ9 = Vector{Int64}(undef,neQ9*QU9n)
        e2gQ9 = Vector{Int64}(undef,neQ9)
        for i = 1:neQ9
          line = split(readline(fhand))
          for j = 1:QU9n
            e2nQ9[j+(i-1)*QU9n] = Meta.parse(line[j])+1
          end
        end
        while (true)
          line = readline(fhand)
          if (contains(line,"# Geometric entity indices"))
            for i = 1:neQ9
              line = readline(fhand)
              e2gQ9[i] = Meta.parse(line)+1
            end
            break
          end
        end
      elseif (contains(line,"tet2 # type name"))
        line = readline(fhand)
        line = readline(fhand)
        line = readline(fhand)
        line = split(readline(fhand))
        neT10 = Meta.parse(line[1])
        line = readline(fhand)
        e2nT10 = Vector{Int64}(undef,neT10*T10n)
        e2gT10 = Vector{Int64}(undef,neT10)
        for i = 1:neT10
          line = split(readline(fhand))
          for j = 1:T10n
            e2nT10[j+(i-1)*T10n] = Meta.parse(line[j])+1
          end
        end
        while (true)
          line = readline(fhand)
          if (contains(line,"# Geometric entity indices"))
            for i = 1:neT10
              line = readline(fhand)
              e2gT10[i] = Meta.parse(line)
            end
            break
          end
        end
      elseif (contains(line,"prism2 # type name"))
        line = readline(fhand)
        line = readline(fhand)
        line = readline(fhand)
        line = split(readline(fhand))
        neP18 = Meta.parse(line[1])
        line = readline(fhand)
        e2nP18 = Vector{Int64}(undef,neP18*P18n)
        e2gP18 = Vector{Int64}(undef,neP18)
        for i = 1:neP18
          line = split(readline(fhand))
          for j = 1:P18n
            e2nP18[j+(i-1)*P18n] = Meta.parse(line[j])+1
          end
        end
        while (true)
          line = readline(fhand)
          if (contains(line,"# Geometric entity indices"))
            for i = 1:neP18
              line = readline(fhand)
              e2gP18[i] = Meta.parse(line)
            end
            break
          end
        end
      elseif (contains(line,"hex2 # type name"))
        line = readline(fhand)
        line = readline(fhand)
        line = readline(fhand)
        line = split(readline(fhand))
        neH27 = Meta.parse(line[1])
        line = readline(fhand)
        e2nH27 = Vector{Int64}(undef,neH27*H27n)
        e2gH27 = Vector{Int64}(undef,neH27)
        for i = 1:neH27
          line = split(readline(fhand))
          for j = 1:H27n
            e2nH27[j+(i-1)*H27n] = Meta.parse(line[j])+1
          end
        end
        while (true)
          line = readline(fhand)
          if (contains(line,"# Geometric entity indices"))
            for i = 1:neH27
              line = readline(fhand)
              e2gH27[i] = Meta.parse(line)
            end
            break
          end
        end
      end
    end
  end
  # if var is not defined and assign default value
  if (!@isdefined e2gT6)
    neT6 = 0
    e2gT6 = [0]
  end
  if (!@isdefined e2gQ9)
    neQ9 = 0
    e2gQ9 = [0]
  end
  if (!@isdefined e2gT10)
    neT10 = 0
    e2gT10 = [0]
  end
  if (!@isdefined e2gP18)
    neP18 = 0
    e2gP18 = [0]
  end
  if (!@isdefined e2gH27)
    neH27 = 0
    e2gH27 = [0]
  end

  nVol = maximum([maximum(e2gT10), maximum(e2gP18), maximum(e2gH27)])
  nSurf = maximum([maximum(e2gT6), maximum(e2gQ9)])

  surf = Array{Int64}(undef,nSurf)
  vol  = Array{Int64}(undef,nVol)
  fill!(surf,-1)
  fill!(vol,-1)

  mesh.nΩ = size(Ω_list)[1]
  mesh.nΓ = size(Γ_list)[1]

  for i = 1:mesh.nΩ
    for j = 1:size(Ω_list[i])[1]
      ith_b = Ω_list[i][j]
      vol[ith_b] = i
    end
  end

  for i = 1:mesh.nΓ
    for j = 1:size(Γ_list[i])[1]
      ith_b = Γ_list[i][j]
      surf[ith_b] = i
    end
  end

  neΓT6  = zeros(Int64,mesh.nΓ)
  neΓQ9  = zeros(Int64,mesh.nΓ)
  neΩT10 = zeros(Int64,mesh.nΩ)
  neΩW18 = zeros(Int64,mesh.nΩ)
  neΩH27 = zeros(Int64,mesh.nΩ)

  # number of nodes for supported surface elements
  Se_nn = [TR3n, TR6n, QU4n, QU8n, QU9n]
  #Set_Γ = [:TR3, :TR6, :QU4, :QU8, :QU9]
  # number of nodes for supported volume elements
  Ve_nn = [TE4n, T10n, HE8n, H20n, H27n, PE6n, P15n, P18n]
  #Set_Ω = [:TE4, :T10, :HE8, :H20, :H27, :PE6, :P15, :P18]
  #
  neΩ = zeros(Int64,(size(Ve_nn)[1], mesh.nΩ ))
  neΓ = zeros(Int64,(size(Se_nn)[1], mesh.nΓ ))

  for i = 1:neT6
    ith_s = e2gT6[i]
    if (surf[ith_s]>0)
      neΓ[2,surf[ith_s]] += 1
    end
  end

  for i = 1:neQ9
    ith_s = e2gQ9[i]
    if (surf[ith_s]>0)
      neΓ[5,surf[ith_s]] += 1
    end
  end

  for i = 1:neT10
    ith_s = e2gT10[i]
    if (vol[ith_s]>0)
      neΩ[2,vol[ith_s]] += 1
    end
  end

  for i = 1:neP18
    ith_s = e2gP18[i]
    if (vol[ith_s]>0)
      neΩ[7,vol[ith_s]] += 1
    end
  end

  for i = 1:neH27
    ith_s = e2gH27[i]
    if (vol[ith_s]>0)
      neΩ[5,vol[ith_s]] += 1
    end
  end
  #
  mesh.Ω = [Domain() for i in 1:mesh.nΩ]
  mesh.Γ = [Boundary() for i = 1:mesh.nΓ]
  # fill mesh data structure
  for b = 1:mesh.nΓ
    iΓ = mesh.Γ[b]
    #
    iΓ.Sen = size(Se_nn)[1]
    iΓ.Set = [:TR3, :TR6, :QU4, :QU8, :QU9]
    iΓ.Senn = Se_nn
    iΓ.ne = neΓ[:,b]
    iΓ.bcdofs = bc_dof[b]
    iΓ.bcvals = bc_vals[b]
    iΓ.eskip = zeros(Int64,iΓ.Sen)
    nn = 0
    for j = 1:iΓ.Sen-1
      nn += neΓ[j,b]*Se_nn[j]
      iΓ.eskip[j+1] = nn
    end
    nn += neΓ[end,b]*Se_nn[end]
    iΓ.e2n = Vector{Int64}(undef,nn)
  end
  #
  for d = 1:mesh.nΩ
    iΩ = mesh.Ω[d]
    #
    iΩ.Sen = size(Ve_nn)[1]
    iΩ.Set = [:TE4, :T10, :HE8, :H20, :H27, :PE6, :P15, :P18]
    iΩ.Senn = Ve_nn
    iΩ.ne = neΩ[:,d]
    iΩ.mat = MORFE_load_material(mat[d])
    iΩ.eskip = zeros(Int64,iΩ.Sen)
    #
    nn = 0
    for j = 1:iΩ.Sen-1
      nn += neΩ[j,d]*Ve_nn[j]
      iΩ.eskip[j+1] = nn
    end
    nn += neΩ[end,d]*Ve_nn[end]
    iΩ.e2n = Vector{Int64}(undef,nn)
  end
  # fill boundaries connectivity
  for b = 1:mesh.nΓ
    iΓ = mesh.Γ[b]
    #
    skip = iΓ.eskip[2]+1
    for e = 1:neT6
      ith_s = e2gT6[e]
      ith_s = surf[ith_s]
      if (ith_s==b)
        iΓ.e2n[skip:skip+TR6n-1] = e2nT6[1+(e-1)*TR6n:e*TR6n]
        skip += TR6n
      end
    end
    #
    skip = iΓ.eskip[5]+1
    for e = 1:neQ9
      ith_s = e2gQ9[e]
      ith_s = surf[ith_s]
      if (ith_s==b)
        iΓ.e2n[skip:skip+QU9n-1] = e2nQ9[1+(e-1)*QU9n:e*QU9n]
        skip += QU9n
      end
    end
    #
  end
  #
  for d = 1:mesh.nΩ
    iΩ = mesh.Ω[d]
    #
    skip = iΩ.eskip[2]
    for e = 1:neT10
      ith_s = e2gT10[e]
      ith_s = vol[ith_s]
      #
      if (ith_s==d)
        iΩ.e2n[skip+1]  =  e2nT10[1+(e-1)*T10n]
        iΩ.e2n[skip+2]  =  e2nT10[2+(e-1)*T10n]
        iΩ.e2n[skip+3]  =  e2nT10[4+(e-1)*T10n]
        iΩ.e2n[skip+4]  =  e2nT10[3+(e-1)*T10n]
        iΩ.e2n[skip+5]  =  e2nT10[5+(e-1)*T10n]
        iΩ.e2n[skip+6]  =  e2nT10[9+(e-1)*T10n]
        iΩ.e2n[skip+7]  =  e2nT10[8+(e-1)*T10n]
        iΩ.e2n[skip+8]  =  e2nT10[7+(e-1)*T10n]
        iΩ.e2n[skip+9]  =  e2nT10[6+(e-1)*T10n]
        iΩ.e2n[skip+10] = e2nT10[10+(e-1)*T10n]
        skip += T10n
      end
      #
    end
    #
    skip = iΩ.eskip[5]+1
    for e = 1:neH27
      ith_s = e2gH27[e]
      ith_s = vol[ith_s]
      if (ith_s==d)
        iΩ.e2n[skip:skip+H27n-1] = e2nH27[1+(e-1)*H27n:e*H27n]
        skip += H27n
      end
    end
    #
    skip = iΩ.eskip[7]
    for e = 1:neP18
      ith_s = e2gP18[e]
      ith_s = vol[ith_s]
      if (ith_s==d)
        iΩ.e2n[skip+1]  =  e2nP18[1+(e-1)*P18n]
        iΩ.e2n[skip+2]  =  e2nP18[2+(e-1)*P18n]
        iΩ.e2n[skip+3]  =  e2nP18[3+(e-1)*P18n]
        iΩ.e2n[skip+4]  =  e2nP18[4+(e-1)*P18n]
        iΩ.e2n[skip+5]  =  e2nP18[5+(e-1)*P18n]
        iΩ.e2n[skip+6]  =  e2nP18[6+(e-1)*P18n]
        iΩ.e2n[skip+7]  =  e2nP18[7+(e-1)*P18n]
        iΩ.e2n[skip+8]  =  e2nP18[9+(e-1)*P18n]
        iΩ.e2n[skip+9]  =  e2nP18[8+(e-1)*P18n]
        iΩ.e2n[skip+10] = e2nP18[16+(e-1)*P18n]
        #
        iΩ.e2n[skip+11] = e2nP18[18+(e-1)*P18n]
        iΩ.e2n[skip+12] = e2nP18[17+(e-1)*P18n]
        iΩ.e2n[skip+13] = e2nP18[10+(e-1)*P18n]
        iΩ.e2n[skip+14] = e2nP18[12+(e-1)*P18n]
        iΩ.e2n[skip+15] = e2nP18[15+(e-1)*P18n]
        #iΩ.e2n[skip+16] = e2nP18[11+(e-1)*P18n]
        #iΩ.e2n[skip+17] = e2nP18[13+(e-1)*P18n]
        #iΩ.e2n[skip+18] = e2nP18[14+(e-1)*P18n]
        #
        skip += P15n
      end
    end
    #
  end
  #
  return nothing
  #
end

"""
> import_mesh!( mesh::Grid, mesh_file::String, Ω_list::Vector{Vector{Int64}}, Γ_list::Vector{Vector{Int64}}, mat::Vector{material}, bc_dof::Vector{Vector{Int64}}, bc_vals::Vector{Vector{Float64}}, ::Type{Val{:MED}})
It reads a MED mesh format
- mesh_file : name of the mesh file
- Ω_list : list of domains in which the structures is organised
- Γ_list : list of boundaries
- mat : list of materials associated to each domain
- bc_dof : list of constrained degrees of freedom of each boundary
- bc_vals : boundary condition values
- :MED : specializes function to read MED mesh format
"""
function import_mesh!(mesh::Grid, mesh_file::String,
                      Ω_list,
                      Γ_list,
                      mat::Vector{String}, bc_dof::Vector{Vector{Int64}},
                      bc_vals::Vector{Vector{Float64}}, ::Type{Val{:MED}})
  #
  pth = joinpath(".","mesh",mesh_file)
  #
  MeshFile = aster_read_mesh(pth)
  # mesh nodes
  nodes = MeshFile["nodes"]
  mesh.nn = length(nodes)
  mesh.n2c = Vector{Float64}(undef,mesh.nn*dim)
  for i = 1:mesh.nn
    for j = 1:dim
      mesh.n2c[j+(i-1)*dim] = nodes[i][j]
    end
  end
  # allocate domains and boundaries
  mesh.nΩ = size(Ω_list)[1]
  mesh.nΓ = size(Γ_list)[1]
  mesh.Ω = [Domain() for i ∈ 1:mesh.nΩ]
  mesh.Γ = [Boundary() for i ∈ 1:mesh.nΓ]
  #
  Elements = MeshFile["elements"]
  ElementTypes = MeshFile["element_types"]
  ElementGroups = MeshFile["element_sets"]
  #
  neΓ = zeros(Int64,(Snse,mesh.nΓ))
  neΩ = zeros(Int64,(Vnse,mesh.nΩ))
  #
  for i = 1:mesh.nΓ
    sname = Γ_list[i]
    selem = ElementGroups[sname]
    for j = 1:size(selem)[1]
      etype = ElementTypes[selem[j]]
      if (etype == :TR3)
        neΓ[1,i] += 1
      elseif (etype == :TR6)
        neΓ[2,i] += 1
      elseif (etype == :QU4)
        neΓ[3,i] += 1
      elseif (etype == :QU8)
        neΓ[4,i] += 1
      elseif (etype == :QU9)
        neΓ[5,i] += 1
      else
        println("Warning, element not supported")
      end
    end
  end
  #
  for i = 1:mesh.nΩ
    vname = Ω_list[i]
    velem = ElementGroups[vname]
    for j = 1:size(velem)[1]
      etype = ElementTypes[velem[j]]
      if (etype == :TE4)
        neΩ[1,i] += 1
      elseif (etype == :T10)
        neΩ[2,i] += 1
      elseif (etype == :HE8)
        neΩ[3,i] += 1
      elseif (etype == :H20)
        neΩ[4,i] += 1
      elseif (etype == :H27)
        neΩ[5,i] += 1
      elseif (etype == :PE6)
        neΩ[6,i] += 1
      elseif (etype == :P15)
        neΩ[7,i] += 1
      elseif (etype == :P18)
        neΩ[8,i] += 1
      else
        println("Warning, element not supported")
      end
    end
  end
  #
  # number of nodes for supported surface elements
  Se_nn = [TR3n, TR6n, QU4n, QU8n, QU9n]
  # number of nodes for supported volume elements
  Ve_nn = [TE4n, T10n, HE8n, H20n, H27n, PE6n, P15n, P18n]
  # allocate connectivity
  for i = 1:mesh.nΩ
    iΩ = mesh.Ω[i]
    #
    iΩ.Sen = size(Ve_nn)[1]
    iΩ.Set = [:TE4, :T10, :HE8, :H20, :H27, :PE6, :P15, :P18]
    iΩ.Senn = Ve_nn
    iΩ.ne = neΩ[:,i]
    iΩ.mat = MORFE_load_material(mat[i])
    iΩ.eskip = zeros(Int64,iΩ.Sen)
    #
    nn = 0
    for j = 1:iΩ.Sen-1
      nn += neΩ[j,i]*Ve_nn[j]
      iΩ.eskip[j+1] = nn
    end
    nn += neΩ[end,i]*Ve_nn[end]
    iΩ.e2n = Vector{Int64}(undef,nn)
  end
  #
  for i = 1:mesh.nΓ
    iΓ = mesh.Γ[i]
    #
    iΓ.Sen = size(Se_nn)[1]
    iΓ.Set = [:TR3, :TR6, :QU4, :QU8, :QU9]
    iΓ.Senn = Se_nn
    iΓ.ne = neΓ[:,i]
    iΓ.bcdofs = bc_dof[i]
    iΓ.bcvals = bc_vals[i]
    iΓ.eskip = zeros(Int64,iΓ.Sen)
    #
    nn = 0
    for j = 1:iΓ.Sen-1
      nn += neΓ[j,i]*Se_nn[j]
      iΓ.eskip[j+1] = nn
    end
    nn += neΓ[end,i]*Se_nn[end]
    iΓ.e2n = Vector{Int64}(undef,nn)
  end
  # fill boundaries with mesh connectivity
  for i = 1:mesh.nΓ
    iΓ = mesh.Γ[i]
    sname = Γ_list[i]
    selem = ElementGroups[sname]
    ecount = zeros(Int64,size(Se_nn)[1] )
    for j = 1:size(selem)[1]
      etype = ElementTypes[selem[j]]
      if (etype == :TR3)
        skip = iΓ.eskip[1]+ecount[1]*Se_nn[1]
        iΓ.e2n[1+skip:skip+Se_nn[1]] = Elements[selem[j]]
        ecount[1] += 1
      elseif (etype == :TR6)
        skip = iΓ.eskip[2]+ecount[2]*Se_nn[2]
        iΓ.e2n[1+skip:skip+Se_nn[2]] = Elements[selem[j]]
        ecount[2] += 1
      elseif (etype == :QU4)
        skip = iΓ.eskip[3]+ecount[3]*Se_nn[3]
        iΓ.e2n[1+skip:skip+Se_nn[3]] = Elements[selem[j]]
        ecount[3] += 1
      elseif (etype == :QU8)
        skip = iΓ.eskip[4]+ecount[4]*Se_nn[4]
        iΓ.e2n[1+skip:skip+Se_nn[4]] = Elements[selem[j]]
        ecount[4] += 1
      elseif (etype == :QU9)
        skip = iΓ.eskip[5]+ecount[5]*Se_nn[5]
        iΓ.e2n[1+skip:skip+Se_nn[5]] = Elements[selem[j]]
        ecount[5] += 1
      else
        println("Warning, element not supported")
      end
    end
    # fill domains with mesh connectivity
    for i = 1:mesh.nΩ
      iΩ = mesh.Ω[i]
      sname = Ω_list[i]
      velem = ElementGroups[sname]
      ecount = zeros(Int64,size(Ve_nn)[1] )
      for j = 1:size(velem)[1]
        etype = ElementTypes[velem[j]]
        if (etype == :TE4)
          skip = iΩ.eskip[1]+ecount[1]*Ve_nn[1]
          iΩ.e2n[1+skip:skip+Ve_nn[1]] = Elements[velem[j]]
          ecount[1] += 1
        elseif (etype == :T10)
          skip = iΩ.eskip[2]+ecount[2]*Ve_nn[2]
          iΩ.e2n[1+skip:skip+Ve_nn[2]] = Elements[velem[j]]
          ecount[2] += 1
        elseif (etype == :HE8)
          skip = iΩ.eskip[3]+ecount[3]*Ve_nn[3]
          iΩ.e2n[1+skip:skip+Ve_nn[3]] = Elements[velem[j]]
          ecount[3] += 1
        elseif (etype == :H20)
          skip = iΩ.eskip[4]+ecount[4]*Ve_nn[4]
          iΩ.e2n[1+skip:skip+Ve_nn[4]] = Elements[velem[j]]
          ecount[4] += 1
        elseif (etype == :H27)
          skip = iΩ.eskip[5]+ecount[5]*Ve_nn[5]
          iΩ.e2n[1+skip:skip+Ve_nn[5]] = Elements[velem[j]]
          ecount[5] += 1
        elseif (etype == :PE6)
          skip = iΩ.eskip[6]+ecount[6]*Ve_nn[6]
          iΩ.e2n[1+skip:skip+Ve_nn[6]] = Elements[velem[j]]
          ecount[6] += 1
        elseif (etype == :P15)
          skip = iΩ.eskip[7]+ecount[7]*Ve_nn[7]
          iΩ.e2n[1+skip:skip+Ve_nn[7]] = Elements[velem[j]]
          ecount[7] += 1
        elseif (etype == :P18)
          skip = iΩ.eskip[8]+ecount[8]*Ve_nn[8]
          iΩ.e2n[1+skip:skip+Ve_nn[8]] = Elements[velem[j]]
          ecount[8] += 1
        else
          println("Warning, element not supported")
        end
      end
    end
  end
  #
  mesh.tcn = 0
  mesh.tnn = 0
  for i = 1:mesh.nΩ
    iΩ = mesh.Ω[i]
    for se = 1:iΩ.Sen
      mesh.tcn += iΩ.ne[se]
      mesh.tnn += iΩ.ne[se]*iΩ.Senn[se]
    end
  end
  #
  return nothing
end


text:
# ---
# title: 1233. Remove Sub-Folders from the Filesystem
# id: problem1233
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, String
# link: <https://leetcode.com/problems/remove-sub-folders-from-the-filesystem/description/>
# hidden: true
# ---
#
# Given a list of folders, remove all sub-folders in those folders and return in
# **any order** the folders after removing.
#
# If a `folder[i]` is located within another `folder[j]`, it is called a sub-
# folder of it.
#
# The format of a path is one or more concatenated strings of the form: `/`
# followed by one or more lowercase English letters. For example, `/leetcode`
# and `/leetcode/problems` are valid paths while an empty string and `/` are
# not.
#
#
#
# **Example 1:**
#
#
#
#     Input: folder = ["/a","/a/b","/c/d","/c/d/e","/c/f"]
#     Output: ["/a","/c/d","/c/f"]
#     Explanation: Folders "/a/b/" is a subfolder of "/a" and "/c/d/e" is inside of folder "/c/d" in our filesystem.
#
#
# **Example 2:**
#
#
#
#     Input: folder = ["/a","/a/b/c","/a/b/d"]
#     Output: ["/a"]
#     Explanation: Folders "/a/b/c" and "/a/b/d/" will be removed because they are subfolders of "/a".
#
#
# **Example 3:**
#
#
#
#     Input: folder = ["/a/b/c","/a/b/ca","/a/b/d"]
#     Output: ["/a/b/c","/a/b/ca","/a/b/d"]
#
#
#
#
# **Constraints:**
#
#   * `1 <= folder.length <= 4 * 10^4`
#   * `2 <= folder[i].length <= 100`
#   * `folder[i]` contains only lowercase letters and '/'
#   * `folder[i]` always starts with character '/'
#   * Each folder name is unique.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

test_suffix = parsed_args["test"] ? "_test" : ""
input = readlines("input/day_16$test_suffix.txt")

mutable struct BitBuffer
    nibbles::Vector{UInt8}
    n::UInt16
    bits_left::UInt8
end

struct Packet
    version:: UInt8
    type:: UInt8
    bitlen::UInt32
    value::Union{UInt64, Vector{Packet}}
end

BitBuffer(nibbles) = BitBuffer(nibbles, 1, 4)

function bitmask(bitcount)
    if bitcount >= 32
        throw(AssertionError("too many bits"))
    end
    (UInt32(1) << bitcount) - 1
end

function get_bits!(buf, bits)
    result::UInt32 = 0
    if buf.bits_left > bits
        shft = buf.bits_left - bits
        result = (buf.nibbles[buf.n] >> shft) & bitmask(bits)
        buf.bits_left = shft
    else
        result = buf.nibbles[buf.n] & bitmask(buf.bits_left)
        bits -= buf.bits_left
        n = buf.n + 1
        while bits >= 4
            result = (result << 4) | buf.nibbles[n]
            bits -= 4
            n += 1
        end
        buf.n = n
        buf.bits_left = 4 - bits
        if bits > 0
            result = (result << bits) | (buf.nibbles[n] >> buf.bits_left)
        end
    end
    result
end

function version_type!(buf)
    version = get_bits!(buf, 3)
    type = get_bits!(buf, 3)
    (version, type)
end

function get_literal!(buf)
    value::UInt64 = 0
    bitlen::UInt32 = 0
    while true
        chunk = get_bits!(buf, 5)
        bitlen += 5
        value = (value << 4) | (chunk & 0x0f)
        if (chunk & 0x10) == 0
            break;
        end
    end
    (bitlen, value)
end

function parse_packet!(buf)
    (v, t) = version_type!(buf)
    if t==4 # literal
        (bitlen, value) = get_literal!(buf)
        packet = Packet(v, t, bitlen + 6, value)
    else
        length_type = get_bits!(buf, 1) # length 6 + 1 = 7
        if length_type == 0
            expected_bits = get_bits!(buf, 15) # length 7 + 15 = 22
            packets = Vector{Packet}()
            bits = 0
            while bits < expected_bits
                p = parse_packet!(buf)
                push!(packets, p)
                bits += p.bitlen
            end
            packet = Packet(v, t, expected_bits + 22, packets)
        else
            packet_count = get_bits!(buf, 11) # length 7 + 11 = 18
            packets = Vector{Packet}()
            while packet_count > 0
                p = parse_packet!(buf)
                push!(packets, p)
                packet_count -= 1
            end
            packet_bits = sum([p.bitlen for p in packets])
            packet = Packet(v, t, packet_bits + 18, packets)
        end
    end
    packet
end

function parse_outer_packet!(buf)
    packet = parse_packet!(buf)
    # tidy up the end of the parse_packet
    if buf.bits_left < 4
        buf.n += 1
        buf.bits_left = 4
    end
    packet
end

function add_versions(packet)
    v::UInt32 = packet.version
    if packet.type != 4
        v += sum([add_versions(p) for p in packet.value])
    end
    v
end

buffer_from_string(str) = BitBuffer([parse(UInt8, c; base=16) for c in str])

function dump_packets(p)
    function do_dump(p, indent)
        println("$(indent)Packet $(p.version), $(p.type)")
        new_ident = "$indent    "
        if p.type == 4
            println("$(new_ident)Value: $(p.value)")
        else
            foreach(ps -> do_dump(ps, new_ident), p.value)
        end
    end
    do_dump(p,"")
end

function evaluate(p)
    value::UInt64 = 0
    if p.type == 4
        value = p.value
    else
        subvalues = [evaluate(sp) for sp in p.value]
        if p.type == 0
            value = sum(subvalues)
        elseif p.type == 1
            value = prod(subvalues)
        elseif p.type == 2
            value = minimum(subvalues)
        elseif p.type == 3
            value = maximum(subvalues)
        elseif p.type == 5
            value = subvalues[1] > subvalues[2] ? 1 : 0
        elseif p.type == 6
            value = subvalues[1] < subvalues[2] ? 1 : 0
        elseif p.type == 7
            value = subvalues[1] == subvalues[2] ? 1 : 0
        end
    end
    value
end

function apply(operator)
    function apply_once(in)
        p = parse_outer_packet!(buffer_from_string(in))
        pre = in[1:6]
        v = operator(p)
        println("$pre, $v")
    end
    foreach(apply_once, input)
end

function part1()
    apply(add_versions)
    "above"
end

function part2()
    apply(evaluate)
    "above"
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
## Element type
struct IIDMVoltMeasure <: AbstractNodeLabel
  V::Complex    # p.u.
  S::Complex       # MW + im*MVar
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
###############################
#
# Retain SVD components
#
###############################

"""
Retain eigenvectors that represent up to `k` percent of power
"""
function retain_svd{T <: AbstractFloat}(A::Array{T, 2}, k::T=0.9)

    ss = svdfact(A)

    pw = cumsum(ss.S ./ sum(ss.S))

    keep = maximum(find(pw .< k)) + 1

    ss = ss.U[:, 1:keep]

    return ss, pw[keep]
end

###############################
#
# Average the epochs
#
###############################

"""
    reduce_epochs(Array{T, 3}, new_num_epochs::Int=30)

Average epochs down to specified number of epochs.
"""
function reduce_epochs{T <: AbstractFloat}(a::Array{T, 3}, new_num_epochs::Int=30)
    if new_num_epochs < size(a, 2)
        ep_per_av = floor(Int, size(a, 2) / new_num_epochs)
        new = zeros(size(a, 1), new_num_epochs, size(a, 3))
        for i in 1:new_num_epochs-1
            new[:, i, :] = mean(a[:, 1+((i-1)*ep_per_av):ep_per_av+((i-1)*ep_per_av), :], 2)
        end
        new[:, new_num_epochs, :] = mean(a[:, 1+((new_num_epochs-1)*ep_per_av):end, :], 2)
        return new
    else
        return a
    end
end

function reduce_epochs(s::SSR,  new_num_epochs::Int=30)

    s.processing["epochs"] = reduce_epochs(s.processing["epochs"], new_num_epochs)
    return s
end

###############################
#
# Correct midline for bilateral
#
###############################

"""
Correct midline estimates from bilateral beamformer

Region supression for bilateral beamformers supresses itself along the midline. This causes innacurate estimates
around 0 on the x axis. To compensate for this take the average of adjacent valid locations.

#### Input

* vi: Volume image
* pm: Locations plus or minus the midline in mm to correct
"""
function correct_midline(v::VolumeImage; pm::Real=0.005, units=Meter)

    Logging.info("Correcting volume image midline errors caused by bilateral region supression")
    Logging.debug("Supressing region ± $(pm * units)")

    x = AbstractFloat[xi / (1 * units) for xi in v.x]
    midline_idxs = abs(x)
    midline_idxs = midline_idxs .<= pm

    valid_idxs = falses(size(midline_idxs))
    valid_idxs[minimum(find(midline_idxs))-1] = true
    valid_idxs[maximum(find(midline_idxs))+1] = true

    new_val = mean(v.data[valid_idxs, :, :], 1)
    for i in find(midline_idxs)
        v.data[i, :, :] = new_val
    end

    return v
end


text:
include("/home/alkan/c/nese/jtb/yeast_dyn.jl")

## ---- find number of coop nodes in a state ---
function ncoop (s,intr)
    n = length(s)
    act = intr*diagm(s) # keep only the columns for active nodes
    nc=0
    for (i=1:n)
#        if (abs(sum(act[i,:])) == sum(abs(act[i,:])))
        if (!in(1,act[i,:]) || !in(-1,act[i,:]))
            nc += 1
        end
    end
    return nc
end

## ---- find coop coef for ALL point attractors of a network ---
function cooppointattr (intr,slf)
    listt = getpointattr(intr,slf)[:,2:end] # first in the list is the zero attr
    ncop=Int64[];
    for (i=1:length(listt[1,:]))
        push!(ncop,ncoop(listt[:,i],intr))
    end
    if (length(ncop)==0)  return -1,-1  # no fixed point other than (0,..0)
    else return (mean(ncop)/n,mean(listt))
    end
end

## ---- find coop coef for sampled cycles of a network ---
function coopattr (intr,slf,ninit)
    n = length(slf)
    (listc,listb,listlen) = getcyclessample(intr,slf,ninit)
    listcop=Float64[];
    listact=Float64[];
    npointattr = 0;
    for (i=1:length(listc[1,:]))        # go over each attractor
        s = listc[:,i]
        if (sum(s) != 0)
            ncp = ncoop(s,intr)
            nac = sum(abs(intr)*s)
            period = 1
            s=upd(s,intr,slf,n)
            while (s != listc[:,i])
                ncp += ncoop(s,intr)
                nac += sum(abs(intr)*s)
                period += 1
                s=upd(s,intr,slf,n)
            end
            push!(listcop,ncp/period)
            push!(listact,nac/period)
            if (period==1) npointattr += 1 end
        end
    end
    if (length(listcop)==0)  return -1,-1,-1,-1 # no fixed point other than (0,..0)
    else                    # also returns the fracion of active nodes
        return (1.0*mean(listcop)/n,mean(listact)/n,length(listcop),npointattr)
    end
end

## ---- find coop coef for sampled cycles weighted by basin size ---
function coopbasinweighted (intr,slf,ninit)
    n = length(slf)
    (listc,listb,listlen) = getcyclessample(intr,slf,ninit)
    listcop=Float64[];
    listact=Float64[];
    totalweight=0;
    for (i=1:length(listc[1,:]))
        s = listc[:,i]
        if (sum(s) != 0)                # exclude the zero attractor
            ncp = ncoop(s,intr)
            nac = sum(abs(intr)*s)
            period = 1
            s=upd(s,intr,slf,n)
            while (s != listc[:,i])
                ncp += ncoop(s,intr)
                nac += sum(abs(intr)*s)
                period += 1
                s=upd(s,intr,slf,n)
            end
            push!(listcop,listb[i]*ncp/period) # weigh with basin weight
            push!(listact,listb[i]*nac/period)
            totalweight += listb[i]            # excludes the zero attractor
        end
    end
    if (length(listcop)==0)  return -1,-1 # no fixed point other than (0,..0)
    else  return (sum(listcop)/(n*totalweight),sum(listact)/(n*totalweight))
    end
end

## ---- find coop coef for ALL cycles of a network ---
function coopattrall (intr,slf,ninit)
    n = length(slf)
    (listc,listb) = getcycles(intr,slf)
    listcop=Float64[];
    listact=Float64[];
    listlen=Float64[];
    for (i=1:length(listc[1,:]))        # go over each attractor
        s = listc[:,i]
        if (sum(s) != 0)
            ncp = ncoop(s,intr)
            nac = sum(abs(intr)*s)
            period = 1
            s=upd(s,intr,slf,n)
            while (s != listc[:,i])
                ncp += ncoop(s,intr)
                nac += sum(abs(intr)*s)
                period += 1
                s=upd(s,intr,slf,n)
            end
            push!(listcop,ncp/period)
            push!(listact,nac/period)
            push!(listlen,period)
        end
    end
    if (length(listcop)==0)  return -1,-1 # no fixed point other than (0,..0)
    else                    # also returns the fracion of active nodes
        return (1.0*mean(listcop)/n,mean(listact)/n)
    end
end

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
module NearestNeighborsBallTree
	using Base.Test
	using NearestNeighbors

	X = readcsv("test/iris.csv")

	t = BallTree(X)

	v = X[:, 21]

	k_nearest(v, t, 1)
	k_nearest(v, t, 2)
	k_nearest(v, t, 3)

	k_nearest(v, t, 1, 21)
	k_nearest(v, t, 2, 21)
	k_nearest(v, t, 3, 21)

	# inball(v, t, 0.5)
	# inball(v, t, 0.5, 21)

	# NearestNeighbors.walk(t)
	@test NearestNeighbors.check(t)
end


text:
using SPH
using Test, Plots

default(;color = :thermal, aspect_ratio = 1, label = "", xaxis=("x", (-0.1, 1.1)), yaxis = ("y", (-0.1, 1.1)))

function M6_hardcoded(q)
    if q < 0.0
        res = 0.0
    elseif q < 1.0
        res = (3 - q)^5 - 6*(2 - q)^5 + 15*(1 - q)^5
    elseif q < 2.0
        res = (3 - q)^5 - 6*(2 - q)^5
    elseif q < 3.0
        res = (3 - q)^5
    else
        res = 0.0
    end
    return res
end

@testset "SPH.jl" begin
   xs = rand(1000) * 6 .- 1
   fs_hardcoded = M6_hardcoded.(xs) .* SPH.σ6[1]
   fs_procedural = M6[1].(xs)
   @test isapprox.(fs_hardcoded, fs_procedural, atol=1e-12) |> all
   #@show xs[a], fs_hardcoded[a], fs_procedural[a]
end

@testset "Density calculations" begin

    particles = SPH.particles_in_box(1)
    @test SPH.compute_distances(particles) ≈ [0.0]

    particles = [SPH.Particle(0.0, 0.0), SPH.Particle(0.0, 1.0)]
    @test SPH.compute_distances(particles) ≈ [0.0 1.0; 1.0 0.0]

    push!(particles, SPH.Particle(1.0, 1.0))
    @test SPH.compute_distances(particles) ≈ [0.0 1.0 √2; 1.0 0.0 1.0; √2 1.0 0.0]

end

N = 100
particles = SPH.particles_in_box(N)

NEdge = floor(Int, sqrt(N))
xs = LinRange(0, 1, NEdge)
ys = LinRange(0, 1, NEdge)

#left = SPH.Particle.(zeros(NEdge), ys)
#bottom = SPH.Particle.(xs, zeros(NEdge))
#top = SPH.Particle.(xs, ones(NEdge))
#right = SPH.Particle.(ones(NEdge), ys)

#append!(particles, left)
#append!(particles, right)
#append!(particles, top)
#append!(particles, bottom)

timestep = 0.01
anim = @animate for i in 1:500
    global particles = SPH.update(particles, dt=timestep);
    if i % 2 == 0
        SPH.plot_particles(particles, clims = (0, 120))
    end
end
gif(anim, "test.gif", fps = 25)


text:
"""
    elapsed_time(res, i)

Elapsed time in seconds for i-th picture
"""
elapsed_time(res::ExpImgsResults, i)::Float64 =
    (Dates.Second(res.time[i] - res.exi.experiment_start)).value

"""
    get_time_series(res::ExpImgsResults, tt)
    get_time_series(res::ExpImgsResults)

tt in seconds since ex.experiment_start.  If no tt then return at times.

Get time series data as matrix:
- mean diameter
- std of diameter
- median diameter
- 10% quantile
- 90% quantile
- scalloping
- symmetry

Example

    dia_mean, dia_std, dia_med, dia_10, dia_90, scalloping, sym = get_time_series(res, tt)

"""
function get_time_series(res::ExpImgsResults, tt)
    @unpack dia_mean, dia_std, dia_quant, scalloping, tb_cor, t = res
    vars = [dia_mean, dia_std, dia_quant[2,:], dia_quant[1,:], dia_quant[3,:], scalloping, tb_cor]
    outs = eltype(vars)[]
    for v in vars
        ipt = IP.interpolate((t,), v, IP.Gridded(IP.Linear()))
        push!(outs, ipt[tt])
    end
    return outs
end
function get_time_series(res::ExpImgsResults)
    @unpack dia_mean, dia_std, dia_quant, scalloping, tb_cor, t = res
    return [dia_mean, dia_std, dia_quant[2,:], dia_quant[1,:], dia_quant[3,:], scalloping, tb_cor]
end

#################
# findpeaks
# put into VAWTools

"""
    findpeaks(vec; min_dominance=0.0, min_prominence=0.0)

Return

- peaks

- index of peaks

- prominence (difference between peak height and the highest low-point to
  traverse to get to a higher point)

- dominance (separation to next higher point, or maximal distance to the edge)

TODO: use troughs vector
"""
function findpeaks(vec; min_dominance=0.0, min_prominence=0.0)
    ## find all peaks and all troughs
    peaks = eltype(vec)[]
    ips = Int[]
    troughs = eltype(vec)[]
    its = Int[]
    for i=2:length(vec)-1
        if vec[i-1]<vec[i] && vec[i+1]<=vec[i]
            push!(peaks, vec[i])
            push!(ips, i)
        elseif vec[i-1]>vec[i] && vec[i+1]>=vec[i]
            push!(troughs, vec[i])
            push!(its, i)
        end
    end

    # calculate dominance and prominence
    doms = Int[]
    proms = eltype(vec)[]
    for (i,ip) in enumerate(ips)
        ele = peaks[i]
        left = @view vec[ip-1:-1:1]
        right = @view vec[ip+1:end]
        dom1 = findfirst(x->x>ele, left)
        dom2 = findfirst(x->x>ele, right)
        if dom1==0 && dom2==0
            # highest peak
            max_dist_to_margin = max(ip, length(vec)-ip)
            push!(doms, max_dist_to_margin)
            # push!(doms, length(vec))
            push!(proms, ele-minimum(vec))
        else
            if dom1==0 # no higher peak on left
                push!(doms,dom2)
                #min1 = minimum(@view vec[ip-1:-1:1])
                min2 = minimum(@view vec[ip+1:ip+dom2])
                min1 = min2
            elseif dom2==0 # no higher peak on right
                push!(doms,dom1)
                min1 = minimum(@view vec[ip-1:-1:ip-dom1])
                #min2 = minimum(@view vec[ip+1:end])
                min2 = min1
            else
                push!(doms,min(dom1,dom2))
                min1 = minimum(@view vec[ip-1:-1:ip-dom1])
                min2 = minimum(@view vec[ip+1:ip+dom2])
            end
            push!(proms, ele - max(min1,min2))
        end
    end

    inds = (doms.>=min_dominance) .& (proms.>=min_prominence)
    return peaks[inds], ips[inds], proms[inds], doms[inds]
end

"""
    findpeaks_onesided(vec, dir::Symbol; min_dominance=0.0, min_prominence=0.0)

As findpeaks but only scans prominence and dominance in direction `dir` (:left,:right)

Return

- peaks

- index of peaks

- prominence (difference between peak height and lowest point to
  traverse to get to a higher point)

- dominance (separation to next higher point, or minimal distance to the edge)

The API inspired by
https://ch.mathworks.com/help/signal/ref/findpeaks.html.
"""
function findpeaks_onesided(vec, dir; min_dominance=0.0, min_prominence=0.0)
    @assert dir in [:left,:right]
    ## find all peaks and all troughs
    peaks = eltype(vec)[]
    ips = Int[]
    troughs = eltype(vec)[]
    its = Int[]
    for i=2:length(vec)-1
        if vec[i-1]<vec[i] && vec[i+1]<=vec[i]
            push!(peaks, vec[i])
            push!(ips, i)
        elseif vec[i-1]>vec[i] && vec[i+1]>=vec[i]
            push!(troughs, vec[i])
            push!(its, i)
        end
    end

    # calculate dominance and prominence
    doms = Int[]
    proms = eltype(vec)[]
    for (i,ip) in enumerate(ips)
        ele = peaks[i]
        side = dir==:left ? @view(vec[ip-1:-1:1]) : @view(vec[ip+1:end])
        dom = findfirst(x->x>ele, side)
        if dom==0
            # highest peak to border
            push!(doms, length(side))
            push!(proms, ele-minimum(side))
        else
            push!(doms,dom)
            side2 = dir==:left ? @view(vec[ip-1:-1:ip-dom]) : @view(vec[ip+1:ip+dom])
            push!(proms, ele - minimum(side2))
        end
    end

    inds = (doms.>=min_dominance) .& (proms.>=min_prominence)
    return peaks[inds], ips[inds], proms[inds], doms[inds]
end

find_first_big_peak(vec, dir=:left) = find_first_big_peak(findpeaks_onesided(vec, dir==:left ? :right : :left)[1:3]..., dir)
function find_first_big_peak(peaks, locs, proms, dir=:left)
    # run with a median filter over it
end

## This does not work well:
# using OnlineStats
# function find_outlier_median(vec, quant=0.1, exp_weight=1/4,
#                              what=[:below,:above,:both][3];
#                              verbose=false)
#     @assert 0<exp_weight<=1
#     # mean for the first 1/0.2 observations, then Exponential.
#     # e-fold of importance given by exp_weight

#     w = ExponentialWeight(exp_weight)

#     # none of them seem quite kosher:
#     Qu = [QuantileMM, QuantileSGD, QuantileMSPI][3]

#     s = Series(w, Qu([quant, 0.5, 1-quant]))

#     outliers = Int[]
#     q1 = Float64[]
#     q2 = Float64[]
#     q3 = Float64[]

#     fit!(s, vec[1], 1.0) # use 1.0 weight, see https://github.com/joshday/OnlineStats.jl/issues/93
#     for i = 2:length(vec)
#         v = vec[i]
#         sv = value(s)[1]
#         push!(q1,sv[1])
#         push!(q2,sv[2])
#         push!(q3,sv[3])
#         if what==:both
#             sv[1]<v<sv[3] || push!(outliers, i)
#         elseif what==:below
#             sv[1]<v || push!(outliers, i)
#         elseif what==:above
#             v<sv[3] || push!(outliers, i)
#         end
#         fit!(s, v)
#     end
#     sv = value(s)[1]
#     push!(q1,sv[1])
#     push!(q2,sv[2])
#     push!(q3,sv[3])

#     if verbose
#         P.figure()
#         P.plot(1:length(vec),vec)
#         map(x->P.plot(2:length(vec)+1,x) ,(q1,q2,q3))
#         P.plot(outliers, vec[outliers], ".")
#     end

#     w = Bounded(EqualWeight(), exp_weight)
#         w = ExponentialWeight(exp_weight)
#     @show Series(vec, w, QuantileMSPI([quant, 0.5, 1-quant]))
#     outliers
# end

"""
    find_spikes_quantile(vec, comp_vec_q=>, quant=0.9,
                          quant_offset=mean(vec)/2, quant_factor=1,
                          pre_window=5, post_window = 2;
                          verbose=false)

Find index of spikes/outliers by comparing each value with a certain quantile
of a running window.

    vec -- input vector
    comp_vec_q=> -- comparison to use: > means remove outliers above
    quant=0.9 -- what quantile to compare against
    quant_offset=median(vec)/2
    quant_factor=1
    pre_window=5, post_window=0 -- window size over which to look at values

The comparison is then `comp_vec_q(v[i], q*quant_factor + quant_offset)`, i.e.
the quantile value `q` is multiplied by the `quant_factor` and `quant_offset` is added.

The window for the i-th value is `i-pre_window:i+post_window-1`, i.e. if
`post_window==0` then the value itself is not considered.

Keyword:
    verbose=false -- if true plot

Return

- index of spikes/outliers
- value of median over window at outliers

"""
function find_spikes_quantile(vec, comp_vec_q= >,
                               quant=0.9, quant_offset=median(vec), quant_factor=1,
                               pre_window=length(vec)÷10, post_window = 0;
                              verbose=false)
    ex = 3
    pre_w = 1./(pre_window+1:-1:2)
    post_w = 1./(1:post_window).^ex
    w = StatsBase.Weights([pre_w; post_w])
    outliers = Int[]
    qs = Float64[]
    med = Float64[]
    for i = pre_window+1:length(vec)-post_window
        vi = vec[i]
        v = @view vec[i-pre_window:i+post_window-1]
        if length(v)>0
            q = quant_factor * quantile(v, w, quant) + quant_offset
            push!(qs, q)
            push!(med, median(v))
            if comp_vec_q(vi, q)
                push!(outliers, i)
            end
        end
    end
    if verbose
        P.plot(1:length(vec),vec)
        map(x->P.plot(pre_window+1:length(vec)-post_window,x) ,(qs,))
        P.plot(outliers, vec[outliers], ".")
    end
    outliers, med[outliers-pre_window]
end

#############
# Saving and Loading
function filename(exi::ExpImgs)
    @unpack dir, thin_num, ns, algo = exi
    "$(dir)-thin$(thin_num)-ns$(length(ns))-$(algo)"
end
filename(res::ExpImgsResults) = filename(res.exi)
"""
    save_result(res::ExpImgsResults; overwrite=false, store_imgs=false)

"""
function save_result(res::ExpImgsResults; overwrite=false, store_imgs=false)
    fln = "output/"*filename(res)*".jld"
    if isfile(fln) && !overwrite
        warn("File exists and overwrite==false, not saving!")
        return nothing
    end
    if store_imgs==false && length(res.imgs)>0
        res = ExpImgsResults(res,
                             imgs=[])
    end
    JLD.jldopen(fln, "w") do file
        file["res"] = res
    end
    nothing
end
load_result(fln) = JLD.load(fln)
load_result(exi::ExpImgs) = load_result(filename(exi)*".jld")

load_legacy(fln, ex) = load_legacy(fln, ex.experiment_start, ex.experiment_running)
function load_legacy(fln, experiment_start, experiment_running)
    @show experiment_start, experiment_running
    res = type2dict(JLD.load(fln)["res"])
    ep = type2dict(pop!(res, :ep))
    ep[:p2m] = RcImages.Pixel2Meter(;type2dict(ep[:p2m])...)
    ep[:extras] = Dict{Symbol,Any}()
    ep[:extras][:fln] = fln
    exi = ExpImgs(;ep...)
    res[:exi] = exi
    res[:extras] = Dict{Symbol,Any}()
    res[:extras][:fln] = fln
    res[:time] = Dates.Time.(pop!(res, :capture_times))
    res[:dia_std] = vec(std(res[:ts_dist]+res[:bs_dist],1))
    # figure out times of flowing water
    res[:t] = time2seconds(res[:time], experiment_start)
    ExpImgsResults(;res...)
end


text:
using SimLynx
using Test

@testset "SimLynx.jl" begin
    @test greet() == "Hello World!"
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

@enum CONTEXT_TYPE CPU=1 GPU=2 CPU_PINNED=3

Base.convert(::Type{CONTEXT_TYPE}, x::Integer) = CONTEXT_TYPE(x)

"""
    Context(dev_type, dev_id)

A context describes the device type and id on which computation should be carried on.
"""
struct Context
  device_type::CONTEXT_TYPE
  device_id::Int

  Context(dev_type::CONTEXT_TYPE, dev_id::Integer = 0) = new(dev_type, dev_id)
end

const _default_ctx = Ref{Context}(Context(CPU, 0))

Context(dev_type::Integer, dev_id::Integer = 0) =
  Context(convert(CONTEXT_TYPE, dev_type), dev_id)

Base.show(io::IO, ctx::Context) =
  print(io, lowercase("$(ctx.device_type)$(ctx.device_id)"))

function _with_context(dev_type::Union{Symbol,Expr}, dev_id, e::Expr)
  global _default_ctx
  quote
    ctx = current_context()
    ctx′ = Context($(esc(dev_type)), $(esc(dev_id)))
    $_default_ctx[] = ctx′
    try
      return $(esc(e))
    finally
      $_default_ctx[] = ctx
    end
  end
end

"""
    @context device_type [device_id] expr

Change the default context in the following expression.

# Examples
```jl-repl
julia> mx.@context mx.GPU begin
         mx.zeros(2, 3)
       end
2×3 NDArray{Float32,2} @ gpu0:
 0.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  0.0f0

julia> @context mx.GPU mx.zeros(3, 2)
3×2 NDArray{Float32,2} @ gpu0:
 0.0f0  0.0f0
 0.0f0  0.0f0
 0.0f0  0.0f0
```
"""
macro context(dev_type, e::Expr)
  _with_context(dev_type, 0, e)
end

macro context(dev_type, dev_id, e::Expr)
  _with_context(dev_type, dev_id, e)
end

for dev ∈ [:cpu, :gpu]
  ctx = QuoteNode(Symbol(uppercase(string(dev))))
  docstring = """
        @$dev [device_id] expr

    A shorthand for `@context mx.GPU`.

    # Examples
    ```jl-repl
    julia> mx.@with_gpu mx.zeros(2, 3)
    2×3 NDArray{Float32,2} @ gpu0:
     0.0f0  0.0f0  0.0f0
     0.0f0  0.0f0  0.0f0
    ```
    """
  @eval begin
    @doc $docstring ->
    macro $dev(e::Expr)
      ctx = $ctx
      quote
        @context $ctx $(esc(e))
      end
    end

    macro $dev(dev_id, e::Expr)
      ctx = $ctx
      quote
        @context $ctx $(esc(dev_id)) $(esc(e))
      end
    end
  end
end  # for dev ∈ [:cpu, :gpu]

"""
    cpu(dev_id)

Get a CPU context with a specific id. `cpu()` is usually the default context for many
operations when no context is specified.

# Arguments
* `dev_id::Integer = 0`: the CPU id.
"""
cpu(dev_id::Integer = 0) = Context(CPU, dev_id)

"""
    gpu(dev_id)

Get a GPU context with a specific id. The K GPUs on a node is typically numbered as 0,...,K-1.

# Arguments
* `dev_id::Integer = 0` the GPU device id.
"""
gpu(dev_id::Integer = 0) = Context(GPU, dev_id)

"""
    num_gpus()

Query CUDA for the number of GPUs present.
"""
function num_gpus()
  n = Ref{Cint}()
  @mxcall :MXGetGPUCount (Ref{Cint},) n
  n[]
end

"""
    empty_cache(ctx::Context = current_context())

Empties the memory cache for the current contexts device.
MXNet utilizes a memory pool to avoid excessive allocations.
Calling empty_cache will empty the memory pool of the contexts
device. This will only free the memory of the unreferenced data.
"""
function empty_cache(ctx::Context = current_context())
  @mxcall :MXStorageEmptyCache (Cint, Cint) ctx.device_type ctx.device_id
  ctx
end

"""
    gpu_memory_info(dev_id = 0)::Tuple{UInt64,UInt64}

Query CUDA for the free and total bytes of GPU global memory.
It returns a tuple of `(free memory, total memory)`.

```julia-repl
julia> mx.gpu_memory_info()
(0x00000003af240000, 0x00000003f9440000)
```
"""
function gpu_memory_info(dev_id = 0)
  free = Ref{UInt64}()
  n = Ref{UInt64}()
  @mxcall :MXGetGPUMemoryInformation64 (Cint, Ref{UInt64}, Ref{UInt64}) dev_id free n
  free[], n[]
end

"""
    current_context()

Return the current context.

By default, `mx.cpu()` is used for all the computations
and it can be overridden by using the `@context` macro.

# Examples
```jl-repl
julia> mx.current_context()
cpu0

julia> mx.@context mx.GPU 1 begin  # Context changed in the following code block
         mx.current_context()
       end
gpu1

julia> mx.current_context()
cpu0
```
"""
current_context() = _default_ctx[]


text:
using Revise
using MendelIHT
using SnpArrays
using Random
using GLM
using DelimitedFiles
using Distributions
using LinearAlgebra
using CSV
using DataFrames
using StatsBase
BLAS.set_num_threads(1)

function simulate_beta(r::Int, p::Int)
    k = [10^i for i in 1:r]
    d = Normal(0, 0.1)
    B = zeros(r, p)
    for i in 1:r
        k = 10^i
        B[i, 1:k] .= rand(d, k)
        shuffle!(@view(B[i, :]))
    end
    return B
end

# function simulate_genotypes(r::Int; seed::Int=2021)
#     # set seed
#     Random.seed!(seed)
#     n, p = 10000, 30000

#     # simulate `.bed` file with no missing data
#     x = simulate_random_snparray(undef, n, p)
#     xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)

#     # intercept is the only nongenetic covariate
#     Z = ones(n, 1)
#     intercepts = zeros(r)' # each trait have 0 intercept

#     # simulate β
#     B = simulate_beta(r, p)
#     writedlm("trueb.txt", B)

#     # between trait covariance matrix
#     Σ = random_covariance_matrix(r)
#     writedlm("true_cov.txt", Σ)

#     # between sample covariance is identity
#     V = Matrix(I, n, n)

#     # simulate using Distributions.jl
#     μ = Z * intercepts + xla * B'
#     Y = rand(MatrixNormal(μ', Σ, V))
#     standardize!(Transpose(Y))

#     return xla, Matrix(Z'), B, Σ, Y
# end

function simulate_genotypes(plinkfile::String, r::Int; seed::Int=2021)
    # set seed
    Random.seed!(seed)

    # simulate `.bed` file with no missing data
    x = SnpArray(plinkfile)
    xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
    n, p = size(x)

    # intercept is the only nongenetic covariate
    Z = ones(n, 1)
    intercepts = zeros(r)' # each trait have 0 intercept

    # simulate β
    B = simulate_beta(r, p)
    writedlm("trueb.txt", B)

    # between trait covariance matrix
    Σ = random_covariance_matrix(r)
    writedlm("true_cov.txt", Σ)

    # between sample covariance is identity
    V = Matrix(I, n, n)

    # simulate using Distributions.jl
    μ = Z * intercepts + xla * B'
    Y = rand(MatrixNormal(μ', Σ, V))
    standardize!(Transpose(Y))

    return xla, Matrix(Z'), B, Σ, Y
end

"""
Simulate 3 phenotypes of different polygenecity
y1: 10 causal SNPs
y2: 100 causal SNPs
y3: 1000 causal SNPs
"""
function one_simulation(seed::Int)
    isdir("sim$seed") || mkdir("sim$seed")
    cd("sim$seed")
    r = 3 # number of phenotypes
    plinkname = "/scratch/users/bbchu/fastphase/10k_K20/ukb.10k.chr10.bed"

    # simulate data
    Random.seed!(seed)
    xla, Z, B, Σ, Y = simulate_genotypes(plinkname, r, seed=seed)
    GC.gc()

    # B[1, findall(!iszero, B[1, :])]
    # B[2, findall(!iszero, B[2, :])]
    # B[3, findall(!iszero, B[3, :])]

    # run multivariate IHT
    path=200:200:2000
    mIHT_time = @elapsed begin
        mses = cv_iht(Y, Transpose(xla), Z, path=path)
        GC.gc()
        k_rough_guess = path[argmin(mses)]
        path = (k_rough_guess - 150):50:(k_rough_guess + 150)
        mses = cv_iht(Y, Transpose(xla), Z, path=path)
        GC.gc()
        k_rough_guess = path[argmin(mses)]
        path = (k_rough_guess - 40):10:(k_rough_guess + 40)
        mses = cv_iht(Y, Transpose(xla), Z, path=path)
        GC.gc()
        k_rough_guess = path[argmin(mses)]
        path = (k_rough_guess - 9):(k_rough_guess + 9)
        mses = cv_iht(Y, Transpose(xla), Z, path=path)
        GC.gc()
        iht_result = fit_iht(Y, Transpose(xla), Z, k=path[argmin(mses)])
        writedlm("iht_beta.txt", iht_result.beta)
    end
    println("mIHT run time = $mIHT_time")

    # check power and FDR
    power, FDR = Float64[], Float64[]
    β1 = findall(!iszero, iht_result.beta[1, :])
    β2 = findall(!iszero, iht_result.beta[2, :])
    β3 = findall(!iszero, iht_result.beta[3, :])
    β1_true = findall(!iszero, B[1, :])
    β2_true = findall(!iszero, B[2, :])
    β3_true = findall(!iszero, B[3, :])
    # power
    push!(power, length(β1 ∩ β1_true) / 10)
    push!(power, length(β2 ∩ β2_true) / 100)
    push!(power, length(β3 ∩ β3_true) / 1000)
    # fdr
    push!(FDR, (count(!iszero, β1) - length(β1 ∩ β1_true)) / max(1, count(!iszero, β1)))
    push!(FDR, (count(!iszero, β2) - length(β2 ∩ β2_true)) / max(1, count(!iszero, β2)))
    push!(FDR, (count(!iszero, β3) - length(β3 ∩ β3_true)) / max(1, count(!iszero, β3)))
    # overall power and FDR
    correct_snps = unique([x[2] for x in findall(!iszero, B)])
    discovered_snps = unique([x[2] for x in findall(!iszero, iht_result.beta)])
    overall_power = length(correct_snps ∩ discovered_snps) / length(correct_snps)
    overall_FDR = (count(!iszero, discovered_snps) - length(discovered_snps ∩ correct_snps)) /
        max(1, count(!iszero, discovered_snps))

    # save summary stats
    open("power_summary.txt", "w") do io
        println(io, "trait1_k10,trait2_k100,trait3_k1000,overall")
        println(io, power[1], ',', power[2], ',', power[3], ',', overall_power)
    end
    open("FDR_summary.txt", "w") do io
        println(io, "trait1_k10,trait2_k100,trait3_k1000,overall")
        println(io, FDR[1], ',', FDR[2], ',', FDR[3], ',', overall_FDR)
    end

    return nothing
end

seed = parse(Int, ARGS[1])
one_simulation(seed)


text:
@testset "Update Function Tests" begin

  @test isdefined(Julz, :update) == true

end


text:
code_0 = 0
code_1 = 1
l_a = [code_0, code_1]
code_a = "a"
code_b = "b"
l_b = [code_a, code_b]
function main()
    for i in l_a
        println(join([i], " "))
    end
    for i in l_b
        println(join([i], " "))
    end
    if "a" in ["a", "b"]
        println(join(["OK"], " "))
    end
end

main()