text:
module HighDimPDE
    using Reexport
    @reexport using DiffEqBase
    using Statistics
    using Flux, Zygote, LinearAlgebra
    using Functors
    # using ProgressMeter: @showprogress
    using CUDA
    using Random
    using SparseArrays

    abstract type HighDimPDEAlgorithm <: DiffEqBase.AbstractODEAlgorithm end

    """
        PIDEProblem(g, f, μ, σ, x, tspan, p = nothing, x=nothing, u_domain=nothing, neumann_bc=nothing)

    Defines a Partial Integro Differential Problem, of the form `du/dt = 1/2 Tr(\\sigma \\sigma^T) Δu(t,x) + μ ∇u(t,x) + \\int f(u,x) dx`; where f is a nonlinear Lipschitz function

    # Arguments
    * `g` : The initial condition g(x, p, t).
    * `f` : The function f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t)
    * `μ` : The drift function of X from Ito's Lemma μ(x, p, t)
    * `σ` : The noise function of X from Ito's Lemma σ(x, p, t)
    * `tspan`: The timespan of the problem.
    * `p`: the parameter
    * `x`: the point of the solution required
    * `u_domain` : if provided, approximating the solution on the hypercube `u_domain[1] × u_domain[2]`.
    * `neumann_bc`: if provided, neumann boundary conditions on the hypercube `neumann_bc[1] × neumann_bc[2]`.
    """
    struct PIDEProblem{uType,G,F,Mu,Sigma,xType,tType,P,UD,NBC,K} <: DiffEqBase.AbstractODEProblem{uType,tType,false}
        u0::uType
        g::G # initial condition
        f::F # nonlinear part
        μ::Mu
        σ::Sigma
        x::xType
        tspan::Tuple{tType,tType}
        p::P
        u_domain::UD # for DeepSplitting only
        neumann_bc::NBC # neumann boundary conditions
        kwargs::K
    end

    function PIDEProblem(g, f, μ, σ, tspan;
                                    p=nothing,
                                    x=nothing,
                                    u_domain=nothing,
                                    neumann_bc=nothing,
                                    kwargs...)

    @assert !isnothing(x) ⊻ !isnothing(u_domain) "Need to provide whether `x` or `u_domain`"
    !isnothing(u_domain) ? x = first(u_domain) : nothing # x is required for even when solving on u_domain
    @assert eltype(x) <: AbstractFloat "`x` should be a Float"
    @assert eltype(tspan) <: AbstractFloat "`tspan` should be a tuple of Float"
    @assert typeof(u_domain) <: Union{Nothing,Tuple} "type of `u_domain` can be whether `Nothing` or Tuple{AbstractVector, AbstractVector}"
    isnothing(u_domain) ? nothing : @assert eltype(eltype(u_domain)) <: AbstractFloat "`tspan` should be a tuple of Float"
    @assert typeof(neumann_bc) <: Union{Nothing,Tuple} "type of `neumann_bc` can be whether `Nothing` or Tuple{AbstractVector, AbstractVector}"
    isnothing(neumann_bc) ? nothing : @assert eltype(eltype(neumann_bc)) <: eltype(x)
    @assert eltype(g(x)) == eltype(x) "Type of `g(x)` must match type of x"
    @assert eltype(f(x, x, g(x), g(x), 0f0, 0f0, p, tspan[1])) == eltype(x) "Type of non linear function `f(x)` must type of x"

    PIDEProblem{typeof(g(x)),
                typeof(g),
                typeof(f),
                typeof(μ),
                typeof(σ),
                typeof(x),
                eltype(tspan),
                typeof(p),
                typeof(u_domain),
                typeof(neumann_bc),
                typeof(kwargs)}(
                g(x), g, f, μ, σ, x, tspan, p, u_domain, neumann_bc, kwargs)
    end

    Base.summary(prob::PIDEProblem) = string(nameof(typeof(prob)))

    function Base.show(io::IO, A::PIDEProblem)
    println(io, summary(A))
    print(io, "timespan: ")
    show(io, A.tspan)
    end

    include("MCSample.jl")
    include("reflect.jl")
    include("DeepSplitting.jl")
    include("MLP.jl")

    export PIDEProblem, DeepSplitting, MLP

    export NormalSampling, UniformSampling, NoSampling, solve
end


text:
module Elf
using Base
using StaticArrays, DataStructures

export iself,
    elfclass,
    endian,
    sections,
    segments,
    symbols,
    elfosabi,
    elftype

include("./constants.jl")
# include("./types.jl")

"""
    iself(bin) -> Bool

Verify ELF Header with Magic numbers

# Examples
```julia-repl
julia> iself(read("/path/to/binary"))
true
```

# Arguments
- `bin::Vector{UInt8}` : binary file to check.

"""
iself(ehdr::Elf64_Ehdr) = ehdr.e_ident[begin:4] == [ELFMAG1, ELFMAG2, ELFMAG3, ELFMAG4]

"""
    elfclass(ehdr::Elf64_Ehdr) -> EHClass

Detrmines ELF class(64bit or 32bit) if invalid class error will occur.
"""
function elfclass(ehdr::Elf64_Ehdr)
    try
        EHClass(ehdr.e_ident[EI_CLASS])
    catch
        error("Invalid class.")
    end
end

"""
    endian(ehdr::Elf64_Ehdr) -> EHEndian

Detrmines ELF binary endian. Returns enum represents endian. If invalid error will occur.
"""

function endian(ehdr::Elf64_Ehdr)
    try
        EHEndian(ehdr.e_ident[EI_DATA])
    catch
        error("Invalid endian.")
    end
end

"""
    elfversion(ehdr::Elf64_Ehdr) -> Symbol

If binary has valid ELF version, else error will occur.
"""

function elfversion(ehdr::Elf64_Ehdr)
    try
        ehdr.e_ident[EI_VERSION]
    catch
        error("Invalid ELF version.")
    end
end

function elfosabi(ehdr::Elf64_Ehdr)
    try
        (arch = EHOsAbi(ehdr.e_ident[EI_OSABI]), version = ehdr.e_ident[EI_OSABI])
    catch
        error("Invalid OS ABI")
    end
end

function elftype(ehdr::Elf64_Ehdr)
    try
        EHType(ehdr.e_type)
    catch
        error("Invalid ELF type.")
    end
end

"""
    sections(bin::Vector{UInt8}) -> OrderedDict{String,Elf64_Shdr}

Section header infomation. Returns dictoinary represents section name and its infomation ordered by section number.
"""
function sections(bin::Vector{UInt8})
    sections_info::OrderedDict{String,Elf64_Shdr} = OrderedDict()
    ehdr = Elf64_Ehdr(bin)
    shstr = Elf64_Shdr(bin, ehdr.e_shoff + ehdr.e_shentsize * ehdr.e_shstrndx)
    for i = 0:(ehdr.e_shnum - 1)
        shdr = Elf64_Shdr(bin, ehdr.e_shoff + ehdr.e_shentsize * i)
        sname = unsafe_string(pointer(bin, shstr.sh_offset + shdr.sh_name + 1))
        push!(sections_info, sname => shdr)
    end

    sections_info
end

"""
    segments(bin::Vector{UInt8} -> OrderedDict{String,Elf64_Phdr}

Program header infomation. Returns dictoinary represents segment name and its infomation ordered by segment number.
"""

function segments(bin::Vector{UInt8})
    ehdr = Elf64_Ehdr(bin)
    shstr = Elf64_Shdr(bin, ehdr.e_shoff + ehdr.e_shentsize * ehdr.e_shstrndx)
    segments_info::OrderedDict{String,Elf64_Phdr} = OrderedDict()

    for i = 0:(ehdr.e_phnum - 1)
        phdr = Elf64_Phdr(bin, ehdr.e_phoff + ehdr.e_phentsize * i)
        for j = 0:(ehdr.e_shnum - 1)
            shdr = Elf64_Shdr(bin, ehdr.e_shoff + ehdr.e_shentsize * j)
            size = shdr.sh_type != UInt32(SHT_NOBITS) ? shdr.sh_size : 0
            if shdr.sh_offset < phdr.p_offset
                continue
            end
            if shdr.sh_offset + size > phdr.p_offset + phdr.p_filesz
                continue
            end
            seg_name = unsafe_string(pointer(bin, shstr.sh_offset + shdr.sh_name + 1))
            push!(segments_info, seg_name => phdr)
        end
    end

    segments_info
end

"""
    symbols(bin::Vector{UInt8}) -> (OrderedDict{String,Elf64_Sym}, OrderedDict{String,Elf64_Rel})

ELF symbols in the file. Returns tuple of ELF symbols and relocatable symbols.

# Examples
```julia-repl
julia> symbols(read("/path/to/binary"))
(OrderedDict{String,Elf64_Sym}(...),OrderedDict{String,Elf64_Rel}())
```

# Arguments
- `bin::Vector{UInt8}`: ELF binary.
"""

function symbols(bin::Vector{UInt8})
    ehdr = Elf64_Ehdr(bin)
    symbols_info::OrderedDict{String,Elf64_Sym} = OrderedDict()
    str = sections(bin)[".strtab"]

    for i = 0:(ehdr.e_shnum - 1)
        shdr = Elf64_Shdr(bin, ehdr.e_shoff + ehdr.e_shentsize * i)
        if shdr.sh_type != UInt32(SHT_SYMTAB)
            continue
        end
        global sym = shdr
        for j = 0:(div(sym.sh_size, sym.sh_entsize) - 1)
            symp = Elf64_Sym(bin, sym.sh_offset + sym.sh_entsize * j)
            if symp.st_name == 0
                continue
            end
            push!(
                symbols_info,
                unsafe_string(pointer(bin, str.sh_offset + symp.st_name + 1)) => symp,
            )
        end
    end

    symbols_rel_info::OrderedDict{String,Elf64_Rel} = OrderedDict()

    for i = 0:(ehdr.e_shnum - 1)
        shdr = Elf64_Shdr(bin, ehdr.e_shoff + ehdr.e_shentsize * i)
        if shdr.sh_type != UInt32(SHT_REL) && shdr.sh_type != UInt32(SHT_RELA)
            continue
        end
        rel = shdr
        for j = 0:(div(rel.sh_size, rel.sh_entsize) - 1)
            relp = Elf64_Rel(bin, rel.sh_offset + rel.sh_entsize * j)
            symp =
                Elf64_Sym(bin, sym.sh_offset + (sym.sh_entsize * elf64_r_sym(relp.r_info)))
            if symp.st_name == 0
                continue
            end
            push!(
                symbols_rel_info,
                unsafe_string(pointer(bin, str.sh_offset + symp.st_name + 1)) => relp,
            )
        end
    end

    symbols_info, symbols_rel_info
end

end # module

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of alignments.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Alignments" begin
    # Left
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │ false  │ 1.0    │ 1      │
│ 2      │ true   │ 2.0    │ 2      │
│ 3      │ false  │ 3.0    │ 3      │
│ 4      │ true   │ 4.0    │ 4      │
│ 5      │ false  │ 5.0    │ 5      │
│ 6      │ true   │ 6.0    │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = :l)
    @test result == expected

    # Center
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│   1    │ false  │  1.0   │   1    │
│   2    │  true  │  2.0   │   2    │
│   3    │ false  │  3.0   │   3    │
│   4    │  true  │  4.0   │   4    │
│   5    │ false  │  5.0   │   5    │
│   6    │  true  │  6.0   │   6    │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = :c)
    @test result == expected

    # Per column configuration
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │      1 │
│ 2      │   true │  2.0   │      2 │
│ 3      │  false │  3.0   │      3 │
│ 4      │   true │  4.0   │      4 │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = [:l,:r,:c,:r])
    @test result == expected

    # Cell override
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │ 1      │
│ 2      │   true │  2.0   │      2 │
│      3 │ false  │  3.0   │   3    │
│ 4      │   true │  4.0   │   4    │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ))
    @test result == expected

    # Headers
    # ==========================================================================

    header = (["A", "B", "C", "D"],
              ["a", "b", "c", "d"])

    expected = """
┌───────────┬───────────┬───────────┬───────────┐
│ A         │     B     │         C │         D │
│ a         │     b     │         c │         d │
├───────────┼───────────┼───────────┼───────────┤
│ 1         │     false │    1.0    │ 1         │
│ 2         │      true │    2.0    │         2 │
│         3 │ false     │    3.0    │     3     │
│ 4         │      true │    4.0    │     4     │
│ 5         │     false │    5.0    │         5 │
│ 6         │      true │    6.0    │ 6         │
└───────────┴───────────┴───────────┴───────────┘
"""

    result = pretty_table(String, data;
                          header = header,
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ),
                          columns_width = 9,
                          header_alignment = [:l,:c,:r,:r])
    @test result == expected

    expected = """
┌───────────┬───────────┬───────────┬───────────┐
│         A │     B     │         C │         D │
│ a         │ b         │         c │     d     │
├───────────┼───────────┼───────────┼───────────┤
│ 1         │     false │    1.0    │ 1         │
│ 2         │      true │    2.0    │         2 │
│         3 │ false     │    3.0    │     3     │
│ 4         │      true │    4.0    │     4     │
│ 5         │     false │    5.0    │         5 │
│ 6         │      true │    6.0    │ 6         │
└───────────┴───────────┴───────────┴───────────┘
"""

    result = pretty_table(String, data;
                          header = header,
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ),
                          columns_width = 9,
                          header_alignment = [:l,:c,:r,:r],
                          header_cell_alignment = Dict( (1,1) => :r,
                                                        (2,2) => :l,
                                                        (2,4) => :c))
    @test result == expected
end

@testset "Column alignment regex" begin
    matrix = [ i == 2 ? missing :
               i == 5 ? nothing : (-1)^j*10.0^(-i+j) for i = 1:7, j = 1:7]

    expected = """
┌────────────┬────────────┬──────────────┬──────────────┬────────────────┬────────────────┬────────────────┐
│     Col. 1 │     Col. 2 │       Col. 3 │       Col. 4 │         Col. 5 │         Col. 6 │         Col. 7 │
├────────────┼────────────┼──────────────┼──────────────┼────────────────┼────────────────┼────────────────┤
│ -1.0       │ 10.0       │ -100.0       │ 1000.0       │ -10000.0       │ 100000.0       │     -1.0e6     │
│    missing │    missing │      missing │      missing │        missing │        missing │        missing │
│ -0.01      │  0.1       │   -1.0       │   10.0       │   -100.0       │   1000.0       │ -10000.0       │
│ -0.001     │  0.01      │   -0.1       │    1.0       │    -10.0       │    100.0       │  -1000.0       │
│    nothing │    nothing │      nothing │      nothing │        nothing │        nothing │        nothing │
│ -1.0e-5    │  0.0001    │   -0.001     │    0.01      │     -0.1       │      1.0       │    -10.0       │
│ -1.0e-6    │  1.0e-5    │   -0.0001    │    0.001     │     -0.01      │      0.1       │     -1.0       │
└────────────┴────────────┴──────────────┴──────────────┴────────────────┴────────────────┴────────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(0 => [r"\."]))

    @test result == expected

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(i => [r"\."] for i = 1:7))

    @test result == expected

    expected = """
┌───────────┬───────────┬─────────────┬─────────────┬───────────────┬───────────────┬───────────────┐
│    Col. 1 │    Col. 2 │      Col. 3 │      Col. 4 │        Col. 5 │        Col. 6 │        Col. 7 │
├───────────┼───────────┼─────────────┼─────────────┼───────────────┼───────────────┼───────────────┤
│ -1.0      │ 10.0      │ -100.0      │ 1000.0      │ -10000.0      │ 100000.0      │     -1.0e6    │
│   missing │   missing │     missing │     missing │       missing │       missing │       missing │
│ -0.01     │  0.1      │   -1.0      │   10.0      │   -100.0      │   1000.0      │ -10000.0      │
│ -0.001    │  0.01     │   -0.1      │    1.0      │    -10.0      │    100.0      │  -1000.0      │
│   nothing │   nothing │     nothing │     nothing │       nothing │       nothing │       nothing │
│ -1.0e-5   │  0.0001   │   -0.001    │    0.01     │     -0.1      │      1.0      │    -10.0      │
│ -1.0e-6   │  1.0e-5   │   -0.0001   │    0.001    │     -0.01     │      0.1      │     -1.0      │
└───────────┴───────────┴─────────────┴─────────────┴───────────────┴───────────────┴───────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(0 => [r"\.", r"^"]))

    @test result == expected

    # Filters
    # ==========================================================================

    expected = """
┌─────────┬─────────────┬──────────┐
│  Col. 1 │      Col. 4 │   Col. 5 │
├─────────┼─────────────┼──────────┤
│    -1.0 │ 1000.0      │ -10000.0 │
│ missing │     missing │  missing │
│   -0.01 │   10.0      │   -100.0 │
│  -0.001 │    1.0      │    -10.0 │
│ nothing │     nothing │  nothing │
│ -1.0e-5 │    0.01     │     -0.1 │
│ -1.0e-6 │    0.001    │    -0.01 │
└─────────┴─────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          filters_col = ((data,j)->j ∈ [1,4,5],))

    @test result == expected

    # Additional rows
    # ==========================================================================

    expected = """
┌─────┬─────────┬─────────────┬──────────┐
│ Row │  Col. 1 │      Col. 4 │   Col. 5 │
├─────┼─────────┼─────────────┼──────────┤
│   1 │    -1.0 │ 1000.0      │ -10000.0 │
│   2 │ missing │     missing │  missing │
│   3 │   -0.01 │   10.0      │   -100.0 │
│   4 │  -0.001 │    1.0      │    -10.0 │
│   5 │ nothing │     nothing │  nothing │
│   6 │ -1.0e-5 │    0.01     │     -0.1 │
│   7 │ -1.0e-6 │    0.001    │    -0.01 │
└─────┴─────────┴─────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          filters_col = ((data,j)->j ∈ [1,4,5],),
                          show_row_number = true)

    @test result == expected

    expected = """
┌─────┬───┬─────────┬─────────────┬──────────┐
│ Row │   │  Col. 1 │      Col. 4 │   Col. 5 │
├─────┼───┼─────────┼─────────────┼──────────┤
│   1 │ a │    -1.0 │ 1000.0      │ -10000.0 │
│   2 │ a │ missing │     missing │  missing │
│   3 │ a │   -0.01 │   10.0      │   -100.0 │
│   4 │ a │  -0.001 │    1.0      │    -10.0 │
│   5 │ a │ nothing │     nothing │  nothing │
│   6 │ a │ -1.0e-5 │    0.01     │     -0.1 │
│   7 │ a │ -1.0e-6 │    0.001    │    -0.01 │
└─────┴───┴─────────┴─────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          filters_col = ((data,j)->j ∈ [1,4,5],),
                          row_names = ["a" for i = 1:7],
                          show_row_number = true)

    @test result == expected

    # Column alignment
    # ==========================================================================

    expected = """
┌─────────┬──────────────────────┬──────────┐
│  Col. 1 │               Col. 4 │   Col. 5 │
├─────────┼──────────────────────┼──────────┤
│    -1.0 │          1000.0      │ -10000.0 │
│ missing │              missing │  missing │
│   -0.01 │            10.0      │   -100.0 │
│  -0.001 │             1.0      │    -10.0 │
│ nothing │              nothing │  nothing │
│ -1.0e-5 │             0.01     │     -0.1 │
│ -1.0e-6 │             0.001    │    -0.01 │
└─────────┴──────────────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          columns_width = [-1, -1, -1, 20, -1, -1, -1],
                          filters_col = ((data,j)->j ∈ [1,4,5],))

    @test result == expected

    expected = """
┌─────────┬──────────────────────┬──────────┐
│ Col. 1  │        Col. 4        │  Col. 5  │
├─────────┼──────────────────────┼──────────┤
│  -1.0   │     1000.0           │ -10000.0 │
│ missing │         missing      │ missing  │
│  -0.01  │       10.0           │  -100.0  │
│ -0.001  │        1.0           │  -10.0   │
│ nothing │         nothing      │ nothing  │
│ -1.0e-5 │        0.01          │   -0.1   │
│ -1.0e-6 │        0.001         │  -0.01   │
└─────────┴──────────────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment = :c,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          columns_width = [-1, -1, -1, 20, -1, -1, -1],
                          filters_col = ((data,j)->j ∈ [1,4,5],))

    @test result == expected

    expected = """
┌─────────┬──────────────────────┬──────────┐
│ Col. 1  │ Col. 4               │ Col. 5   │
├─────────┼──────────────────────┼──────────┤
│ -1.0    │ 1000.0               │ -10000.0 │
│ missing │     missing          │ missing  │
│ -0.01   │   10.0               │ -100.0   │
│ -0.001  │    1.0               │ -10.0    │
│ nothing │     nothing          │ nothing  │
│ -1.0e-5 │    0.01              │ -0.1     │
│ -1.0e-6 │    0.001             │ -0.01    │
└─────────┴──────────────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment = :l,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          columns_width = [-1, -1, -1, 20, -1, -1, -1],
                          filters_col = ((data,j)->j ∈ [1,4,5],))

    @test result == expected

    # Minimum and maximum column width
    # ==========================================================================

    expected = """
┌─────────────────┬─────────────────┬─────────────────┐
│          Col. 1 │          Col. 4 │          Col. 5 │
├─────────────────┼─────────────────┼─────────────────┤
│            -1.0 │     1000.0      │        -10000.0 │
│         missing │         missing │         missing │
│           -0.01 │       10.0      │          -100.0 │
│          -0.001 │        1.0      │           -10.0 │
│         nothing │         nothing │         nothing │
│         -1.0e-5 │        0.01     │            -0.1 │
│         -1.0e-6 │        0.001    │           -0.01 │
└─────────────────┴─────────────────┴─────────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          minimum_columns_width = 15,
                          filters_col = ((data,j)->j ∈ [1,4,5],))

    @test result == expected

    expected = """
┌─────────┬────────────┬──────────┐
│  Col. 1 │     Col. 4 │   Col. 5 │
├─────────┼────────────┼──────────┤
│    -1.0 │ 1000.0   … │ -10000.0 │
│ missing │     missi… │  missing │
│   -0.01 │   10.0   … │   -100.0 │
│  -0.001 │    1.0   … │    -10.0 │
│ nothing │     nothi… │  nothing │
│ -1.0e-5 │    0.01  … │     -0.1 │
│ -1.0e-6 │    0.001 … │    -0.01 │
└─────────┴────────────┴──────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(4 => [r"\.", r"^"]),
                          maximum_columns_width = 10,
                          filters_col = ((data,j)->j ∈ [1,4,5],))

    @test result == expected

    # Fallback alignment
    # ==========================================================================

    expected = """
┌──────────┬──────────┬───────────┬──────────┬────────────┬────────────┬────────────┐
│   Col. 1 │   Col. 2 │    Col. 3 │   Col. 4 │     Col. 5 │     Col. 6 │     Col. 7 │
├──────────┼──────────┼───────────┼──────────┼────────────┼────────────┼────────────┤
│  -1.0    │  10.0    │ -100.0    │ 1000.0   │ -10000.0   │ 100000.0   │     -1.0e6 │
│ missing  │ missing  │  missing  │  missing │    missing │    missing │    missing │
│  -0.01   │   0.1    │   -1.0    │   10.0   │   -100.0   │   1000.0   │ -10000.0   │
│  -0.001  │   0.01   │   -0.1    │    1.0   │    -10.0   │    100.0   │  -1000.0   │
│ nothing  │ nothing  │  nothing  │  nothing │    nothing │    nothing │    nothing │
│  -1.0e-5 │   0.0001 │   -0.001  │    0.01  │     -0.1   │      1.0   │    -10.0   │
│  -1.0e-6 │   1.0e-5 │   -0.0001 │    0.001 │     -0.01  │      0.1   │     -1.0   │
└──────────┴──────────┴───────────┴──────────┴────────────┴────────────┴────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(i => [r"\."] for i = 1:7),
                          alignment_anchor_fallback = :c)

    @test result == expected

    expected = """
┌──────────────┬──────────────┬──────────────┬─────────────┬────────────┬───────────┬─────────────┐
│       Col. 1 │       Col. 2 │       Col. 3 │      Col. 4 │     Col. 5 │    Col. 6 │      Col. 7 │
├──────────────┼──────────────┼──────────────┼─────────────┼────────────┼───────────┼─────────────┤
│      -1.0    │      10.0    │    -100.0    │    1000.0   │  -10000.0  │  100000.0 │      -1.0e6 │
│ missing      │ missing      │ missing      │ missing     │ missing    │ missing   │ missing     │
│      -0.01   │       0.1    │      -1.0    │      10.0   │    -100.0  │    1000.0 │  -10000.0   │
│      -0.001  │       0.01   │      -0.1    │       1.0   │     -10.0  │     100.0 │   -1000.0   │
│ nothing      │ nothing      │ nothing      │ nothing     │ nothing    │ nothing   │ nothing     │
│      -1.0e-5 │       0.0001 │      -0.001  │       0.01  │      -0.1  │       1.0 │     -10.0   │
│      -1.0e-6 │       1.0e-5 │      -0.0001 │       0.001 │      -0.01 │       0.1 │      -1.0   │
└──────────────┴──────────────┴──────────────┴─────────────┴────────────┴───────────┴─────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(i => [r"\."] for i = 1:7),
                          alignment_anchor_fallback = :r)

    @test result == expected

    expected = """
┌────────────┬──────────┬──────────────┬─────────────┬────────────────┬────────────────┬────────────────┐
│     Col. 1 │   Col. 2 │       Col. 3 │      Col. 4 │         Col. 5 │         Col. 6 │         Col. 7 │
├────────────┼──────────┼──────────────┼─────────────┼────────────────┼────────────────┼────────────────┤
│ -1.0       │  10.0    │ -100.0       │    1000.0   │ -10000.0       │ 100000.0       │     -1.0e6     │
│    missing │ missing  │      missing │ missing     │        missing │        missing │        missing │
│ -0.01      │   0.1    │   -1.0       │      10.0   │   -100.0       │   1000.0       │ -10000.0       │
│ -0.001     │   0.01   │   -0.1       │       1.0   │    -10.0       │    100.0       │  -1000.0       │
│    nothing │ nothing  │      nothing │ nothing     │        nothing │        nothing │        nothing │
│ -1.0e-5    │   0.0001 │   -0.001     │       0.01  │     -0.1       │      1.0       │    -10.0       │
│ -1.0e-6    │   1.0e-5 │   -0.0001    │       0.001 │     -0.01      │      0.1       │     -1.0       │
└────────────┴──────────┴──────────────┴─────────────┴────────────────┴────────────────┴────────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(0 => [r"\."]),
                          alignment_anchor_fallback_override = Dict(2 => :c, 4 => :r, 5 => :l))

    @test result == expected

    expected = """
┌──────────┬──────────┬───────────┬─────────────┬────────────────┬────────────┬────────────┐
│   Col. 1 │   Col. 2 │    Col. 3 │      Col. 4 │         Col. 5 │     Col. 6 │     Col. 7 │
├──────────┼──────────┼───────────┼─────────────┼────────────────┼────────────┼────────────┤
│  -1.0    │  10.0    │ -100.0    │    1000.0   │ -10000.0       │ 100000.0   │     -1.0e6 │
│ missing  │ missing  │  missing  │ missing     │        missing │    missing │    missing │
│  -0.01   │   0.1    │   -1.0    │      10.0   │   -100.0       │   1000.0   │ -10000.0   │
│  -0.001  │   0.01   │   -0.1    │       1.0   │    -10.0       │    100.0   │  -1000.0   │
│ nothing  │ nothing  │  nothing  │ nothing     │        nothing │    nothing │    nothing │
│  -1.0e-5 │   0.0001 │   -0.001  │       0.01  │     -0.1       │      1.0   │    -10.0   │
│  -1.0e-6 │   1.0e-5 │   -0.0001 │       0.001 │     -0.01      │      0.1   │     -1.0   │
└──────────┴──────────┴───────────┴─────────────┴────────────────┴────────────┴────────────┘
"""

    result = pretty_table(String, matrix,
                          alignment_anchor_regex = Dict(0 => [r"\."]),
                          alignment_anchor_fallback = :c,
                          alignment_anchor_fallback_override = Dict(2 => :c, 4 => :r, 5 => :l))

    @test result == expected
end


text:
""" Representing an #undef value """
struct JLUNDEF
end

"""
Creates a new array that has `JLUNDEF()` objects in places of undefined values.
"""
function undefs_replaced(arr::AbstractArray)
   ret = Array{Any}(undef, size(arr))
   @inbounds for i in eachindex(arr)
      if isassigned(arr, i)
         ret[i] = arr[i]
      else
         ret[i] = JLUNDEF()
      end
   end
   ret
end

"""
Copies all elements that are not `JLUNDEF()` from `src`
to the corresponding index in `dst`
"""
function copy_defined!(dst::Array, src::Array)
   @inbounds for i in eachindex(src)
      if src[i] != JLUNDEF()
         dst[i] = src[i]
      end
   end
   dst
end


text:
using QuadGK
using SpecialFunctions
using Roots
using RCall

function Gammadist(gamma)
	return ((adap.be^adap.al)/SpecialFunctions.gamma(adap.al))*(gamma^(adap.al-1))*exp(-adap.be*gamma)
end

function pip0(gamma)
	U = 4*adap.theta_f*adap.Lf/(2.0*adap.NN)
	R = 2*adap.Lf*adap.rho/(2.0*adap.NN)
	return Gammadist(gamma)*exp(-(Gammadist(gamma)*U/(2.0*adap.NN))/(gamma/(adap.NN+0.0)+R/(2.0*adap.NN)))
end

function intpip0()

	f(gam) = pip0(gam)
	return QuadGK.quadgk(f,0,1000)[1]
end

function calcBGam(L,alpha,beta,theta)
	# not sure whats going on here, seems to overestimate or under sometimes
	u = 2*theta/(2.0*adap.NN)
	r = adap.rho/(2.0*adap.NN)

	a = -1.0*(2^(alpha+1))*exp(L*adap.NN*r*beta)*L*adap.N*((1.0/(L*adap.N*r))^(1-alpha))*u*(beta^alpha)
	b =  convert(Float64,SpecialFunctions.gamma_inc(1-alpha,L*adap.NN*r*beta,1)[2])
	#fudge = 1-gamInt.cdf(0.2,a=adap.al2,scale=1./adap.be2)
	#fudge = 1.
	fudge = 0.25

	c = exp(a*b*fudge)

	return c
end

function set_theta_f_gam()
	i(theta) = calcBGam(adap.Lf,adap.al2,adap.be2,theta)-adap.B
	# theta_f  = fsolve(lambda theta ,0.0000001)
	theta_f  = Roots.find_zero(i,0.0000001)
	adap.theta_f = theta_f
end

function calcB(L,theta)
	t = -1.0*adap.gam_neg/(adap.NN+0.)
	u = 2.0*theta/(2.0*adap.NN)
	r = adap.rho/(2.0*adap.NN)

	#return u*t/((t+r*L)^2)
	#return u*t/(2*(t+(1.-np.exp(-2*r*L))/2.)^2)

	# Nordborg models
	#####return u/(2*r*t*(1+((1.-np.exp(-2*r*L))/2.)*(1-t)/t)^2)
	#return u/(t*(1+((1-np.exp(-2.0*r*L))/2.)*(1-t)/t)^2)
	#return u/(t*(1+r*L*(1-t)/t)^2)
	#####return u/(t*(1+(np.exp(-2*r*L)/2)/t)^2)
end

# Br() needs L and theta, not working right now
function get_B_vals()
	ret = Array{Float64}(undef,30)
	for i in 20:50
		L = convert(Int64,round(1.3^i,digits=0))
		ret[i] = (Br(L),L)
	end
	return ret
end

"""
	Estimating and plotting MAP using locfit and ggplot2 in R. It assume your folder contains the posterior estimated through ABCreg
"""
function plotMap(;analysis::String,output::String)

	try
		@eval using RCall
		@eval R"""library(ggplot2);library(abc)"""

		out = filter(x -> occursin("post",x), readdir(analysis,join=true))
		out = filter(x -> !occursin(".1.",x),out)

		open(x) = Array(CSV.read(GZip.open(x),DataFrame,header=false))
		posteriors = open.(out)

		maxp = DataFrame(Array{Float64,2}(undef,size(posteriors,1),5),[:aw,:as,:a,:gamNeg,:shape])
		R"""getmap <- function(df){
				temp = as.data.frame(df)
			    d <-locfit(~temp[,1],temp);
			    map<-temp[,1][which.max(predict(d,newdata=temp))]
			}"""
		getmap(x) = rcopy(R"""matrix(apply($x,2,getmap),nrow=1)""")
		tmp = getmap.(posteriors)
		maxp = DataFrame(vcat(tmp...),[:aw,:as,:a,:gamNeg,:shape])

		al  = maxp[:,1:3]
		gam  = maxp[:,4:end]
		p = rcopy(R"""al = $al
			dal = melt(al)
			pal = ggplot(dal) + geom_density(aes(x=value,fill=variable),alpha=0.5) + scale_fill_manual(values=c('#30504f', '#e2bd9a', '#ab2710'))
			ggsave(pal,filename=paste0($output))
			""")
		return(maxp)
	catch
		println("Please install R, ggplot2 and abc in your system before execute this function")
	end
end

function set_ppos!(param::parameters)

	function f!(F,x,param=param)
		F[1] = alpha_exp_sim_tot(param,x[1],x[2])-param.al_tot
		F[2] = alpha_exp_sim_low(param,x[1],x[2])-param.al_low
	end

	ppos_l,ppos_h = nlsolve(f!,[0.0; 0.0]).zero

	param.ppos_l,param.ppos_h = ppos_l, ppos_h

end

function poisson_fixation(;observed_values::Vector{S}, λds::Vector{Float64}, λdn::Vector{Float64},λweak::Vector{Float64},λstrong::Vector{Float64}) where S <: Union{Int64,Float64}

		ds = @. λds / (λds + λdn) * observed_values
		dn = @. λdn / (λds + λdn) * observed_values
		dweak = @. λweak / (λds + λdn) * observed_values
		dstrong = @. λstrong / (λds + λdn) * observed_values

		sampled_ds     = pois_rand.(ds)
		sampled_dn     = pois_rand.(dn)
		sampled_weak   = pois_rand.(dweak)

		if all(dstrong .> 0)
			sampled_strong = pois_rand.(dstrong)
		else
			sampled_strong = .-pois_rand.(.-dstrong)
		end

		if all(sampled_weak .≈ 0)
			sampled_weak .= 0
		elseif all(sampled_strong .≈ 0)
			sampled_strong .= 0
		end

		alphas = @. [sampled_weak/sampled_dn sampled_strong/sampled_dn (sampled_weak+sampled_strong)/sampled_dn]

		out = alphas,sampled_dn, sampled_ds
		return out
end


text:
module ICNF

    using
        CUDA,
        ComputationalResources,
        DataFrames,
        DiffEqSensitivity,
        DifferentialEquations,
        Distributions,
        DistributionsAD,
        Flux,
        ForwardDiff,
        IterTools,
        MLJBase,
        MLJModelInterface,
        MLUtils,
        NNlib,
        Optim,
        Optimisers,
        Optimization,
        OptimizationFlux,
        OptimizationOptimJL,
        OptimizationOptimisers,
        OrdinaryDiffEq,
        SciMLBase,
        ScientificTypes,
        Zygote,
        LinearAlgebra,
        Random,
        Statistics

    include("types.jl")
    include("defaults.jl")
    include("core_icnf.jl")
    include("core_cond_icnf.jl")

    include("ffjord.jl")
    include("rnode.jl")
    include("planar.jl")

    include("cond_ffjord.jl")
    include("cond_rnode.jl")
    include("cond_planar.jl")

    include("metrics.jl")

    include("utils.jl")

end


text:
"""
    logseriescdf(p::T) where T <: Real

Returns the vector of samples of the discrete cdf of logarithmic distribution
"""
function logseriescdf(p::T) where T <: Real
  cdfs = T.([0.])
  for i in 1:100000000
    @inbounds push!(cdfs, cdfs[i]-(p^i)/(i*log(1-p)))
    # one can change atol here
    if isapprox(cdfs[i],  1.)
      return cdfs
    end
  end
  cdfs
end

"""
    levyel(θ::Real, u1::Real, u2::Real)

u1 and u2 are supposed to be random numbers form [0,1]
An element from Levy stable distribution with parameters α = 1/θ, β = 1,
γ = (cos(pi/(2*θ)))^θ and δ = 0.
Return Float, given parameter ϴ of dostribution
"""
function levyel(θ::Real, u1::Real, u2::Real)
  ϕ = pi*u1-pi/2
  v = quantile.(Exponential(1.), u2)
  γ = (cos(pi/(2*θ)))^θ
  v = ((cos(pi/(2*θ)+(1/θ-1)*ϕ))/v)^(θ-1)
  γ*v*sin(1/θ*(pi/2+ϕ))*(cos(pi/(2*θ))*cos(ϕ))^(-θ)
end

"""
    tiltedlevygen(V0::Real, α::Real; rng::AbstractRNG)

Returns a Real sampled  from the expotencialy tilted Levy stable pdf
f(x; V0, α) = exp(-V0^α) g(x; α)/exp(-V0), where g(x; α) is a stable Levy pdf
with parameters α = 1/θ, β = 1, γ = (cos(pi/(2*θ)))^θ and δ = 0.
"""
function tiltedlevygen(V0::T, α::T; rng::AbstractRNG) where T <: Real
  x = levyel(α, rand(rng, T), rand(rng, T))
  u = rand(rng, T)
  while exp(-V0^α*x)/(1500*exp(-V0)) < u
    x = levyel(α, rand(rng, T), rand(rng, T))
    u = rand(rng, T)
  end
  return x.*V0.^α
end

"""
    Ginv(y::Real, α::Real)

Returns Real, helper for the joe/frank nested copula generator
"""
Ginv(y::Real, α::Real) = ((1-y)*SpecialFunctions.gamma(1-α))^(-1/α)

"""
    InvlaJ(n::Int, α::Real)

Returns Real from the inverse laplacea transform of generator of Joe nested copula
parametered by α
"""
InvlaJ(n::Int, α::T) where T <: Real = 1-1/(n*beta(T(n), 1-α))

"""
  sampleInvlaJ(α::Real, v::Real)

Returns Int, a sample of inverce laplacea transform of generator of Joe nested copula,
given a parameter α and a random numver v ∈ [0,1]
"""
function sampleInvlaJ(α::Real, v::Real)
  if v <= α
    return 1
  else
    G = Ginv(v, α)
    if G > 2^62
      return 2^62
    else
      return (InvlaJ(floor(Int, G), α) < v) ? ceil(Int, G) : floor(Int, G)
    end
  end
end

"""
  elInvlaF(θ₁::Real, θ₀::Real, logseriescdf::Vector{Real}; rng::AbstractRNG)

Returns Int, a single sample of the inverse laplacea transform of the generator
of nested Frank copula
"""
function elInvlaF(θ₁::T, θ₀::T, logseriescdf::Vector{T}; rng::AbstractRNG) where T <: Real
  c1 = 1-exp(-θ₁)
  α = θ₀/θ₁
  if θ₀ <= 1
    v = rand(rng, T)
    X = findlast(logseriescdf .< rand(rng, T))[1]
    while v > 1/((X-α)*beta(X, 1-α))
      v = rand(rng, T)
      X = findlast(logseriescdf .< rand(rng, T))[1]
    end
    return X
  else
    v = rand(rng, T)
    X = sampleInvlaJ(α, rand(rng, T))
    while v > c1^(X-1)
      X = sampleInvlaJ(α, rand(rng, T))
      v = rand(rng, T)
    end
    return X
  end
end

"""
  nestedfrankgen(θ₁::T, θ₀::T, V0::Int, logseriescdf::Vector{T}; rng::AbstractRNG) where T <: Real

Return int, sampled from the Inverse Laplace trensform of nested
Frank copula given parametes θ₁ θ₀ (child and parent)
and an element of V0 - vector of samples of invlaplace of the parents copula
"""
function nestedfrankgen(θ₁::T, θ₀::T, V0::Int, logseriescdf::Vector{T}; rng::AbstractRNG) where T <: Real
    return sum([elInvlaF(θ₁, θ₀, logseriescdf; rng = rng) for j in 1:V0])
end


text:
module TestUtils

using BroadcastableStructs: foldlargs
using Test

@testset "foldlargs" begin
    @testset for n in 1:20
        @test foldl(string, 1:n, init = 0) == foldlargs(string, 0, 1:n...)
    end
end

end  # module


text:
__precompile__()

module WriteVTK

# All the code is based on the VTK file specification [1], plus some
# undocumented stuff found around the internet...
# [1] http://www.vtk.org/VTK/img/file-formats.pdf

export VTKCellTypes, VTKCellType
export MeshCell
export vtk_grid, vtk_save, vtk_point_data, vtk_cell_data
export vtk_multiblock
export paraview_collection, collection_add_timestep, paraview_collection_load
export vtk_write_array

import CodecZlib: ZlibCompressorStream
import TranscodingStreams

using LightXML

using Base64: base64encode

import Base: close, isopen, show

# Cell type definitions as in vtkCellType.h
include("VTKCellTypes.jl")

## Constants ##
const DEFAULT_COMPRESSION_LEVEL = 6
const IS_LITTLE_ENDIAN = ENDIAN_BOM == 0x04030201

## Types ##
abstract type VTKFile end

_compression_level(x::Bool) = x ? DEFAULT_COMPRESSION_LEVEL : 0
_compression_level(x) = Int(x)

struct DatasetFile <: VTKFile
    xdoc::XMLDocument
    path::String
    grid_type::String
    Npts::Int           # Number of grid points.
    Ncls::Int           # Number of cells.
    compression_level::Int  # Compression level for zlib (if 0, compression is disabled)
    appended::Bool      # Data is appended? (otherwise it's written inline, base64-encoded)
    buf::IOBuffer       # Buffer with appended data.
    function DatasetFile(xdoc, path, grid_type, Npts, Ncls, compression,
                         appended)
        buf = IOBuffer()
        clevel = _compression_level(compression)
        if !(0 ≤ clevel ≤ 9)
            error("Unexpected value of `compress` argument: $compression.\n",
                  "It must be a `Bool` or a value between 0 and 9.")
        end
        if !appended  # in this case we don't need a buffer
            close(buf)
        end
        new(xdoc, path, grid_type, Npts, Ncls, clevel, appended, buf)
    end
end

function show(io::IO, vtk::DatasetFile)
    open_str = isopen(vtk) ? "open" : "closed"
    print(io, "VTK file '$(vtk.path)' ($(vtk.grid_type) file, $open_str)")
end

struct MultiblockFile <: VTKFile
    xdoc::XMLDocument
    path::String
    blocks::Vector{VTKFile}
    # Constructor.
    MultiblockFile(xdoc, path) = new(xdoc, path, VTKFile[])
end

struct CollectionFile <: VTKFile
    xdoc::XMLDocument
    path::String
    timeSteps::Vector{String}
    # Constructor.
    CollectionFile(xdoc, path) = new(xdoc, path, VTKFile[])
end

struct MeshCell{V <: AbstractVector{<:Integer}}
    ctype::VTKCellTypes.VTKCellType  # cell type identifier (see VTKCellTypes.jl)
    connectivity::V      # indices of points (one-based, following the convention in Julia)
    function MeshCell{V}(ctype::VTKCellTypes.VTKCellType, conn::V) where V
        if ctype.nodes ∉ (length(conn), -1)
            error("Wrong number of nodes in connectivity vector.")
        end
        new(ctype, conn)
    end
end

MeshCell(ctype, conn::V) where V = MeshCell{V}(ctype, conn)

close(vtk::VTKFile) = free(vtk.xdoc)
isopen(vtk::VTKFile) = (vtk.xdoc.ptr != C_NULL)

# Add a default extension to the filename,
# unless the user have already given the correct one
function add_extension(filename, default_extension) :: String
    path, ext = splitext(filename)
    if ext == default_extension
        return filename
    end
    if !isempty(ext)
        @warn("detected extension '$(ext)' does not correspond to " *
              "dataset type.\nAppending '$(default_extension)' to filename.")
    end
    filename * default_extension
end

# Multiblock-specific functions and types.
include("gridtypes/multiblock.jl")
include("gridtypes/ParaviewCollection.jl")

# Grid-specific functions and types.
include("gridtypes/structured.jl")
include("gridtypes/unstructured.jl")
include("gridtypes/rectilinear.jl")
include("gridtypes/imagedata.jl")
include("gridtypes/array.jl")

# Common functions.
include("gridtypes/common.jl")

# This allows using do-block syntax for generation of VTK files.
for func in (:vtk_grid, :vtk_multiblock, :paraview_collection,
             :paraview_collection_load)
    @eval begin
        """
            $($func)(f::Function, args...; kwargs...)

        Create VTK file and apply `f` to it.
        The file is automatically closed by the end of the call.

        This allows to use the do-block syntax for creating VTK files:

        ```julia
        saved_files = $($func)(args...; kwargs...) do vtk
            # do stuff with the `vtk` file handler
        end
        ```
        """
        function ($func)(f::Function, args...; kwargs...)
            vtk = ($func)(args...; kwargs...)
            local outfiles
            try
                f(vtk)
            finally
                outfiles = vtk_save(vtk)
            end
            outfiles :: Vector{String}
        end
    end
end

end


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000), camera=campixel!);
    screen = display(scene)

    nrows = 2
    ncols = 2

    maingl = GridLayout(scene, 1, 1, alignmode=Outside(30))

    gridgl = maingl[1, 1] = GridLayout(
        nrows, ncols)

    r1 = gridgl[1, 1] = LRect(scene, color=:green, strokewidth=10)
    gridgl[1, 2] = LRect(scene, color=:blue)
    gridgl[2, 1] = LRect(scene, color=:red)

    gridgl[2, 2] = LAxis(scene)
    gridgl[2, 2, Top()] = LRect(scene, height=40, strokevisible=false)
    gridgl[2, 2, Left()] = LRect(scene, strokevisible=false)
    gridgl[2, 2, Bottom()] = LRect(scene, strokevisible=false)

    slgl = gridgl[3, :] = GridLayout(1, 1)
    sl = slgl[1, 1] = LSlider(scene, height = 40, range = 0:100)
    but = slgl[1, 2] = LButton(scene, label="Hello")

    sl2 = slgl[2, 1] = LSlider(scene, height = 40, range = 0:100)
    but2 = slgl[2, 2] = LButton(scene, label="Hello Hello")

    on(sl.value) do val
        r1.attributes.strokewidth = val
    end
    nothing
end

# save("layout.png", scene)


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT
module SchellingExample

export run_schelling_model!

using Agents
using Random
using InteractiveDynamics
using CairoMakie
#using GLMakie # fails the CI tests (on linux container)
using Statistics: mean
using Dates

mutable struct SchellingAgent <: AbstractAgent
    id::Int
    pos::NTuple{2,Int}
    mood::Bool
    group::Int
end

# This was causing some weird problem when I try to use
# it as part of a module.
# @agent SchellingAgent GridAgent{2} begin
#     mood::Bool
#     group::Int
# end

function initialize(;
    numagents = 320,
    griddims = (20,20),
    min_to_be_happy = 3,
    seed = 125
    )
    # define all parameters for the model
    space = GridSpace(griddims, periodic=false)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    rng = Random.MersenneTwister(seed)
    # create the model
    model = ABM(
        SchellingAgent,
        space;
        properties,
        rng,
        scheduler = Schedulers.randomly
    )
    # populate the model with Agents
    for n in 1:numagents
        agent = SchellingAgent(
            n, # id
            (1,1), # does not matter see add_agent_single!
            false, # mood
            n < numagents / 2 ? 1 : 2 # group-no
            )
        # adds the agent in a random empty location in the grid
        add_agent_single!(agent, model)
    end
    return model
end

function agent_step!(agent, model)
    minhappy = model.min_to_be_happy
    count_neighbors_same_group = 0
    # count the number of nearby neighbors belonging to same group
    for neighbor in nearby_agents(agent, model)
        if agent.group == neighbor.group
            count_neighbors_same_group += 1
        end
    end
    if count_neighbors_same_group ≥ minhappy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
    return
end

function video_filename(file_prefix)
    file_prefix * Dates.format(Dates.now(),"yymmddHHMMSS") * ".mp4"
end

function run_schelling_model!(total_ticks, file_prefix)
    # The agent data to collect each step
    x(agent) = agent.pos[1]
    adata = [(:mood, sum), (x, mean)]
    alabels = ["happy", "avg. x"]

    model = initialize(; numagents = 300)
    video_model = deepcopy(model)

    groupcolor(a) = a.group == 1 ? :blue : :orange
    groupmarker(a) = a.group == 1 ? :circle : :rect
    # A bit of a hack!
    # Unfortunately the abm_video does not allow us to record data
    # so, we are creating two copies of the same model and running
    # them one after the other. The same number of timesteps
    # as the model is deterministic because of the random-seed
    # repeated runs should be exactly same and produce same results
    data, _ = run!(model, agent_step!, total_ticks; adata)
    filename = video_filename(file_prefix)
    abm_video(
        filename,
        video_model,
        agent_step!;
        ac = groupcolor,
        am = groupmarker,
        as = 10,
        framerate = 4,
        frames = total_ticks + 1, # plus one to account for 0th frame
        title = "Schelling's segregation model",
        adata = adata,
        alabels = alabels,
        # There is a SPF (steps per frame param)
        # That defaults to 1 and should never be changed for this hack to make sense
    )
    return (data, filename)
end

# (data, _) = run_model!(20,"schelling")
# data

# To create interactive application
# parange = Dict(:min_to_be_happy => 0:8)

# figure, adf, mdf = abm_data_exploration(
#     model, agent_step!, dummystep, parange;
#     ac = groupcolor, am = groupmarker, as = 10,
#     adata, alabels
# )

# To generate a static plot of the model

# figure, _ = abm_plot(
#     model;
#     ac = groupcolor,
#     am = groupmarker,
#     as = 10
# )
# figure

# To generate a video (post processed) of the ABM

end


text:
roomList = generateTreeView(
    ("Room", "X", "Y", "Width", "Height"),
    getTreeData(loadedState.map),
    resize=true,
    sortable=false,
    editable=Bool[
        true,
        true,
        true,
        false,
        false
    ],
    callbacks=listViewCallbackUnion[
        function(store, col, row, v)
            if v != ""
                room = Maple.getRoomByName(loadedState.map, store[row, 1])
                wantedRoom = Maple.getRoomByName(loadedState.map, v)

                if isa(wantedRoom, Maple.Room) && room != wantedRoom
                    info_dialog("The selected room name is already in use.", window)

                else
                    deleteDrawableRoomCache(loadedState.map, room)

                    store[row, col] = v

                    room.name = v

                    loadedState.room = room
                    loadedState.roomName = v

                    persistence["files_lastroom"] = v
                end
            end
        end,
        function(store, col, row, v)
            try
                newX = parse(Int, v)
                room = Maple.getRoomByName(loadedState.map, store[row, 1])
                x, y = room.position

                simple = get(config, "use_simple_room_values", true)
                multiplier = simple ? 8 : 1

                room.position = (newX * multiplier, y)
                store[row, col] = newX

                draw(canvas)

            catch

            end
        end,
        function(store, col, row, v)
            try
                newY = parse(Int, v)
                room = Maple.getRoomByName(loadedState.map, store[row, 1])
                x, y = room.position

                simple = get(config, "use_simple_room_values", true)
                multiplier = simple ? 8 : 1

                room.position = (x, newY * multiplier)
                store[row, col] = newY

                draw(canvas)

            catch

            end
        end,
        false,
        false
    ]
)

function roomListRowHandler(list::ListContainer, row)
    selected = row[1]

    if selected != loadedState.roomName
        previousRoom = loadedState.room

        loadedState.roomName = selected
        loadedState.room = getRoomByName(loadedState.map, loadedState.roomName)
        persistence["files_lastroom"] = loadedState.roomName

        reset!(camera)
        setCamera!(camera, loadedState.room.position...)
        updateDrawingLayers!(loadedState.map, loadedState.room)

        handleRoomChanged(loadedState.map, loadedState.room, previousRoom)

        draw(canvas)
    end
end

connectChanged(roomListRowHandler, roomList)

function selectLoadedRoom!(roomList::ListContainer)
    # Select the specified room or the first one
    if loadedState.room !== nothing
        selectRow!(roomList, r -> r[1] == loadedState.roomName)

    else
        selectRow!(roomList)
    end
end

function showAllRoomListColumns()
    set_gtk_property!.(roomList.cols, :visible, true)

    persistence["room_list_column_visibility"] = "all"
end

function showRoomListNameColumn()
    set_gtk_property!(roomList.cols[1], :visible, true)
    set_gtk_property!.(roomList.cols[2:end], :visible, false)

    persistence["room_list_column_visibility"] = "roomOnly"
end

function hideAllRoomListColumns()
    set_gtk_property!.(roomList.cols, :visible, false)

    persistence["room_list_column_visibility"] = "none"
end

const roomListVisiblityIndices = [
    "all",
    "roomOnly",
    "none"
]

const roomListVisiblityFunctions = Dict{String, Function}(
    "all" => showAllRoomListColumns,
    "roomOnly" => showRoomListNameColumn,
    "none" => hideAllRoomListColumns
)

function getRoomListVisibilityRadioIndex()
    return something(findfirst(isequal(get(persistence, "room_list_column_visibility", "all")), roomListVisiblityIndices), 1)
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
using EndowmentsLH
using Random, Statistics, Test

# function marginals_test()
# 	@testset "Marginals" begin
# 	end
# end

function endowment_test(e)
    @testset "Endowment $e" begin
        rng = MersenneTwister(43);
        n = 50_000;
        m = marginal(e);
        # println(m);
        @test isa(m, AbstractMarginal);
        if EndowmentsLH.isbounded(m)
            @test all(EndowmentsLH.ub(m) .> EndowmentsLH.lb(m))
        end

        # println(e);
        @test isa(name(e), Symbol);
        x = EndowmentsLH.draw_test_endowments(e, n, rng);
        @test validate_draws(e, x)
        @test eltype(x) == eltype(e);
        if applicable(EndowmentsLH.mean, m)
            @test isapprox(Statistics.mean(x), EndowmentsLH.mean(m), atol = 0.1);
            @test isapprox(Statistics.std(x), EndowmentsLH.std(m), atol = 0.1);
        end

        # Quantiles
        pctV = collect(0.0 : 0.2 : 1.0);
        quantileV = marginal_quantile(m, pctV);
        if !isnothing(quantileV)
            @test size(pctV) == size(quantileV)
            @test all(diff(quantileV) .>= 0.0)
            if isbounded(m)
                @test isapprox(quantileV[1], lb(m))
                @test isapprox(quantileV[end], ub(m))
            end
        end
    end
end

function endowment_draws_test()
    rng = MersenneTwister(45);
    @testset "EndowmentDraws" begin
        # Empty
        # eV = EndowmentsLH.make_test_endowment_vector();
        ed = EndowmentDraws();
        @test isempty(ed)
        @test length(ed) == 0
        @test isnothing(names(ed))
        @test isnothing(get_meta(ed, :test))

        n = 50;
        ed = EndowmentsLH.make_test_endowment_draws(n);
        @test validate_draws(ed)
        @test length(ed) == n
        @test !isempty(ed)
        @test length(names(ed)) > 1
        @test !has_endowment(ed, :test)
        @test isnothing(get_meta(ed, :test))

        for dName in names(ed)
            @test has_endowment(ed, dName)
            endow = get_meta(ed, dName);
            @test name(endow) == dName
            @test isequal(get_label(ed, dName), label(endow))

            drawV = get_draws(ed, dName);
            @test length(drawV) == length(ed)

            draw23V = get_draws(ed, dName, 2:3);
            @test draw23V == drawV[2:3]

            # if eltype(drawV) <: Real
            #     qV = endow_quantiles(ed, dName, [0.3, 0.7]);
            #     @test eltype(qV) == eltype(drawV)
            #     @test size(qV) == (2,)
            #     @test all(diff(qV) .> 0.0)

            #     q = endow_quantiles(ed, dName, 0.7);
            #     @test isapprox(qV[2], q)
            #     nLess = sum(get_draws(ed, dName) .<= q);
            #     @test abs(nLess / n - 0.7) < 1.5 / n
            # end

            newDrawV = EndowmentsLH.draw_test_endowments(endow, length(ed), rng);
            replace_draws!(ed, dName, newDrawV);
            @test get_draws(ed, dName) == newDrawV

            replace_draws!(ed, dName, newDrawV[2]);
            @test all(x -> newDrawV[2] == x,  get_draws(ed, dName))
        end

        for j = 1 : n
            tpE = type_endowments(ed, j);
            @test size(tpE) == size(names(ed))
        end

        # Select rows
        dName = names(ed)[2];
        idxV = 2 : 3;
        ed2 = select_rows(ed, idxV);
        @test length(ed2) == length(idxV)
        for (j, idx) in enumerate(idxV)
            @test isequal(get_draws(ed, dName, idx),  get_draws(ed2, dName, j))
        end
    end
end

function endowment_corr_test()
    rng = MersenneTwister(45);
    @testset "Correlation matrix" begin
        n = 15;
        ed = EndowmentsLH.make_test_endowment_draws(n);
        nameV = names(ed);
        nVars = length(nameV);
        corrM = corr_matrix(ed, nameV);
        @test size(corrM) == (nVars, nVars)
        fcM = formatted_corr_matrix(ed, nameV);
        @test size(fcM) == (nVars + 1, nVars + 1)
        for j = 1 : nVars
            @test fcM[j+1, j+1] == "1.0"
        end

        @test isapprox(corrM[1, 3],
            cor(get_draws(ed, nameV[1]), get_draws(ed, nameV[3])))

        corr3M = corr_matrix(ed);
        @test isapprox(corrM, corr3M; nans = true)

        corr2M = corr_matrix(ed, nameV[[1,3,5]]);
        @test isapprox(corr2M[1,2],  corrM[1,3])
        @test isapprox(corr2M[1,3],  corrM[1,5])
    end
end

## -----------  Custom endowment type, non-scalar values

# This is the custom type to be drawn.
struct CustomEndow{T1}
    x :: Vector{T1}
end

# Custom marginal that can handle a custom type
# # could be built in!
# struct CustomMarginal{T1} <: AbstractMarginal{T1} end
# EndowmentsLH.isbounded(::CustomMarginal{T1}) where T1 = false;
# Need to define `validate_draws` for each marginal
# EndowmentsLH.validate_draws(endow1 :: CustomMarginal{T1}, drawV) where T1 = true;

# Drawing endowments happens outside of the Endowments package.
function draw_endowments(T1, n, rng :: AbstractRNG)
    x = [CustomEndow(rand(rng, T1, 3))  for j = 1 : n];
    return x
end

function custom_type_test()
    rng = MersenneTwister(123);
    n = 7;

    @testset "Custom type" begin
        ed = EndowmentDraws();

        # Add a regular set of draws
        endow2 = Endowment(:endow2, "Endow2", PercentileMarginal{Float32}());
        draw2V = EndowmentsLH.draw_test_endowments(endow2, n, rng);
        add_draws!(ed, endow2, draw2V);
        @test validate_draws(ed)

        # Custom endowment
        T1 = Float16;
        drawV = draw_endowments(T1, n, rng);
        marginal1 = UnknownMarginal{CustomEndow{T1}}();
        endow1 = Endowment(:endow1, "Endow1", marginal1);
        add_draws!(ed, endow1, drawV);
        @test length(ed) == n
        @test length(names(ed)) == 2
        @test validate_draws(ed)
    end
end

@testset "All" begin
    ev = EndowmentsLH.make_test_endowment_vector();
    for e in ev
        endowment_test(e)
    end
    endowment_draws_test()
    endowment_corr_test()
    custom_type_test()
    include("fixed_percentiles.jl");
    include("quantile_test.jl");
end

# -----------


text:
import FileIO
import CSV
import GeometryBasics
import Observables
import AbstractPlotting
import Makie
import DataFrames

using Porta

startframe = 1
frames = 90
resolution = (360, 360)
maxsamples = 90
segments = 90
speed = 1
α = 40 / 180 * pi
solidtop = U1(pi - 2α)
solidbottom = U1(-pi)
ghosttop = U1(pi)
ghostbottom = solidtop
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[1]
modelname = "planethopf"
FPS = 24

# Map from S² into itself
fmap(b::S²) = b

"""
sample(dataframe, part, max)

Samples points from a dataframe with the given `dataframe`, `part` id, and the `max` number
of samples limit.
"""
function sample(dataframe, part, max)
    groupdataframe = DataFrames.groupby(dataframe, :partid)
    subdataframe = groupdataframe[(partid=part,)]
    ϕ = subdataframe.x ./ 180 .* pi
    θ = subdataframe.y ./ 180 .* pi
    coordinates = map(x -> Geographic(1, x[1], x[2]), eachrow([ϕ θ]))
    sampledpoints = Array{Geographic,1}(undef, max)
    count = length(coordinates)
    if count > max
        indices = convert(Array{Int64}, floor.(range(1, stop=count, length=max)))
        sampledpoints = coordinates[indices]
    else
        sampledpoints = coordinates
    end
    sampledpoints
end

# Use QGIS to design a geo map
colorref = FileIO.load("data/basemap_color.png")
color = FileIO.load("data/basemap_mask.png")
# Made with Natural Earth.
# Free vector and raster map data @ naturalearthdata.com.
countries = Dict("iran" => Geographic(1, 0.29826 * pi, 0.36031 * pi / 2),
                 "us" => Geographic(1, -0.53173 * pi, 0.41211 * pi / 2),
                 "china" => Geographic(1, 0.57886 * pi, 0.39846 * pi / 2),
                 "ukraine" => Geographic(1, 0.17314 * pi, 0.53754 * pi / 2),
                 "australia" => Geographic(1, 0.74319 * pi, -0.28082 * pi / 2),
                 "germany" => Geographic(1, 0.05806 * pi, 0.56850 * pi / 2),
                 "israel" => Geographic(1, 0.19362 * pi, 0.34495 * pi / 2),
                 "canada" => Geographic(1, -0.59081 * pi, 0.62367 * pi / 2),
                 "india" => Geographic(1, 0.43868 * pi, 0.22881 * pi / 2),
                 "southkorea" => Geographic(1, 0.70981 * pi, 0.39897 * pi / 2),
                 "france" => Geographic(1, 0.01229 * pi, 0.51364 * pi / 2),
                 "uganda" => Geographic(1, 0.1794 * pi , 0.0153 * pi / 2),
                 "brazil" => Geographic(1, 0.2884 * -pi, 0.1581 * -pi / 2),
                 "antarctica" => Geographic(1, 0.75 * pi, -0.92069 * pi / 2))
highlighted = ["us", "china", "iran"]
# The path to the dataset
path = "data/natural_earth_vector"
# The scene object that contains other visual objects
# AbstractPlotting.reasonable_resolution() = (800, 800)
AbstractPlotting.primary_resolution() = resolution
scene = AbstractPlotting.Scene(backgroundcolor = :black,
                               show_axis = false,
                               resolution = resolution)

s3rotation = Quaternion(1, 0, 0, 0)
config = Biquaternion(ℝ³(0, 0, 0))
solidwhirls = []
ghostwhirls = []
for country in countries
    countryname = country[1]
    dataframe = DataFrames.DataFrame(CSV.File(joinpath(path, "$(countryname)-nodes.csv")))
    center = country[2]
    x = Int(floor((center.ϕ + pi) / 2pi * size(color, 2)))
    y = size(color, 1) - Int(floor((center.θ + pi / 2) / pi * size(color, 1)))
    r, g, b, a = colorref[y, x].r, colorref[y, x].g, colorref[y, x].b, colorref[y, x].alpha
    solidcolor = AbstractPlotting.RGBAf0(r, g, b, 0.9)
    ghostcolor = AbstractPlotting.RGBAf0(r, g, b, 0.1)
    transparency = false
    if countryname == "antarctica"
        transparency = true
        solidcolor = AbstractPlotting.RGBAf0(r, g, b, 0.1)
    end
    parts = countryname == "us" ? 2 : 1 # countryname in highlighted ? 12 : 1
    for part in 0:parts-1
        # Sample a random subset of the points
        points = sample(dataframe, part, maxsamples)
        solidwhirl = Whirl(scene,
                           points,
                           τmap,
                           fmap,
                           top = solidtop,
                           bottom = solidbottom,
                           s3rotation = s3rotation,
                           config = config,
                           segments = segments,
                           color = solidcolor,
                           transparency = transparency)
        push!(solidwhirls, solidwhirl)
        if countryname in highlighted
            ghostwhirl = Whirl(scene,
                               points,
                               τmap,
                               fmap,
                               top = ghosttop,
                               bottom = ghostbottom,
                               s3rotation = s3rotation,
                               config = config,
                               segments = segments,
                               color = ghostcolor,
                               transparency = true)
            push!(ghostwhirls, ghostwhirl)
        end
    end
end

framesprites = []
framesprite1 = Frame(scene,
                     ghostbottom,
                     τmap,
                     fmap,
                     color,
                     s3rotation = s3rotation,
                     config = config,
                     segments = segments,
                     transparency = true)
framesprite2 = Frame(scene,
                     ghosttop,
                     τmap,
                     fmap,
                     color,
                     s3rotation = s3rotation,
                     config = config,
                     segments = segments,
                     transparency = true)
push!(framesprites, framesprite1)
push!(framesprites, framesprite2)

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i)
    step = (i - 1) / frames
    u = ℝ³(0, 0, 1)
    τ = step * speed * 2pi
    q = Quaternion(τ, u)
    #β = step * 2π
    #α = 80 / 180 * pi
    #solidtop = U1(pi - α) * U1(β)
    #solidbottom = U1(-pi) * U1(β)
    #ghosttop = U1(pi) * U1(β)
    #ghostbottom = solidtop
    #f(b::S²) = Cartesian(rotate(ℝ³(Cartesian(b)), q))

    for item in solidwhirls
        update(item, q)
        #update(item, σmap, f)
        #update(item, solidtop, solidbottom)
    end

    for item in ghostwhirls
        update(item, q)
        #update(item, σmap, f)
        #update(item, ghosttop, ghostbottom)
    end
    update(framesprite1, q)
    update(framesprite2, q)
end

n = ℝ³(0, 0, 1)
v = ℝ³(-2, 2, 1) * 0.9
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display

if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                 framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    if !isdir(directory)
        mkdir(directory)
    end
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        println("Generating frame $i took $elapsed (s).")

        start = time()
        imagename = "drorbarnatan2010_$(resolution[2])_$i.jpeg"
        filepath = joinpath(directory, imagename)
        FileIO.save(filepath,
                    scene;
                    resolution = resolution,
                    pt_per_unit = 100.0,
                    px_per_unit = 100.0)
        elapsed = time() - start
        println("Saving file $filepath took $elapsed (s).")

        step = (i - 1) / frames
        println("Completed step $(100step).\n")
        sleep(1)
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using LinearRegressionKit
using Test, DataFrames, StatsModels

include("test_sweep_operator.jl")
include("test_utilities.jl")
include("test_LinearRegression.jl")
include("test_cooksd.jl")
include("test_lessthanfullrank.jl")
include("test_noint.jl")
include("test_heteroscedasticity.jl")
include("test_kfold.jl")
include("test_ridge.jl")


text:
module TensND

import Base: @pure, eltype
import LinearAlgebra: normalize, dot, tr

using LinearAlgebra, SymPy, Tensors, OMEinsum, Rotations

include("array_utils.jl")
include("bases.jl")
include("tens.jl")
include("tens_isotropic.jl")
include("special_tens.jl")
include("coorsystems.jl")

# function __init__()
# end

end # module


text:
struct StiffnessTensor{T}
    η11::Matrix{T}
    η33::Matrix{T}
    η55::Matrix{T}
    η13::Matrix{T}
    η15::Matrix{T}
    η35::Matrix{T}
end

"""
    Initialise 2D isotropic viscous tensor
"""
𝓒init(nel,nip) = StiffnessTensor(fill(4/3,nel,nip),
                                 fill(4/3,nel,nip),
                                 fill(2.0,nel,nip),
                                 fill(-2/3,nel,nip),
                                 fill(0.0,nel,nip),
                                 fill(0.0,nel,nip))

"""
    𝓒 = Tensor: [η11 η22 η33 η44 η55 η66 η12 η13 η23]
    a = Semi-axes: [a₁ a₂ a₃]
    ϕ = Volume fraction
    w = viscosity cut-off
"""
struct DEM
    𝓒::Matrix{Float64}
    a::Matrix{Float64}
    ϕ::Float64
    w::Float64
    a1a2_blk::Vector{Float64}
    a2a3_blk::Vector{Float64}
    permutation_blk::Vector{Int64}
    sblk::Int64
    nblk::Int64
    parameterization::ParameterizationWeak{Float64}
end

"""
    DEMloader(fname::String)
        Load DEM hdf5 file
        fname -> file name, including path and extension
"""
function DEMloader(fname::String)
    fid = h5open(fname,"r")
    tensor = read(fid,"C")
    axes = read(fid,"ax")
    volume = read(fid,"vol")
    return tensor,axes,volume
end

"""
    getDEM(fname::String)
        Initialise DEM structure reading .h5 from fname
"""
function getDEM(fname, Δη, ϕ)
    tensor,axis,volume = DEMloader(fname)
    # fixtensor!(tensor)

    η_cutoff = weakening(Δη, ϕ)

    # -- Get and process axes from the DEM
    ax1   = view(axis,:,1)
    ax2   = 1.0
    ax3   = view(axis,:,3)
    a1a2  = @. log10(ax1/ax2)
    a2a3  = @. log10(ax2/ax3)
    a2a3_blk, permutation_blk, sblk, nblk = sort_axes(a2a3, a1a2)
    idx = 1:sblk:length(a2a3)
    # idx = 1:sblk
    a1a2_blk = view(a1a2,idx)

    DEM(
        tensor,
        axis,
        # dropdims(volume,dims=2),
        volume,
        η_cutoff,
        Array(a1a2_blk),
        Array(a2a3_blk),
        permutation_blk,
        sblk,
        nblk,
        parameterization(Δη, ϕ)
    )

end

function get_stride(v::Vector)
    i=1
    @inbounds while v[i]==v[i+1]
        i+=1
    end
    return i
end

function sort_axes(a2a3, a1a2)
    sblk = get_stride(a1a2) # block size
    nblk = div(length(a2a3),sblk)+1 # number of blocksa1
    blk = 1:sblk
    # blk  = 1:sblk:length(a2a3)
    sorted_blk = sort(view(a2a3,blk)) # sorted block
    permutation_blk = sortperm(view(a2a3,blk)) # sorting permutations
    return sorted_blk, permutation_blk, sblk, nblk
end

function argminsorted(v)
    m = length(v)
    for i in 1:m-1
        if @inbounds v[i] < v[i+1]
            return i
        end
    end
    return m
end

function fixtensor!(C)
    dummy = zeros(2)
    i1 = [1,3]
    i2 = [7,9]
    i3 = [4,6]
    @inbounds for i in axes(C,1)
        dummy[1], dummy[2] = abs(C[i,2]-C[i,1]), abs(C[i,2]-C[i,3])
        idx = i1[argmin(dummy)]
        C[i,idx] = C[i,2] = (C[i,idx]+C[i,2])*0.5

        dummy[1], dummy[2] = abs(C[i,8]-C[i,7]), abs(C[i,8]-C[i,9])
        idx = i2[argmin(dummy)]
        C[i,8] = C[i,idx] = (C[i,8]+C[i,idx])*0.5

        dummy[1], dummy[2] = abs(C[i,4]-C[i,5]), abs(C[i,6]-C[i,5])
        idx = i3[argmin(dummy)]
        C[i,idx] = C[i,5] = (C[i,idx]+C[i,5])*0.5
    end
end

function anisotropic_tensor(FSE::Array{FiniteStrainEllipsoid{Float64},2}, D::DEM, ::Val{T}) where {T}
    # Allocate arrays
    nu_11 = Matrix{Float64}(undef, size(FSE, 1), 7)
    nu_33 = similar(nu_11)
    nu_55 = similar(nu_11)
    nu_13 = similar(nu_11)
    nu_15 = similar(nu_11)
    nu_35 = similar(nu_11)
    v₁ = [similar(D.a1a2_blk) for _ in 1:Threads.nthreads()]
    v₂ = [similar(D.a2a3_blk) for _ in 1:Threads.nthreads()]
    # max_a1a2, max_a2a3 = maximum(D.a1a2_blk), maximum(D.a2a3_blk)

    # Get η from data base and rotate it
    @batch for i in eachindex(FSE)
        get_tensor_and_rotate!( nu_11, nu_33, nu_55, nu_13, nu_15, nu_35,
                                FSE[i], D, i, v₁, v₂)
    end

    # 7th ip
    one_sixth = 1/6
    @tturbo for i in 1:size(nu_11,1)
        nu_11[i,7] = (nu_11[i,1] + nu_11[i,2] + nu_11[i,3] + nu_11[i,4] + nu_11[i,5] + nu_11[i,6]) * one_sixth
        nu_33[i,7] = (nu_33[i,1] + nu_33[i,2] + nu_33[i,3] + nu_33[i,4] + nu_33[i,5] + nu_33[i,6]) * one_sixth
        nu_55[i,7] = (nu_55[i,1] + nu_55[i,2] + nu_55[i,3] + nu_55[i,4] + nu_55[i,5] + nu_55[i,6]) * one_sixth
        nu_13[i,7] = (nu_13[i,1] + nu_13[i,2] + nu_13[i,3] + nu_13[i,4] + nu_13[i,5] + nu_13[i,6]) * one_sixth
        nu_15[i,7] = (nu_15[i,1] + nu_15[i,2] + nu_15[i,3] + nu_15[i,4] + nu_15[i,5] + nu_15[i,6]) * one_sixth
        nu_35[i,7] = (nu_35[i,1] + nu_35[i,2] + nu_35[i,3] + nu_35[i,4] + nu_35[i,5] + nu_35[i,6]) * one_sixth
    end

    StiffnessTensor(nu_11, nu_33, nu_55, nu_13, nu_15, nu_35)

end

function argminsortedsecond(v)
    imin = 0

    for i in 1:length(v)-1
        if @inbounds v[i] < v[i+1]
            imin = i
        end
    end

    if imin != 0

        if  @inbounds v[imin+1] < v[imin-1]
            return imin, imin+1
        else
            return imin, imin-1
        end

    elseif imin == 0
        return length(v), length(v)

    elseif imin == 1
        return 1, 1
    end
end

function get_tensor_and_rotate!(nu_11, nu_33, nu_55, nu_13, nu_15, nu_35,
                                FSEᵢ, D, i, v₁, v₂)

    r₁, r₂ = fabric_shape(FSEᵢ.a1, 1.0, D.parameterization)
    nt = Threads.threadid()

    @inbounds for j in eachindex(D.a1a2_blk)
        v₁[nt][j] = abs(r₁-D.a1a2_blk[j])
    end
    r₁_imin = argminsorted(v₁[nt])

    @inbounds for j in eachindex(D.a2a3_blk)
        v₂[nt][j] = abs(r₂-D.a2a3_blk[j])
    end
    r₂_imin = D.permutation_blk[argminsorted(v₂[nt])]

    im = D.sblk*(r₁_imin-1) + r₂_imin
    # Allocate stiffness tensor
    C = @SMatrix [D.𝓒[im, 1]  D.𝓒[im, 7]  D.𝓒[im, 8]    0                    0                      0
                  D.𝓒[im, 7]  D.𝓒[im, 2]  D.𝓒[im, 9]    0                    0                      0
                  D.𝓒[im, 8]  D.𝓒[im, 9]  D.𝓒[im, 3]    0                    0                      0
                  0           0           0             max(D.𝓒[im, 4],D.w)  0                      0
                  0           0           0             0                    max(D.𝓒[im, 5],D.w)    0
                  0           0           0             0                    0                      max(D.𝓒[im, 6],D.w)]

    # Rotation matrix
    R = rotmatrix(FSEᵢ.x1, FSEᵢ.y1)
    # Rotate tensor
    η11,η33,η55,η13,η15,η35 = directRotation2D(R,C)
    # Fill array
    nu_11[i] = η11
    nu_33[i] = η33
    nu_55[i] = max(η55, D.w)
    nu_13[i] = η13
    nu_15[i] = η15
    nu_35[i] = η35
end

function rotmatrix(x, y)
    h = @muladd 1/√(y*y + x*x)
    sina = y*h
    cosa = x*h
    return (cosa, sina, -sina)
end

function unrotate_anisotropic_tensor(𝓒,vx1,vx2,vy1,vy2)
    η11 = 𝓒.η11[:,1:6]
    η33 = 𝓒.η33[:,1:6]
    η55 = 𝓒.η55[:,1:6]
    η13 = 𝓒.η13[:,1:6]
    η15 = 𝓒.η15[:,1:6]
    η35 = 𝓒.η35[:,1:6]

    # -- Get η from data base and rotate it
    for i in CartesianIndices(η11)
        # Allocate stiffness tensor
        C = @SMatrix [ η11[i]  0  η13[i]  0  η15[i]  0
                       0       0  0       0  0       0
                       η13[i]  0  η33[i]  0  η35[i]  0
                       0       0  0       0  0       0
                       η15[i]  0  η35[i]  0  η55[i]  0
                       0       0  0       0  0       0]
        ## =================================================================

        ## ROTATE VISCOUS TENSOR ===========================================
        R = @SMatrix [vx1[i]  0.0   vx2[i]
                      0.0     1.0   0.0
                      vy1[i]  0.0   vy2[i]]
        # ----- Rotate tensor (fast version derived from symbolic calculus)
        # nu_11[i], nu_33[i], nu_55[i], nu_13[i], nu_15[i], nu_35[i] = directRotation2D(R,C)
        n11,n33,n55,n13,n15,n35 = directRotation2D(inv(R),C)
        𝓒.η11[i] = n11
        𝓒.η33[i] = n33
        𝓒.η55[i] = n55
        𝓒.η13[i] = n13
        𝓒.η15[i] = n15
        𝓒.η35[i] = n35
    end

    return 𝓒

end ### END rotate_tensor FUNCTION #############################################

function anisotropic_tensor_blks(FSE::Array{FiniteStrainEllipsoid{Float64},2}, D::DEM)
    # -- Allocate arrays
    n = size(FSE, 1)
    nu0 = Vector{Float64}(undef, 9)
    nu_interp = fill(0.0,6,6)
    nu_11 = Matrix{Float64}(undef, n, 6)
    nu_33 = similar(nu_11)
    nu_55 = similar(nu_11)
    nu_13 = similar(nu_11)
    nu_15 = similar(nu_11)
    nu_35 = similar(nu_11)

    # -- Fitting coefficients of the axes parameterisation
    R1,R2 = fittingcoefficients()

    # -- threading blocks
    nblk = Threads.nthreads()
    sblk = div(length(FSE),nblk)
    blk = 1:sblk

    # -- Get η from data base and rotate it
    @sync begin
        for i in 1:nblk-1
            Threads.@spawn get_tensor_and_rotate_blk!( nu_11, nu_33, nu_55, nu_13, nu_15, nu_35,
                                FSE,R1,R2,D, blk.*(1 + (i-1)))
        end
        get_tensor_and_rotate_blk!( nu_11, nu_33, nu_55, nu_13, nu_15, nu_35,
                                    FSE,R1,R2,D, sblk*(nblk-1):length(FSE))
    end

    return StiffnessTensor(nu_11, nu_33, nu_55, nu_13, nu_15, nu_35)

end

function get_tensor_and_rotate_blk!(nu_11, nu_33, nu_55, nu_13, nu_15, nu_35,
    FSE,R1,R2,D,iblk)
    @inbounds for i in iblk
        ## GET VISCOUS TENSOR ==============================================
        # Average fabric -> r₁ = log10(a1/a2) and r₂ = log10(a2/a3)
        F = FSE[i]
        r₁, r₂ = fabric_parametrisation(F.a1,
                                        1.0,
                                        F.a2,
                                        R1,
                                        R2)

        v₁ = abs.(r₁.-D.a1a2_blk)
        v₂ = abs.(r₂.-D.a2a3_blk)
        r₁_imin = argminsorted(v₁)
        r₂_imin = D.permutation_blk[argminsorted(v₂)]
        im = D.sblk*(r₁_imin-1) + r₂_imin
        𝓒 = view(D.𝓒,im,:)
        # Allocate stiffness tensor
        C = @SMatrix [𝓒[1]  𝓒[9]  𝓒[8]  0     0    0
                      𝓒[9]  𝓒[2]  𝓒[7]  0     0    0
                      𝓒[8]  𝓒[7]  𝓒[3]  0     0    0
                      0     0     0     𝓒[4] 0     0
                      0     0     0     0    𝓒[5]  0
                      0     0     0     0    0     𝓒[6]]
        ## =================================================================

        ## ROTATE VISCOUS TENSOR ===========================================
        R = @SMatrix [F.x1  0   F.x2
                      0          1   0
                      F.y1  0   F.y2]
        # ----- Rotate tensor
        # Cr = TensorRotation(R,C)
        nu_11[i], nu_33[i], nu_55[i], nu_13[i], nu_15[i], nu_35[i] = directRotation2D(R,C)
        # ==================================================================
    end
end

function custommin(a::AbstractArray, b::Float64)
    imin = 1 # index of min value
    vmin = (a[1] - b)^2 # min value

    @inbounds for i = 2:length(a)
        d = (a[i] - b)^2
        if vmin > d
            imin = i
            vmin = d
        end

    end
    return imin
end

"""
    Allocate K matrix (based on Bowers 'Applied Mechanics of Solids', Chapter 3]
"""
function getK(R)
    # 11 block
    K11 = R[1,1]*R[1,1]
    K12 = R[1,2]*R[1,2]
    K13 = R[1,3]*R[1,3]
    K21 = R[2,1]*R[2,1]
    K22 = R[2,2]*R[2,2]
    K23 = R[2,3]*R[2,3]
    K31 = R[3,1]*R[3,1]
    K32 = R[3,2]*R[3,2]
    K33 = R[3,3]*R[3,3]
    # 12 block
    K14 = 2*R[1,2]*R[1,3]
    K15 = 2*R[1,3]*R[1,1]
    K16 = 2*R[1,1]*R[1,2]
    K24 = 2*R[2,2]*R[2,3]
    K25 = 2*R[2,3]*R[2,1]
    K26 = 2*R[2,1]*R[2,2]
    K34 = 2*R[3,2]*R[3,3]
    K35 = 2*R[3,3]*R[3,1]
    K36 = 2*R[3,1]*R[3,2]
    # 21 block
    K41 =  R[2,1]*R[3,1]
    K42 =  R[2,2]*R[3,2]
    K43 =  R[2,3]*R[3,3]
    K51 =  R[3,1]*R[1,1]
    K52 =  R[3,2]*R[1,2]
    K53 =  R[3,3]*R[1,3]
    K61 =  R[1,1]*R[2,1]
    K62 =  R[1,2]*R[2,2]
    K63 =  R[1,3]*R[2,3]
    # 22 block
    K44 =  R[2,2]*R[3,3]+R[2,3]*R[3,2]
    K54 =  R[2,3]*R[3,1]+R[2,1]*R[3,3]
    K64 =  R[2,1]*R[3,2]+R[2,2]*R[3,1]
    K45 =  R[3,2]*R[1,3]+R[3,3]*R[1,2]
    K55 =  R[3,3]*R[1,1]+R[3,1]*R[1,3]
    K65 =  R[3,1]*R[1,2]+R[3,2]*R[1,1]
    K46 =  R[1,2]*R[2,3]+R[1,3]*R[2,2]
    K56 =  R[1,3]*R[2,1]+R[1,1]*R[2,3]
    K66 =  R[1,1]*R[2,2]+R[1,2]*R[2,1]

    K = @SMatrix [K11 K12 K13 K14 K15 K16
                  K21 K22 K23 K24 K25 K26
                  K31 K32 K33 K34 K35 K36
                  K41 K42 K43 K44 K45 K46
                  K51 K52 K53 K54 K55 K56
                  K61 K62 K63 K64 K65 K66]

end

function TensorRotation(R,C)
    K = getK(R) # -> rotation 6x6 matrix
    K*C*K' # -> rotated 4th rank tensor
end

function directRotation2D(R,C)
    ## Tensor coefficients
    a = C[1,1]
    # b = C[2,2]
    c = C[3,3]
    # d = C[4,4]
    e = C[5,5]
    # f = C[6,6]
    # g = C[1,2]
    h = C[1,3]
    # l = C[2,3]
    ## Rotation matrix
    x1 = y2 = R[1]
    x2 = R[2]
    y1 = R[3]
    ## Rotated components
    x1x1 = x1*x1
    x2x2 = x2*x2
    y1y1 = y1*y1
    y2y2 = y2*y2
    n11 = @muladd (c * x2x2 + h * x1x1 + 4 * e * x1x1) * x2x2 + (a * x1x1 + h * x2x2) * x1x1
    n33 = @muladd (c * y2y2 + h * y1y1 + 4 * e * y1y1) * y2y2 + (a * y1y1 + h * y2y2) * y1y1
    n55 = @muladd (c * x2 * y2 + h * x1 * y1) * x2 * y2 + (x1 * y2 + x2 * y1) *(x1 * y2 + x2 * y1) * e + (a * x1 * y1 + h * x2 * y2) * x1 * y1
    n13 = @muladd ((c * x2x2 + h * x1x1) * y2 + 4 * e * x1 * x2 * y1) * y2 + (a * x1x1 + h * x2x2) * y1y1
    n15 = @muladd ((c * x2x2 + h * x1x1) * y2 + 2 * (x1 * y2 + x2 * y1) * e * x1) * x2 + (a * x1x1 + h * x2x2) * x1 * y1
    n35 = @muladd ((c * y2y2 + h * y1y1) * x2 + 2 * (x1 * y2 + x2 * y1) * e * y1) * y2 + (a * y1y1 + h * y2y2) * x1 * y1

    return n11, n33, n55, n13, n15, n35
end

function directRotation(R,C)
    @fastmath begin
        ## Tensor coefficients
        a = C[1,1]
        b = C[2,2]
        c = C[3,3]
        d = C[4,4]
        e = C[5,5]
        f = C[6,6]
        g = C[1,2]
        h = C[1,3]
        l = C[2,3]
        ## Rotation matrix
        x1 = R[1,1]
        x2 = R[1,2]
        x3 = R[1,3]
        y1 = R[2,1]
        y2 = R[2,2]
        y3 = R[2,3]
        z1 = R[3,1]
        z2 = R[3,2]
        z3 = R[3,3]
        ## Rotated components
        n11 = 4 * ((e * x3 ^ 2 + f * x2 ^ 2) * x1*x1 + d * x2 ^ 2 * x3 ^ 2) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * x3 ^ 2 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * x1*x1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * x2 ^ 2
        n13 = 4 * ((e * x3 * z3 + f * x2 * z2) * x1 * z1 + d * x2 * x3 * z2 * z3) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * z3 ^ 2 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * z1 ^ 2 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * z2 ^ 2
        n15 = 2 * (((x1 * z3 + x3 * z1) * e * x1 + (y1 * z3 + y3 * z1) * d * x2) * x3 + (x1 * y3 + x3 * y1) * f * x1 * x2) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * x3 * z3 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * x1 * z1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * x2 * z2
        n33 = 4 * ((e * z3 ^ 2 + f * z2 ^ 2) * z1 ^ 2 + d * z2 ^ 2 * z3 ^ 2) + (h * z1 ^ 2 + l * z2 ^ 2 + c * z3 ^ 2) * z3 ^ 2 + (g * z2 ^ 2 + h * z3 ^ 2 + a * z1 ^ 2) * z1 ^ 2 + (g * z1 ^ 2 + l * z3 ^ 2 + b * z2 ^ 2) * z2 ^ 2
        n35 = 2 * (((x1 * z3 + x3 * z1) * e * z1 + (y1 * z3 + y3 * z1) * d * z2) * z3 + (x1 * y3 + x3 * y1) * f * z1 * z2) + (h * z1 ^ 2 + l * z2 ^ 2 + c * z3 ^ 2) * x3 * z3 + (g * z2 ^ 2 + h * z3 ^ 2 + a * z1 ^ 2) * x1 * z1 + (g * z1 ^ 2 + l * z3 ^ 2 + b * z2 ^ 2) * x2 * z2
        n55 = (x1 * z3 + x3 * z1) ^ 2 * e + (y1 * z3 + y3 * z1) ^ 2 * d + (x1 * y3 + x3 * y1) ^ 2 * f + (h * x1 * z1 + l * x2 * z2 + c * x3 * z3) * x3 * z3 + (g * x2 * z2 + h * x3 * z3 + a * x1 * z1) * x1 * z1 + (g * x1 * z1 + l * x3 * z3 + b * x2 * z2) * x2 * z2
    end

    return n11,n33,n55,n13,n15,n35

end

function directRotation!(n, R, C)
    ## Tensor
    a = C[1,1]
    b = C[2,2]
    c = C[3,3]
    d = C[4,4]
    e = C[5,5]
    f = C[6,6]
    g = C[1,2]
    h = C[1,3]
    l = C[2,3]
    ## Rotation matrix
    x1 = R[1,1]
    x2 = R[1,2]
    x3 = R[1,3]
    y1 = R[2,1]
    y2 = R[2,2]
    y3 = R[2,3]
    z1 = R[3,1]
    z2 = R[3,2]
    z3 = R[3,3]
    ## Rotated components
    n[1,1] = 4 * ((e * x3 ^ 2 + f * x2 ^ 2) * x1*x1 + d * x2 ^ 2 * x3 ^ 2) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * x3 ^ 2 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * x1*x1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * x2 ^ 2
    n[1,2] = n[1,2] = 4 * ((e * x3 * y3 + f * x2 * y2) * x1 * y1 + d * x2 * x3 * y2 * y3) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * y3 ^ 2 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * y1*y1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * y2*y2
    n[2,2] = 4 * ((e * y3 ^ 2 + f * y2*y2) * y1*y1 + d * y2*y2 * y3 ^ 2) + (h * y1*y1 + l * y2*y2 + c * y3 ^ 2) * y3 ^ 2 + (g * y2*y2 + h * y3 ^ 2 + a * y1*y1) * y1*y1 + (g * y1*y1 + l * y3 ^ 2 + b * y2*y2) * y2*y2
    n[1,3] = n[1,3] = 4 * ((e * x3 * z3 + f * x2 * z2) * x1 * z1 + d * x2 * x3 * z2 * z3) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * z3 ^ 2 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * z1 ^ 2 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * z2 ^ 2
    n[2,3] = n[3,2] = 4 * ((e * y3 * z3 + f * y2 * z2) * y1 * z1 + d * y2 * y3 * z2 * z3) + (h * y1*y1 + l * y2*y2 + c * y3 ^ 2) * z3 ^ 2 + (g * y2*y2 + h * y3 ^ 2 + a * y1*y1) * z1 ^ 2 + (g * y1*y1 + l * y3 ^ 2 + b * y2*y2) * z2 ^ 2
    n[3,3] = 4 * ((e * z3 ^ 2 + f * z2 ^ 2) * z1 ^ 2 + d * z2 ^ 2 * z3 ^ 2) + (h * z1 ^ 2 + l * z2 ^ 2 + c * z3 ^ 2) * z3 ^ 2 + (g * z2 ^ 2 + h * z3 ^ 2 + a * z1 ^ 2) * z1 ^ 2 + (g * z1 ^ 2 + l * z3 ^ 2 + b * z2 ^ 2) * z2 ^ 2
    n[1,4] = n[1,4] = 2 * (((x2 * z3 + x3 * z2) * e * x1 + (y2 * z3 + y3 * z2) * d * x2) * x3 + (x2 * y3 + x3 * y2) * f * x1 * x2) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * y3 * z3 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * y1 * z1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * y2 * z2
    n[2,4] = n[4,2] = 2 * (((x2 * z3 + x3 * z2) * e * y1 + (y2 * z3 + y3 * z2) * d * y2) * y3 + (x2 * y3 + x3 * y2) * f * y1 * y2) + (h * y1*y1 + l * y2*y2 + c * y3 ^ 2) * y3 * z3 + (g * y2*y2 + h * y3 ^ 2 + a * y1*y1) * y1 * z1 + (g * y1*y1 + l * y3 ^ 2 + b * y2*y2) * y2 * z2
    n[3,4] = n[4,3] = 2 * (((x2 * z3 + x3 * z2) * e * z1 + (y2 * z3 + y3 * z2) * d * z2) * z3 + (x2 * y3 + x3 * y2) * f * z1 * z2) + (h * z1 ^ 2 + l * z2 ^ 2 + c * z3 ^ 2) * y3 * z3 + (g * z2 ^ 2 + h * z3 ^ 2 + a * z1 ^ 2) * y1 * z1 + (g * z1 ^ 2 + l * z3 ^ 2 + b * z2 ^ 2) * y2 * z2
    n[4,4] = (x2 * z3 + x3 * z2) ^ 2 * e + (y2 * z3 + y3 * z2) ^ 2 * d + (x2 * y3 + x3 * y2) ^ 2 * f + (h * y1 * z1 + l * y2 * z2 + c * y3 * z3) * y3 * z3 + (g * y2 * z2 + h * y3 * z3 + a * y1 * z1) * y1 * z1 + (g * y1 * z1 + l * y3 * z3 + b * y2 * z2) * y2 * z2
    n[1,5] = n[5,1] = 2 * (((x1 * z3 + x3 * z1) * e * x1 + (y1 * z3 + y3 * z1) * d * x2) * x3 + (x1 * y3 + x3 * y1) * f * x1 * x2) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * x3 * z3 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * x1 * z1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * x2 * z2
    n[2,5] = n[5,2] = 2 * (((x1 * z3 + x3 * z1) * e * y1 + (y1 * z3 + y3 * z1) * d * y2) * y3 + (x1 * y3 + x3 * y1) * f * y1 * y2) + (h * y1*y1 + l * y2*y2 + c * y3 ^ 2) * x3 * z3 + (g * y2*y2 + h * y3 ^ 2 + a * y1*y1) * x1 * z1 + (g * y1*y1 + l * y3 ^ 2 + b * y2*y2) * x2 * z2
    n[3,5] = n[5,3] = 2 * (((x1 * z3 + x3 * z1) * e * z1 + (y1 * z3 + y3 * z1) * d * z2) * z3 + (x1 * y3 + x3 * y1) * f * z1 * z2) + (h * z1 ^ 2 + l * z2 ^ 2 + c * z3 ^ 2) * x3 * z3 + (g * z2 ^ 2 + h * z3 ^ 2 + a * z1 ^ 2) * x1 * z1 + (g * z1 ^ 2 + l * z3 ^ 2 + b * z2 ^ 2) * x2 * z2
    n[4,5] = n[5,4] = (x1 * z3 + x3 * z1) * (x2 * z3 + x3 * z2) * e + (y1 * z3 + y3 * z1) * (y2 * z3 + y3 * z2) * d + (x1 * y3 + x3 * y1) * (x2 * y3 + x3 * y2) * f + (h * y1 * z1 + l * y2 * z2 + c * y3 * z3) * x3 * z3 + (g * y2 * z2 + h * y3 * z3 + a * y1 * z1) * x1 * z1 + (g * y1 * z1 + l * y3 * z3 + b * y2 * z2) * x2 * z2
    n[5,5] = (x1 * z3 + x3 * z1) ^ 2 * e + (y1 * z3 + y3 * z1) ^ 2 * d + (x1 * y3 + x3 * y1) ^ 2 * f + (h * x1 * z1 + l * x2 * z2 + c * x3 * z3) * x3 * z3 + (g * x2 * z2 + h * x3 * z3 + a * x1 * z1) * x1 * z1 + (g * x1 * z1 + l * x3 * z3 + b * x2 * z2) * x2 * z2
    n[1,6] = n[6,1] = 2 * (((x1 * z2 + x2 * z1) * e * x1 + (y1 * z2 + y2 * z1) * d * x2) * x3 + (x1 * y2 + x2 * y1) * f * x1 * x2) + (h * x1*x1 + l * x2 ^ 2 + c * x3 ^ 2) * x3 * y3 + (g * x2 ^ 2 + h * x3 ^ 2 + a * x1*x1) * x1 * y1 + (g * x1*x1 + l * x3 ^ 2 + b * x2 ^ 2) * x2 * y2
    n[2,6] = n[6,2] = 2 * (((x1 * z2 + x2 * z1) * e * y1 + (y1 * z2 + y2 * z1) * d * y2) * y3 + (x1 * y2 + x2 * y1) * f * y1 * y2) + (h * y1*y1 + l * y2*y2 + c * y3 ^ 2) * x3 * y3 + (g * y2*y2 + h * y3 ^ 2 + a * y1*y1) * x1 * y1 + (g * y1*y1 + l * y3 ^ 2 + b * y2*y2) * x2 * y2
    n[3,6] = n[6,3] = 2 * (((x1 * z2 + x2 * z1) * e * z1 + (y1 * z2 + y2 * z1) * d * z2) * z3 + (x1 * y2 + x2 * y1) * f * z1 * z2) + (h * z1 ^ 2 + l * z2 ^ 2 + c * z3 ^ 2) * x3 * y3 + (g * z2 ^ 2 + h * z3 ^ 2 + a * z1 ^ 2) * x1 * y1 + (g * z1 ^ 2 + l * z3 ^ 2 + b * z2 ^ 2) * x2 * y2
    n[4,6] = n[6,4] = (x1 * z2 + x2 * z1) * (x2 * z3 + x3 * z2) * e + (y1 * z2 + y2 * z1) * (y2 * z3 + y3 * z2) * d + (x1 * y2 + x2 * y1) * (x2 * y3 + x3 * y2) * f + (h * y1 * z1 + l * y2 * z2 + c * y3 * z3) * x3 * y3 + (g * y2 * z2 + h * y3 * z3 + a * y1 * z1) * x1 * y1 + (g * y1 * z1 + l * y3 * z3 + b * y2 * z2) * x2 * y2
    n[5,6] = n[6,5] = (x1 * z2 + x2 * z1) * (x1 * z3 + x3 * z1) * e + (y1 * z2 + y2 * z1) * (y1 * z3 + y3 * z1) * d + (x1 * y2 + x2 * y1) * (x1 * y3 + x3 * y1) * f + (h * x1 * z1 + l * x2 * z2 + c * x3 * z3) * x3 * y3 + (g * x2 * z2 + h * x3 * z3 + a * x1 * z1) * x1 * y1 + (g * x1 * z1 + l * x3 * z3 + b * x2 * z2) * x2 * y2
    n[6,6] = (x1 * z2 + x2 * z1) ^ 2 * e + (y1 * z2 + y2 * z1) ^ 2 * d + (x1 * y2 + x2 * y1) ^ 2 * f + (h * x1 * y1 + l * x2 * y2 + c * x3 * y3) * x3 * y3 + (g * x2 * y2 + h * x3 * y3 + a * x1 * y1) * x1 * y1 + (g * x1 * y1 + l * x3 * y3 + b * x2 * y2) * x2 * y2
end

tensor2voigt(𝓒) = [[𝓒.η11[i] 𝓒.η13[i] 𝓒.η15[i]
                    𝓒.η13[i] 𝓒.η33[i] 𝓒.η35[i]
                    𝓒.η15[i] 𝓒.η35[i] 𝓒.η55[i]] for i in CartesianIndices(𝓒.η11)]

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
using Test
using ServiceMesh
using ReinforcementLearningBase

@testset "ReinforcementLearning interface test" begin
    env = ServiceMeshEnv(
        microservices = 3,
        maxinstances = 100,
        maxqueue = 10,
        steplength = 1,
        boottime = 3,
        instancecost = 0.05,
        jobtypes = [
            JobParams(1, [1, 2, 3], [2, 3, 5], ones(3), 100, 10.0, SingleArrival(5))
            JobParams(2, [1, 3], [5, 1], ones(2), 100, 10.0, SingleArrival(15))
        ],
        seed = 37,
    )

    @test RLBase.test_interfaces!(env)
end

@testset "Test Propagation" begin
    env = ServiceMeshEnv(
        microservices = 3,
        maxinstances = 100,
        maxqueue = 10,
        steplength = 1,
        boottime = 3,
        instancecost = 0.05,
        jobtypes = [
            JobParams(1, [1, 2, 3], [2, 3, 5], ones(3), 100, 10.0, SingleArrival(5))
            JobParams(2, [1, 3], [5, 1], ones(2), 100, 10.0, SingleArrival(15))
        ],
        seed = 37,
    )

    action = fill(5, 3)

    @test state(env) = []
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
struct SmoothAnalysis
    options::Dict{Symbol, Any}
end

SmoothAnalysis(; kwargs...) = SmoothAnalysis(Dict{Symbol, Any}(kwargs))

function (l::SmoothAnalysis)(le::Entry)
    return splitapply(le) do entry
        options = copy(l.options)
        npoints = pop!(l.options, :npoints, 200)
        labels, mappings = map(getlabel, entry.mappings), map(getvalue, entry.mappings)
        x, y = mappings.positional
        min, max = extrema(x)
        min < max || return Entry[]
        model = Loess.loess(Float64.(x), Float64.(y); options...)
        x̂ = collect(range(min, max, length=npoints))
        ŷ = Loess.predict(model, x̂)
        labeled_result = map(Labeled, labels.positional, [x̂, ŷ])
        default_plottype = Lines
        return Entry(
            AbstractPlotting.plottype(entry.plottype, default_plottype),
            Arguments(labeled_result),
            entry.attributes
        )
    end
end

"""
    smooth(span=0.75, degreee=2)

Fit a loess model. `span` is the degree of smoothing, typically in `[0,1]`.
Smaller values result in smaller local context in fitting.
`degree` is the polynomial degree used in the loess model.
"""
smooth(; kwargs...) = Layer((SmoothAnalysis(; kwargs...),))


text:
@testset "Testing LRU Cache" begin
    # Test `Priority`
    x = CachedArrays.Priority(0, 10)
    y = CachedArrays.Priority(1, 20)
    z = CachedArrays.Priority(0, 30)

    @test (CachedArrays.unwrap).((x, y, z)) == (10, 20, 30)
    @test x < y
    @test z < y
    @test !(z < x)
    @test !(x < z)
    @test z == x
    @test z != y
    @test hash(z) == hash(x)
    @test hash(z) != hash(y)

    # Test `LRU`
    lru = CachedArrays.LRU{Int}()
    @test isempty(lru)

    # LRU is a minheap, so should always return the item with the lowest priority.
    push!(lru, 1, 10)
    push!(lru, 2, 20)
    @test length(lru) == 2
    @test pop!(lru) == 1
    @test length(lru) == 1
    @test first(lru) == 2
    push!(lru, 3, 10)
    @test first(lru) == 3
    @test length(lru) == 2

    # Now, update "3" to have a higher priority.
    CachedArrays.update!(lru, 3, 30)
    @test length(lru) == 2
    @test first(lru) == 2
    @test pop!(lru) == 2
    @test pop!(lru) == 3
    @test isempty(lru)

    # Finally, test that deletion works.
    push!(lru, 3, 30)
    push!(lru, 1, 10)
    push!(lru, 2, 20)
    @test length(lru) == 3
    delete!(lru, 2)
    @test length(lru) == 2
    @test pop!(lru) == 1
    @test pop!(lru) == 3
    @test isempty(lru)
end


text:
using CompScienceMeshes, BEAST

Γ = readmesh(Pkg.dir("BEAST", "examples", "sphere.in"))
X, Y = raviartthomas(Γ), buffachristiansen(Γ)

κ = ω = 1.0; γ = κ*im
T = MWSingleLayer3D(γ)
N = NCross()

Txx = sdata(assemble(T,X,X))
Tyy = sdata(assemble(T,Y,Y))
Nxy = sdata(assemble(N,X,Y))

iNxy = inv(Nxy)
A = iNxy' * Tyy * iNxy * Txx

@show cond(Txx)
@show cond(Tyy)
@show cond(A)

wx, wy, wp = eigvals(Txx), eigvals(Tyy), eigvals(A)


text:
function read_fa_wms_tr(path::String)
    wms=Vector{Matrix{Float64}}()
    wm=zeros(1,4)
    f=open(path)
    for line in eachline(f)
        prefix=line[1:2]
        prefix == "01" && (wm=transpose([parse(Float64,i) for i in split(line)[2:end]]))
        prefix != "01" && prefix != "NA" && prefix != "PO" && prefix != "//" && (wm=vcat(wm, transpose([parse(Float64,i) for i in split(line)[2:end]])))
        prefix == "//" && push!(wms, wm)
    end
    return wms
end

#wm_samples are in decimal probability space, not log space
function assemble_source_priors(no_sources::Integer, wm_samples::AbstractVector{<:AbstractMatrix{<:AbstractFloat}}, prior_wt::AbstractFloat=PRIOR_WT) #estimate a dirichlet prior on wm_samples inputs; if the number of samples is lower than the number of sources, return a false bool for init and permutation functions
    source_priors = Vector{Union{Vector{Dirichlet{Float64}},Bool}}()
    for source in 1:no_sources
        if source <= length(wm_samples)
            push!(source_priors, estimate_dirichlet_prior_on_wm(wm_samples[source], prior_wt))
        else
            push!(source_priors, false)
        end
    end
    return source_priors
end

            function estimate_dirichlet_prior_on_wm(wm::AbstractMatrix{<:AbstractFloat}, wt::AbstractFloat=PRIOR_WT)
                for i in 1:size(wm)[1]
                    !(isprobvec(wm[i,:])) && throw(DomainError("Bad weight vec supplied to estimate_dirichlet_prior_on_wm! $(wm[i,:])"))
                end
                prior = Vector{Dirichlet{Float64}}()
                for position in 1:size(wm)[1]
                    normvec=wm[position,:]
                    zero_idxs=findall(isequal(0.),wm[position,:])
                    normvec[zero_idxs].+=10^-99
                    push!(prior, Dirichlet(normvec.*wt))
                end
                return prior
            end

function cluster_mix_prior!(df::DataFrame, wms::AbstractVector{<:AbstractMatrix{<:AbstractFloat}})
    mix=falses(size(df,1),length(wms))
    for (o, row) in enumerate(eachrow(df))
        row.cluster != 0 && (mix[o,row.cluster]=true)
    end

    represented_sources=unique(df.cluster)
    return mix[:,represented_sources]
end

function infocenter_wms_trim(wm::AbstractMatrix{<:AbstractFloat}, trimsize::Integer)
    !(size(wm,2)==4) && throw(DomainError("Bad wm! 2nd dimension should be size 4"))
    infovec=get_pwm_info(wm, logsw=false)
    maxval, maxidx=findmax(infovec)
    upstream_extension=Int(floor((trimsize-1)/2))
    downstream_extension=Int(ceil((trimsize-1)/2))
    1+upstream_extension+downstream_extension > size(wm,1) && throw(DomainError("Src too short for trim! $upstream_extension $downstream_extension"))
    return wm[max(1,maxidx-upstream_extension):min(maxidx+downstream_extension,size(wm,1)),:]
end

function filter_priors(target_src_no::Integer, target_src_size::Integer, prior_wms::AbstractVector{<:AbstractMatrix{<:AbstractFloat}}, prior_mix::BitMatrix)
    wms=Vector{Matrix{Float64}}(undef, target_src_no)
    freqsort_idxs=sortperm([sum(prior_mix[:,s]) for s in 1:length(prior_wms)])
    for i in 1:target_src_no
        target_src_idx=freqsort_idxs[i]
        wms[i]=infocenter_wms_trim(prior_wms[target_src_idx], target_src_size)
    end
    return wms
end

function combine_filter_priors(target_src_no::Integer, target_src_size::Integer, prior_wms::Tuple{<:AbstractVector{<:AbstractMatrix{<:AbstractFloat}},AbstractVector{<:AbstractMatrix{<:AbstractFloat}}}, prior_mix::Tuple{BitMatrix,BitMatrix})
    wms=Vector{Matrix{Float64}}(undef, target_src_no)
    cat_wms=vcat(prior_wms[1],prior_wms[2])
    first_freq=[sum(prior_mix[1][:,s]) for s in 1:length(prior_wms[1])]
    second_freq=[sum(prior_mix[2][:,s]) for s in 1:length(prior_wms[2])]
    freqsort_idxs=sortperm(vcat(first_freq,second_freq))
    for i in 1:target_src_no
        target_src_idx=freqsort_idxs[i]
        wms[i]=infocenter_wms_trim(cat_wms[target_src_idx], target_src_size)
    end
    return wms
end