text:
using Base.Iterators
using Scapin.Brick

function expected_avg_gradient_operator(h::NTuple{d,T}) where {d,T<:Number}
    ℒ = cell_vertices(d)
    D = zeros(T, d, fill(2, d)...)
    for i = 1:d, p in ℒ
        D[i, p] = (-1)^p[i] / 2^(d - 1) / h[i]
    end
    return D
end

function expected_avg_strain_displacement_operator(h::NTuple{d,T}) where {d,T<:Number}
    ℒ = cell_vertices(d)
    D = expected_avg_gradient_operator(h)
    B = zeros(T, d, d, d, fill(2, d)...)
    for i = 1:d, j = 1:d, k = 1:d, p in ℒ
        B[i, j, k, p] = (T(j == k) * D[i, p] + T(i == k) * D[j, p]) / 2
    end
    return B
end

function expected_K_λ(h::NTuple{d,T}) where {d,T<:Number}
    K_λ = zeros(Float64, d, fill(2, d)..., d, fill(2, d)...)
    if d == 2
        value = -1 / 4
        K_λ[1, 1, 1, 2, 1, 2] = value
        K_λ[1, 1, 1, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 1, 1] = value
        K_λ[1, 2, 2, 2, 2, 1] = value
        K_λ[2, 1, 1, 1, 2, 1] = value
        K_λ[2, 1, 1, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 1, 1] = value
        K_λ[2, 1, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 2, 2] = value
        K_λ[2, 2, 2, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 1, 2] = value

        value = 1 / 4
        K_λ[1, 1, 1, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 1, 1] = value
        K_λ[1, 1, 2, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 1, 2] = value
        K_λ[1, 2, 1, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 2, 2] = value
        K_λ[2, 1, 1, 1, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 1] = value
        K_λ[2, 1, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 1, 2] = value
        K_λ[2, 2, 2, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 2] = value

        value = -h[1] / (3 * h[2])
        K_λ[2, 1, 1, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 1] = value

        value = -h[1] / (6 * h[2])
        K_λ[2, 1, 1, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 1, 1] = value

        value = h[1] / (6 * h[2])
        K_λ[2, 1, 1, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 1, 2] = value

        value = h[1] / (3 * h[2])
        K_λ[2, 1, 1, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 2, 2] = value

        value = -h[2] / (3 * h[1])
        K_λ[1, 1, 1, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 1, 2] = value

        value = -h[2] / (6 * h[1])
        K_λ[1, 1, 1, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 1, 1] = value

        value = h[2] / (6 * h[1])
        K_λ[1, 1, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 2, 1] = value

        value = h[2] / (3 * h[1])
        K_λ[1, 1, 1, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 2, 2] = value
    elseif d == 3
        value = -h[1] / 12
        K_λ[2, 1, 1, 1, 3, 1, 1, 2] = value
        K_λ[2, 1, 1, 1, 3, 1, 2, 2] = value
        K_λ[2, 1, 1, 2, 3, 1, 1, 2] = value
        K_λ[2, 1, 1, 2, 3, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 3, 1, 1, 1] = value
        K_λ[2, 1, 2, 1, 3, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 3, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 3, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 3, 2, 1, 2] = value
        K_λ[2, 2, 1, 1, 3, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 3, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 3, 2, 2, 2] = value
        K_λ[2, 2, 2, 1, 3, 2, 1, 1] = value
        K_λ[2, 2, 2, 1, 3, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 3, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 3, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 1, 2, 2] = value
        K_λ[3, 1, 1, 2, 2, 1, 1, 1] = value
        K_λ[3, 1, 1, 2, 2, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 2, 1, 2, 1] = value
        K_λ[3, 1, 2, 1, 2, 1, 2, 2] = value
        K_λ[3, 1, 2, 2, 2, 1, 1, 1] = value
        K_λ[3, 1, 2, 2, 2, 1, 1, 2] = value
        K_λ[3, 2, 1, 1, 2, 2, 2, 1] = value
        K_λ[3, 2, 1, 1, 2, 2, 2, 2] = value
        K_λ[3, 2, 1, 2, 2, 2, 1, 1] = value
        K_λ[3, 2, 1, 2, 2, 2, 1, 2] = value
        K_λ[3, 2, 2, 1, 2, 2, 2, 1] = value
        K_λ[3, 2, 2, 1, 2, 2, 2, 2] = value
        K_λ[3, 2, 2, 2, 2, 2, 1, 1] = value
        K_λ[3, 2, 2, 2, 2, 2, 1, 2] = value

        value = -h[1] / 24
        K_λ[2, 1, 1, 1, 3, 2, 1, 2] = value
        K_λ[2, 1, 1, 1, 3, 2, 2, 2] = value
        K_λ[2, 1, 1, 2, 3, 2, 1, 2] = value
        K_λ[2, 1, 1, 2, 3, 2, 2, 2] = value
        K_λ[2, 1, 2, 1, 3, 2, 1, 1] = value
        K_λ[2, 1, 2, 1, 3, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 3, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 3, 2, 2, 1] = value
        K_λ[2, 2, 1, 1, 3, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 3, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 3, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 3, 1, 2, 2] = value
        K_λ[2, 2, 2, 1, 3, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 3, 1, 2, 1] = value
        K_λ[2, 2, 2, 2, 3, 1, 1, 1] = value
        K_λ[2, 2, 2, 2, 3, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 2, 2, 2] = value
        K_λ[3, 1, 1, 2, 2, 2, 1, 1] = value
        K_λ[3, 1, 1, 2, 2, 2, 1, 2] = value
        K_λ[3, 1, 2, 1, 2, 2, 2, 1] = value
        K_λ[3, 1, 2, 1, 2, 2, 2, 2] = value
        K_λ[3, 1, 2, 2, 2, 2, 1, 1] = value
        K_λ[3, 1, 2, 2, 2, 2, 1, 2] = value
        K_λ[3, 2, 1, 1, 2, 1, 2, 1] = value
        K_λ[3, 2, 1, 1, 2, 1, 2, 2] = value
        K_λ[3, 2, 1, 2, 2, 1, 1, 1] = value
        K_λ[3, 2, 1, 2, 2, 1, 1, 2] = value
        K_λ[3, 2, 2, 1, 2, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 2, 1, 2, 2] = value
        K_λ[3, 2, 2, 2, 2, 1, 1, 1] = value
        K_λ[3, 2, 2, 2, 2, 1, 1, 2] = value

        value = h[1] / 24
        K_λ[2, 1, 1, 1, 3, 2, 1, 1] = value
        K_λ[2, 1, 1, 1, 3, 2, 2, 1] = value
        K_λ[2, 1, 1, 2, 3, 2, 1, 1] = value
        K_λ[2, 1, 1, 2, 3, 2, 2, 1] = value
        K_λ[2, 1, 2, 1, 3, 2, 1, 2] = value
        K_λ[2, 1, 2, 1, 3, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 3, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 3, 2, 2, 2] = value
        K_λ[2, 2, 1, 1, 3, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 3, 1, 2, 1] = value
        K_λ[2, 2, 1, 2, 3, 1, 1, 1] = value
        K_λ[2, 2, 1, 2, 3, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 3, 1, 1, 2] = value
        K_λ[2, 2, 2, 1, 3, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 3, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 3, 1, 2, 2] = value
        K_λ[3, 1, 1, 1, 2, 2, 1, 1] = value
        K_λ[3, 1, 1, 1, 2, 2, 1, 2] = value
        K_λ[3, 1, 1, 2, 2, 2, 2, 1] = value
        K_λ[3, 1, 1, 2, 2, 2, 2, 2] = value
        K_λ[3, 1, 2, 1, 2, 2, 1, 1] = value
        K_λ[3, 1, 2, 1, 2, 2, 1, 2] = value
        K_λ[3, 1, 2, 2, 2, 2, 2, 1] = value
        K_λ[3, 1, 2, 2, 2, 2, 2, 2] = value
        K_λ[3, 2, 1, 1, 2, 1, 1, 1] = value
        K_λ[3, 2, 1, 1, 2, 1, 1, 2] = value
        K_λ[3, 2, 1, 2, 2, 1, 2, 1] = value
        K_λ[3, 2, 1, 2, 2, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 2, 1, 1, 1] = value
        K_λ[3, 2, 2, 1, 2, 1, 1, 2] = value
        K_λ[3, 2, 2, 2, 2, 1, 2, 1] = value
        K_λ[3, 2, 2, 2, 2, 1, 2, 2] = value

        value = h[1] / 12
        K_λ[2, 1, 1, 1, 3, 1, 1, 1] = value
        K_λ[2, 1, 1, 1, 3, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 3, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 3, 1, 2, 1] = value
        K_λ[2, 1, 2, 1, 3, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 3, 1, 2, 2] = value
        K_λ[2, 1, 2, 2, 3, 1, 1, 2] = value
        K_λ[2, 1, 2, 2, 3, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 3, 2, 1, 1] = value
        K_λ[2, 2, 1, 1, 3, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 3, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 3, 2, 2, 1] = value
        K_λ[2, 2, 2, 1, 3, 2, 1, 2] = value
        K_λ[2, 2, 2, 1, 3, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 3, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 3, 2, 2, 2] = value
        K_λ[3, 1, 1, 1, 2, 1, 1, 1] = value
        K_λ[3, 1, 1, 1, 2, 1, 1, 2] = value
        K_λ[3, 1, 1, 2, 2, 1, 2, 1] = value
        K_λ[3, 1, 1, 2, 2, 1, 2, 2] = value
        K_λ[3, 1, 2, 1, 2, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 2, 1, 1, 2] = value
        K_λ[3, 1, 2, 2, 2, 1, 2, 1] = value
        K_λ[3, 1, 2, 2, 2, 1, 2, 2] = value
        K_λ[3, 2, 1, 1, 2, 2, 1, 1] = value
        K_λ[3, 2, 1, 1, 2, 2, 1, 2] = value
        K_λ[3, 2, 1, 2, 2, 2, 2, 1] = value
        K_λ[3, 2, 1, 2, 2, 2, 2, 2] = value
        K_λ[3, 2, 2, 1, 2, 2, 1, 1] = value
        K_λ[3, 2, 2, 1, 2, 2, 1, 2] = value
        K_λ[3, 2, 2, 2, 2, 2, 2, 1] = value
        K_λ[3, 2, 2, 2, 2, 2, 2, 2] = value

        value = -h[2] / 12
        K_λ[1, 1, 1, 1, 3, 1, 1, 2] = value
        K_λ[1, 1, 1, 1, 3, 2, 1, 2] = value
        K_λ[1, 1, 1, 2, 3, 1, 1, 2] = value
        K_λ[1, 1, 1, 2, 3, 2, 1, 2] = value
        K_λ[1, 1, 2, 1, 3, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 3, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 3, 1, 2, 2] = value
        K_λ[1, 1, 2, 2, 3, 2, 2, 2] = value
        K_λ[1, 2, 1, 1, 3, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 3, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 3, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 3, 2, 1, 1] = value
        K_λ[1, 2, 2, 1, 3, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 3, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 3, 1, 2, 1] = value
        K_λ[1, 2, 2, 2, 3, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 1, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 1, 2] = value
        K_λ[3, 1, 1, 2, 1, 1, 1, 1] = value
        K_λ[3, 1, 1, 2, 1, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 1, 2, 2, 1] = value
        K_λ[3, 1, 2, 1, 1, 2, 2, 2] = value
        K_λ[3, 1, 2, 2, 1, 1, 2, 1] = value
        K_λ[3, 1, 2, 2, 1, 1, 2, 2] = value
        K_λ[3, 2, 1, 1, 1, 2, 1, 1] = value
        K_λ[3, 2, 1, 1, 1, 2, 1, 2] = value
        K_λ[3, 2, 1, 2, 1, 1, 1, 1] = value
        K_λ[3, 2, 1, 2, 1, 1, 1, 2] = value
        K_λ[3, 2, 2, 1, 1, 2, 2, 1] = value
        K_λ[3, 2, 2, 1, 1, 2, 2, 2] = value
        K_λ[3, 2, 2, 2, 1, 1, 2, 1] = value
        K_λ[3, 2, 2, 2, 1, 1, 2, 2] = value

        value = -h[2] / 24
        K_λ[1, 1, 1, 1, 3, 1, 2, 2] = value
        K_λ[1, 1, 1, 1, 3, 2, 2, 2] = value
        K_λ[1, 1, 1, 2, 3, 1, 2, 2] = value
        K_λ[1, 1, 1, 2, 3, 2, 2, 2] = value
        K_λ[1, 1, 2, 1, 3, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 3, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 3, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 3, 2, 1, 2] = value
        K_λ[1, 2, 1, 1, 3, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 3, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 3, 1, 2, 1] = value
        K_λ[1, 2, 1, 2, 3, 2, 2, 1] = value
        K_λ[1, 2, 2, 1, 3, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 3, 2, 1, 1] = value
        K_λ[1, 2, 2, 2, 3, 1, 1, 1] = value
        K_λ[1, 2, 2, 2, 3, 2, 1, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 2, 2] = value
        K_λ[3, 1, 1, 2, 1, 1, 2, 1] = value
        K_λ[3, 1, 1, 2, 1, 1, 2, 2] = value
        K_λ[3, 1, 2, 1, 1, 2, 1, 1] = value
        K_λ[3, 1, 2, 1, 1, 2, 1, 2] = value
        K_λ[3, 1, 2, 2, 1, 1, 1, 1] = value
        K_λ[3, 1, 2, 2, 1, 1, 1, 2] = value
        K_λ[3, 2, 1, 1, 1, 2, 2, 1] = value
        K_λ[3, 2, 1, 1, 1, 2, 2, 2] = value
        K_λ[3, 2, 1, 2, 1, 1, 2, 1] = value
        K_λ[3, 2, 1, 2, 1, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 1, 2, 1, 1] = value
        K_λ[3, 2, 2, 1, 1, 2, 1, 2] = value
        K_λ[3, 2, 2, 2, 1, 1, 1, 1] = value
        K_λ[3, 2, 2, 2, 1, 1, 1, 2] = value

        value = h[2] / 24
        K_λ[1, 1, 1, 1, 3, 1, 2, 1] = value
        K_λ[1, 1, 1, 1, 3, 2, 2, 1] = value
        K_λ[1, 1, 1, 2, 3, 1, 2, 1] = value
        K_λ[1, 1, 1, 2, 3, 2, 2, 1] = value
        K_λ[1, 1, 2, 1, 3, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 3, 2, 1, 1] = value
        K_λ[1, 1, 2, 2, 3, 1, 1, 1] = value
        K_λ[1, 1, 2, 2, 3, 2, 1, 1] = value
        K_λ[1, 2, 1, 1, 3, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 3, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 3, 1, 2, 2] = value
        K_λ[1, 2, 1, 2, 3, 2, 2, 2] = value
        K_λ[1, 2, 2, 1, 3, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 3, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 3, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 3, 2, 1, 2] = value
        K_λ[3, 1, 1, 1, 1, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 1, 1, 2, 2] = value
        K_λ[3, 1, 1, 2, 1, 2, 2, 1] = value
        K_λ[3, 1, 1, 2, 1, 2, 2, 2] = value
        K_λ[3, 1, 2, 1, 1, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 1, 1, 1, 2] = value
        K_λ[3, 1, 2, 2, 1, 2, 1, 1] = value
        K_λ[3, 1, 2, 2, 1, 2, 1, 2] = value
        K_λ[3, 2, 1, 1, 1, 1, 2, 1] = value
        K_λ[3, 2, 1, 1, 1, 1, 2, 2] = value
        K_λ[3, 2, 1, 2, 1, 2, 2, 1] = value
        K_λ[3, 2, 1, 2, 1, 2, 2, 2] = value
        K_λ[3, 2, 2, 1, 1, 1, 1, 1] = value
        K_λ[3, 2, 2, 1, 1, 1, 1, 2] = value
        K_λ[3, 2, 2, 2, 1, 2, 1, 1] = value
        K_λ[3, 2, 2, 2, 1, 2, 1, 2] = value

        value = h[2] / 12
        K_λ[1, 1, 1, 1, 3, 1, 1, 1] = value
        K_λ[1, 1, 1, 1, 3, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 3, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 3, 2, 1, 1] = value
        K_λ[1, 1, 2, 1, 3, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 3, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 3, 1, 2, 1] = value
        K_λ[1, 1, 2, 2, 3, 2, 2, 1] = value
        K_λ[1, 2, 1, 1, 3, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 3, 2, 1, 2] = value
        K_λ[1, 2, 1, 2, 3, 1, 1, 2] = value
        K_λ[1, 2, 1, 2, 3, 2, 1, 2] = value
        K_λ[1, 2, 2, 1, 3, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 3, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 3, 1, 2, 2] = value
        K_λ[1, 2, 2, 2, 3, 2, 2, 2] = value
        K_λ[3, 1, 1, 1, 1, 1, 1, 1] = value
        K_λ[3, 1, 1, 1, 1, 1, 1, 2] = value
        K_λ[3, 1, 1, 2, 1, 2, 1, 1] = value
        K_λ[3, 1, 1, 2, 1, 2, 1, 2] = value
        K_λ[3, 1, 2, 1, 1, 1, 2, 1] = value
        K_λ[3, 1, 2, 1, 1, 1, 2, 2] = value
        K_λ[3, 1, 2, 2, 1, 2, 2, 1] = value
        K_λ[3, 1, 2, 2, 1, 2, 2, 2] = value
        K_λ[3, 2, 1, 1, 1, 1, 1, 1] = value
        K_λ[3, 2, 1, 1, 1, 1, 1, 2] = value
        K_λ[3, 2, 1, 2, 1, 2, 1, 1] = value
        K_λ[3, 2, 1, 2, 1, 2, 1, 2] = value
        K_λ[3, 2, 2, 1, 1, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 1, 1, 2, 2] = value
        K_λ[3, 2, 2, 2, 1, 2, 2, 1] = value
        K_λ[3, 2, 2, 2, 1, 2, 2, 2] = value

        value = -(h[1] * h[2]) / (9 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 1, 2] = value
        K_λ[3, 1, 1, 2, 3, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 3, 1, 2, 2] = value
        K_λ[3, 1, 2, 2, 3, 1, 2, 1] = value
        K_λ[3, 2, 1, 1, 3, 2, 1, 2] = value
        K_λ[3, 2, 1, 2, 3, 2, 1, 1] = value
        K_λ[3, 2, 2, 1, 3, 2, 2, 2] = value
        K_λ[3, 2, 2, 2, 3, 2, 2, 1] = value

        value = -(h[1] * h[2]) / (18 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 2, 2] = value
        K_λ[3, 1, 1, 1, 3, 2, 1, 2] = value
        K_λ[3, 1, 1, 2, 3, 1, 2, 1] = value
        K_λ[3, 1, 1, 2, 3, 2, 1, 1] = value
        K_λ[3, 1, 2, 1, 3, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 3, 2, 2, 2] = value
        K_λ[3, 1, 2, 2, 3, 1, 1, 1] = value
        K_λ[3, 1, 2, 2, 3, 2, 2, 1] = value
        K_λ[3, 2, 1, 1, 3, 1, 1, 2] = value
        K_λ[3, 2, 1, 1, 3, 2, 2, 2] = value
        K_λ[3, 2, 1, 2, 3, 1, 1, 1] = value
        K_λ[3, 2, 1, 2, 3, 2, 2, 1] = value
        K_λ[3, 2, 2, 1, 3, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 3, 2, 1, 2] = value
        K_λ[3, 2, 2, 2, 3, 1, 2, 1] = value
        K_λ[3, 2, 2, 2, 3, 2, 1, 1] = value

        value = -(h[1] * h[2]) / (36 * h[3])
        K_λ[3, 1, 1, 1, 3, 2, 2, 2] = value
        K_λ[3, 1, 1, 2, 3, 2, 2, 1] = value
        K_λ[3, 1, 2, 1, 3, 2, 1, 2] = value
        K_λ[3, 1, 2, 2, 3, 2, 1, 1] = value
        K_λ[3, 2, 1, 1, 3, 1, 2, 2] = value
        K_λ[3, 2, 1, 2, 3, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 3, 1, 1, 2] = value
        K_λ[3, 2, 2, 2, 3, 1, 1, 1] = value

        value = (h[1] * h[2]) / (36 * h[3])
        K_λ[3, 1, 1, 1, 3, 2, 2, 1] = value
        K_λ[3, 1, 1, 2, 3, 2, 2, 2] = value
        K_λ[3, 1, 2, 1, 3, 2, 1, 1] = value
        K_λ[3, 1, 2, 2, 3, 2, 1, 2] = value
        K_λ[3, 2, 1, 1, 3, 1, 2, 1] = value
        K_λ[3, 2, 1, 2, 3, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 3, 1, 1, 1] = value
        K_λ[3, 2, 2, 2, 3, 1, 1, 2] = value

        value = (h[1] * h[2]) / (18 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 3, 2, 1, 1] = value
        K_λ[3, 1, 1, 2, 3, 1, 2, 2] = value
        K_λ[3, 1, 1, 2, 3, 2, 1, 2] = value
        K_λ[3, 1, 2, 1, 3, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 3, 2, 2, 1] = value
        K_λ[3, 1, 2, 2, 3, 1, 1, 2] = value
        K_λ[3, 1, 2, 2, 3, 2, 2, 2] = value
        K_λ[3, 2, 1, 1, 3, 1, 1, 1] = value
        K_λ[3, 2, 1, 1, 3, 2, 2, 1] = value
        K_λ[3, 2, 1, 2, 3, 1, 1, 2] = value
        K_λ[3, 2, 1, 2, 3, 2, 2, 2] = value
        K_λ[3, 2, 2, 1, 3, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 3, 2, 1, 1] = value
        K_λ[3, 2, 2, 2, 3, 1, 2, 2] = value
        K_λ[3, 2, 2, 2, 3, 2, 1, 2] = value

        value = (h[1] * h[2]) / (9 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 1, 1] = value
        K_λ[3, 1, 1, 2, 3, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 3, 1, 2, 1] = value
        K_λ[3, 1, 2, 2, 3, 1, 2, 2] = value
        K_λ[3, 2, 1, 1, 3, 2, 1, 1] = value
        K_λ[3, 2, 1, 2, 3, 2, 1, 2] = value
        K_λ[3, 2, 2, 1, 3, 2, 2, 1] = value
        K_λ[3, 2, 2, 2, 3, 2, 2, 2] = value

        value = -h[3] / 12
        K_λ[1, 1, 1, 1, 2, 1, 2, 1] = value
        K_λ[1, 1, 1, 1, 2, 2, 2, 1] = value
        K_λ[1, 1, 1, 2, 2, 1, 2, 2] = value
        K_λ[1, 1, 1, 2, 2, 2, 2, 2] = value
        K_λ[1, 1, 2, 1, 2, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 2, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 2, 1, 2, 2] = value
        K_λ[1, 1, 2, 2, 2, 2, 2, 2] = value
        K_λ[1, 2, 1, 1, 2, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 2, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 2, 1, 1, 2] = value
        K_λ[1, 2, 1, 2, 2, 2, 1, 2] = value
        K_λ[1, 2, 2, 1, 2, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 2, 2, 1, 1] = value
        K_λ[1, 2, 2, 2, 2, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 2, 2, 1, 2] = value
        K_λ[2, 1, 1, 1, 1, 2, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 2, 2, 1] = value
        K_λ[2, 1, 1, 2, 1, 2, 1, 2] = value
        K_λ[2, 1, 1, 2, 1, 2, 2, 2] = value
        K_λ[2, 1, 2, 1, 1, 1, 1, 1] = value
        K_λ[2, 1, 2, 1, 1, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 1, 1, 1, 2] = value
        K_λ[2, 1, 2, 2, 1, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 1, 2, 1, 1] = value
        K_λ[2, 2, 1, 1, 1, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 1, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 1, 2, 2, 2] = value
        K_λ[2, 2, 2, 1, 1, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 1, 1, 2, 1] = value
        K_λ[2, 2, 2, 2, 1, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 1, 1, 2, 2] = value

        value = -h[3] / 24
        K_λ[1, 1, 1, 1, 2, 1, 2, 2] = value
        K_λ[1, 1, 1, 1, 2, 2, 2, 2] = value
        K_λ[1, 1, 1, 2, 2, 1, 2, 1] = value
        K_λ[1, 1, 1, 2, 2, 2, 2, 1] = value
        K_λ[1, 1, 2, 1, 2, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 2, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 2, 1, 2, 1] = value
        K_λ[1, 1, 2, 2, 2, 2, 2, 1] = value
        K_λ[1, 2, 1, 1, 2, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 2, 2, 1, 2] = value
        K_λ[1, 2, 1, 2, 2, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 2, 2, 1, 1] = value
        K_λ[1, 2, 2, 1, 2, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 2, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 2, 1, 1, 1] = value
        K_λ[1, 2, 2, 2, 2, 2, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 2, 1, 2] = value
        K_λ[2, 1, 1, 1, 1, 2, 2, 2] = value
        K_λ[2, 1, 1, 2, 1, 2, 1, 1] = value
        K_λ[2, 1, 1, 2, 1, 2, 2, 1] = value
        K_λ[2, 1, 2, 1, 1, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 1, 1, 2, 2] = value
        K_λ[2, 1, 2, 2, 1, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 1, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 1, 2, 1, 2] = value
        K_λ[2, 2, 1, 1, 1, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 1, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 1, 2, 2, 1] = value
        K_λ[2, 2, 2, 1, 1, 1, 1, 2] = value
        K_λ[2, 2, 2, 1, 1, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 1, 1, 1, 1] = value
        K_λ[2, 2, 2, 2, 1, 1, 2, 1] = value

        value = h[3] / 24
        K_λ[1, 1, 1, 1, 2, 1, 1, 2] = value
        K_λ[1, 1, 1, 1, 2, 2, 1, 2] = value
        K_λ[1, 1, 1, 2, 2, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 2, 2, 1, 1] = value
        K_λ[1, 1, 2, 1, 2, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 2, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 2, 1, 1, 1] = value
        K_λ[1, 1, 2, 2, 2, 2, 1, 1] = value
        K_λ[1, 2, 1, 1, 2, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 2, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 2, 1, 2, 1] = value
        K_λ[1, 2, 1, 2, 2, 2, 2, 1] = value
        K_λ[1, 2, 2, 1, 2, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 2, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 2, 1, 2, 1] = value
        K_λ[1, 2, 2, 2, 2, 2, 2, 1] = value
        K_λ[2, 1, 1, 1, 1, 1, 1, 2] = value
        K_λ[2, 1, 1, 1, 1, 1, 2, 2] = value
        K_λ[2, 1, 1, 2, 1, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 1, 1, 2, 1] = value
        K_λ[2, 1, 2, 1, 1, 2, 1, 2] = value
        K_λ[2, 1, 2, 1, 1, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 1, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 1, 2, 2, 1] = value
        K_λ[2, 2, 1, 1, 1, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 1, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 1, 1, 1, 1] = value
        K_λ[2, 2, 1, 2, 1, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 1, 2, 1, 2] = value
        K_λ[2, 2, 2, 1, 1, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 1, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 1, 2, 2, 1] = value

        value = h[3] / 12
        K_λ[1, 1, 1, 1, 2, 1, 1, 1] = value
        K_λ[1, 1, 1, 1, 2, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 2, 1, 1, 2] = value
        K_λ[1, 1, 1, 2, 2, 2, 1, 2] = value
        K_λ[1, 1, 2, 1, 2, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 2, 2, 1, 1] = value
        K_λ[1, 1, 2, 2, 2, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 2, 2, 1, 2] = value
        K_λ[1, 2, 1, 1, 2, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 2, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 2, 1, 2, 2] = value
        K_λ[1, 2, 1, 2, 2, 2, 2, 2] = value
        K_λ[1, 2, 2, 1, 2, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 2, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 2, 1, 2, 2] = value
        K_λ[1, 2, 2, 2, 2, 2, 2, 2] = value
        K_λ[2, 1, 1, 1, 1, 1, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 1, 1, 1, 2] = value
        K_λ[2, 1, 1, 2, 1, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 1, 2, 1, 1] = value
        K_λ[2, 1, 2, 1, 1, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 1, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 1, 2, 2, 2] = value
        K_λ[2, 2, 1, 1, 1, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 1, 1, 2, 1] = value
        K_λ[2, 2, 1, 2, 1, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 1, 1, 2, 2] = value
        K_λ[2, 2, 2, 1, 1, 2, 1, 1] = value
        K_λ[2, 2, 2, 1, 1, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 1, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 1, 2, 2, 2] = value

        value = -(h[1] * h[3]) / (9 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 2, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 2, 2, 2, 2] = value
        K_λ[2, 2, 2, 1, 2, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 2, 2, 1, 2] = value

        value = -(h[1] * h[3]) / (18 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 2, 2] = value
        K_λ[2, 1, 1, 1, 2, 2, 2, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 2, 2, 2, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 2, 2, 1, 2] = value
        K_λ[2, 2, 1, 1, 2, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 2, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 2, 2, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 2, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 2, 2, 1, 1] = value

        value = -(h[1] * h[3]) / (36 * h[2])
        K_λ[2, 1, 1, 1, 2, 2, 2, 2] = value
        K_λ[2, 1, 1, 2, 2, 2, 2, 1] = value
        K_λ[2, 1, 2, 1, 2, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 2, 2, 1, 1] = value
        K_λ[2, 2, 1, 1, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 1, 1] = value

        value = (h[1] * h[3]) / (36 * h[2])
        K_λ[2, 1, 1, 1, 2, 2, 1, 2] = value
        K_λ[2, 1, 1, 2, 2, 2, 1, 1] = value
        K_λ[2, 1, 2, 1, 2, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 2, 2, 2, 1] = value
        K_λ[2, 2, 1, 1, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 2, 1] = value

        value = (h[1] * h[3]) / (18 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 1, 2] = value
        K_λ[2, 1, 1, 1, 2, 2, 1, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 2, 2, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 2, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 2, 2, 2] = value
        K_λ[2, 2, 1, 1, 2, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 2, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 2, 1, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 2, 2, 1] = value

        value = (h[1] * h[3]) / (9 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 2, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 2, 2, 1, 2] = value
        K_λ[2, 2, 2, 1, 2, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 2, 2, 2, 2] = value

        value = -(h[2] * h[3]) / (9 * h[1])
        K_λ[1, 1, 1, 1, 1, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 2, 1, 2] = value
        K_λ[1, 1, 2, 1, 1, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 1, 2, 2, 2] = value
        K_λ[1, 2, 1, 1, 1, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 1, 1, 2, 1] = value
        K_λ[1, 2, 2, 2, 1, 1, 2, 2] = value

        value = -(h[2] * h[3]) / (18 * h[1])
        K_λ[1, 1, 1, 1, 1, 2, 1, 2] = value
        K_λ[1, 1, 1, 1, 1, 2, 2, 1] = value
        K_λ[1, 1, 1, 2, 1, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 2, 2, 2] = value
        K_λ[1, 1, 2, 1, 1, 2, 1, 1] = value
        K_λ[1, 1, 2, 1, 1, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 1, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 2, 2, 1] = value
        K_λ[1, 2, 1, 1, 1, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 1, 1, 2, 1] = value
        K_λ[1, 2, 1, 2, 1, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 1, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 1, 1, 2, 2] = value
        K_λ[1, 2, 2, 2, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 1, 2, 1] = value

        value = -(h[2] * h[3]) / (36 * h[1])
        K_λ[1, 1, 1, 1, 1, 2, 2, 2] = value
        K_λ[1, 1, 1, 2, 1, 2, 2, 1] = value
        K_λ[1, 1, 2, 1, 1, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 1, 1, 1, 2, 2] = value
        K_λ[1, 2, 1, 2, 1, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 1, 1, 1] = value

        value = (h[2] * h[3]) / (36 * h[1])
        K_λ[1, 1, 1, 1, 1, 1, 2, 2] = value
        K_λ[1, 1, 1, 2, 1, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 1, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 1, 2, 2, 1] = value
        K_λ[1, 2, 2, 1, 1, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 2, 1, 1] = value

        value = (h[2] * h[3]) / (18 * h[1])
        K_λ[1, 1, 1, 1, 1, 1, 1, 2] = value
        K_λ[1, 1, 1, 1, 1, 1, 2, 1] = value
        K_λ[1, 1, 1, 2, 1, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 1, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 1, 1, 2, 2] = value
        K_λ[1, 1, 2, 2, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 1, 2, 1, 2] = value
        K_λ[1, 2, 1, 1, 1, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 2, 2, 2] = value
        K_λ[1, 2, 2, 1, 1, 2, 1, 1] = value
        K_λ[1, 2, 2, 1, 1, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 1, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 2, 2, 1] = value

        value = (h[2] * h[3]) / (9 * h[1])
        K_λ[1, 1, 1, 1, 1, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 1, 1, 2, 1] = value
        K_λ[1, 1, 2, 2, 1, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 2, 1, 2] = value
        K_λ[1, 2, 2, 1, 1, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 1, 2, 2, 2] = value
    else
    end
    return K_λ
end

function expected_K_μ(h::NTuple{d,T}) where {d,T<:Number}
    K_μ = zeros(Float64, d, fill(2, d)..., d, fill(2, d)...)
    if d == 2
        value = -1 / 4
        K_μ[1, 1, 1, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 1, 2] = value
        K_μ[2, 1, 1, 1, 1, 2] = value
        K_μ[2, 1, 1, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 2, 1] = value

        value = 1 / 4
        K_μ[1, 1, 1, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 1] = value
        K_μ[2, 1, 1, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 2, 2] = value

        value = h[1] / (6 * h[2]) - (2 * h[2]) / (3 * h[1])
        K_μ[1, 1, 1, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 1, 2] = value

        value = (-h[2] / (3 * h[1])) - h[1] / (6 * h[2])
        K_μ[1, 1, 1, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 1, 1] = value

        value = h[1] / (3 * h[2]) - h[2] / (3 * h[1])
        K_μ[2, 1, 1, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 1, 2] = value

        value = (-h[2] / (6 * h[1])) - h[1] / (3 * h[2])
        K_μ[2, 1, 1, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 1, 1] = value

        value = h[2] / (6 * h[1]) - (2 * h[1]) / (3 * h[2])
        K_μ[2, 1, 1, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 2, 1] = value

        value = h[2] / (3 * h[1]) - h[1] / (3 * h[2])
        K_μ[1, 1, 1, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 2, 1] = value

        value = h[2] / (3 * h[1]) + (2 * h[1]) / (3 * h[2])
        K_μ[2, 1, 1, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 2, 2] = value

        value = (2 * h[2]) / (3 * h[1]) + h[1] / (3 * h[2])
        K_μ[1, 1, 1, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 2, 2] = value
    elseif d == 3
        value = -h[1] / 12
        K_μ[2, 1, 1, 1, 3, 1, 2, 1] = value
        K_μ[2, 1, 1, 1, 3, 1, 2, 2] = value
        K_μ[2, 1, 1, 2, 3, 1, 1, 1] = value
        K_μ[2, 1, 1, 2, 3, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 3, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 3, 1, 2, 2] = value
        K_μ[2, 1, 2, 2, 3, 1, 1, 1] = value
        K_μ[2, 1, 2, 2, 3, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 3, 2, 2, 1] = value
        K_μ[2, 2, 1, 1, 3, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 3, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 3, 2, 1, 2] = value
        K_μ[2, 2, 2, 1, 3, 2, 2, 1] = value
        K_μ[2, 2, 2, 1, 3, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 3, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 3, 2, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 1, 2, 2] = value
        K_μ[3, 1, 1, 2, 2, 1, 1, 2] = value
        K_μ[3, 1, 1, 2, 2, 1, 2, 2] = value
        K_μ[3, 1, 2, 1, 2, 1, 1, 1] = value
        K_μ[3, 1, 2, 1, 2, 1, 2, 1] = value
        K_μ[3, 1, 2, 2, 2, 1, 1, 1] = value
        K_μ[3, 1, 2, 2, 2, 1, 2, 1] = value
        K_μ[3, 2, 1, 1, 2, 2, 1, 2] = value
        K_μ[3, 2, 1, 1, 2, 2, 2, 2] = value
        K_μ[3, 2, 1, 2, 2, 2, 1, 2] = value
        K_μ[3, 2, 1, 2, 2, 2, 2, 2] = value
        K_μ[3, 2, 2, 1, 2, 2, 1, 1] = value
        K_μ[3, 2, 2, 1, 2, 2, 2, 1] = value
        K_μ[3, 2, 2, 2, 2, 2, 1, 1] = value
        K_μ[3, 2, 2, 2, 2, 2, 2, 1] = value

        value = -h[1] / 24
        K_μ[2, 1, 1, 1, 3, 2, 2, 1] = value
        K_μ[2, 1, 1, 1, 3, 2, 2, 2] = value
        K_μ[2, 1, 1, 2, 3, 2, 1, 1] = value
        K_μ[2, 1, 1, 2, 3, 2, 1, 2] = value
        K_μ[2, 1, 2, 1, 3, 2, 2, 1] = value
        K_μ[2, 1, 2, 1, 3, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 3, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 3, 2, 1, 2] = value
        K_μ[2, 2, 1, 1, 3, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 3, 1, 2, 2] = value
        K_μ[2, 2, 1, 2, 3, 1, 1, 1] = value
        K_μ[2, 2, 1, 2, 3, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 3, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 3, 1, 2, 2] = value
        K_μ[2, 2, 2, 2, 3, 1, 1, 1] = value
        K_μ[2, 2, 2, 2, 3, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 2, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 2, 2, 2] = value
        K_μ[3, 1, 1, 2, 2, 2, 1, 2] = value
        K_μ[3, 1, 1, 2, 2, 2, 2, 2] = value
        K_μ[3, 1, 2, 1, 2, 2, 1, 1] = value
        K_μ[3, 1, 2, 1, 2, 2, 2, 1] = value
        K_μ[3, 1, 2, 2, 2, 2, 1, 1] = value
        K_μ[3, 1, 2, 2, 2, 2, 2, 1] = value
        K_μ[3, 2, 1, 1, 2, 1, 1, 2] = value
        K_μ[3, 2, 1, 1, 2, 1, 2, 2] = value
        K_μ[3, 2, 1, 2, 2, 1, 1, 2] = value
        K_μ[3, 2, 1, 2, 2, 1, 2, 2] = value
        K_μ[3, 2, 2, 1, 2, 1, 1, 1] = value
        K_μ[3, 2, 2, 1, 2, 1, 2, 1] = value
        K_μ[3, 2, 2, 2, 2, 1, 1, 1] = value
        K_μ[3, 2, 2, 2, 2, 1, 2, 1] = value

        value = h[1] / 24
        K_μ[2, 1, 1, 1, 3, 2, 1, 1] = value
        K_μ[2, 1, 1, 1, 3, 2, 1, 2] = value
        K_μ[2, 1, 1, 2, 3, 2, 2, 1] = value
        K_μ[2, 1, 1, 2, 3, 2, 2, 2] = value
        K_μ[2, 1, 2, 1, 3, 2, 1, 1] = value
        K_μ[2, 1, 2, 1, 3, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 3, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 3, 2, 2, 2] = value
        K_μ[2, 2, 1, 1, 3, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 3, 1, 1, 2] = value
        K_μ[2, 2, 1, 2, 3, 1, 2, 1] = value
        K_μ[2, 2, 1, 2, 3, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 3, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 3, 1, 1, 2] = value
        K_μ[2, 2, 2, 2, 3, 1, 2, 1] = value
        K_μ[2, 2, 2, 2, 3, 1, 2, 2] = value
        K_μ[3, 1, 1, 1, 2, 2, 1, 1] = value
        K_μ[3, 1, 1, 1, 2, 2, 2, 1] = value
        K_μ[3, 1, 1, 2, 2, 2, 1, 1] = value
        K_μ[3, 1, 1, 2, 2, 2, 2, 1] = value
        K_μ[3, 1, 2, 1, 2, 2, 1, 2] = value
        K_μ[3, 1, 2, 1, 2, 2, 2, 2] = value
        K_μ[3, 1, 2, 2, 2, 2, 1, 2] = value
        K_μ[3, 1, 2, 2, 2, 2, 2, 2] = value
        K_μ[3, 2, 1, 1, 2, 1, 1, 1] = value
        K_μ[3, 2, 1, 1, 2, 1, 2, 1] = value
        K_μ[3, 2, 1, 2, 2, 1, 1, 1] = value
        K_μ[3, 2, 1, 2, 2, 1, 2, 1] = value
        K_μ[3, 2, 2, 1, 2, 1, 1, 2] = value
        K_μ[3, 2, 2, 1, 2, 1, 2, 2] = value
        K_μ[3, 2, 2, 2, 2, 1, 1, 2] = value
        K_μ[3, 2, 2, 2, 2, 1, 2, 2] = value

        value = h[1] / 12
        K_μ[2, 1, 1, 1, 3, 1, 1, 1] = value
        K_μ[2, 1, 1, 1, 3, 1, 1, 2] = value
        K_μ[2, 1, 1, 2, 3, 1, 2, 1] = value
        K_μ[2, 1, 1, 2, 3, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 3, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 3, 1, 1, 2] = value
        K_μ[2, 1, 2, 2, 3, 1, 2, 1] = value
        K_μ[2, 1, 2, 2, 3, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 3, 2, 1, 1] = value
        K_μ[2, 2, 1, 1, 3, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 3, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 3, 2, 2, 2] = value
        K_μ[2, 2, 2, 1, 3, 2, 1, 1] = value
        K_μ[2, 2, 2, 1, 3, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 3, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 3, 2, 2, 2] = value
        K_μ[3, 1, 1, 1, 2, 1, 1, 1] = value
        K_μ[3, 1, 1, 1, 2, 1, 2, 1] = value
        K_μ[3, 1, 1, 2, 2, 1, 1, 1] = value
        K_μ[3, 1, 1, 2, 2, 1, 2, 1] = value
        K_μ[3, 1, 2, 1, 2, 1, 1, 2] = value
        K_μ[3, 1, 2, 1, 2, 1, 2, 2] = value
        K_μ[3, 1, 2, 2, 2, 1, 1, 2] = value
        K_μ[3, 1, 2, 2, 2, 1, 2, 2] = value
        K_μ[3, 2, 1, 1, 2, 2, 1, 1] = value
        K_μ[3, 2, 1, 1, 2, 2, 2, 1] = value
        K_μ[3, 2, 1, 2, 2, 2, 1, 1] = value
        K_μ[3, 2, 1, 2, 2, 2, 2, 1] = value
        K_μ[3, 2, 2, 1, 2, 2, 1, 2] = value
        K_μ[3, 2, 2, 1, 2, 2, 2, 2] = value
        K_μ[3, 2, 2, 2, 2, 2, 1, 2] = value
        K_μ[3, 2, 2, 2, 2, 2, 2, 2] = value

        value = -h[2] / 12
        K_μ[1, 1, 1, 1, 3, 2, 1, 1] = value
        K_μ[1, 1, 1, 1, 3, 2, 1, 2] = value
        K_μ[1, 1, 1, 2, 3, 1, 1, 1] = value
        K_μ[1, 1, 1, 2, 3, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 3, 2, 2, 1] = value
        K_μ[1, 1, 2, 1, 3, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 3, 1, 2, 1] = value
        K_μ[1, 1, 2, 2, 3, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 3, 2, 1, 1] = value
        K_μ[1, 2, 1, 1, 3, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 3, 1, 1, 1] = value
        K_μ[1, 2, 1, 2, 3, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 3, 2, 2, 1] = value
        K_μ[1, 2, 2, 1, 3, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 3, 1, 2, 1] = value
        K_μ[1, 2, 2, 2, 3, 1, 2, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 1, 2, 1, 2] = value
        K_μ[3, 1, 1, 2, 1, 1, 1, 2] = value
        K_μ[3, 1, 1, 2, 1, 2, 1, 2] = value
        K_μ[3, 1, 2, 1, 1, 1, 2, 2] = value
        K_μ[3, 1, 2, 1, 1, 2, 2, 2] = value
        K_μ[3, 1, 2, 2, 1, 1, 2, 2] = value
        K_μ[3, 1, 2, 2, 1, 2, 2, 2] = value
        K_μ[3, 2, 1, 1, 1, 1, 1, 1] = value
        K_μ[3, 2, 1, 1, 1, 2, 1, 1] = value
        K_μ[3, 2, 1, 2, 1, 1, 1, 1] = value
        K_μ[3, 2, 1, 2, 1, 2, 1, 1] = value
        K_μ[3, 2, 2, 1, 1, 1, 2, 1] = value
        K_μ[3, 2, 2, 1, 1, 2, 2, 1] = value
        K_μ[3, 2, 2, 2, 1, 1, 2, 1] = value
        K_μ[3, 2, 2, 2, 1, 2, 2, 1] = value

        value = -h[2] / 24
        K_μ[1, 1, 1, 1, 3, 2, 2, 1] = value
        K_μ[1, 1, 1, 1, 3, 2, 2, 2] = value
        K_μ[1, 1, 1, 2, 3, 1, 2, 1] = value
        K_μ[1, 1, 1, 2, 3, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 3, 2, 1, 1] = value
        K_μ[1, 1, 2, 1, 3, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 3, 1, 1, 1] = value
        K_μ[1, 1, 2, 2, 3, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 3, 2, 2, 1] = value
        K_μ[1, 2, 1, 1, 3, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 3, 1, 2, 1] = value
        K_μ[1, 2, 1, 2, 3, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 3, 2, 1, 1] = value
        K_μ[1, 2, 2, 1, 3, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 3, 1, 1, 1] = value
        K_μ[1, 2, 2, 2, 3, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 2, 2] = value
        K_μ[3, 1, 1, 1, 1, 2, 2, 2] = value
        K_μ[3, 1, 1, 2, 1, 1, 2, 2] = value
        K_μ[3, 1, 1, 2, 1, 2, 2, 2] = value
        K_μ[3, 1, 2, 1, 1, 1, 1, 2] = value
        K_μ[3, 1, 2, 1, 1, 2, 1, 2] = value
        K_μ[3, 1, 2, 2, 1, 1, 1, 2] = value
        K_μ[3, 1, 2, 2, 1, 2, 1, 2] = value
        K_μ[3, 2, 1, 1, 1, 1, 2, 1] = value
        K_μ[3, 2, 1, 1, 1, 2, 2, 1] = value
        K_μ[3, 2, 1, 2, 1, 1, 2, 1] = value
        K_μ[3, 2, 1, 2, 1, 2, 2, 1] = value
        K_μ[3, 2, 2, 1, 1, 1, 1, 1] = value
        K_μ[3, 2, 2, 1, 1, 2, 1, 1] = value
        K_μ[3, 2, 2, 2, 1, 1, 1, 1] = value
        K_μ[3, 2, 2, 2, 1, 2, 1, 1] = value

        value = h[2] / 24
        K_μ[1, 1, 1, 1, 3, 1, 2, 1] = value
        K_μ[1, 1, 1, 1, 3, 1, 2, 2] = value
        K_μ[1, 1, 1, 2, 3, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 3, 2, 2, 2] = value
        K_μ[1, 1, 2, 1, 3, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 3, 1, 1, 2] = value
        K_μ[1, 1, 2, 2, 3, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 3, 2, 1, 2] = value
        K_μ[1, 2, 1, 1, 3, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 3, 1, 2, 2] = value
        K_μ[1, 2, 1, 2, 3, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 3, 2, 2, 2] = value
        K_μ[1, 2, 2, 1, 3, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 3, 1, 1, 2] = value
        K_μ[1, 2, 2, 2, 3, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 3, 2, 1, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 2, 1] = value
        K_μ[3, 1, 1, 1, 1, 2, 2, 1] = value
        K_μ[3, 1, 1, 2, 1, 1, 2, 1] = value
        K_μ[3, 1, 1, 2, 1, 2, 2, 1] = value
        K_μ[3, 1, 2, 1, 1, 1, 1, 1] = value
        K_μ[3, 1, 2, 1, 1, 2, 1, 1] = value
        K_μ[3, 1, 2, 2, 1, 1, 1, 1] = value
        K_μ[3, 1, 2, 2, 1, 2, 1, 1] = value
        K_μ[3, 2, 1, 1, 1, 1, 2, 2] = value
        K_μ[3, 2, 1, 1, 1, 2, 2, 2] = value
        K_μ[3, 2, 1, 2, 1, 1, 2, 2] = value
        K_μ[3, 2, 1, 2, 1, 2, 2, 2] = value
        K_μ[3, 2, 2, 1, 1, 1, 1, 2] = value
        K_μ[3, 2, 2, 1, 1, 2, 1, 2] = value
        K_μ[3, 2, 2, 2, 1, 1, 1, 2] = value
        K_μ[3, 2, 2, 2, 1, 2, 1, 2] = value

        value = h[2] / 12
        K_μ[1, 1, 1, 1, 3, 1, 1, 1] = value
        K_μ[1, 1, 1, 1, 3, 1, 1, 2] = value
        K_μ[1, 1, 1, 2, 3, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 3, 2, 1, 2] = value
        K_μ[1, 1, 2, 1, 3, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 3, 1, 2, 2] = value
        K_μ[1, 1, 2, 2, 3, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 3, 2, 2, 2] = value
        K_μ[1, 2, 1, 1, 3, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 3, 1, 1, 2] = value
        K_μ[1, 2, 1, 2, 3, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 3, 2, 1, 2] = value
        K_μ[1, 2, 2, 1, 3, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 3, 1, 2, 2] = value
        K_μ[1, 2, 2, 2, 3, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 3, 2, 2, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 1, 1] = value
        K_μ[3, 1, 1, 1, 1, 2, 1, 1] = value
        K_μ[3, 1, 1, 2, 1, 1, 1, 1] = value
        K_μ[3, 1, 1, 2, 1, 2, 1, 1] = value
        K_μ[3, 1, 2, 1, 1, 1, 2, 1] = value
        K_μ[3, 1, 2, 1, 1, 2, 2, 1] = value
        K_μ[3, 1, 2, 2, 1, 1, 2, 1] = value
        K_μ[3, 1, 2, 2, 1, 2, 2, 1] = value
        K_μ[3, 2, 1, 1, 1, 1, 1, 2] = value
        K_μ[3, 2, 1, 1, 1, 2, 1, 2] = value
        K_μ[3, 2, 1, 2, 1, 1, 1, 2] = value
        K_μ[3, 2, 1, 2, 1, 2, 1, 2] = value
        K_μ[3, 2, 2, 1, 1, 1, 2, 2] = value
        K_μ[3, 2, 2, 1, 1, 2, 2, 2] = value
        K_μ[3, 2, 2, 2, 1, 1, 2, 2] = value
        K_μ[3, 2, 2, 2, 1, 2, 2, 2] = value

        value = -h[3] / 12
        K_μ[1, 1, 1, 1, 2, 2, 1, 1] = value
        K_μ[1, 1, 1, 1, 2, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 2, 2, 1, 2] = value
        K_μ[1, 1, 1, 2, 2, 2, 2, 2] = value
        K_μ[1, 1, 2, 1, 2, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 2, 1, 2, 1] = value
        K_μ[1, 1, 2, 2, 2, 1, 1, 2] = value
        K_μ[1, 1, 2, 2, 2, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 2, 2, 1, 1] = value
        K_μ[1, 2, 1, 1, 2, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 2, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 2, 2, 2, 2] = value
        K_μ[1, 2, 2, 1, 2, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 2, 1, 2, 1] = value
        K_μ[1, 2, 2, 2, 2, 1, 1, 2] = value
        K_μ[1, 2, 2, 2, 2, 1, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 1, 2, 1] = value
        K_μ[2, 1, 1, 1, 1, 2, 2, 1] = value
        K_μ[2, 1, 1, 2, 1, 1, 2, 2] = value
        K_μ[2, 1, 1, 2, 1, 2, 2, 2] = value
        K_μ[2, 1, 2, 1, 1, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 1, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 1, 1, 2, 2] = value
        K_μ[2, 1, 2, 2, 1, 2, 2, 2] = value
        K_μ[2, 2, 1, 1, 1, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 1, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 1, 1, 1, 2] = value
        K_μ[2, 2, 1, 2, 1, 2, 1, 2] = value
        K_μ[2, 2, 2, 1, 1, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 1, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 1, 1, 1, 2] = value
        K_μ[2, 2, 2, 2, 1, 2, 1, 2] = value

        value = -h[3] / 24
        K_μ[1, 1, 1, 1, 2, 2, 1, 2] = value
        K_μ[1, 1, 1, 1, 2, 2, 2, 2] = value
        K_μ[1, 1, 1, 2, 2, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 2, 2, 2, 1] = value
        K_μ[1, 1, 2, 1, 2, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 2, 1, 2, 2] = value
        K_μ[1, 1, 2, 2, 2, 1, 1, 1] = value
        K_μ[1, 1, 2, 2, 2, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 2, 2, 1, 2] = value
        K_μ[1, 2, 1, 1, 2, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 2, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 2, 2, 2, 1] = value
        K_μ[1, 2, 2, 1, 2, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 2, 1, 2, 2] = value
        K_μ[1, 2, 2, 2, 2, 1, 1, 1] = value
        K_μ[1, 2, 2, 2, 2, 1, 2, 1] = value
        K_μ[2, 1, 1, 1, 1, 1, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 2, 2, 2] = value
        K_μ[2, 1, 1, 2, 1, 1, 2, 1] = value
        K_μ[2, 1, 1, 2, 1, 2, 2, 1] = value
        K_μ[2, 1, 2, 1, 1, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 1, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 1, 1, 2, 1] = value
        K_μ[2, 1, 2, 2, 1, 2, 2, 1] = value
        K_μ[2, 2, 1, 1, 1, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 1, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 1, 1, 1, 1] = value
        K_μ[2, 2, 1, 2, 1, 2, 1, 1] = value
        K_μ[2, 2, 2, 1, 1, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 1, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 1, 1, 1, 1] = value
        K_μ[2, 2, 2, 2, 1, 2, 1, 1] = value

        value = h[3] / 24
        K_μ[1, 1, 1, 1, 2, 1, 1, 2] = value
        K_μ[1, 1, 1, 1, 2, 1, 2, 2] = value
        K_μ[1, 1, 1, 2, 2, 1, 1, 1] = value
        K_μ[1, 1, 1, 2, 2, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 2, 2, 1, 2] = value
        K_μ[1, 1, 2, 1, 2, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 2, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 2, 2, 2, 1] = value
        K_μ[1, 2, 1, 1, 2, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 2, 1, 2, 2] = value
        K_μ[1, 2, 1, 2, 2, 1, 1, 1] = value
        K_μ[1, 2, 1, 2, 2, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 2, 2, 1, 2] = value
        K_μ[1, 2, 2, 1, 2, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 2, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 2, 2, 2, 1] = value
        K_μ[2, 1, 1, 1, 1, 1, 1, 2] = value
        K_μ[2, 1, 1, 1, 1, 2, 1, 2] = value
        K_μ[2, 1, 1, 2, 1, 1, 1, 1] = value
        K_μ[2, 1, 1, 2, 1, 2, 1, 1] = value
        K_μ[2, 1, 2, 1, 1, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 1, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 1, 1, 1, 1] = value
        K_μ[2, 1, 2, 2, 1, 2, 1, 1] = value
        K_μ[2, 2, 1, 1, 1, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 1, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 1, 1, 2, 1] = value
        K_μ[2, 2, 1, 2, 1, 2, 2, 1] = value
        K_μ[2, 2, 2, 1, 1, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 1, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 1, 1, 2, 1] = value
        K_μ[2, 2, 2, 2, 1, 2, 2, 1] = value

        value = h[3] / 12
        K_μ[1, 1, 1, 1, 2, 1, 1, 1] = value
        K_μ[1, 1, 1, 1, 2, 1, 2, 1] = value
        K_μ[1, 1, 1, 2, 2, 1, 1, 2] = value
        K_μ[1, 1, 1, 2, 2, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 2, 2, 1, 1] = value
        K_μ[1, 1, 2, 1, 2, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 2, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 2, 2, 2, 2] = value
        K_μ[1, 2, 1, 1, 2, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 2, 1, 2, 1] = value
        K_μ[1, 2, 1, 2, 2, 1, 1, 2] = value
        K_μ[1, 2, 1, 2, 2, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 2, 2, 1, 1] = value
        K_μ[1, 2, 2, 1, 2, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 2, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 2, 2, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 1, 1, 1] = value
        K_μ[2, 1, 1, 1, 1, 2, 1, 1] = value
        K_μ[2, 1, 1, 2, 1, 1, 1, 2] = value
        K_μ[2, 1, 1, 2, 1, 2, 1, 2] = value
        K_μ[2, 1, 2, 1, 1, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 1, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 1, 1, 1, 2] = value
        K_μ[2, 1, 2, 2, 1, 2, 1, 2] = value
        K_μ[2, 2, 1, 1, 1, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 1, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 1, 1, 2, 2] = value
        K_μ[2, 2, 1, 2, 1, 2, 2, 2] = value
        K_μ[2, 2, 2, 1, 1, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 1, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 1, 1, 2, 2] = value
        K_μ[2, 2, 2, 2, 1, 2, 2, 2] = value

        value =
            (-(2 * h[2] * h[3]) / (9 * h[1])) +
            (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 1, 2, 1, 2] = value
        K_μ[1, 1, 2, 1, 1, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 1, 2, 2, 2] = value
        K_μ[1, 2, 1, 1, 1, 1, 1, 1] = value
        K_μ[1, 2, 1, 2, 1, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 1, 1, 2, 1] = value
        K_μ[1, 2, 2, 2, 1, 1, 2, 2] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) - (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (36 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 1, 2, 2, 2] = value
        K_μ[1, 1, 2, 1, 1, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 1, 2, 1, 2] = value
        K_μ[1, 2, 1, 1, 1, 1, 2, 1] = value
        K_μ[1, 2, 1, 2, 1, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 1, 1, 1, 1] = value
        K_μ[1, 2, 2, 2, 1, 1, 1, 2] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) + (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 1, 2] = value
        K_μ[1, 1, 1, 2, 1, 2, 1, 1] = value
        K_μ[1, 1, 2, 1, 1, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 1, 2, 2, 1] = value
        K_μ[1, 2, 1, 1, 1, 1, 1, 2] = value
        K_μ[1, 2, 1, 2, 1, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 1, 1, 2, 2] = value
        K_μ[1, 2, 2, 2, 1, 1, 2, 1] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) +
            (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 1, 1] = value
        K_μ[3, 1, 1, 2, 3, 2, 1, 2] = value
        K_μ[3, 1, 2, 1, 3, 2, 2, 1] = value
        K_μ[3, 1, 2, 2, 3, 2, 2, 2] = value
        K_μ[3, 2, 1, 1, 3, 1, 1, 1] = value
        K_μ[3, 2, 1, 2, 3, 1, 1, 2] = value
        K_μ[3, 2, 2, 1, 3, 1, 2, 1] = value
        K_μ[3, 2, 2, 2, 3, 1, 2, 2] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) +
            (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 1, 1] = value
        K_μ[2, 1, 1, 2, 2, 2, 1, 2] = value
        K_μ[2, 1, 2, 1, 2, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 2, 2, 2, 2] = value
        K_μ[2, 2, 1, 1, 2, 1, 1, 1] = value
        K_μ[2, 2, 1, 2, 2, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 2, 1, 2, 1] = value
        K_μ[2, 2, 2, 2, 2, 1, 2, 2] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) - (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (36 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 2, 1] = value
        K_μ[2, 1, 1, 2, 2, 2, 2, 2] = value
        K_μ[2, 1, 2, 1, 2, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 2, 2, 1, 2] = value
        K_μ[2, 2, 1, 1, 2, 1, 2, 1] = value
        K_μ[2, 2, 1, 2, 2, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 2, 1, 1, 1] = value
        K_μ[2, 2, 2, 2, 2, 1, 1, 2] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) - (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 2, 1] = value
        K_μ[3, 1, 1, 2, 3, 2, 2, 2] = value
        K_μ[3, 1, 2, 1, 3, 2, 1, 1] = value
        K_μ[3, 1, 2, 2, 3, 2, 1, 2] = value
        K_μ[3, 2, 1, 1, 3, 1, 2, 1] = value
        K_μ[3, 2, 1, 2, 3, 1, 2, 2] = value
        K_μ[3, 2, 2, 1, 3, 1, 1, 1] = value
        K_μ[3, 2, 2, 2, 3, 1, 1, 2] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) - (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (36 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 2, 2] = value
        K_μ[1, 1, 1, 2, 1, 2, 2, 1] = value
        K_μ[1, 1, 2, 1, 1, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 1, 2, 1, 1] = value
        K_μ[1, 2, 1, 1, 1, 1, 2, 2] = value
        K_μ[1, 2, 1, 2, 1, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 1, 1, 1, 2] = value
        K_μ[1, 2, 2, 2, 1, 1, 1, 1] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) + (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 1, 2] = value
        K_μ[3, 1, 1, 2, 3, 2, 1, 1] = value
        K_μ[3, 1, 2, 1, 3, 2, 2, 2] = value
        K_μ[3, 1, 2, 2, 3, 2, 2, 1] = value
        K_μ[3, 2, 1, 1, 3, 1, 1, 2] = value
        K_μ[3, 2, 1, 2, 3, 1, 1, 1] = value
        K_μ[3, 2, 2, 1, 3, 1, 2, 2] = value
        K_μ[3, 2, 2, 2, 3, 1, 2, 1] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) + (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 1, 2] = value
        K_μ[2, 1, 1, 2, 2, 2, 1, 1] = value
        K_μ[2, 1, 2, 1, 2, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 2, 2, 2, 1] = value
        K_μ[2, 2, 1, 1, 2, 1, 1, 2] = value
        K_μ[2, 2, 1, 2, 2, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 2, 1, 2, 2] = value
        K_μ[2, 2, 2, 2, 2, 1, 2, 1] = value

        value =
            (-(h[2] * h[3]) / (36 * h[1])) - (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (36 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 2, 2] = value
        K_μ[2, 1, 1, 2, 2, 2, 2, 1] = value
        K_μ[2, 1, 2, 1, 2, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 2, 2, 1, 1] = value
        K_μ[2, 2, 1, 1, 2, 1, 2, 2] = value
        K_μ[2, 2, 1, 2, 2, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 2, 1, 1, 2] = value
        K_μ[2, 2, 2, 2, 2, 1, 1, 1] = value

        value =
            (-(h[2] * h[3]) / (36 * h[1])) - (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 2, 2] = value
        K_μ[3, 1, 1, 2, 3, 2, 2, 1] = value
        K_μ[3, 1, 2, 1, 3, 2, 1, 2] = value
        K_μ[3, 1, 2, 2, 3, 2, 1, 1] = value
        K_μ[3, 2, 1, 1, 3, 1, 2, 2] = value
        K_μ[3, 2, 1, 2, 3, 1, 2, 1] = value
        K_μ[3, 2, 2, 1, 3, 1, 1, 2] = value
        K_μ[3, 2, 2, 2, 3, 1, 1, 1] = value

        value =
            (h[2] * h[3]) / (36 * h[1]) - (h[1] * h[3]) / (9 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 2, 2] = value
        K_μ[2, 1, 1, 2, 2, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 2, 1, 1, 2] = value
        K_μ[2, 1, 2, 2, 2, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 2, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 2, 2, 2, 1] = value
        K_μ[2, 2, 2, 1, 2, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 2, 2, 1, 1] = value

        value =
            (h[2] * h[3]) / (36 * h[1]) - (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 2, 2] = value
        K_μ[3, 1, 1, 2, 3, 1, 2, 1] = value
        K_μ[3, 1, 2, 1, 3, 1, 1, 2] = value
        K_μ[3, 1, 2, 2, 3, 1, 1, 1] = value
        K_μ[3, 2, 1, 1, 3, 2, 2, 2] = value
        K_μ[3, 2, 1, 2, 3, 2, 2, 1] = value
        K_μ[3, 2, 2, 1, 3, 2, 1, 2] = value
        K_μ[3, 2, 2, 2, 3, 2, 1, 1] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) - (2 * h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 2, 1] = value
        K_μ[2, 1, 1, 2, 2, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 2, 1, 1, 1] = value
        K_μ[2, 1, 2, 2, 2, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 2, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 2, 2, 2, 2] = value
        K_μ[2, 2, 2, 1, 2, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 2, 2, 1, 2] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) - (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 2, 1] = value
        K_μ[3, 1, 1, 2, 3, 1, 2, 2] = value
        K_μ[3, 1, 2, 1, 3, 1, 1, 1] = value
        K_μ[3, 1, 2, 2, 3, 1, 1, 2] = value
        K_μ[3, 2, 1, 1, 3, 2, 2, 1] = value
        K_μ[3, 2, 1, 2, 3, 2, 2, 2] = value
        K_μ[3, 2, 2, 1, 3, 2, 1, 1] = value
        K_μ[3, 2, 2, 2, 3, 2, 1, 2] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) - (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 2, 2] = value
        K_μ[1, 1, 1, 2, 1, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 1, 1, 1, 2] = value
        K_μ[1, 1, 2, 2, 1, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 1, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 1, 2, 2, 1] = value
        K_μ[1, 2, 2, 1, 1, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 1, 2, 1, 1] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) + (h[1] * h[3]) / (18 * h[2]) -
            (2 * h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 1, 2] = value
        K_μ[3, 1, 1, 2, 3, 1, 1, 1] = value
        K_μ[3, 1, 2, 1, 3, 1, 2, 2] = value
        K_μ[3, 1, 2, 2, 3, 1, 2, 1] = value
        K_μ[3, 2, 1, 1, 3, 2, 1, 2] = value
        K_μ[3, 2, 1, 2, 3, 2, 1, 1] = value
        K_μ[3, 2, 2, 1, 3, 2, 2, 2] = value
        K_μ[3, 2, 2, 2, 3, 2, 2, 1] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) + (h[1] * h[3]) / (9 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 1, 2] = value
        K_μ[2, 1, 1, 2, 2, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 2, 1, 2, 2] = value
        K_μ[2, 1, 2, 2, 2, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 2, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 2, 2, 1, 1] = value
        K_μ[2, 2, 2, 1, 2, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 2, 2, 2, 1] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) - (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 2, 1] = value
        K_μ[1, 1, 1, 2, 1, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 1, 1, 1, 1] = value
        K_μ[1, 1, 2, 2, 1, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 1, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 1, 2, 2, 2] = value
        K_μ[1, 2, 2, 1, 1, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 1, 2, 1, 2] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) + (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 1, 2] = value
        K_μ[1, 1, 1, 2, 1, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 1, 1, 2, 2] = value
        K_μ[1, 1, 2, 2, 1, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 1, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 1, 2, 1, 1] = value
        K_μ[1, 2, 2, 1, 1, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 1, 2, 2, 1] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) +
            (h[1] * h[3]) / (9 * h[2]) +
            (2 * h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 1, 1] = value
        K_μ[3, 1, 1, 2, 3, 1, 1, 2] = value
        K_μ[3, 1, 2, 1, 3, 1, 2, 1] = value
        K_μ[3, 1, 2, 2, 3, 1, 2, 2] = value
        K_μ[3, 2, 1, 1, 3, 2, 1, 1] = value
        K_μ[3, 2, 1, 2, 3, 2, 1, 2] = value
        K_μ[3, 2, 2, 1, 3, 2, 2, 1] = value
        K_μ[3, 2, 2, 2, 3, 2, 2, 2] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) +
            (2 * h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 1, 1] = value
        K_μ[2, 1, 1, 2, 2, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 2, 1, 2, 1] = value
        K_μ[2, 1, 2, 2, 2, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 2, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 2, 2, 1, 2] = value
        K_μ[2, 2, 2, 1, 2, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 2, 2, 2, 2] = value

        value =
            (2 * h[2] * h[3]) / (9 * h[1]) +
            (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 1, 1] = value
        K_μ[1, 1, 1, 2, 1, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 1, 1, 2, 1] = value
        K_μ[1, 1, 2, 2, 1, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 1, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 1, 2, 1, 2] = value
        K_μ[1, 2, 2, 1, 1, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 1, 2, 2, 2] = value
    else
    end
    return K_μ
end

@testset "Brick module" begin

    @testset "Gauss–Legendre integration" begin
        order = 2
        for d = 1:3
            @testset "Gauss-Legendre integration, $(d)d" begin
                for n in map(Tuple, CartesianIndices(Tuple(fill(1:(2*order-1), d))))
                    h = Tuple([1.1, 1.2, 1.3][1:d])
                    actual = integrate(x -> prod(x .^ n), h)
                    n1 = n .+ 1
                    expected = prod(((h ./ 2) .^ n1 .- (h ./ (-2)) .^ n1) ./ n1)
                    @test actual ≈ expected rtol = 1e-15 atol = 1e-15
                end
            end
        end
    end

    @testset "Node numbering" begin
        for d = 1:3
            h = Tuple([1.1, 1.2, 1.3][1:d])
            ℒ = cell_vertices(d)
            for p ∈ ℒ
                x = h .* (-1) .^ Tuple(p) ./ 2
                N = shape(x, h)
                for q ∈ ℒ
                    @test N[q] == Float64(p == q)
                end
            end
        end
    end

    @testset "Brick element operators" begin
        for d = 2:3
            @testset "Brick element operators, $(d)d" begin
                # Note: reference values where computed with maxima
                # (see `ref/nd_brick_elasticity.mac`).
                h = Tuple([1.1, 1.2, 1.3][1:d])
                C = Hooke{d,Float64}(5.6, 0.3)

                @testset "Average gradient_operator, $(d)d" begin
                    D_act = avg_gradient_operator(h)
                    D_exp = expected_avg_gradient_operator(h)
                    @test D_act ≈ D_exp rtol = 1e-15 atol = 1e-15
                end

                @testset "Average strain-displacement matrix, $(d)d" begin
                    B_act = avg_strain_displacement_operator(h)
                    B_exp = expected_avg_strain_displacement_operator(h)
                    @test B_act ≈ B_exp rtol = 1e-15 atol = 1e-15
                end

                @testset "Stiffness operator, $(d)d" begin
                    K_act = stiffness_operator(h, C)
                    K_λ = expected_K_λ(h)
                    K_μ = expected_K_μ(h)
                    K_exp = C.λ * K_λ + C.μ * K_μ
                    @test K_act ≈ K_exp rtol = 1e-15 atol = 1e-15
                end
            end
        end
    end
end


text:
import Pkg

Pkg.activate(".")
Pkg.add("Statistics")
Pkg.add("YAML")
Pkg.add("JLD2")
Pkg.add(url="https://github.com/felixhorger/DICOM.jl")
Pkg.add(url="https://github.com/felixhorger/DICOMTools.jl")
Pkg.add(url="https://github.com/felixhorger/MRIRelax.jl")

ext: jl


text:
"""
Module Help

Provides functions for displaying help about the commands during speech recognition.

# Functions
- [`help`](@ref)

To see a description of a function type `?<functionname>`.
"""
module Help

import ..JustSayIt: @voiceargs, TYPE_MODEL_NAME, command, command_names, next_token

const COMMANDS_KEYWORD = "commands"

"Show help for your commands or a spoken command or module."
function help()
    valid_input = [COMMANDS_KEYWORD, command_names()...]
    keyword = next_token(valid_input)
    if keyword == COMMANDS_KEYWORD
        cmd_length_max = maximum(length.(command_names()))
        @info join(["", "Your commands:",
                    map(sort([command_names()...])) do x
                       join((x, command(x)), " "^(cmd_length_max+1-length(x)) * "=> ")
                    end...
                    ], "\n")
    elseif keyword in command_names()
        @info "Command $keyword" ""=Base.Docs.doc(command(keyword))
    else
        @info "Keyword not recognized."
    end
    return
end

end # module Help


text:
# Test basic methods
@testset "Basics" begin
    # test _macro_error
    @testset "_macro_error" begin
        source = LineNumberNode(42, :bob)
        @test_throws ErrorException InfiniteOpt._macro_error(:test, (:x,), source, "hello")
    end
    # test _esc_non_constant
    @testset "_esc_non_constant" begin
        @test InfiniteOpt._esc_non_constant(42) == 42
        @test InfiniteOpt._esc_non_constant(:(quote 42 end)) == :(quote 42 end)
        @test InfiniteOpt._esc_non_constant(:(x[i])) == esc(:(x[i]))
        @test InfiniteOpt._esc_non_constant(:x) == esc(:x)
    end
    # test _get_name
    @testset "_get_name" begin
        @test InfiniteOpt._get_name(:a) == :a
        @test InfiniteOpt._get_name(nothing) == ()
        @test InfiniteOpt._get_name("bob") == "bob"
        @test InfiniteOpt._get_name(:("a$(x)")) == :("a$(x)")
        @test InfiniteOpt._get_name(:(x[1:2])) == :x
    end
    # test _name_call
    @testset "_name_call" begin
        @test InfiniteOpt._name_call("bob", []) == "bob"
        @test InfiniteOpt._name_call("", []) == ""
        expected = :(string("x", "[", string($(Expr(:escape, "2"))), ",", string($(Expr(:escape, "1"))), "]"))
        @test InfiniteOpt._name_call("x", ["2", "1"]) == expected
    end
    # test _extract_kwargs
    @testset "_extract_kwargs" begin
        # test normal
        args = (:x, :(d=42), :(container=Cat))
        @test InfiniteOpt._extract_kwargs(args)[1] == [:x]
        @test InfiniteOpt._extract_kwargs(args)[2] == [:(d=42)]
        @test InfiniteOpt._extract_kwargs(args)[3] == :(Cat)
        @test InfiniteOpt._extract_kwargs(args)[4] === nothing
        # test with semi-colon keywords
        args = (Expr(:parameters, :(d=42), :(base_name = "cat")), :x, :(container=Cat))
        @test InfiniteOpt._extract_kwargs(args)[1] == [:x]
        @test InfiniteOpt._extract_kwargs(args)[2] == [:(d=42)]
        @test InfiniteOpt._extract_kwargs(args)[3] == :(Cat)
        @test InfiniteOpt._extract_kwargs(args)[4] == esc("cat")
    end
    # test _add_kw_args
    @testset "_add_kwargs" begin
        esc_kw = esc(Expr(:kw, :d, 2))
        call = :(f(a))
        @test InfiniteOpt._add_kwargs(call, [:(d = 2)]) isa Nothing
        @test call == :(f(a, $esc_kw))
    end
     # test _has_idxvars
     @testset "_has_idxvars" begin
        vars = Any[:i, :j]
        @test !InfiniteOpt._has_idxvars(:w, vars)
        @test InfiniteOpt._has_idxvars(:i, vars)
        @test InfiniteOpt._has_idxvars(:(a[i]), vars)
        @test !InfiniteOpt._has_idxvars(:(ones(a)), vars)
    end
    # test _valid_model
    @testset "_valid_model" begin
        m = InfiniteModel()
        m2 = 2
        @test InfiniteOpt._valid_model(error, m, :m) isa Nothing
        @test_throws ErrorException InfiniteOpt._valid_model(error, m2, :m2)
    end
    # test _error_if_cannot_register
    @testset "_error_if_cannot_register" begin
        m = InfiniteModel()
        @test InfiniteOpt._error_if_cannot_register(error, m, :x) isa Nothing
        object_dictionary(m)[:x] = 2
        @test_throws ErrorException InfiniteOpt._error_if_cannot_register(error, m, :x)
        @test_throws ErrorException InfiniteOpt._error_if_cannot_register(error, m, :((x, 2)))
    end
    # test _macro_assign_and_return
    @testset "_macro_assign_and_return" begin
        m = InfiniteModel()
        esc_m = esc(:m)
        code = :(GeneralVariableRef($esc_m, 1, FiniteVariableIndex))
        @test InfiniteOpt._macro_assign_and_return(error, code, :a, esc_m) isa Expr
    end
    # test _finalize_macro
    @testset "_finalize_macro" begin
        m = InfiniteModel()
        esc_m = esc(:m)
        code = :(GeneralVariableRef($esc_m, 1, FiniteVariableIndex))
        source = LineNumberNode(42, :bob)
        @test InfiniteOpt._finalize_macro(error, esc_m, code, source) isa Expr
    end
end


text:
module EMLPIDIncr
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("PIDIncr/PIDIncr.jl")
end


text:
using BinaryBuilder

name = "glslang"
version = v"11.7.0"

source = "https://github.com/KhronosGroup/glslang.git"
commit = "925503088e2bcd76921b1e102c37fc320bace254" # Nov 11, 2021 (tag 11.7.0)

sources = [
    GitSource(source, commit)
]

script = raw"""
cd glslang

CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Explicitly use our cmake toolchain file
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# The products that we will ensure are always built
products = [
    ExecutableProduct("glslangValidator", :glslangValidator),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("SPIRV_Tools_jll")
]

platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
function DataDictionary(time_start,time_stop,time_step)

	# load stoichiometry
	stoichiometric_matrix = readdlm("./stoichiometry.dat")
	# Setup default flux bounds array
	default_bounds_array = [
		0 100.0; # 1 1.0*m_A_e<uptake:>1.0*m_A_c
		0 1.0; # 2 1.0*m_A_c<catalyze:>1.0*m_B_c
		0 1.0; # 3 1.0*m_B_c<catalyze:>1.0*m_A_c
		0 2.0; # 4 1.0*m_A_c<catalyze:>1.0*m_C_c
		0 100.0; # 5 1.0*m_C_c<catalyze:>1.0*m_B_c
		0 100.0; # 6 1.0*m_B_c<secrete:>1.0*m_B_e
		0 100.0; # 7 1.0*m_C_c<secrete:>1.0*m_C_e
	]

	# Setup default species bounds array
	species_bounds_array = [
		-10.0 10.0; # 1 m_A_e
		-10.0 10.0; # 2 m_B_e
		-10.0 10.0; # 3 m_C_e
		 0.0 0.0; # 4 m_A_c
		 0.0 0.0; # 5 m_B_c
		 0.0 0.0; # 6 m_C_c
	]

	# Setup the objective coefficient array
	objective_coefficient_array = [
		0.0; # 1 1.0*m_A_e<uptake:>1.0*m_A_c
		0.0; # 2 1.0*m_A_c<catalyze:>1.0*m_B_c
		0.0; # 3 1.0*m_B_c<catalyze:>1.0*m_A_c
		0.0; # 4 1.0*m_A_c<catalyze:>1.0*m_C_c
		0.0; # 5 1.0*m_C_c<catalyze:>1.0*m_B_c
		0.0; # 6 1.0*m_B_c<secrete:>1.0*m_B_e
		0.0; # 7 1.0*m_C_c<secrete:>1.0*m_C_e
	]

	# Min/Max flag - default is minimum -
	is_minimum_flag = true

	# List of reation strings - used to write flux report
	list_of_reaction_strings = [
		"1.0*m_A_e<uptake:>1.0*m_A_c"
		"1.0*m_A_c<catalyze:>1.0*m_B_c"
		"1.0*m_B_c<catalyze:>1.0*m_A_c"
		"1.0*m_A_c<catalyze:>1.0*m_C_c"
		"1.0*m_C_c<catalyze:>1.0*m_B_c"
		"1.0*m_B_c<secrete:>1.0*m_B_e"
		"1.0*m_C_c<secrete:>1.0*m_C_e"
	]

	# List of metabolite strings - used to write flux report
	list_of_metabolite_symbols = [
		"m_A_e"
		"m_B_e"
		"m_C_e"
		"m_A_c"
		"m_B_c"
		"m_C_c"
	]

	# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{AbstractString,Any}()
	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
	data_dictionary["objective_coefficient_array"] = objective_coefficient_array
	data_dictionary["default_flux_bounds_array"] = default_bounds_array;
	data_dictionary["species_bounds_array"] = species_bounds_array
	data_dictionary["list_of_reaction_strings"] = list_of_reaction_strings
	data_dictionary["list_of_metabolite_symbols"] = list_of_metabolite_symbols
	data_dictionary["is_minimum_flag"] = is_minimum_flag
	# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
	return data_dictionary
end

# maximize_product_dictionary
function maximize_product_dictionary(time_start,time_stop,time_step)
	# load the data dictionary -
	data_dictionary = DataDictionary(time_start,time_stop,time_step)
	# Modify the data dictionary -
	objective_coefficient_array = data_dictionary["objective_coefficient_array"]
	objective_coefficient_array[6] = -1
	objective_coefficient_array[7] = -1
	# return -
	return data_dictionary
end


text:
"""
    choose_last_parameter!(taller, arity, n_way)

Given a test set where the first k-1 parameters are chosen and the last parameter
has not been chosen, this fills in the last parameters for each test case.
"""
function choose_last_parameter!(taller, allc, matcher = case_partial_cover)
    param_idx  = size(taller, 1)
    for set_col_idx in 1:size(taller, 2)
        if any(taller[:, set_col_idx] .== 0)
            match_hist = matches_from_missing(allc, taller[:, set_col_idx], param_idx, matcher)
            if (any(match_hist .> 0))
                # The argmax tie-breaks in a consistent manner.
                taller[param_idx, set_col_idx] = argmax(match_hist)
                add_coverage!(allc, taller[:, set_col_idx])
            end  # else don't set this entry by leaving it zero.
        end
    end
end

function choose_last_parameter_filter!(taller, allc, disallow)
    param_idx  = size(taller, 1)
    putative = similar(taller[:, 1])
    for set_col_idx in 1:size(taller, 2)
        if any(taller[:, set_col_idx] .== 0)
            match_hist = matches_from_missing(
                allc, taller[:, set_col_idx], param_idx, case_partial_cover)
            for max_idx in (midx for midx in sortperm(match_hist, rev = true)
                    if match_hist[midx] > 0)
                # The argmax tie-breaks in a consistent manner.
                putative .= taller[:, set_col_idx]
                putative[param_idx] = max_idx
                if !disallow(putative)
                    taller[param_idx, set_col_idx] = max_idx
                    add_coverage!(allc, taller[:, set_col_idx])
                    break
                end
            end  # else don't set this entry by leaving it zero.
        end
    end
end

"""
The given test cases have missing values, which are zeros.
This fills in missing values anywhere they cover tuples.
"""
function fill_missing_test_set_values!(taller, allc, matcher = case_compatible_with_tuple)
    param_idx  = size(taller, 1)
    for missing_col_idx in 1:size(taller, 2)
        nonzero = sum(taller[1:(param_idx - 1), missing_col_idx] .> 0)
        if nonzero < param_idx - 1
            # The found_values has what format?
            found_entry = fill_consistent_matches(allc, taller[:, missing_col_idx], matcher)
            remain_zero = sum(found_entry .> 0)
            if remain_zero < nonzero
                taller[:, missing_col_idx] .= found_entry
                add_coverage!(allc, found_entry)
            end  # else nothing found for this row.
        end
    end
end

"""
If not all tuples are covered by a test set, this creates new test cases
to cover all remaining tuples.
"""
function cover_remaining_by_creating_cases(allc, matcher = case_compatible_with_tuple)
    param_idx = parameter_cnt(allc)
    add_entries = Array{Array{eltype(allc),1},1}()
    while remaining_uncovered(allc) > 0
        # add a new row. Fill with necessary tuples.
        entry = first_match_for_parameter(allc, param_idx)
        filled = fill_consistent_matches(allc, entry, matcher)
        add_coverage!(allc, filled)
        push!(add_entries, filled)
    end
    add_entries
end

function put_tuple_in_case(tuple, case)
    for i in 1:length(tuple)
        if tuple[i] != 0
            if case[i] == 0
                case[i] = tuple[i]
            elseif case[i] != tuple[i]
                error("tuple doesn't match case $(tuple), $(case)")
            end  # else the values match.
        end
    end
    case
end

"""
Loop over remaining tuples instead of looping over each test.
For each tuple, loop over tests where that tuple could go. If that fails,
add the tuple at the end as its own test.
"""
function insert_tuple_into_tests(test_set, allc, matcher = case_compatible_with_tuple)
    add_tests = Array{eltype(allc), 1}[]
    for find_cover_idx in allc.remain:-1:1
        tuple = allc.allc[:, find_cover_idx]
        unmatched = true
        for test_idx in 1:size(test_set, 2)
            test_case = test_set[:, test_idx]
            matches = matcher(test_case, tuple)
            if matches
                test_set[:, test_idx] = put_tuple_in_case(tuple, test_case)
                unmatched = false
                break
            end
        end
        if unmatched
            for tc_idx in 1:length(add_tests)
                test_case = add_tests[tc_idx]
                if matcher(test_case, tuple)
                    add_tests[tc_idx] = put_tuple_in_case(tuple, test_case)
                    unmatched = false
                    break
                end
            end
        end
        if unmatched
            push!(add_tests, tuple)
        end
    end
    allc.remain = 0
    hcat(test_set, add_tests...)
end

function putative_allowed(buffer, base, replace, disallow)
    buffer .= base
    buffer[replace .!= 0] .= replace[replace .!= 0]
    !disallow(buffer)
end

function insert_tuple_into_tests_filter(test_set, allc, disallow)
    add_tests = Array{eltype(allc), 1}[]
    putative = similar(test_set[:, 1])
    for find_cover_idx in allc.remain:-1:1
        tuple = allc.allc[:, find_cover_idx]
        unmatched = true
        for test_idx in 1:size(test_set, 2)
            test_case = test_set[:, test_idx]
            matches = case_compatible_with_tuple(test_case, tuple)
            if matches && putative_allowed(putative, test_case, tuple, disallow)
                test_set[:, test_idx] = put_tuple_in_case(tuple, test_case)
                unmatched = false
                break
            end
        end
        if unmatched
            for tc_idx in 1:length(add_tests)
                test_case = add_tests[tc_idx]
                if case_compatible_with_tuple(test_case, tuple) &&
                        putative_allowed(putative, test_case, tuple, disallow)
                    add_tests[tc_idx] = put_tuple_in_case(tuple, test_case)
                    unmatched = false
                    break
                end
            end
        end
        if unmatched
            push!(add_tests, tuple)
        end
    end
    allc.remain = 0
    hcat(test_set, add_tests...)
end

"""
As a finishing step on a set of test cases, this fills in missing values
which aren't needed to cover the tuples but can increase coverage
of higher tuples.
"""
function fill_remaining_missing_values!(test_set, arity)
    param_cnt = size(test_set, 1)
    # We could have zero values at the end, so fill them in with the
    # least-used values.
    hist = zeros(Int, param_cnt, maximum(arity))
    for hist_entry in 1:size(test_set, 2)
        for hist_param in 1:size(test_set, 1)
            if test_set[hist_param, hist_entry] > 0
                hist[hist_param, test_set[hist_param, hist_entry]] += 1
            end
        end
    end
    for fill_col in 1:size(test_set, 2)
        for fill_param in 1:size(test_set, 1)
            if test_set[fill_param, fill_col] == 0
                fill_val = argmin(hist[fill_param, 1:arity[fill_param]])
                test_set[fill_param, fill_col] = fill_val
                hist[fill_param, fill_val] += 1
            end
        end
    end
end

function fill_remaining_missing_values_filter!(test_set, arity, disallow)
    param_cnt = size(test_set, 1)
    # We could have zero values at the end, so fill them in with the
    # least-used values.
    hist = zeros(Int, param_cnt, maximum(arity))
    for hist_entry in 1:size(test_set, 2)
        for hist_param in 1:size(test_set, 1)
            if test_set[hist_param, hist_entry] > 0
                hist[hist_param, test_set[hist_param, hist_entry]] += 1
            end
        end
    end
    putative = similar(test_set[:, 1])
    for fill_col in 1:size(test_set, 2)
        for fill_param in 1:size(test_set, 1)
            if test_set[fill_param, fill_col] == 0
                putative .= test_set[:, fill_col]
                for fill_val in sortperm(hist[fill_param, 1:arity[fill_param]])
                    putative[fill_param] = fill_val
                    if !disallow(putative)
                        test_set[fill_param, fill_col] = fill_val
                        hist[fill_param, fill_val] += 1
                        break
                    end
                end
            end
        end
    end
end

"""
    ipog(arity, n_way)

The `arity` is an integer array of the number of possible values each
parameter can take. The `n_way` is whether each parameter must appear
once in the test suite, or whether each pair of parameters must appear
together, or whether each triple must appear together. The wayness
can be set as high as the length of the arity.

This uses the in-parameter-order-general algorithm. It will always return
the same set of values. It takes less time and memory than most other
approaches.

This function represents a test set as a two-dimensional array of
the same integer type as the input arity. Each value of the array
is either an integer number, from 1 to the arity of that parameter,
or it is 0 for what the paper calls a don't-care value.

Lei, Yu, Raghu Kacker, D. Richard Kuhn, Vadim Okun, and James Lawrence.
2008. “IPOG/IPOG-D: Efficient Test Generation for Multi-Way Combinatorial
Testing.” Software Testing, Verification & Reliability 18 (3): 125–48.
"""
function ipog(arity, n_way)
    nonincreasing = sortperm(arity, rev = true)
    original_arity = arity
    arity = arity[nonincreasing]
    original_order = sortperm(nonincreasing)

    param_cnt = length(arity)
    # Setup by taking first n_way parameters.
    # This is a 2D array.
    test_set = all_combinations(arity[1:n_way], n_way)

    for param_idx in (n_way + 1):param_cnt
        taller = zeros(eltype(arity), param_idx, size(test_set, 2))
        taller[1:(param_idx - 1), :] .= test_set

        # Seed test cases by adding them once params are covered and not double-covering.
        # Make mixed strength here, once all params at a strength are covered.
        allc = one_parameter_combinations_matrix(arity[1:param_idx], n_way)

        choose_last_parameter!(taller, allc)

        test_set = insert_tuple_into_tests(taller, allc)
    end

    fill_remaining_missing_values!(test_set, arity)

    # reorder test columns with `original_order`.
    test_set[original_order, :]
end

function keep_allowed(test_set, disallow)
    allow_cnt = 0
    allowed = zeros(Int, size(test_set, 2))
    for idx in 1:size(test_set, 2)
        if !disallow(test_set[:, idx])
            allow_cnt += 1
            allowed[allow_cnt] = idx
        end  # else not allowed so not put into the list.
    end
    test_set[:, allowed[1:allow_cnt]]
end

function reorder_disallow(disallow, original_order)
    let oo = original_order, dis = disallow
        choices -> begin
            v = zeros(eltype(oo), size(oo)...)
            v[1:length(choices)] .= choices
            dis(v[oo])
        end
    end
end

function ipog_multi(arity, n_way, disallow, seed)
    nonincreasing = sortperm(arity, rev = true)
    original_arity = arity
    arity = arity[nonincreasing]
    original_order = sortperm(nonincreasing)
    forbid = reorder_disallow(disallow, original_order)

    param_cnt = length(arity)
    if seed !== nothing
        seed = seed[nonincreasing, :]
        seed_tests = seed[1:n_way, :]
    else
        seed = zeros(Int, n_way, 0)
        seed_tests = zeros(Int, n_way, 0)
    end
    # Setup by taking first n_way parameters.
    # This is a 2D array.
    combo_tests = keep_allowed(all_combinations(arity[1:n_way], n_way), forbid)
    test_set = add_tests_to_seeds(seed[1:n_way, :], combo_tests)

    for param_idx in (n_way + 1):param_cnt
        taller = zeros(eltype(arity), param_idx, size(test_set, 2))
        taller[1:(param_idx - 1), :] .= test_set
        # Seed test cases by adding them once params are covered and not double-covering.
        # Make mixed strength here, once all params at a strength are covered.
        allc = one_parameter_combinations_matrix(arity[1:param_idx], n_way)
        remove_combinations!(allc, forbid)

        if size(seed, 2) > 0
            taller[param_idx, 1:size(seed, 2)] = seed[param_idx, :]
            for seed_cover_idx in 1:size(seed, 2)
                add_coverage!(allc, taller[:, seed_cover_idx])
            end
        end
        choose_last_parameter_filter!(taller, allc, forbid)

        test_set = insert_tuple_into_tests_filter(taller, allc, forbid)
    end

    fill_remaining_missing_values_filter!(test_set, arity, forbid)

    # reorder test columns with `original_order`.
    test_set[original_order, :]
end

"""
The seeds may not be unique, so we can't combine tests and seeds
with the unique function. The problem is that seeds may not be
unique for the few parameters in the combination test.
"""
function add_tests_to_seeds(seeds, tests)
    append_cnt = 0
    to_append = zeros(Int, size(tests, 2))
    for tidx in 1:size(tests, 2)
        found = false
        test = tests[:, tidx]
        for sidx in 1:size(seeds, 2)
            if case_compatible_with_tuple(test, seeds[:, sidx])
                for pidx in 1:length(test)
                    if test[pidx] != 0
                        seeds[pidx, sidx] = test[pidx]
                    end
                end
                found = true
                break
            end
        end
        if !found
            append_cnt += 1
            to_append[append_cnt] = tidx
        end
    end
    if append_cnt > 0
        extracted = tests[:, to_append[1:append_cnt]]
        hcat(seeds, extracted)
    else
        seeds
    end
end

function ipog_inner(arity, n_way, forbid, seed)
    param_cnt = length(arity)
    # Setup by taking first n_way parameters.
    # This is a 2D array.
    combo_tests = keep_allowed(all_combinations(arity[1:n_way], n_way), forbid)
    test_set = add_tests_to_seeds(seed[1:n_way, :], combo_tests)

    for param_idx in (n_way + 1):param_cnt
        taller = zeros(eltype(arity), param_idx, size(test_set, 2))
        taller[1:(param_idx - 1), :] .= test_set
        # Seed test cases by adding them once params are covered and not double-covering.
        # Make mixed strength here, once all params at a strength are covered.
        allc = one_parameter_combinations_matrix(arity[1:param_idx], n_way)
        remove_combinations!(allc, forbid)

        if size(seed, 2) > 0
            taller[param_idx, 1:size(seed, 2)] = seed[param_idx, :]
            for seed_cover_idx in 1:size(seed, 2)
                add_coverage!(allc, taller[:, seed_cover_idx])
            end
        end
        choose_last_parameter_filter!(taller, allc, forbid)

        test_set = insert_tuple_into_tests_filter(taller, allc, forbid)
    end
    test_set
end

struct WayWork
    indices
    arity
    n_way
    combo_cnt
end

function ipog_multi_way(arity, n_way, levels, disallow, seed)
    param_cnt = length(arity)
    levels[n_way] = [collect(1:param_cnt)]
    waynesses = sort(collect(keys(levels)), rev = true)
    work = WayWork[]
    for wayness in waynesses
        for indices_idx = 1:length(levels[wayness])
            indices = levels[wayness][indices_idx]
            an_arity = arity[indices]
            # We sort arity so that the in-parameter-order tackles harder parms first.
            arity_order = sortperm(an_arity, rev = true)
            indices_sorted = indices[arity_order]
            arity_sorted = an_arity[arity_order]
            combo_cnt = total_combinations(arity_sorted, wayness)
            push!(work, WayWork(indices_sorted, arity_sorted, wayness, combo_cnt))
        end
    end
    # I guess this will work better if we put the stiffer wayness first.
    sort!(work; by = (x -> (x.n_way, x.combo_cnt)), rev = true)

    if seed === nothing
        test_cases = zeros(Int, param_cnt, 0)
    else
        test_cases = seed
    end
    for round in work
        round_seed = test_cases[round.indices, :]
        whole_order = vcat(round.indices, [ind for ind in 1:param_cnt if ind ∉ round.indices])
        forbid = reorder_disallow(disallow, sortperm(whole_order))
        rotated_cases = test_cases[round.indices, :]
        tests_rotated = ipog_inner(round.arity, round.n_way, forbid, rotated_cases)
        expand_cases = zeros(Int, param_cnt, size(tests_rotated, 2))
        expand_cases[:, 1:size(test_cases, 2)] .= test_cases
        # reverse ordering?
        expand_cases[round.indices, :] .= tests_rotated
        test_cases = expand_cases
    end

    fill_remaining_missing_values_filter!(test_cases, arity, disallow)
    test_cases
end

function ipog_bytuple_instrumented(arity, n_way, strategy)
    nonincreasing = sortperm(arity, rev = true)
    original_arity = arity
    arity = arity[nonincreasing]
    original_order = sortperm(nonincreasing)

    param_cnt = length(arity)
    widen = zeros(Int, param_cnt)
    fillz = zeros(Int, param_cnt)
    cover = zeros(Int, param_cnt)
    # Setup by taking first n_way parameters.
    # This is a 2D array.
    test_set = all_combinations(arity[1:n_way], n_way)

    for param_idx in (n_way + 1):param_cnt
        taller = zeros(eltype(arity), param_idx, size(test_set, 2))
        taller[1:(param_idx - 1), :] .= test_set

        allc = one_parameter_combinations_matrix(arity[1:param_idx], n_way)

        choose_last_parameter!(taller, allc, strategy[:lastparam])

        widen[param_idx] = size(allc.allc, 2) - allc.remain

        test_set = insert_tuple_into_tests(taller, allc, strategy[:expand])
    end

    # remove the part to fill missing values at the end.

    # reorder test columns with `original_order`.
    (test_set[original_order, :], widen, fillz, cover)
end


text:
using BSON
using RCall

R"load('Ricker_sigle_lognormal.rda')"
R"ls()"

@rget y X y_valid X_valid Y_valid tprime

bson("Ricker_sigle_lognormal.bson", y = deepcopy(y),
    X = deepcopy(X), y_valid = deepcopy(y_valid),
    X_valid = deepcopy(X_valid), Y_valid = deepcopy(Y_valid),
    tprime = deepcopy(tprime))

R"rm(list=ls())"

R"load('Ricker_twolag_lognormal.rda')"
R"ls()"

@rget y X y_valid X_valid Y_valid tprime

bson("Ricker_twolag_lognormal.bson", y = deepcopy(y),
    X = deepcopy(X), y_valid = deepcopy(y_valid),
    X_valid = deepcopy(X_valid), Y_valid = deepcopy(Y_valid),
    tprime = deepcopy(tprime))

R"rm(list=ls())"

R"load('Ricker_single.rda')"
R"ls()"

@rget y X y_valid X_valid Y_valid tprime

bson("Ricker_single.bson", y = deepcopy(y),
    X = deepcopy(X), y_valid = deepcopy(y_valid),
    X_valid = deepcopy(X_valid), Y_valid = deepcopy(Y_valid),
    tprime = deepcopy(tprime))

### plotting

# using Plotly
# nplot = 500

# trace0 = scatter3d(Dict(
#     :x => X[1:nplot,1],
#     :y => X[1:nplot,2],
#     :z => y[1:nplot],
#     :marker => Dict(
#         :color => "rgba(0,0,0,0.70)",
#         :size => 2
#     ),
#     :mode => "markers"
# ))
# data = [trace0]
# maximum(X[1:nplot,1:2])

# layout = Layout(Dict(
#   :hovermode => "closest",
#   :margin => Dict(
#     :r => 10,
#     :t => 25,
#     :b => 40,
#     :l => 60
#   ),
#   :scene => Dict(
#     :xaxis => Dict(:title => "y[t-1]",
#                    :range => [0.0, 12.0]),
#     :yaxis => Dict(:title => "y[t-2]",
#                     :range => [0.0, 12.0]),
#     :zaxis => Dict(:title => "y[t]",
#                     :range => [0.0, 12.0]),
#   :showlegend => false
# )))

# p = plot(data, layout)

# nplot = 250
# Rplot = R"plot"
#
# R"pdf('appdata2_plots/appdata2_lag2only_hsced.pdf', width=4.5, height=3.5);
#   par(mar=c(4,4,1,1)+0.1)"
# Rplot(X[1:nplot,2], y[1:nplot], ylab="y[t]", xlab="y[t-2]", main="")
# R"dev.off()"
#
# R"pdf('appdata2_plots/appdata2_lag2only.pdf', width=4.5, height=3.5);
#   par(mar=c(4,4,1,1)+0.1)"
# Rplot(X[1:nplot,2], y[1:nplot], ylab="y[t]", xlab="y[t-2]", main="")
# R"dev.off()"


text:
using Flux, Random

struct ImmutableChain{T<:Tuple}
    layers::T
end

ImmutableChain(x, y, xs...) = ImmutableChain((x, y, xs...))

Flux.children(c::ImmutableChain) = c.layers
Flux.mapchildren(f, c::ImmutableChain) = ImmutableChain(f.(c.layers))

applyall(::Tuple{}, x) = x
applyall(fs::Tuple, x) = applyall(Base.tail(fs), first(fs)(x))

(c::ImmutableChain)(x) = applyall(c.layers, x)

# Helper function to take a recursive data structure (Such as an ImmutableChain)
# and flatten it out into a tuple
@generated function flatten_tuple(t)
    stmts = Any[]
    tuple_args = Any[]
    function collect_args(t, sym)
        for i = 1:fieldcount(t)
            fT = fieldtype(t, i)
            if fT === Nothing
                continue
            end
            g = gensym()
            push!(stmts, :($g = getfield($sym, $i)))
            if !(fT <: AbstractArray)
                collect_args(fT, g)
            else
                push!(tuple_args, g)
            end
        end
    end
    collect_args(t, :t)
    return quote
        $(stmts...)
        tuple($(tuple_args...))
    end
end

# The inverse of flatten_tuple()
@generated function unflatten_tuple(pattern, t)
    function expand_arg(pattern, idx)
        fields = Any[]
        for i = 1:fieldcount(pattern)
            fT = fieldtype(pattern, i)
            if fT === Nothing
                push!(fields, nothing)
                continue
            end
            if fT <: Tuple || fT <: NamedTuple
                arg, idx = expand_arg(fT, idx)
                push!(fields, arg)
            else
                push!(fields, Expr(:call, :getfield, :t, idx))
                idx += 1
            end
        end
        Expr(:new, pattern, fields...), idx
    end
    res = expand_arg(pattern, 1)[1]
    return res
end

# Include our TPUBatchNorm object
include("tpu_batchnorm.jl")
export map_to_tpu, map_to_cpu, TPUBatchNorm, tpu
function map_to_tpu end
const tpu = map_to_tpu

## Mapping utilities
# Convert scalars to single-element XRTArrays with eltype Float32:
map_to_tpu(x::Real) = XRTArray(convert(Float32, x))

# Convert arrays to XRTArrays with eltype Float32
map_to_tpu(x::AbstractArray) = XRTArray(Float32.(x))

# Strip off the TrackedArray coating to get at the data underneath
map_to_tpu(x::TrackedArray) = map_to_tpu(Flux.data(x))

# Turn Chain objects into ImmutableChain objects which store the computation within their type signature
map_to_tpu(x::Chain) = ImmutableChain(tuple(map(map_to_tpu, x.layers)...))

# Convert BatchNorm layers into TPUBatchNorm layers, passing all children straight through,
# except for the "active" child, which is not used by the TPUBatchNorm
map_to_tpu(x::Flux.BatchNorm) = TPUBatchNorm(map(map_to_tpu, Flux.children(x))[1:end-1]...)

# We deal manually with recurrence, so Recur{T} values get turned into just T values
map_to_tpu(x::Flux.Recur) = Flux.mapchildren(map_to_tpu, x.cell)

# For all other objects, just map the children through `map_to_tpu`.
map_to_tpu(x) = Flux.mapchildren(map_to_tpu, x)

map_to_cpu(x) = Flux.mapchildren(map_to_cpu, x)
map_to_cpu(x::XRTArray) = convert(Array, x)
map_to_cpu(x::ImmutableChain) = Chain(map(map_to_cpu, x.layers)...)
map_to_cpu(x::TPUBatchNorm) = Flux.BatchNorm(map(map_to_cpu, Flux.children(x))..., false)
map_to_cpu(x::Flux.LSTMCell) = Flux.Recur(Flux.mapchildren(map_to_cpu, x))

## Layer reimplementations that need to know about XRTArray
import NNlib: softmax, logsoftmax, ∇softmax, ∇logsoftmax
import Flux: logitcrossentropy
export crossentropy, logitcrossentropy, softmax, ∇softmax, logsoftmax, ∇logsoftmax

# Work around inference limitations to dynamically find the number of
# batches held within a tensor.
neg_nbatch(x) = XRTArray(Float32(-size(x, ndims(x))))
Zygote.@nograd neg_nbatch

function crossentropy(ŷ::XRTArray, y::XRTArray)
    return sum(y .* log.(ŷ)) / neg_nbatch(y)
end

function logitcrossentropy(logŷ::XRTArray, y::XRTArray)
    return sum(y .* logsoftmax(logŷ)) / neg_nbatch(y)
end

# Define softmax and its gradient for XRTArrays
function softmax(xs::XRTArray)
    ys = xs .- maximum(xs, dims=1)
    zs = exp.(ys)
    return zs ./ sum(zs, dims=1)
end
function ∇softmax(Δ, xs::XRTArray)
    sf = softmax(xs)
    return sf .* (Δ .- sum(Δ .* sf, dims = 1))
end
Zygote.@adjoint softmax(xs::XRTArray) = softmax(xs), Δ -> (∇softmax(Δ, xs),)

# Define logsoftmax and its gradient for XRTArrays
function logsoftmax(xs::XRTArray)
    ys = xs .- maximum(xs, dims=1)
    zs = sum(exp.(ys), dims=1)
    return ys .- log.(zs)
end
∇logsoftmax(Δ, xs::XRTArray) = ∇softmax(Δ ./ max.(eps(eltype(xs)),softmax(xs)), xs)
Zygote.@adjoint logsoftmax(xs::XRTArray) = logsoftmax(xs), Δ -> (∇logsoftmax(Δ, xs),)

function Base.sort(x::XRTArray{T,<:Any,1}; keys=x) where {T}
    return XLA.HloSort{typeof(<)}(0, false)(<, keys, x)[2]
end

Base.rand(::Type{XRTArray{T}}, shape::Int...) where {T} = XLA.HloRng(T, shape, 1)(XRTArray(typemin(T)),XRTArray(typemax(T)))
Base.rand(::Type{XRTArray{F}}, shape::Int...) where {F <: AbstractFloat} = XLA.HloRng(F, shape, 1)(XRTArray(0.0f0),XRTArray(1.0f0))
Base.randn(::Type{XRTArray{F}}, shape::Int...) where {F <: AbstractFloat} = XLA.HloRng(F, shape, 2)(XRTArray(0.0f0),XRTArray(1.0f0))

# Shuffle by repeatedly sorting x with a random vector as the sorting key.
# see https://github.com/tensorflow/tensorflow/blob/85deaa11cae878ba2c0e5284085956f75434b5b2/tensorflow/compiler/tf2xla/kernels/random_ops.cc#L83-L143 for more
Base.@pure calc_rounds(L; exponent=3) = ceil(UInt32, exponent * log(L) / log(typemax(UInt32)))
function shuffle(x::XRTArray)
    round_idx = XRTArray(calc_rounds(length(x)))
    while round_idx > XRTArray(0x00000000)
        keys = rand(XRTArray{UInt32}, length(x))
        x = sort(x, keys=keys)
        round_idx -= XRTArray(0x00000001)
    end
    return x
end