text:
################################################################################
# Implement abstract type of CutPruner
################################################################################
export AbstractCutPruningAlgo
export CutPruner, AbstractCutPruner
# High-level functions
export addcuts!, ncuts
# Low-level functions
export appendcuts!, replacecuts!, keeponlycuts!, removecuts!

abstract type AbstractCutPruningAlgo end

"""
A cut pruner maintains a matrix `A` and a vector `b` such that
represents `size(A, 1)` (` == length(b)`) cuts.
Let `a_i` be `-A[i,:]` if `lazy_minus` and `A[i,:]` otherwise and `β_i` be `b[i]`, the meaning of the cut depends on the sense.
If `sense` is
* `:Min`, then the cut pruner represents the concave polyhedral function `min ⟨a_i, x⟩ + β_i`;
* `:Max`, then the cut pruner represents the convex polyhedral function `max ⟨a_i, x⟩ + β_i`;
* `:≤`, then the cut pruner represents the polyhedra defined by the intersection of the half-space `⟨a_i, x⟩ ≤ β_i`;
* `:≥`, then the cut pruner represents the polyhedra defined by the intersection of the half-space `⟨a_i, x⟩ ≥ β_i`.

Internally, instead of `sense`, the booleans `isfun` and `islb` are stored.
The mapping between `sense` and these two booleans is given by the following table

| `sense` | `isfun` | `islb` |
| ------- | ------- | ------ |
| Min     | true    | false  |
| Max     | true    | true   |
| ≤       | false   | false  |
| ≥       | false   | true   |

"""
abstract type AbstractCutPruner{N, T} end

function gettype(sense::Symbol)
    if sense == :Min
        true, false
    elseif sense == :Max
        true, true
    elseif sense == :≤
        false, false
    elseif sense == :≥
        false, true
    else
        throw(ArgumentError("Invalid value `$sense' for sense. It should be :Min, :Max, :≤ or :≥."))
    end
end
function getsense(isfun::Bool, islb::Bool)
    if isfun
        islb ? :Max : :Min
    else
        islb ? :≥ : :≤
    end
end
getsense(pruner::AbstractCutPruner) = getsense(isfun(pruner), islb(pruner))
isfun(pruner::AbstractCutPruner) = pruner.isfun
islb(pruner::AbstractCutPruner) = pruner.islb

struct CutPruner{N, T} end

"""Return whether the CutPruner `man` has any cut."""
function Base.isempty(man::AbstractCutPruner)
    return isempty(man.b)
end

"""Return number of cuts in CutPruner `man`."""
function ncuts(man::AbstractCutPruner)
    return length(man.b)
end

# COMPARISON
hastrust(man::AbstractCutPruner) = true
hasterritories(man::AbstractCutPruner) = false
"""Get current `trust` of CutPruner `man`."""
gettrust(man::AbstractCutPruner) = man.trust

function _indmin(a::Vector, tiebreaker::Vector)
    imin = 1
    for i in 2:length(a)
        if a[i] < a[imin] || (a[i] == a[imin] && tiebreaker[i] < tiebreaker[imin])
            imin = i
        end
    end
    imin
end

"""Remove `num` cuts in CutPruner `man`."""
function choosecutstoremove(man::AbstractCutPruner, num::Int)
    # MergeSort is stable so in case of equality, the oldest cut loose
    # However PartialQuickSort is a lot faster

    trust = gettrust(man)
    if num == 1
        [_indmin(trust, man.ids)]
    else
        # /!\ PartialQuickSort is unstable, here it does not matter
        function _lt(i, j)
            # If cuts have same trust, remove oldest cut
            if trust[i] == trust[j]
                man.ids[i] < man.ids[j]
            # Else, remove cuts with lowest trust
            else
                trust[i] < trust[j]
            end
        end
        # Return index of `num` cuts with lowest trusts
        sort(1:length(trust), alg=PartialQuickSort(num), lt=_lt)[1:num]
    end
end

"""Test if cut `i` is better than `newcuttrust`."""
isbetter(man::AbstractCutPruner, i::Int, mycut::Bool) = gettrust(man)[i] > initialtrust(man, mycut)

# CHANGE

function getnreplaced(man::AbstractCutPruner, R, ncur, nnew, nmycut)
    # Check if some new cuts should be ignored
    take = man.maxncuts - ncur
    # get number of cuts to remove
    nreplaced = length(R)
    # Start:
    # |      | R      j|
    # | take |         |
    # |  mycut  |      |
    # End:
    # |      | R   j   |
    # | take       |   |
    # |  mycut  |      |
    while take + length(R) - nreplaced < nnew
        # I first try to see if the nmycut cuts generated by me can be taken
        # because they have better trust than the nnew-nmycut others
        if isbetter(man, R[nreplaced], take < nmycut)
            nreplaced -= 1
        else
            take += 1
        end
    end
    take, nreplaced
end

# Add cuts Ax >= b
# If mycut then the cut has been added because of one of my trials
function addcuts!(man::AbstractCutPruner{N, T},
                  A::AbstractMatrix{T},
                  b::AbstractVector{T},
                  mycut::AbstractVector{Bool}) where {N, T}
    # get current number of cuts:
    ncur = ncuts(man)
    ncutinitial = size(A, 1)

    # First pass: remove redundancy among new cuts
    # check redundancy among new cuts (require more than one cut in Anew)
    redundant1 = checkredundancy(A, b, man.isfun, man.islb, man.TOL_EPS)
    tokeep1 = setdiff(collect(1:ncutinitial), redundant1)
    if length(tokeep1) < size(A, 1)
        A = A[tokeep1, :]
        b = b[tokeep1]
        mycut = mycut[tokeep1]
    end

    # Second pass: remove redundancy between old and new cuts
    # if specify, check if new cuts are redundants with old cuts
    nincumbents = size(A, 1)
    if man.excheck
        redundants2 = checkredundancy(man.A, man.b, A, b, man.isfun, man.islb, man.TOL_EPS)
        tokeep2 = setdiff(collect(1:nincumbents), redundants2)
        if !isempty(redundants2)

            # if all cuts are redundants, then do nothing:
            if length(tokeep2) == 0
                return zeros(Int, nincumbents)
            end
            A = A[tokeep2, :]
            b = b[tokeep2]
            mycut = mycut[tokeep2]
        end

        redundants = vcat(redundant1, tokeep1[redundants2])
        tokeep = tokeep1[tokeep2]
    else
        redundants = redundant1
        tokeep = tokeep1
    end

    # get number of new cuts in A:
    nnew = size(A, 1)
    @assert length(mycut) == length(b) == nnew
    @assert nnew > 0

    if man.maxncuts == -1 || ncur + nnew <= man.maxncuts
        # If enough room, just append cuts
        status = (VERSION < v"0.7-" ? ncur + (1:nnew) : ncur .+ (1:nnew))
        appendcuts!(man, A, b, mycut)
    else
        # Otherwise, we need need to remove some cuts

        # get indexes of cuts with lowest trusts:
        R = choosecutstoremove(man, ncur + nnew - man.maxncuts)

        nmycut = sum(mycut)
        take, nreplaced = getnreplaced(man, R, ncur, nnew, nmycut)

        # Cuts that will be replaced
        R = @view R[1:nreplaced]
        # Nowe we split A, b into
        # * A, b   : cuts to be pushed
        # * Ar, br : cuts to replace old cuts
        # * _, _   : cuts removed
        if nreplaced == nnew
            status = R
            Ar = A
            br = b
            mycutr = mycut
            A = similar(A, 0, size(A, 2))
            b = similar(b, 0)
            mycut = similar(mycut, 0)
        else
            status = zeros(Int, nnew)
            if take < nnew
                # Remove ignored cuts
                takemycut = min(take, nmycut)
                takenotmycut = take - takemycut
                takeit = zeros(Bool, nnew)
                for i in 1:nnew
                    ok = false
                    if mycut[i]
                        if takemycut > 0
                            takemycut -= 1
                            ok = true
                        end
                    else
                        if takenotmycut > 0
                            takenotmycut -= 1
                            ok = true
                        end
                    end
                    if ok
                        takeit[i] = true
                    end
                end

                takeit = findall(takeit)
            else
                takeit = collect(1:nnew)
            end
            @assert take == length(takeit)
            replaced = takeit[1:nreplaced]
            pushed = takeit[(nreplaced+1):end]

            status[replaced] = R
            Ar = @view A[replaced,:]
            br = @view b[replaced]
            mycutr = @view mycut[replaced]
            status[pushed] = (VERSION < v"0.7-" ? ncur + (1:length(pushed)) : ncur .+ (1:length(pushed)))
            status[pushed] = ncur .+ (1:length(pushed))
            appendcuts!(man, (@view A[pushed,:]), (@view b[pushed]), (@view mycut[pushed]))
        end
        if nreplaced > 0
            man.A[R, :] .= Ar
            man.b[R] .= br
            replacecuts!(man, R, Ar, br, mycutr)
        end
    end

    if isempty(redundants)
        fullstatus = status
    else
        fullstatus = Array{Int}(undef, nnew + length(redundants))
        fullstatus[tokeep] .= status
        fullstatus[redundants] .= 0
    end
    fullstatus
end

# Low-level functions
function _keeponlycuts!(man::AbstractCutPruner, K::AbstractVector{Int})
    man.A = man.A[K, :]
    man.b = man.b[K]
    man.ids = man.ids[K]
    if hastrust(man)
        man.trust = gettrust(man)[K]
    end
    if hasterritories(man)
        man.territories = man.territories[K]
    end
end

"""Keep only cuts whose indexes are in Vector `K` in CutPruner `man`. If `K` is not sorted, the cuts will change their order accordingly."""
function keeponlycuts!(man::AbstractCutPruner, K::AbstractVector{Int})
    _keeponlycuts!(man, K)
end

"""Remove cuts whose indexes are in Vector `K` in CutPruner `man`."""
function removecuts!(man::AbstractCutPruner, K::AbstractVector{Int})
    keeponlycuts!(man, setdiff(1:ncuts(man), K))
end

function _replacecuts!(man::AbstractCutPruner, K::AbstractVector{Int}, A, b)
    man.A[K, :] = A
    man.b[K] = b
    man.ids[K] = newids(man, length(K))
end

"""Replace cuts at indexes in `K` by cuts in (A, b, mycut) in CutPruner `man`."""
function replacecuts!(man::AbstractCutPruner, K::AbstractVector{Int}, A, b, mycut::AbstractVector{Bool})
    _replacecuts!(man, K, A, b)
    man.trust[K] = initialtrusts(man, mycut)
end

function _appendcuts!(man::AbstractCutPruner, A, b)
    man.A = [man.A; A]
    man.b = [man.b; b]
    append!(man.ids, newids(man, length(b)))
end

"""Append cuts (A, b, mycut) in CutPruner `man`."""
function appendcuts!(man::AbstractCutPruner, A, b, mycut::AbstractVector{Bool})
    _appendcuts!(man, A, b)
    append!(man.trust, initialtrusts(man, mycut))
end

# Update trust
export addusage!, addposition!
function addusage!(pruner::AbstractCutPruner, position)
    # do nothing by default
end
function addposition!(pruner::AbstractCutPruner, position)
    # do nothing by default
end

# Unexported utilities

function newids(man::AbstractCutPruner, n::Int)
    (man.id+1):(man.id += n)
end

"""Get a Vector of Float64 specifying the initial trusts of `mycut`."""
function initialtrusts(man::AbstractCutPruner, mycut::AbstractVector{Bool})
    Float64[initialtrust(man, mc) for mc in mycut]
end


text:
@testset "Test LogisticFunction"  begin
    L = LogisticFunction(3)

    x = [1, 2, 3]
    z = z_buffer(L)
    dx = dx_buffer(L)

    evaluate!(L, x, z, dx=dx)

    @test z ≈ [0.7310585786, 0.8807970779, 0.9525741268]
    @test dx ≈ [0.1966119332 0 0;
                0 0.1049935854 0;
                0 0 0.0451766597]
end

@testset "Test HyperbolicTangent"  begin
    H = HyperbolicTangent(3)

    x = [1, 2, 3]
    z = z_buffer(H)
    dx = dx_buffer(H)

    evaluate!(H, x, z, dx=dx)

    @test z ≈ [0.7615941559, 0.9640275801, 0.9950547537]
    @test dx ≈ [0.4199743416 0 0;
                0 0.0706508249 0;
                0 0 0.0098660372]
end

@testset "Test AffineMap"  begin
    @testset "Full Input"  begin
        A = AffineMap([2, 2, 4], [1, 2, 3])

        x = [1, 2, 3]
        z = z_buffer(A)
        dx = dx_buffer(A)

        evaluate!(A, x, z, dx=dx)

        @test z ≈ [3.0, 6.0, 15.0]
        @test dx ≈ [2.0 0 0;
                    0 2.0 0;
                    0 0 4.0]
    end
    @testset "Default Bias"  begin
        A = AffineMap([2, 2, 4])

        x = [1, 2, 3]
        z = z_buffer(A)
        dx = dx_buffer(A)

        evaluate!(A, x, z, dx=dx)

        @test z ≈ [2.0, 4.0, 12.0]
        @test dx ≈ [2.0 0 0;
                    0 2.0 0;
                    0 0 4.0]
    end

    @testset "Argument Error"  begin
        @test_throws ArgumentError AffineMap([2, 2, 4],
                                             [1, 1, 1, 1])
    end
end

@testset "Test Identity"  begin
    H = Identity(3)

    x = [1, 2, 3]
    z = z_buffer(H)
    dx = dx_buffer(H)

    evaluate!(H, x, z, dx=dx)

    @test z ≈ [1.0, 2.0, 3.0]
    @test dx ≈ [1.0 0 0;
                0 1.0 0;
                0 0 1.0]
end

ext: jl


text:
"""
    D1matrix(nvertices)
Returns the first difference of dimension `nvercies-1 x nvertices`.
It is important to note that the numerical first derivative of a
vector produces a vector with 1 component less than the input vector.

#Examples
```julia julia-repl
julia> D1matrix(3)
2×3 Matrix{Float64}:
 -1.0   1.0  0.0
  0.0  -1.0  1.0
```
"""
function D1matrix(nvertices::Int)
    D1matrix = zeros(Float64,nvertices-1,nvertices)
    for i=1:nvertices-1
        D1matrix[i,i] = -1.
        D1matrix[i,i+1] = 1.
    end
    return D1matrix
end

"""
    D2matrix(nvertices)
Returns the second difference of dimension `nvercies-2 x nvertices`.
It is important to note that the numerical second derivative of a
vector produces a vector with 2 component less than the input vector.

#Examples
```julia julia-repl
julia> D2matrix(4)
2×4 Matrix{Float64}:
 1.0  -2.0   1.0  0.0
 0.0   1.0  -2.0  1.0
```
"""
function D2matrix(nvertices::Int)

    D2matrix = zeros(Float64,nvertices-2,nvertices)

    for i=1:nvertices-2
        D2matrix[i,i] = 1.
        D2matrix[i,i+1] = -2.
        D2matrix[i,i+2] = 1.
    end

    return D2matrix
end

"""
    PD2matrices2dim(nverticesRows, nverticesCols)
Given as inputs the number of vertices associated to the rows and
the number of vertices associated to the columns, the function
returns the penalization matrix for rows, PRows and the penalization
matrix for columns, PCols. The penalization with this matrices would
be: λrows*PRows + λcols*Pcols

#Examples
```julia julia-repl
julia> PRows, PCols = PD2matrices2dim(10,10);
julia> PRows
100×100 Matrix{Float64}:
  1.0  -2.0   1.0   0.0   0.0  0.0   0.0   0.0   0.0   0.0  …   0.0  -0.0   0.0   0.0   0.0  0.0   0.0   0.0   0.0   0.0
 -2.0   5.0  -4.0   1.0   0.0  0.0   0.0   0.0   0.0   0.0     -0.0   0.0  -0.0   0.0   0.0  0.0   0.0   0.0   0.0   0.0
  1.0  -4.0   6.0  -4.0   1.0  0.0   0.0   0.0   0.0   0.0      0.0  -0.0   0.0  -0.0   0.0  0.0   0.0   0.0   0.0   0.0
  0.0   1.0  -4.0   6.0  -4.0  1.0   0.0   0.0   0.0   0.0      0.0   0.0  -0.0   0.0  -0.0  0.0   0.0   0.0   0.0   0.0
  ⋮                            ⋮                            ⋱   ⋮                            ⋮
  0.0   0.0   0.0   0.0   0.0  0.0  -0.0   0.0  -0.0   0.0      0.0   0.0   0.0   0.0   0.0  1.0  -4.0   6.0  -4.0   1.0
  0.0   0.0   0.0   0.0   0.0  0.0   0.0  -0.0   0.0  -0.0      0.0   0.0   0.0   0.0   0.0  0.0   1.0  -4.0   5.0  -2.0
  0.0   0.0   0.0   0.0   0.0  0.0   0.0   0.0  -0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0   0.0   1.0  -2.0   1.0

julia> PCols
100×100 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  -2.0  -0.0  …   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  -0.0  -2.0      0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  -0.0  -0.0      0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  -0.0  -0.0      0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮                         ⋮          ⋱        ⋮                        ⋮
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0     -0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0     -0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0     -2.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
```
"""
function PD2matrices2dim(nverticesRows::Int, nverticesCols::Int)

    D2Rows = D2matrix(nverticesRows)
    D2Cols = D2matrix(nverticesCols)

    IRows = (1.0*I)(nverticesRows)
    ICols = (1.0*I)(nverticesCols)

    PRows = kronecker(ICols, D2Rows)'* kronecker(ICols, D2Rows)
    PCols = kronecker(D2Cols,IRows)'* kronecker(D2Cols,IRows)
    return Matrix(PRows), Matrix(PCols)
end

"""
    D3matrix(nvertices)
Returns the second difference of dimension `nvercies-3 x nvertices`.
It is important to note that the numerical second derivative of a
vector produces a vector with 2 component less than the input vector.

#Examples
```julia julia-repl
julia> D3matrix(5)
2×5 Matrix{Float64}:
 -1.0   3.0  -3.0   1.0  0.0
  0.0  -1.0   3.0  -3.0  1.0
```
"""
function D3matrix(nvertices::Int)

    D3matrix = zeros(Float64,nvertices-3,nvertices)

    for i=1:nvertices-3
        D3matrix[i,i] = -1.
        D3matrix[i,i+1] = 3.
        D3matrix[i,i+2] = -3.
        D3matrix[i,i+3] = 1.
    end

    return D3matrix
end

"""
    Ω1matrix(nvertices)
Returns a matrix full of zeros with the appropiate dimensions to apply
a selective penalization to the first differences vector. The default
0 value, implies no penalization on the first derivative.

#Examples
```julia julia-repl
julia> Ω1matrix(5)
4×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
```
"""
function Ω1matrix(nvertices)
    return zeros(Float64, nvertices - 1, nvertices - 1)
end

"""
    Ω2matrix(nvertices)
Returns a matrix full of zeros with the appropiate dimensions to apply
a selective penalization to the second differences vector. The default
0 value, implies no penalization on the second derivative.

#Examples
```julia julia-repl
julia> Ω2matrix(5)
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
```
"""
function Ω2matrix(nvertices)
    return zeros(Float64, nvertices - 2, nvertices - 2)
end

"""
    Ω3matrix(nvertices)
Returns a matrix full of zeros with the appropiate dimensions to apply
a selective penalization to the third differences vector. The default
0 value, implies no penalization on the third derivative.

#Examples
```julia julia-repl
julia> Ω3matrix(5)
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0
```
"""
function Ω3matrix(nvertices)
    return zeros(Float64, nvertices - 3, nvertices - 3)
end

"""
    Wmatrix(nvertices)
Given the total number of experimental points,
build the weight matrix for that number of points.
Note that by default all the weights are equal to 1.
#Examples
```julia julia-repl
julia> Wmatrix(5)
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0
```
"""
function Wmatrix(nexp)
    W = zeros(Float64, nexp, nexp)
    W[diagind(W)] .= 1.
    return W
end


text:
@testset "InferenceData" begin
    data = load_arviz_data("centered_eight")

    @testset "construction" begin
        pydata = PyObject(data)
        @test InferenceData(pydata) isa InferenceData
        @test PyObject(InferenceData(pydata)) === pydata
        @test InferenceData(data) === data
        data2 = InferenceData(; posterior = data.posterior)
        @test data2 isa InferenceData
        @test :posterior in propertynames(data2)
    end

    @testset "properties" begin
        @test :posterior in propertynames(data)
        @test length(propertynames(data)) > 1
        data3 = InferenceData(; posterior = data.posterior, prior = data.prior)
        @test :prior in propertynames(data3)
        delete!(data3, :prior)
        @test :prior ∉ propertynames(data3)
    end

    @testset "conversion" begin
        @test pyisinstance(PyObject(data), ArviZ.arviz.InferenceData)
        data4 = convert(InferenceData, PyObject(data))
        @test data4 isa InferenceData
        @test PyObject(data4) === PyObject(data)
    end
end

@testset "+(::InferenceData, ::InferenceData)" begin
    rng = MersenneTwister(42)
    idata1 = from_dict(posterior = Dict(
        "A" => randn(rng, 2, 10, 2),
        "B" => randn(rng, 2, 10, 5, 2),
    ))
    idata2 = from_dict(prior = Dict(
        "C" => randn(rng, 2, 10, 2),
        "D" => randn(rng, 2, 10, 5, 2),
    ))

    new_idata = idata1 + idata2
    @test new_idata isa InferenceData
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]),
        new_idata,
    ) == []
end

@testset "concat" begin
    rng = MersenneTwister(42)
    idata1 = from_dict(posterior = Dict(
        "A" => randn(rng, 2, 10, 2),
        "B" => randn(rng, 2, 10, 5, 2),
    ))
    idata2 = from_dict(prior = Dict(
        "C" => randn(rng, 2, 10, 2),
        "D" => randn(rng, 2, 10, 5, 2),
    ))

    new_idata = concat(idata1, idata2)
    @test new_idata isa InferenceData
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]),
        new_idata,
    ) == []

    new_idata = concat!(idata1, idata2)
    @test new_idata === idata1
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]),
        idata1,
    ) == []
end

@testset "convert_to_inference_data" begin
    rng = MersenneTwister(42)

    @testset "convert_to_inference_data(::Dict)" begin
        dataset = Dict("A" => randn(rng, 2, 10, 2), "B" => randn(rng, 2, 10, 5, 2))
        idata1 = convert_to_inference_data(dataset)
        @test idata1 isa InferenceData
        @test check_multiple_attrs(Dict(:posterior => ["A", "B"]), idata1) == []
    end

    @testset "convert_to_inference_data(::Array)" begin
        arr = randn(rng, 2, 10, 2)
        idata2 = convert_to_inference_data(arr)
        @test check_multiple_attrs(Dict(:posterior => ["x"]), idata2) == []
    end
end

@testset "from_dict" begin
    rng = MersenneTwister(42)

    posterior = Dict("A" => randn(rng, 2, 10, 2), "B" => randn(rng, 2, 10, 5, 2))
    prior = Dict("C" => randn(rng, 2, 10, 2), "D" => randn(rng, 2, 10, 5, 2))

    idata = from_dict(posterior; prior = prior)
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]),
        idata,
    ) == []

    idata2 = from_dict(; prior = prior)
    @test check_multiple_attrs(Dict(:prior => ["C", "D"]), idata2) == []
end

@testset "netcdf roundtrip" begin
    data = load_arviz_data("centered_eight")
    mktemp() do path, _
        to_netcdf(data, path)
        data2 = from_netcdf(path)
        @test data2 isa InferenceData
        @test propertynames(data) == propertynames(data2)
    end
end

@testset "load_arviz_data" begin
    data = load_arviz_data("centered_eight")
    @test data isa InferenceData

    datasets = load_arviz_data()
    @test datasets isa Dict
end


text:
#Reference:
#须寅, 龙驭球, 龙志飞等. 引入泡状位移含旋转自由度的广义协调三角形膜元[J]. 工程力学, 2000, 17(3): 1-9.
function K_GT9(elm::Tria)::Matrix{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    K=Matrix{Float64}(undef,12,12)
    A=elm.A
    xᵢ=[x₁,x₂,x₃]
    yᵢ=[y₁,y₂,y₃]
    a₁,a₂,a₃=[xᵢ[j]*yᵢ[m]-xᵢ[m]*yᵢ[j] for (j,m) in zip([2,3,1],[3,1,2])]
    b₁,b₂,b₃=[yᵢ[j]-yᵢ[m] for (j,m) in zip([2,3,1],[3,1,2])]
    c₁,c₂,c₃=[-xᵢ[j]+xᵢ[m] for (j,m) in zip([2,3,1],[3,1,2])]
    function BtDB(x)
        x,y=x[1],x[2]
        K[1,1]=0.25*b₁^2*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₁^2*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[1,2]=0.25*b₁*c₁*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₁*c₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[1,3]=0.125*b₁^2*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*c₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₁*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2))
        K[1,4]=0.25*b₂*b₁*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₁*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[1,5]=0.25*b₁*c₂*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*c₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[1,6]=0.125*b₂*b₁*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*c₁*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₁*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[1,7]=0.25*b₁*b₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₁*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[1,8]=0.25*b₁*c₃*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₃*c₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[1,9]=0.125*b₁*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*c₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₁*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[2,1]=0.25*b₁*c₁*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₁*c₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[2,2]=0.25*c₁^2*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₁^2*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[2,3]=0.125*c₁^2*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*b₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₁*c₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2))
        K[2,4]=0.125*b₁*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2)) + 0.25*b₂*c₁*E₀*ν₀/(A^2*(1 - ν₀^2))
        K[2,5]=0.25*c₁*c₂*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*b₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[2,6]=0.125*c₁*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*b₁*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₂*c₁*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[2,7]=0.125*b₁*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2)) + 0.25*b₃*c₁*E₀*ν₀/(A^2*(1 - ν₀^2))
        K[2,8]=0.25*c₁*c₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₁*b₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[2,9]=0.125*c₁*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*b₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₃*c₁*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[3,1]=0.5*b₁*(0.25*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*c₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[3,2]=0.5*c₁*(0.25*c₁*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*b₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*b₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[3,3]=0.25*b₁*(0.25*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)*(0.25*c₁*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*b₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + (1/2)*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
        K[3,4]=0.5*b₂*(0.25*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*c₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[3,5]=0.5*c₂*(0.25*c₁*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*b₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*b₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[3,6]=0.25*b₂*(0.25*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)*(0.25*c₁*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*b₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + (1/2)*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(1 - ν₀^2)
        K[3,7]=0.5*b₃*(0.25*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*c₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[3,8]=0.5*c₃*(0.25*c₁*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*b₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*b₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[3,9]=0.25*b₃*(0.25*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)*(0.25*c₁*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)) + 0.25*b₁*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A*(1 - ν₀^2)))/A + (1/2)*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*E₀*(1 - ν₀)/(1 - ν₀^2)
        K[4,1]=0.25*b₂*b₁*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₁*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[4,2]=0.125*b₁*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2)) + 0.25*b₂*c₁*E₀*ν₀/(A^2*(1 - ν₀^2))
        K[4,3]=0.125*b₂*b₁*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*c₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₂*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2))
        K[4,4]=0.25*b₂^2*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₂^2*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[4,5]=0.25*b₂*c₂*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[4,6]=0.125*b₂^2*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*c₂*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₂*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[4,7]=0.25*b₂*b₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₂*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[4,8]=0.25*b₂*c₃*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₃*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[4,9]=0.125*b₂*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*c₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₂*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[5,1]=0.25*b₁*c₂*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*c₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[5,2]=0.25*c₁*c₂*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*b₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[5,3]=0.125*c₁*c₂*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*b₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₁*c₂*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2))
        K[5,4]=0.25*b₂*c₂*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[5,5]=0.25*c₂^2*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₂^2*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[5,6]=0.125*c₂^2*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*b₂*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₂*c₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[5,7]=0.125*b₂*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2)) + 0.25*b₃*c₂*E₀*ν₀/(A^2*(1 - ν₀^2))
        K[5,8]=0.25*c₂*c₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*b₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[5,9]=0.125*c₂*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*b₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₃*c₂*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[6,1]=0.5*b₁*(0.25*b₂*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*c₁*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2))
        K[6,2]=0.5*c₁*(0.25*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*b₁*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2))
        K[6,3]=0.25*b₁*(0.25*b₂*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(0.25*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + (1/2)*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(1 - ν₀^2)
        K[6,4]=0.5*b₂*(0.25*b₂*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*c₂*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2))
        K[6,5]=0.5*c₂*(0.25*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*b₂*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2))
        K[6,6]=0.25*b₂*(0.25*b₂*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.25*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + (1/2)*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)^2*(1 - ν₀)/(1 - ν₀^2)
        K[6,7]=0.5*b₃*(0.25*b₂*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*c₃*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2))
        K[6,8]=0.5*c₃*(0.25*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*b₃*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2))
        K[6,9]=0.25*b₃*(0.25*b₂*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(0.25*c₂*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₂*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + (1/2)*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(1 - ν₀^2)
        K[7,1]=0.25*b₁*b₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₁*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[7,2]=0.125*b₁*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2)) + 0.25*b₃*c₁*E₀*ν₀/(A^2*(1 - ν₀^2))
        K[7,3]=0.125*b₁*b₃*E₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*c₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₃*c₁*E₀*ν₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2))
        K[7,4]=0.25*b₂*b₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₂*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[7,5]=0.125*b₂*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2)) + 0.25*b₃*c₂*E₀*ν₀/(A^2*(1 - ν₀^2))
        K[7,6]=0.125*b₂*b₃*E₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*c₃*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₃*c₂*E₀*ν₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[7,7]=0.25*b₃^2*E₀/(A^2*(1 - ν₀^2)) + 0.125*c₃^2*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[7,8]=0.25*b₃*c₃*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₃*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[7,9]=0.125*b₃^2*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*c₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₃*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[8,1]=0.25*b₁*c₃*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₃*c₁*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[8,2]=0.25*c₁*c₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₁*b₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[8,3]=0.125*c₁*c₃*E₀*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*b₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₁*c₃*E₀*ν₀*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/(A^2*(1 - ν₀^2))
        K[8,4]=0.25*b₂*c₃*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₃*c₂*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[8,5]=0.25*c₂*c₃*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₂*b₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[8,6]=0.125*c₂*c₃*E₀*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*b₃*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₂*c₃*E₀*ν₀*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[8,7]=0.25*b₃*c₃*E₀*ν₀/(A^2*(1 - ν₀^2)) + 0.125*b₃*c₃*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[8,8]=0.25*c₃^2*E₀/(A^2*(1 - ν₀^2)) + 0.125*b₃^2*E₀*(1 - ν₀)/(A^2*(1 - ν₀^2))
        K[8,9]=0.125*c₃^2*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2)) + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*b₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2)) + 0.125*b₃*c₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A^2*(1 - ν₀^2))
        K[9,1]=0.5*b₁*(0.25*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*c₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[9,2]=0.5*c₁*(0.25*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*b₁*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[9,3]=0.25*b₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)*(0.25*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*c₁*(0.25*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + (1/2)*(0.25*b₁*(-0.5*c₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*c₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*c₁*(-0.5*b₂*(a₃ + x*b₃ + y*c₃)/A + 0.5*b₃*(a₂ + x*b₂ + y*c₂)/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(0.5*b₂*c₃/A - 0.5*b₃*c₂/A)/A + 0.25*(a₁ + x*b₁ + y*c₁)*(-0.5*b₂*c₃/A + 0.5*b₃*c₂/A)/A)*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*E₀*(1 - ν₀)/(1 - ν₀^2)
        K[9,4]=0.5*b₂*(0.25*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*c₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[9,5]=0.5*c₂*(0.25*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*b₂*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[9,6]=0.25*b₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)*(0.25*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*c₂*(0.25*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + (1/2)*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*E₀*(0.25*b₂*(0.5*c₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*c₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₂*(0.5*b₁*(a₃ + x*b₃ + y*c₃)/A - 0.5*b₃*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(-0.5*b₁*c₃/A + 0.5*b₃*c₁/A)/A + 0.25*(a₂ + x*b₂ + y*c₂)*(0.5*b₁*c₃/A - 0.5*b₃*c₁/A)/A)*(1 - ν₀)/(1 - ν₀^2)
        K[9,7]=0.5*b₃*(0.25*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*c₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[9,8]=0.5*c₃*(0.25*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)*b₃*E₀*(1 - ν₀)/(A*(1 - ν₀^2))
        K[9,9]=0.25*b₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)*(0.25*b₃*E₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*c₃*E₀*ν₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))/A + 0.25*c₃*(0.25*c₃*E₀*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)) + 0.25*b₃*E₀*ν₀*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/(A*(1 - ν₀^2)))*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + (1/2)*(0.25*b₃*(-0.5*c₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*c₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*c₃*(-0.5*b₁*(a₂ + x*b₂ + y*c₂)/A + 0.5*b₂*(a₁ + x*b₁ + y*c₁)/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(-0.5*b₁*c₂/A + 0.5*b₂*c₁/A)/A + 0.25*(a₃ + x*b₃ + y*c₃)*(0.5*b₁*c₂/A - 0.5*b₂*c₁/A)/A)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
        return K
    end
    K=hcubature(BtDB,[-1,-1],[1,1])[1]
    #9x9 to 18x18
    I=1:9
    J=[1,2,6,7,8,12,13,14,18]
    L=sparse(I,J,1.,9,18)
    return L'*K*L
end


text:
#WARNING- look, think about market caps and scaling v log
#treat thsi table as a one-off since it has no common elements w/ other tables
function maketable10(cumex::AbstractDataFrame)

  T10_ROWKEYS::Vector{Symbol} = [
    :all1232,
    :lowyield,
    :medyield,
    :highyield,
    :taxdiv,
    :nontaxdiv,
    :monthlydiv,
    :quarterlydiv,
    :semiandannualdiv,
    :specialdiv,
    :nonspecialdiv,
    :lastdiv40,
    :lastdiv4080,
    :lastdiv80,
    :compustat,
    :noncompustat,
    :bkliab13,
    :bkliab1323,
    :bkliab23,
    :assets100,
    :assets1003000,
    :assets3000]

  T10_ROWLABEL::Dict = Dict(
    :all1232 => "All (Code1232)",
    :lowyield => "\\midrule Low Yield \$(\\delta < 0.75\\%)\$",
    :medyield => "Mid Yield \$(0.75\\% < \\delta < 1.5\\%)\$",
    :highyield => "High Yield \$(1.5\\%) < \\delta\$",
    :taxdiv => "\\midrule Taxable Dividend",
    :nontaxdiv => "Non-Taxable",
    :monthlydiv => "\\midrule Monthly Dividend",
    :quarterlydiv => "Quarterly",
    :semiandannualdiv => "Semi- and Annual",
    :specialdiv => "\\midrule Special Dividend",
    :nonspecialdiv => "Not Special",
    :lastdiv40 => "\\midrule Last Dividend <40days",
    :lastdiv4080 => "40-80days",
    :lastdiv80 => ">80days",
    :compustat => "\\midrule w/ Compustat (t-1)",
    :noncompustat => "w/o Compustat",
    :bkliab13 => "\\midrule bkliab < 1/3",
    :bkliab1323 => "bkliab 1/3 - 2/3",
    :bkliab23 => "bkliab > 2/3",
    :assets100 => "\\midrule Small (Assets < \\\$100m)",
    :assets1003000 => "Medium (Assets \\\$100m-\\\$3,000m)",
    :assets3000 => "Large (Assets < \\\$3000m)"
    )

  T10_COLKEYS::Vector{Symbol} = [
    :in,
    :volcoef,
    :voltstat,
    :retcoef,
    :rettstat,
    :N]

  T10_COLLABEL::Dict = Dict(
    :in=>"in \\%",
    :volcoef=>"Coef",
    :voltstat=>"T-stat",
    :retcoef=>"Coef",
    :rettstat=>"T-stat",
    :N=>"N")

  T10_COLNAMES::Vector{Vector{String}} =[
    ["", "Cum-Ex Change from Dividend", ""],
    ["Avg \$\\delta\$", "Volatility", "Avg Return", ""],
    (s::Symbol->T10_COLLABEL[s]).(T10_COLKEYS)]

  T10_WIDTHCOLNAMES::Vector{Vector{Int}} = [
    [1,4,1],
    [1,2,2,1],
    ones(Int,length(T10_COLKEYS))]

  T10_ALIGNMENTCOLNAMES =   [
    ["c", "c", "c"],
    ["c", "c", "c", "c"],
    ["r" for i in 1:length(T10_COLKEYS)]]

  T10_DECIMALS = 1
  T10_COEF_DECIMALS = 2

  T10_ABS_YSYM = :Dabs
  T10_RET_YSYM = :Dret
  T10_NAME::String = "t10-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"

  table10(cumex,
    Yabssym = T10_ABS_YSYM,
    Yretsym = T10_RET_YSYM,
    tablename = T10_NAME,
    rowkeys = T10_ROWKEYS,
    rowlabel = T10_ROWLABEL,
    colkeys = T10_COLKEYS,
    colnames = T10_COLNAMES,
    widthcolnames = T10_WIDTHCOLNAMES,
    decimals = T10_DECIMALS,
    coefdecimals = T10_COEF_DECIMALS,
    alignmentcolnames=T10_ALIGNMENTCOLNAMES
    )

end

function table10(cumex::AbstractDataFrame;
  Yabssym::Symbol = error("Yvolsym is required"),
  Yretsym::Symbol = error("Yretsym is required"),
  tablename::String = error("tablename is required"),
  rowkeys::Vector{Symbol} = error("rowkeys is required"),
  rowlabel::Dict = error("rowlabel is required"),
  colkeys::Vector{Symbol} = error("colkeys is required"),
  colnames::Vector{Vector{String}} = error("colnames is required"),
  widthcolnames::Vector{Vector{Int}} = error("widthcolnames is required"),
  tablepath::String = TABLE_PATH,
  decimals::Int = error("decimals is required"),
  coefdecimals::Int = error("coefdecimals is required"),
  alignmentcolnames::Vector{Vector{String}} = error("alignmentcolnames is required")
  )

  #some of the regression constants and parameters should overlap
  local Xexpr = :yield
  local Xnames= (string).([:intercept, :yield])

  #housekeeping and pre-allocation
  colidx::Dict = Dict(colkey=>i for (i,colkey) ∈ enumerate(colkeys))
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, sum(widthcolnames[1]))).(1:length(rowkeys))
  rowidx::Dict = Dict(rowkey=>desccontent[i] for (i,rowkey) ∈ enumerate(rowkeys))
  descrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys)

  #apply the conditions that will apply to all rows
  #NOTE: Already filtered in formation process
  #scumex = view(cumex, dist.match .* dist.oneex .* dist.dclrexok, :)
  scumex = view(cumex, completecases(cumex[!, [:yield, Yretsym, Yabssym]]), :)

  #constructs a row for the table from the regression object and the name of the coefficient
  @inline function regrow!(rdf::AbstractDataFrame,  descrow::Vector{String})
    #start by running the regressions
    retreg = FMLM(rdf,  Xexpr, Yretsym, Xnames=Xnames, Yname=string(Yretsym), containsmissings=false)
    retΣ::Matrix{Float64} = homoskedasticΣ!(retreg)

    absreg = FMLM(rdf,  Xexpr, Yabssym, Xnames=Xnames, Yname=string(Yabssym), containsmissings=false)
    absΣ::Matrix{Float64} = homoskedasticΣ!(absreg)

    #make sure the names align
    @assert retreg.Xnames[2] == "yield"
    @assert absreg.Xnames[2] == "yield"

    #get the values for the table row
    in::Float64 = mean(rdf.yield)
    volcoef::Float64 = absreg.β[2]
    voltstat::Float64 = absreg.β[2] / absΣ[2,2]^0.5
    retcoef::Float64  = retreg.β[2]
    rettstat::Float64 = retreg.β[2] / retΣ[2,2]^0.5
    N::Float64 = size(rdf,1)

    #format and write the output
    descrow[colidx[:in]] = num2str(in, decimals, scalefactor=100.)
    descrow[colidx[:volcoef]] = brnum(volcoef, coefdecimals, scalefactor=100.)
    descrow[colidx[:voltstat]] = brnum(voltstat, decimals, scalefactor=1.)
    descrow[colidx[:retcoef]] = brnum(retcoef, coefdecimals, scalefactor=100.)
    descrow[colidx[:rettstat]] = brnum(rettstat, decimals, scalefactor=1.)
    descrow[colidx[:N]] = num2str(N, 0, Ints=true)
  end

  #first row is all rows (same result as before)
  scumex1232::SubDataFrame = view(scumex, scumex.ordinary, :)
  regrow!(scumex1232, rowidx[:all1232])

  #performs the regressions within each group and builds the corresponding table rows
  @inline function reggroupby!(df::AbstractDataFrame, rowidx::Dict, Fgroup::Symbol)
    local sdf::SubDataFrame = view(df, (!ismissing).(df[!,Fgroup]), :)

    (size(sdf,1) > 10) || error(
      "sdf below minimum size for effective groupings of column $Fgroup")

    #for each group classification
    for ssdf::SubDataFrame ∈ groupby(sdf, Fgroup)
      local grp::Symbol = ssdf[1,Fgroup]
      haskey(rowidx, grp) || error(
        "Group classification $grp from column $Fgroup not found as key in rowidx")

      regrow!(ssdf, rowidx[grp])
    end
  end

  #build the table rows
  reggroupby!(scumex1232, rowidx, :Gyield)
  reggroupby!(scumex, rowidx, :Gtaxdiv)
  reggroupby!(scumex, rowidx, :Gfreqdiv)
  reggroupby!(scumex, rowidx, :Gspecialdiv)
  reggroupby!(scumex, rowidx, :Glastdiv)
  reggroupby!(scumex1232, rowidx, :Gcompustat)
  reggroupby!(scumex1232, rowidx, :Gbkliab)
  reggroupby!(scumex1232, rowidx, :Gassets)

  #println(desccontent)

  tt::String = textable(colnames=colnames,
    descrownames=descrownames,
    desccontent=desccontent,
    alignmentcolnames=alignmentcolnames,
    widthcolnames=widthcolnames)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Base: @deprecate

@deprecate anynull(x) any(isnull, x)
@deprecate allnull(x) all(isnull, x)


text:
data_union = Union{Nothing, AbstractString, Dict, Array{T, 1} where T}
init_union = Union{Nothing, StanBase.Init, AbstractString, Dict, Array{T, 1} where T}

"""
$(SIGNATURES)

Default `output_base` data files, in the same directory as the model. Internal, not exported.
"""
data_file_path(output_base::AbstractString, id::Int) = output_base * "_data_$(id).R"

"""
$(SIGNATURES)

Default `output_base` init files, in the same directory as the model. Internal, not exported.
"""
init_file_path(output_base::AbstractString, id::Int) = output_base * "_init_$(id).R"

"""
$(SIGNATURES)

Sample `model.n_chains` from `model` using `data` and `init`.

Return the full paths of the sample files and logs as pairs.

In case of an error with a chain, the first value is `nothing`.

`output_base` is used to write the data file (using `StanDump.stan_dump`) and to determine
the resulting names for the sampler output. It defaults to the source file name without the
extension.

When `data` or `init` are provided as a `NamedTuple`or a `Dict`, it is written using
`StanDump.stan_dump` first. If an AbstractString is specified, it is used as a path
to an existing file which will be copied to the output directory unless the length == 0.

When `rm_samples` (default: `true`), remove potential pre-existing sample files after
compiling the model.
"""
function stan_sample(model::T; kwargs...) where {T <: CmdStanModels}
  n_chains = 4
  rm_samples = true
  diagnostics = false

  if :n_chains in keys(kwargs)
    n_chains = kwargs[:n_chains]
    set_n_chains(model, n_chains)
  end
  :init in keys(kwargs) && update_R_files(model, kwargs[:init], n_chains, "init")
  :data in keys(kwargs) && update_R_files(model, kwargs[:data], n_chains, "data")
  if :diagnostics in keys(kwargs)
    diagnostics = kwargs[:diagnostics]
    setup_diagnostics(model, get_n_chains(model))
  end

  _stan_sample(model;  rm_samples = rm_samples)

end

function _stan_sample(model::T; rm_samples = true) where {T <: CmdStanModels}

    rm_samples && rm.(StanRun.find_samples(model.sm))
    cmds_and_paths = [stan_cmd_and_paths(model, id)
                      for id in 1:get_n_chains(model)]
    pmap(cmds_and_paths) do cmd_and_path
        cmd, (sample_path, log_path) = cmd_and_path
        success(cmd) ? sample_path : nothing, log_path
    end

end

"""
$(SIGNATURES)

Run a Stan command. Internal, not exported.
"""
function stan_cmd_and_paths(model::T, id::Integer) where {T <: CmdStanModels}

    append!(model.sample_file, [StanRun.sample_file_path(model.output_base, id)])
    model.output.file = model.sample_file[id]
    if length(model.diagnostic_file) > 0
      model.output.diagnostic_file = model.diagnostic_file[id]
    end
    append!(model.log_file, [StanRun.log_file_path(model.output_base, id)])
    append!(model.cmds, [cmdline(model, id)])
    pipeline(model.cmds[id]; stdout=model.log_file[id]), (model.sample_file[id], model.log_file[id])

end

function update_R_files(model, input, n_chains, fname_part="data")

  model_field = fname_part == "data" ? model.data_file : model.init_file
  if typeof(input) <: NamedTuple || typeof(input) <: Dict
    for i in 1:n_chains
      stan_dump(model.output_base*"_$(fname_part)_$i.R", input, force=true)
      append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
    end
  elseif  typeof(input) <: Array
    if length(input) == n_chains
      for (i, d) in enumerate(input)
        stan_dump(model.output_base*"_$(fname_part)_$i.R", d, force=true)
        append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
      end
    else
      @info "Data vector length does not match number of chains,"
      @info "only first element in data vector will be used,"
      for i in 1:n_chains
        stan_dump(model.output_base*"_$(fname_part)_$i.R", input[1], force=true)
        append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
      end
    end
  elseif typeof(input) <: AbstractString && length(input) > 0
    for i in 1:n_chains
      cp(input, "$(model.output_base)_$(fname_part)_$i.R", force=true)
      append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
    end
  else
    error("\nUnrecognized input argument: $(typeof(input))\n")
  end

end

function setup_diagnostics(model, n_chains)

  for i in 1:n_chains
    append!(model.diagnostic_file, [model.output_base*"_diagnostic_$i.log"])
  end

end


text:
struct DGMesh{dim,T}
    ncells::Int
    domain::Vector{IntervalBox{dim,T}}
    connectivity::Vector{Vector{Tuple{Int,Int}}}
    isinteriorcell::Vector{Bool}
    cellsign::Vector{Int}
    facemaps::Vector{LineMap{dim,T}}
    facescale::Vector{T}
    cellmap::CellMap{dim,T}
    function DGMesh(domain::Vector{IntervalBox{dim,T}},
        connectivity::Vector{Vector{Tuple{Int,Int}}},
        isinteriorcell::Vector{Bool},
        cellsign::Vector{Int},facemaps::Vector{LineMap{dim,T}},
        facescale::Vector{T},cellmap::CellMap{dim,T}) where {dim,T}

        @assert dim == 2
        nfaces = number_of_faces(dim)

        ncells = length(domain)
        @assert length(connectivity) == ncells
        @assert length(cellsign) == ncells
        @assert length(facemaps) == nfaces
        @assert length(isinteriorcell) == ncells

        new{dim,T}(ncells,domain,connectivity,isinteriorcell,
            cellsign,facemaps,facescale,cellmap)

    end
end

function DGMesh(mesh::UniformMesh,coeffs,poly)
    dim = dimension(mesh)
    domain = cell_domain(mesh)
    @assert length(domain) > 0
    cellmap = CellMap(domain[1])
    facescale = face_determinant_jacobian(cellmap)
    connectivity = cell_connectivity(mesh)
    isinteriorcell = interior_cells(connectivity)
    cellsign = cell_signatures(coeffs,poly)
    facemaps = reference_cell_facemaps(dim)
    return DGMesh(domain,connectivity,isinteriorcell,
        cellsign,facemaps,facescale,cellmap)
end

function cell_domain(mesh)
    dim = dimension(mesh)
    T = default_float_type()
    domain = Vector{IntervalBox{dim,T}}(undef,mesh.total_number_of_elements)
    for idx in 1:mesh.total_number_of_elements
        xL,xR = CartesianMesh.element(mesh,idx)
        domain[idx] = IntervalBox(xL,xR)
    end
    return domain
end

function cell_connectivity(mesh)
    nfaces = faces_per_cell(mesh)
    ncells = number_of_elements(mesh)
    connectivity = [Vector{Tuple{Int,Int}}(undef,nfaces) for i = 1:ncells]

    for cellid in 1:ncells
        nbrcellids = neighbors(mesh,cellid)
        nbrfaceids = [nc == 0 ? 0 : neighbor_faceid(faceid) for (faceid,nc) in enumerate(nbrcellids)]
        nbrcellandface = collect(zip(nbrcellids,nbrfaceids))
        connectivity[cellid][:] .= nbrcellandface
    end
    return connectivity
end

function is_interior_cell(cellconnectivity)
    return all([c[1] != 0 for c in cellconnectivity])
end

function interior_cells(connectivity)
    ncells = length(connectivity)
    isinterior = zeros(Bool,ncells)
    for cellid in 1:ncells
        isinterior[cellid] = is_interior_cell(connectivity[cellid])
    end
    return isinterior
end

function cell_signatures!(cellsign,coeffs,poly)

    nf,ncells = size(coeffs)

    @assert length(cellsign) == ncells

    dim = dimension(poly)
    box = reference_cell(dim)
    fill!(cellsign,0)

    for idx in 1:ncells
        update!(poly,coeffs[:,idx])
        cellsign[idx] = sign(poly,box)
    end
end

function cell_signatures(coeffs,poly)
    nf,ncells = size(coeffs)
    cellsign = zeros(Int,ncells)
    cell_signatures!(cellsign,coeffs,poly)
    return cellsign
end

function dimension(dgmesh::DGMesh{dim}) where dim
    return dim
end

function number_of_cells(dgmesh::DGMesh)
    return dgmesh.ncells
end


text:
using SparseArrays
using Serialization
include("constants.jl")
include("data.jl")
include("agent_profile.jl")
include("additional_activity.jl")


text:
# Autogenerated wrapper script for OpenSSL_jll for aarch64-linux-gnu
export libcrypto, libssl, openssl

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libcrypto`
const libcrypto_splitpath = ["lib", "libcrypto.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libcrypto_path = ""

# libcrypto-specific global declaration
# This will be filled out by __init__()
libcrypto_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libcrypto = "libcrypto.so.1.1"

# Relative path to `libssl`
const libssl_splitpath = ["lib", "libssl.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libssl_path = ""

# libssl-specific global declaration
# This will be filled out by __init__()
libssl_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libssl = "libssl.so.1.1"

# Relative path to `openssl`
const openssl_splitpath = ["bin", "openssl"]

# This will be filled out by __init__() for all products, as it must be done at runtime
openssl_path = ""

# openssl-specific global declaration
function openssl(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(openssl_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"OpenSSL")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global libcrypto_path = normpath(joinpath(artifact_dir, libcrypto_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libcrypto_handle = dlopen(libcrypto_path)
    push!(LIBPATH_list, dirname(libcrypto_path))

    global libssl_path = normpath(joinpath(artifact_dir, libssl_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libssl_handle = dlopen(libssl_path)
    push!(LIBPATH_list, dirname(libssl_path))

    global openssl_path = normpath(joinpath(artifact_dir, openssl_splitpath...))

    push!(PATH_list, dirname(openssl_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
#!/usr/bin/env julia

using Luxor

function dot(pos)
    gsave()
    sethue("red")
    circle(pos, 5, :fill)
    grestore()
end

function showt(c, p, ha, va, n)
  text(c, p, halign=ha, valign=va)
  gsave()
  setopacity(0.1)
  text(string(n), p)
  grestore()
end

function text_alignment_tests(fname)
    legend = String[]
    Drawing(1400,1400, fname)
    origin()
    setopacity(0.8)
    sethue("black")
    fontsize(65)
    tiles = Tiler(1000, 1000, 4, 3, margin=50)
    haligns = (:left, :center, :right)
    valigns = (:baseline, :top, :middle, :bottom)
    current_h = 1
    current_v = 1
    for (pos, n) in tiles
      gsave()
      h = haligns[current_h]
      v = valigns[current_v]
      dot(Point(pos))
      showt("Å˰̀Ά", pos, h, v, n)
      grestore()
      push!(legend, "$n h: $h v: $v")
      current_h += 1
      if current_h > 3
        current_v += 1
        current_h = 1
      end
      if current_v > 4
        break
      end
    end
    fontsize(8)
    text(join(legend, "; "), 0, 600, halign=:center)
    finish()
    println("finished test: output in $(fname)")
end

text_alignment_tests("/tmp/text-alignment-tests.pdf")


text:
import CLBLAS
import OpenCL

    const clblas = CLBLAS
    clblas.setup()
    const cl = OpenCL
    alpha = cl.cl_float(5)
    beta = cl.cl_float(2)
    A = float32([1.0 2.0 3.0; 4.0 5.0 6.0])  #rand(Float32, 2, 3)
    B = float32([1.0 2.0; 3.0 4.0; 5.0 6.0]) #rand(Float32, 3, 2)
    C = float32([1.0 2.0; 3.0 4.0])#rand(Float32, 2, 2)

    data = (5.0*A*B) + (2.0*C)

    future = clblas.clblasSgemm(UInt32(0), UInt32(0), alpha, A, B, beta, C)

    result = fetch(future)
    println(result)
    println("------------")
    println(data)
    @assert(result == data)
    info("success!")

    clblas.teardown()


text:
macro proto( expr )
    if expr.head != Symbol("struct")
        throw(ArgumentError("Expected expression to be a type definition."))
    end

    name = expr.args[2]
    type_parameters = nothing
    type_parameter_names = []
    type_parameter_types = []
    if !(name isa Symbol)
        type_parameters = name.args[2:end]
        name = name.args[1]
        type_parameter_names = map(type_parameters) do t
                                    if t isa Symbol
                                        t
                                    else
                                        t.args[1]
                                    end
                                end
        type_parameter_types = map(type_parameters) do t
                                    if t isa Symbol
                                        :Any
                                    else
                                        t.args[2]
                                    end
                                end
        type_parameter_types = Dict( type_parameter_names[i] => type_parameter_types[i] for i in eachindex(type_parameters))
    end

    fields = map(expr.args[3].args[2:2:length(expr.args[3].args)]) do field
                    if field isa Symbol
                        return field
                    end
                    if field.head == :(=)
                        field.args[1]
                    else
                        field
                    end
                end
    field_info = map(fields) do field
                        return if field isa Symbol
                            (field, Any)
#                        elseif field.head == :(=) && !(field.args[1] isa Symbol)
#                            (field.args[1].args[1], field.args[1].args[2])
                        else
                            (field.args[1], field.args[2])
                        end
                    end
    field_names = Tuple(getindex.(field_info, 1))
    field_types = quote Tuple{$(getindex.(field_info, 2)...)} end
    field_subtype_info = map(getindex.(field_info, 2)) do ft
        if ft in type_parameter_names
            return type_parameter_types[ft]
        else
            return ft
        end
    end

    params_ex = Expr(:parameters)
    call_args = Any[]

    Base._kwdef!(expr.args[3], params_ex.args, call_args)

    # remove escapes
    params_ex.args = map(params_ex.args) do ex
        if ex isa Symbol return ex end
        ex.args[2] = ex.args[2].args[1]
        ex
    end

    ex = quote
            if !@isdefined $name
                struct $name{NT<:NamedTuple}
                    properties::NT
                end # struct
            else
                the_methods = collect(methods($name))
                Base.delete_method(the_methods[1])
                Base.delete_method(the_methods[3])
            end # if

            $(
                if type_parameters === nothing
                    :( $name(args...) = $name(NamedTuple{$field_names, $field_types}(args)) )
                else
                    :( $name($(fields...)) where {$(type_parameters...)} = $name(NamedTuple{$field_names, $field_types}(($(field_names...),))) )
                end
            )

            function $name($params_ex)
                $name($(call_args...))
            end

            function Base.getproperty( o::$name, s::Symbol )
                return getproperty( getfield(o, :properties), s )
            end # function

            function Base.propertynames( o::$name )
                return propertynames( getfield(o, :properties) )
            end # function
        end # quote
    ex |> esc
end # macro


text:
using TermWin

TermWin.initsession()
#s = "a test"
v = newTwEntry( rootTwScreen, Rational{Int}; width=25, posy=:random, posx=:random, showHelp=true, box=true )
v.title = "Input: "
v.data.tickSize = 1
activateTwObj( rootTwScreen )
ret = v.value
#unregisterTwObj( rootTwScreen, v )
TermWin.endsession()
println( "You entered ", string( ret ) )


text:
using ColorPlus
using Test

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


text:
using MOCNeutronTransport

# I can't remember if the bottom left corner of the CMFD grid needs to be (0, 0, 0)
# To compensate, we do some extra math to ensure that happens.
grid_offset = 0.1 # Buffer between boundary of the model and the edge of the CMFD grid (cm)
# Bounding box for the CMFD grid
# xmin, xmax, ymin, ymax
bb = (0.0, 80.62054 + 2*grid_offset,
      0.0, 46.2645 + 2*grid_offset)

# Model setup
# ----------------------------------------------------------------------------------------
gmsh.initialize()
gmsh.option.set_number("Geometry.OCCParallel", 1)
gmsh.option.set_number("General.NumThreads", 32)
# Import CAD file
gmsh.merge("fnr_core_slice.step")
# Visualize the model
# gmsh.fltk.run()

# Shift the model to the origin of the xy-plane
# Get the dim tags of all dimension 2 entites
dim_tags = gmsh.model.get_entities(2)
# Translate it so that the corner is at the origin
gmsh.model.occ.translate(dim_tags, 407.51125 + 10*grid_offset,
                                   307.848 + 10*grid_offset,
                                   -300)
# Model is in mm. Convert to cm by shrinking everything by 1/10
gmsh.model.occ.dilate(dim_tags, 0, 0, 0, 1//10, 1//10, 0)
# Synchronize the CAD and Gmsh models
gmsh.model.occ.synchronize()
# Verify results
gmsh.fltk.run()

# Assign materials
# ------------------------------------------------------------------------------------------
all_tags = getindex.(gmsh.model.get_entities(2), 2)
# Assign the materials as physical groups
p = gmsh.model.add_physical_group(2, all_tags)
gmsh.model.set_physical_name(2, p, "MATERIAL_UO2")

# Overlay CMFD/hierarchical grid
# -------------------------------------------------------------------------------------------
nx = [10, 1]
ny = [6, 1]
grid_tags = gmsh_overlay_rectangular_grid(bb, "MATERIAL_MODERATOR", nx, ny)
# Verify results
gmsh.fltk.run()

# Mesh
# ------------------------------------------------------------------------------------------------
# Set the characteristic edge length of the mesh cells
lc = 0.3 # cm
gmsh.model.mesh.set_size(gmsh.model.get_entities(0), lc)
# Optional mesh optimization:
niter = 2 # The optimization iterations

# 2nd order triangles
gmsh.model.mesh.generate(2) # Triangles first for high order meshes.
gmsh.model.mesh.set_order(2)
for () in 1:niter
    gmsh.model.mesh.optimize("HighOrderElastic")
    gmsh.model.mesh.optimize("Relocate2D")
    gmsh.model.mesh.optimize("HighOrderElastic")
end

# gmsh.fltk.run()

# Mesh conversion
# ----------------------------------------------------------------------------------------------------
# We want to write an Abaqus file from Gmsh, read the Abaqus file into Julia, and convert that into
# a rectangularly partitioned hierarchical XDMF file, the final input to MPACT
#
# Write Abaqus file
gmsh.write("fnr_slice.inp")
gmsh.finalize() # done with Gmsh. Finalize
mesh = read_abaqus_2d("fnr_slice.inp")
# Convert mesh into Hierarchical Rectangularly Partitioned Mesh
HRPM = partition_rectangularly(mesh)
write_xdmf_2d("fnr_slice.xdmf", HRPM)
#
## Mass conservation
## ----------------------------------------------------------------------------------------------------
## See how well your mesh conserves mass. Analytic solution not available, so fine tri6 mesh used as
## reference
#fuel_area = area(mesh, "MATERIAL_UO2")
#fuel_area_ref = 6.119632103711788
#println("Fuel area (reference): $fuel_area_ref")
#println("Fuel area ( computed): $fuel_area")
#err = 100*(fuel_area - fuel_area_ref)/fuel_area_ref
#println("Error (relative): $err %")
#println("")
#
#clad_area = area(mesh, "MATERIAL_CLAD")
#clad_area_ref = 19.884367501087866
#println("Clad area (reference): $clad_area_ref")
#println("Clad area ( computed): $clad_area")
#err = 100*(clad_area - clad_area_ref)/clad_area_ref
#println("Error (relative): $err %")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
function load(f::File{format"OMETIFF"}; dropunused=true, inmemory=true)
    open(f) do s
        ret = load(s; dropunused=dropunused, inmemory=inmemory)
    end
end

"""
    load(io; dropunused, inmemory) -> ImageMetadata.ImageMeta

Load an OMETIFF file using the stream `io`.

**Arguments**
- `dropunused::Bool`: controls whether dimensions of length 1 are dropped
  automatically (default) or not.
- `inmemory::Bool`: controls whether arrays are fully loaded into memory
  (default) or left on disk and specific parts only loaded when accessed.

!!! tip
    The `inmemory=false` flag currently returns a read-only view of the data on
    the disk for data integrity reasons. In order to modify the contents, you
    must copy the data into an in-memory container--at least until
    [#52](https://github.com/tlnagy/OMETIFF.jl/issues/52) is fixed--like so:

    ```
    copy(arr)
    ```
"""
function load(io::Stream{format"OMETIFF"}; dropunused=true, inmemory=true)
    if io.filename != nothing && !occursin(".ome.tif", io.filename)
        throw(FileIO.LoaderError("OMETIFF", "Not an OME TIFF file!"))
    end

    orig_file = TiffFile(io)
    summary = load_comments(orig_file)

    # load master OME-XML that contains all information about this dataset
    omexml = load_master_xml(orig_file)

    # find all images in this dataset, can either have the Image or Pixel tag
    containers = findall("//*[@DimensionOrder]", omexml, ["ns" => namespace(omexml)])

    pos_names = nodecontent.(findall("/ns:OME/ns:Image/ns:StageLabel[@Name]/@Name", omexml, ["ns"=>namespace(omexml)]))
    # if all position names aren't unique then substitute names
    if length(pos_names) == 0 || !allunique(pos_names)
        pos_names = ["Pos$i" for i in 1:length(containers)]
    end

    master_rawtype = nothing
    mappedtype = Int64
    dimlist = []
    axeslist = []

    for (idx, container) in enumerate(containers)
        dims, axes_info = build_axes(container)
        push!(dimlist, dims)
        push!(axeslist, axes_info)
        rawtype, mappedtype = type_mapping[container["Type"]]

        if master_rawtype === nothing
            master_rawtype = rawtype
        elseif master_rawtype != rawtype
            throw(FileIO.LoaderError("OMETIFF", "Multiple different storage types are not yet support in a multi position image"))
        end
    end

    ifd_indices = OrderedDict{Int, NTuple{4, Int}}()
    ifd_files = OrderedDict{Int, Tuple{String, String}}()
    obs_filepaths = Dict{String, Int}()
    for (idx, container) in enumerate(containers)
        OMETIFF.ifdindex!(ifd_indices, ifd_files, obs_filepaths, container, dimlist[idx], orig_file, idx)
    end

    files, ifds = get_ifds(orig_file, ifd_indices, ifd_files)

    # determine size of all the dims from the ifds in the tiff file instead
    # of the sizes embedded in the Pixel node
    true_dims = [Set{Int}() for i in 1:4]
    for ifd_index in keys(ifds), dim in 1:4
        push!(true_dims[dim], ifd_index[dim])
    end

    # generate new master dim list with the true dims from above
    master_dims = dimlist[1]
    dimnames = keys(master_dims)[3:6]
    new_data = [dimnames[i]=>length(true_dims[i]) for i in 1:4]
    master_dims = merge(master_dims, new_data)

    masteraxis = copy(axeslist[1])
    masteraxis[6] = Axis{:position}(OMETIFF.to_symbol.(pos_names))

    # check if the axes computed earlier are the correct length, if not,
    # we're not sure about the information through the whole movie so lets
    # strip the units
    err_axes = findall(length.(masteraxis) .!== values(master_dims))
    for ax in err_axes
        masteraxis[ax] = masteraxis[ax](1:master_dims[ax])
    end

    elapsed_times = get_elapsed_times(containers, master_dims, masteraxis)

    if inmemory
        img = inmemoryarray(ifds, master_dims, master_rawtype, mappedtype)
    else
        img = ReadonlyTiffDiskArray(Gray{mappedtype}, master_rawtype, ifds, values(master_dims));
    end

    # find dimensions of length 1 and remove them
    if dropunused
        unused_dims = findall(values(master_dims) .== 1)
        img = dropdims(img; dims=tuple(unused_dims...))
        deleteat!(masteraxis, unused_dims)
    end

    ImageMeta(AxisArray(img, masteraxis...),
              Description=summary,
              Elapsed_Times=elapsed_times)
end

"""
    inmemoryarray(ifds, dims, rawtype, mappedtype) -> Array

Builds an in-memory high-dimensional image using the mapping provided by `ifds`
from indices to [`OMETIFF.IFD`](@ref) objects. The IFD objects store handles to
the file objects and the offsets for the data. `dims` stores the size of each
named dimension. The `rawtype` parameter describes the storage layout of each
element on disk and `mappedtype` is the corresponding fixed or floating point
type.
"""
function inmemoryarray(ifds::OrderedDict{NTuple{4, Int}, IFD},
                       dims::NamedTuple,
                       rawtype::Type,
                       mappedtype::Type)

    data = Array{rawtype, length(dims)}(undef, dims...)

    width, height = dims[1], dims[2]
    # assume no strips
    tmp = Array{rawtype}(undef, height, width)

    # iterate over each IFD
    for (indices, ifd) in ifds
        n_strips = length(ifd.strip_offsets)
        strip_len = floor(Int, (width * height) / n_strips)

        # if the data is stripped and we haven't fix tmp's layout then lets make
        # tmp equal to one strip. This'll be fixed in Julia 1.4
        if n_strips > 1 && size(tmp) != (strip_len, )
            tmp = Array{rawtype}(undef, strip_len)
        end

        target = view(data, :, :, indices...)
        _read_ifd_data!(ifd, target, tmp)

        # transposition must happen here since the on-disk variant does this on access
        if ndims(tmp) == 2
            target .= tmp'
        end
    end

    reinterpret(Gray{mappedtype}, data)
end


text:
abstract AbstractBoundaryCondition

type Value{V} <: AbstractBoundaryCondition end
Value(V) = Value{Val{V}}
function apply_value{V}(BC::Type{Value{V}},I::Symbol, fi::Int, x::Symbol, xi::Int)
    quote
        xmin, xmax = extrema(getfield($I,$fi))
        if $x[$xi] < xmin || $x[$xi] > xmax
            return $V.parameters[1]
        end
    end
end

type Nearest <: AbstractBoundaryCondition end
function apply_nearest(I::Symbol, fi::Int, x::Symbol, xi::Int)
    :($x[$xi] = clamp($x[$xi], extrema(getfield($I,$fi))...))
end

type Reflect <: AbstractBoundaryCondition end
function apply_reflect(I::Symbol, fi::Int, x::Symbol, xi::Int)
    quote
        xmin, xmax = extrema(getfield($I,$fi))
        if $x[$xi] < xmin
            delX = xmax - xmin
            n = floor((xmin - $x[$xi])/delX)
            $x[$xi] = xmax - ((xmin - $x[$xi]) - n*delX)
        elseif $x[$xi] > xmax
            delX = xmax - xmin
            n = floor(($x[$xi] - xmax)/delX)
            $x[$xi] = xmin + ($x[$xi] - xmax - n*delX)
        end
    end
end

type None <: AbstractBoundaryCondition end
function apply_none()
    :()
end

type Error <: AbstractBoundaryCondition end
function apply_error(I::Symbol,fi::Int,x::Symbol,xi::Int)
    quote
        xmin, xmax = extrema(getfield($I,$fi))
        if $x[$xi] < xmin || $x[$xi] > xmax
            return throw(ArgumentError("Out of range"))
        end
    end
end


text:
"""
    SubAction <: AbstractAction

A SubAction can be used in the keyword arguments of an [`Action`](@ref) to define small
sub actions on the action function, such as [`appear`](@ref).

A SubAction should not be created by hand but instead by using one of the constructors.

# Fields
- `frames::Frames`: the frames relative to the parent [`Action`](@ref)
- `anim::Animation`: defines the interpolation function for the transitions
- `func::Function`: the function that gets called in each of those frames.
    Takes the following arguments: `video, action, subaction, rel_frame`
- `transitions::Vector{Transition}`: A list of transitions like [`Translation`](@ref)
- `internal_transitions::Vector{InternalTransition}`:
    A list of internal transitions which store the current transition for a specific frame.
- `defs::Dict{Symbol, Any}` any kind of definitions that are relevant for the subaction.
"""
mutable struct SubAction <: AbstractAction
    frames::Frames
    anim::Animation
    func::Function
    transitions::Vector{Transition}
    internal_transitions::Vector{InternalTransition}
    defs::Dict{Symbol,Any}
end

SubAction(transitions::Transition...) = SubAction(:same, transitions...)
SubAction(func::Function) = SubAction(:same, func)

"""
    SubAction(frames, easing::Union{ReversedEasing, Easing}, args...)

A `SubAction` can be defined with frames an
[easing function](https://jkrumbiegel.github.io/Animations.jl/stable/#Easings-1) and either
a function or transformation(s).

# Example
In the following example a filled circle with radius 50 appears in the first 20 frames,
which means the opacity is increased from 0 to 1.0. The interpolation function used here is
sineio from [`Animations.jl`](https://github.com/jkrumbiegel/Animations.jl).
Then it stays at full opacity and disappears the same way in the last 20 frames using a
linear decay.

```julia
javis(demo, [
    BackgroundAction(1:100, ground),
    Action((args...)->circle(O, 50, :fill); subactions = [
        SubAction(1:20, sineio(), appear(:fade)),
        SubAction(81:100, disappear(:fade))
    ])
])
```

# Arguments
- `frames`: A list of frames for which the function should be called.
    - The frame numbers are relative to the parent [`Action`](@ref).
- `easing::Union{ReversedEasing, Easing}`: The easing function for `args...`
- `args...`: Either a function like [`appear`](@ref) or a Transformation
    like [`Translation`](@ref)
"""
SubAction(frames, easing::Union{ReversedEasing,Easing}, args...) =
    SubAction(frames, easing_to_animation(easing), args...)

SubAction(frames, anim::Animation, transition::Transition...) =
    SubAction(frames, anim, (args...) -> 1, transition...)

SubAction(easing::Union{ReversedEasing,Easing}, args...) =
    SubAction(:same, easing_to_animation(easing), args...)

SubAction(anim::Animation, args...) = SubAction(:same, anim, args...)
"""
    SubAction(frames, func::Function)

A `SubAction` can be defined with frames and a function
inside the `subactions` kwarg of an [`Action`](@ref).
In the following example a filled circle with radius 50 appears in the first 20 frames,
which means the opacity is increased from 0 to 1.0.
Then it stays at full opacity and disappears the same way in the last 20 frames.

# Example
javis(demo, [
    BackgroundAction(1:100, ground),
    Action((args...)->circle(O, 50, :fill); subactions = [
        SubAction(1:20, appear(:fade)),
        SubAction(81:100, disappear(:fade))
    ])
])

# Arguments
- `frames`: A list of frames for which the function should be called.
    - The frame numbers are relative to the parent [`Action`](@ref).
- `func::Function`: The function that gets called for the frames.
    - Needs to have four arguments: `video, action, subaction, rel_frame`
    - For [`appear`](@ref) and [`disappear`](@ref) a closure exists,
      such that `appear(:fade)` works.
"""
SubAction(frames, func::Function) = SubAction(frames, easing_to_animation(linear()), func)

"""
    SubAction(frames, trans::Transition...)

A `SubAction` can also be defined this way with having a list of transitions.
This is similar to defining transitions inside [`Action`](@ref)

In the following example a circle is faded in during the first 25 frames then moves to
- `Point(100, 20)` then to `Point(120, -20)` (the translations are added)
- and then back to the origin
In the last 25 frames it disappears from the world.

# Example
```
javis(demo, [
        BackgroundAction(1:200, ground_opacity),
        Action((args...)->circle(O, 50, :fill); subactions = [
            SubAction(1:25, appear(:fade)),
            SubAction(26:75, Translation(Point(100, 20))),
            SubAction(76:100, Translation(Point(20, -40))),
            SubAction(101:175, Translation(Point(-120, 20))),
            SubAction(176:200, disappear(:fade))
        ]),
    ], tempdirectory="current/images", pathname="current/circle_square.gif")
```

# Arguments
- `frames`: A list of frames for which the function should be called.
    - The frame numbers are relative to the parent [`Action`](@ref).
- `trans::Transition...`: A list of transitions that shall be performed.
"""
SubAction(frames, trans::Transition...) =
    SubAction(frames, easing_to_animation(linear()), (args...) -> 1, trans...)

SubAction(frames, anim::Animation, func::Function, transitions::Transition...) =
    SubAction(frames, anim, func::Function, collect(transitions), [], Dict{Symbol,Any}())


text:
module PICDataStructures

export AbstractPICDataStructure,
    ScalarField, ScalarVariable, VectorField, VectorVariable,
    AbstractGrid, AbstractAxisGrid, AxisGrid, SparseAxisGrid, ParticlePositions,
    build_vector, getdomain, mapgrid, axisnames,
    scalarfield, scalarvariable, vectorfield,
    hasunit, unitname,
    downsample,
    scalarness, domain_type, domain_discretization,
    ScalarQuantity, VectorQuantity, LatticeGrid, ParticleGrid,
    plotdata

# Arrays
using StaticArrays
using StructArrays
using RecursiveArrayTools
using BangBang
using ArrayInterface: parameterless_type
using StructArrays: components, component, similar_structarray
using Base.Broadcast: Broadcasted, ArrayStyle
using Base: @propagate_inbounds
# Math
using LinearAlgebra
using CoordinateTransformations
using ImageTransformations
using IntervalSets
# Units
using Unitful
using Unitful: Units
# Plotting
# Track what remains to be moved over to MakieCore
using Makie: @lift, lift, Arrows, Contour, Figure, Colorbar, Axis,
    AxisAspect, DataAspect, replace_automatic!,
    Point, Vec, LScene, SceneSpace,
    arrows!, contour!,
    cam3d!
using Observables
using MakieCore
using MakieCore: convert_arguments, @recipe, theme, Attributes,
    automatic,
    PointBased, DiscreteSurface, VolumeLike
import RecipesBase, UnitfulRecipes

abstract type AbstractPICDataStructure{T,N,G} <: AbstractArray{T,N} end

include("abstractarray.jl")
include("grids/grids.jl")
include("scalar.jl")
include("vector.jl")
include("traits.jl")
include("units.jl")
include("algebra.jl")
# include("transformations.jl")
include("subset/subset.jl")
include("recipes/plot.jl")
include("recipes/recipesbase.jl")
include("utils.jl")
include("show.jl")

end # module PICDataStructures