text:
#=
    Variable bounds
        SingleVariable -in- LessThan
        SingleVariable -in- GreaterThan
        SingleVariable -in- EqualTo
        SingleVariable -in- Interval

        SingleVariable -in- ZeroOne
        SingleVariable -in- Integer
        SingleVariable -in- Semiinteger
        SingleVariable -in- Semicontinuous
=#
constrdict(m::LinQuadOptimizer, ::SVCI{LE}) = cmap(m).upper_bound
constrdict(m::LinQuadOptimizer, ::SVCI{GE}) = cmap(m).lower_bound
constrdict(m::LinQuadOptimizer, ::SVCI{EQ}) = cmap(m).fixed_bound
constrdict(m::LinQuadOptimizer, ::SVCI{IV}) = cmap(m).interval_bound

constrdict(m::LinQuadOptimizer, ::SVCI{MOI.ZeroOne}) = cmap(m).binary
constrdict(m::LinQuadOptimizer, ::SVCI{MOI.Integer}) = cmap(m).integer

constrdict(m::LinQuadOptimizer, ::SVCI{MOI.Semicontinuous{Float64}}) = cmap(m).semicontinuous
constrdict(m::LinQuadOptimizer, ::SVCI{MOI.Semiinteger{Float64}}) = cmap(m).semiinteger

function setvariablebound!(m::LinQuadOptimizer, col::Int, bound::Float64, sense::Cchar)
    change_variable_bounds!(m, [col], [bound], [sense])
end

function setvariablebound!(m::LinQuadOptimizer, v::SinVar, set::LE)
    setvariablebound!(m, getcol(m, v), set.upper, backend_type(m, Val{:Upperbound}()))
end
function setvariablebound!(m::LinQuadOptimizer, v::SinVar, set::GE)
    setvariablebound!(m, getcol(m, v), set.lower, backend_type(m, Val{:Lowerbound}()))
end
function setvariablebound!(m::LinQuadOptimizer, v::SinVar, set::EQ)
    setvariablebound!(m, getcol(m, v), set.value, backend_type(m, Val{:Upperbound}()))
    setvariablebound!(m, getcol(m, v), set.value, backend_type(m, Val{:Lowerbound}()))
end
function setvariablebound!(m::LinQuadOptimizer, v::SinVar, set::IV)
    setvariablebound!(m, getcol(m, v), set.upper, backend_type(m, Val{:Upperbound}()))
    setvariablebound!(m, getcol(m, v), set.lower, backend_type(m, Val{:Lowerbound}()))
end

# add constraint
function MOI.addconstraint!(m::LinQuadOptimizer, v::SinVar, set::S) where S <: LinSets
    setvariablebound!(m, v, set)
    m.last_constraint_reference += 1
    ref = SVCI{S}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    dict[ref] = v.variable
    ref
end

# delete constraint
MOI.candelete(m::LinQuadOptimizer, c::SVCI{S}) where S <: LinSets = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::SVCI{S}) where S <: LinSets
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    vref = dict[c]
    setvariablebound!(m, SinVar(vref), IV(-Inf, Inf))
    delete!(dict, c)
end

# constraint set
MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{S}}) where S <: LinSets = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{LE})
    MOI.LessThan{Float64}(
        get_variable_upperbound(m, getcol(m, m[c]))
    )
end
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{GE})
    MOI.GreaterThan{Float64}(
        get_variable_lowerbound(m, getcol(m, m[c]))
    )
end
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{EQ})
    MOI.EqualTo{Float64}(
        get_variable_lowerbound(m, getcol(m, m[c]))
    )
end
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{IV})
    lb = get_variable_lowerbound(m, getcol(m, m[c]))
    ub = get_variable_upperbound(m, getcol(m, m[c]))
    return MOI.Interval{Float64}(lb, ub)
end

# constraint function
MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{S}}) where S <: LinSets = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::SVCI{<: LinSets})
    return SinVar(m[c])
end

# modify
MOI.canmodifyconstraint(::LinQuadOptimizer, ::SVCI{S}, ::Type{S}) where S <: LinSets = true
function MOI.modifyconstraint!(m::LinQuadOptimizer, c::SVCI{S}, newset::S) where S <: LinSets
    setvariablebound!(m, SinVar(m[c]), newset)
end

#=
    Binary constraints

For some reason CPLEX doesn't respect bounds on a binary variable, so we
should store the previous bounds so that if we delete the binary constraint
we can revert to the old bounds

Xpress is worse, once binary, the bounds are changed independently of what the user does
=#
function MOI.addconstraint!(m::LinQuadOptimizer, v::SinVar, set::MOI.ZeroOne)
    m.last_constraint_reference += 1
    ref = SVCI{MOI.ZeroOne}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    ub = get_variable_upperbound(m, getcol(m, v))
    lb = get_variable_lowerbound(m, getcol(m, v))
    dict[ref] = (v.variable, lb, ub)
    change_variable_types!(m, [getcol(m, v)], [backend_type(m, set)])
    setvariablebound!(m, getcol(m, v), 1.0, backend_type(m, Val{:Upperbound}()))
    setvariablebound!(m, getcol(m, v), 0.0, backend_type(m, Val{:Lowerbound}()))
    make_problem_type_integer(m)
    ref
end

MOI.candelete(m::LinQuadOptimizer, c::SVCI{MOI.ZeroOne}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::SVCI{MOI.ZeroOne})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    (v, lb, ub) = dict[c]
    change_variable_types!(m, [getcol(m, v)], [backend_type(m, Val{:Continuous}())])
    setvariablebound!(m, getcol(m, v), ub, backend_type(m, Val{:Upperbound}()))
    setvariablebound!(m, getcol(m, v), lb, backend_type(m, Val{:Lowerbound}()))
    delete!(dict, c)
    if !hasinteger(m)
        make_problem_type_continuous(m)
    end
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{MOI.ZeroOne}}) = true
MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{MOI.ZeroOne}) = MOI.ZeroOne()

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{MOI.ZeroOne}}) = true
MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::SVCI{MOI.ZeroOne}) = SinVar(m[c][1])

#=
    Integer constraints
=#

function MOI.addconstraint!(m::LinQuadOptimizer, v::SinVar, set::MOI.Integer)
    change_variable_types!(m, [getcol(m, v)], [backend_type(m, set)])
    m.last_constraint_reference += 1
    ref = SVCI{MOI.Integer}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    dict[ref] = v.variable
    make_problem_type_integer(m)
    ref
end

MOI.candelete(m::LinQuadOptimizer, c::SVCI{MOI.Integer}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::SVCI{MOI.Integer})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    v = dict[c]
    change_variable_types!(m, [getcol(m, v)], [backend_type(m, Val{:Continuous}())])
    delete!(dict, c)
    if !hasinteger(m)
        make_problem_type_continuous(m)
    end
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{MOI.Integer}}) = true
MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{MOI.Integer}) = MOI.Integer()

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{MOI.Integer}}) = true
MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::SVCI{MOI.Integer}) = SinVar(m[c])

#=
    Semicontinuous / Semiinteger constraints
=#
const SEMI_TYPES = Union{MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64}}
function MOI.addconstraint!(m::LinQuadOptimizer, v::SinVar, set::SEMI_TYPES)
    change_variable_types!(m, [getcol(m, v)], [backend_type(m, set)])
    setvariablebound!(m, getcol(m, v), set.upper, backend_type(m, Val{:Upperbound}()))
    setvariablebound!(m, getcol(m, v), set.lower, backend_type(m, Val{:Lowerbound}()))

    m.last_constraint_reference += 1
    ref = SVCI{typeof(set)}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    dict[ref] = v.variable
    make_problem_type_integer(m)
    ref
end

MOI.candelete(m::LinQuadOptimizer, c::SVCI{<:SEMI_TYPES}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::SVCI{<:SEMI_TYPES})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    v = dict[c]
    change_variable_types!(m, [getcol(m, v)], [backend_type(m, Val{:Continuous}())])
    setvariablebound!(m, getcol(m, v), Inf, backend_type(m, Val{:Upperbound}()))
    setvariablebound!(m, getcol(m, v), -Inf, backend_type(m, Val{:Lowerbound}()))
    delete!(dict, c)
    if !hasinteger(m)
        make_problem_type_continuous(m)
    end
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{S}}) where S <:SEMI_TYPES = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::SVCI{S}) where S <: SEMI_TYPES
    dict = constrdict(m, c)
    v = dict[c]
    lb = get_variable_lowerbound(m, getcol(m, v))
    ub = get_variable_upperbound(m, getcol(m, v))
    return S(lb, ub)
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{SVCI{S}}) where S <:SEMI_TYPES = true
MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::SVCI{<:SEMI_TYPES}) = SinVar(m[c])


text:
module PlusOneRefillFrostHelper

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/PlusOneRefill" PlusOneRefill(x::Integer, y::Integer, oneUse::Bool=false, directory::String="objects/FrostHelper/plusOneRefill", dashCount::Integer=1, respawnTime::Number=2.5, particleColor::String="White", recoverStamina::Bool=true)

const placements = Ahorn.PlacementDict(
    "Plus One Refill (Frost Helper)" => Ahorn.EntityPlacement(
        PlusOneRefill
    )
)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

Ahorn.editingOptions(entity::PlusOneRefill) = Dict{String, Any}(
    "particleColor" => colors
)

function Ahorn.selection(entity::PlusOneRefill)
    x, y = Ahorn.position(entity)
	texture = get(entity.data, "directory", "objects/FrostHelper/plusOneRefill")
	realTexture = "$texture/idle00.png"
    return Ahorn.getSpriteRectangle(realTexture, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PlusOneRefill, room::Maple.Room)
	texture = get(entity.data, "directory", "objects/FrostHelper/plusOneRefill")
	realTexture = "$texture/idle00.png"
	Ahorn.drawSprite(ctx, realTexture, 0, 0)
end

end

ext: jl

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
@testset "ECC Tests" begin

    @testset "DER Signature Parsing and Serialization" begin
        @testset "der2sig" begin
            der = hex2bytes("304402201f62993ee03fca342fcb45929993fa6ee885e00ddad8de154f268d98f083991402201e1ca12ad140c04e0e022c38f7ce31da426b8009d02832f0b44f39a6b178b7a1")
            sig = Signature(parse(BigInt, "1f62993ee03fca342fcb45929993fa6ee885e00ddad8de154f268d98f0839914", base=16),
                            parse(BigInt, "1e1ca12ad140c04e0e022c38f7ce31da426b8009d02832f0b44f39a6b178b7a1", base=16))
            @test der2sig(der) == sig
        end
        @testset "sig2der" begin
            testcases = (
                (1, 2),
                (rand(big.(0:big(2)^255)), rand(big.(0:big(2)^255))),
                (rand(big.(0:big(2)^255)), rand(big.(0:big(2)^255))),
                (rand(big.(0:big(2)^255)), rand(big.(0:big(2)^223))))
            for x in testcases
                sig = Signature(x[1], x[2])
                der = sig2der(sig)
                sig2 = der2sig(der)
                @test sig2 == sig
            end
        end
    end

    @testset "Signature Verification" begin
        pk = PrivateKey(rand(big.(0:big(2)^256)))
        𝑧 = rand(big.(0:big(2)^256))
        𝑠 = pksign(pk, 𝑧)
        @test verify(pk.𝑃, 𝑧, 𝑠)
    end

    @testset "scep256k1" begin

        @testset "Order" begin
            point = ECC.N * ECC.G
            @test typeof(point) == S256Point{ECC.Infinity}
        end

        @testset "Public Point" begin
            points = (
                # secret, x, y
                (7, big"0x5cbdf0646e5db4eaa398f365f2ea7a0e3d419b7e0330e39ce92bddedcac4f9bc", big"0x6aebca40ba255960a3178d6d861a54dba813d0b813fde7b5a5082628087264da"),
                (1485, big"0xc982196a7466fbbbb0e27a940b6af926c1a74d5ad07128c82824a11b5398afda", big"0x7a91f9eae64438afb9ce6448a1c133db2d8fb9254e4546b6f001637d50901f55"),
                (big(2)^128, big"0x8f68b9d2f63b5f339239c1ad981f162ee88c5678723ea3351b7b444c9ec4c0da", big"0x662a9f2dba063986de1d90c2b6be215dbbea2cfe95510bfdf23cbf79501fff82"),
                (big(2)^240 + 2^31, big"0x9577ff57c8234558f293df502ca4f09cbc65a6572c842b39b366f21717945116", big"0x10b49c67fa9365ad7b90dab070be339a1daf9052373ec30ffae4f72d5e66d053"),
            )

            for n ∈ points
                point = S256Point(n[2], n[3])
                @test n[1] * ECC.G == point
            end
        end

        @testset "point2sec" begin
            coefficient = 999^3
            uncompressed = "049d5ca49670cbe4c3bfa84c96a8c87df086c6ea6a24ba6b809c9de234496808d56fa15cc7f3d38cda98dee2419f415b7513dde1301f8643cd9245aea7f3f911f9"
            compressed = "039d5ca49670cbe4c3bfa84c96a8c87df086c6ea6a24ba6b809c9de234496808d5"
            point = coefficient * ECC.G
            @test point2sec(point,false) == hex2bytes(uncompressed)
            @test point2sec(point,true) == hex2bytes(compressed)
            coefficient = 123
            uncompressed = "04a598a8030da6d86c6bc7f2f5144ea549d28211ea58faa70ebf4c1e665c1fe9b5204b5d6f84822c307e4b4a7140737aec23fc63b65b35f86a10026dbd2d864e6b"
            compressed = "03a598a8030da6d86c6bc7f2f5144ea549d28211ea58faa70ebf4c1e665c1fe9b5"
            point = coefficient * ECC.G
            @test point2sec(point,false) == hex2bytes(uncompressed)
            @test point2sec(point,true) == hex2bytes(compressed)
            coefficient = 42424242
            uncompressed = "04aee2e7d843f7430097859e2bc603abcc3274ff8169c1a469fee0f20614066f8e21ec53f40efac47ac1c5211b2123527e0e9b57ede790c4da1e72c91fb7da54a3"
            compressed = "03aee2e7d843f7430097859e2bc603abcc3274ff8169c1a469fee0f20614066f8e"
            point = coefficient * ECC.G
            @test point2sec(point,false) == hex2bytes(uncompressed)
            @test point2sec(point,true) == hex2bytes(compressed)
        end

        @testset "sec2point" begin
            sec_bin = hex2bytes("0349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278a")
            point = sec2point(sec_bin)
            want = big"0xa56c896489c71dfc65701ce25050f542f336893fb8cd15f4e8e5c124dbf58e47"
            @test point.𝑦.𝑛 == want
            sec_bin = hex2bytes("049d5ca49670cbe4c3bfa84c96a8c87df086c6ea6a24ba6b809c9de234496808d56fa15cc7f3d38cda98dee2419f415b7513dde1301f8643cd9245aea7f3f911f9")
            point = sec2point(sec_bin)
            want = parse(BigInt, "6fa15cc7f3d38cda98dee2419f415b7513dde1301f8643cd9245aea7f3f911f9", base=16)
            @test point.𝑦.𝑛 == want
        end

        @testset "Signature Verification" begin
            point = S256Point(
                big"0x887387e452b8eacc4acfde10d9aaf7f6d9a0f975aabb10d006e4da568744d06c",
                big"0x61de6d95231cd89026e286df3b6ae4a894a3378e393e93a0f45b666329a0ae34")
            z = big"0xec208baa0fc1c19f708a9ca96fdeff3ac3f230bb4a7ba4aede4942ad003c0f60"
            r = big"0xac8d1c87e51d0d441be8b3dd5b05c8795b48875dffe00b7ffcfac23010d3a395"
            s = big"0x68342ceff8935ededd102dd876ffd6ba72d6a427a3edb13d26eb0781cb423c4"
            @test verify(point, z, Signature(r, s))
            z = big"0x7c076ff316692a3d7eb3c3bb0f8b1488cf72e1afcd929e29307032997a838a3d"
            r = big"0xeff69ef2b1bd93a66ed5219add4fb51e11a840f404876325a1e8ffe0529a2c"
            s = big"0xc7207fee197d27c618aea621406f6bf5ef6fca38681d82b2f06fddbdce6feab6"
            @test verify(point, z, Signature(r, s))
        end
    end

end


text:
"""
    sixel_encode([io], src, [encoder]; kwargs...)

Encode colorant sequence `src` as sixel sequence and write it into a writable io-like object `io`.

# Arguments

- `io`: the output io-like object, which is expected to be writable. Typical choices are: `stdout`,
  `IOBuffer` and `IOStream`. The default value is `stdout`.
- `src`: generic container object(e.g., `AbstractArray`) that contains the colorant sequence.
- `encoder::AbstractSixelEncoder`: the sixel encoder. Currently, only
  [`LibSixelEncoder`](@ref Sixel.LibSixel.LibSixelEncoder) is available.

!!! warning
    `sixel_encode` for tiny images (including vectors) is undefined behavior; sixel format requires
    at least `6` pixels in the row direction. For better visualization result, it's recommended to
    `repeat`(inner) the data so that its size is larger than `(6, 6)` in both dimensions.

!!! note
    For array with `ndims(src) >= 3`, it will be concatenated into a 2d array in row order before
    encoding.

# Parameters

- `transpose::Bool`: whether we need to permute the image's width and height dimension before encoding.
  The default value is `false`.

# Examples

```julia
using Sixel, TestImages

# 2d RGB image
img = testimage("mandril_color")
sixel_encode(img)

# 3d RGB image will be concatenated in row order.
img3d = testimage("mri-stack")
sixel_encode(img3d)

# You can also write the encoded sixel sequences into an `IOBuffer` or `IOStream`(via `open`)
io = IOBuffer()
sixel_encode(io, img)
println(String(take!(io)))

open("out.sixel", "w") do io
    sixel_encode(io, img)
end
```

# References

- [1] VT330/VT340 Programmer Reference Manual, Volume 1: Text Programming
- [2] VT330/VT340 Programmer Reference Manual, Volume 2: Graphics Programming
- [3] https://github.com/saitoha/libsixel
"""
function sixel_encode end

sixel_encode(img::AbstractArray, enc::SEC=default_encoder(img); kwargs...) =
    sixel_encode(stdout, img, enc; kwargs...)

sixel_encode(io::IO, img::AbstractVector, enc::SEC=default_encoder(img); kwargs...) =
    sixel_encode(io, reshape(img, :, 1), enc; kwargs...)

function sixel_encode(io::IO, img::AbstractArray, enc::SEC=default_encoder(img); transpose=false, kwargs...)
    # Conversion from OffsetArrays to Array is not very well supported, so we have to de-offset first.
    img = OffsetArrays.no_offset_view(img)

    # make sure it always tiles along row order
    @assert ndims(img) >= 3
    nrow = transpose ? prod(size(img)[3:end]) : 1
    sixel_encode(io, mosaicview(img; nrow=nrow, rowmajor=true), enc; transpose=transpose, kwargs...)
end

# This is expected to be the only caller of `enc(io, bytes)`; all other methods should eventually
# call this one.
function sixel_encode(
        io::IO,
        img::AbstractMatrix,
        enc::SEC=default_encoder(img);
        transpose=false, kwargs...)
    # Conversion from OffsetArrays to Array is not very well supported, so we have to de-offset first.
    img = OffsetArrays.no_offset_view(img)
    if enc isa LibSixelEncoder
        T = canonical_sixel_eltype(enc, eltype(img))
        AT = Array{T, ndims(img)}
        # libsixel is a C library and assumes row-major memory order, thus `collect` the data into
        # contiguous memeory layout already makes a transpose.
        bytes = transpose ? convert(AT, img) : convert(AT, PermutedDimsArray(img, (2, 1)))
    else
        throw(ArgumentError("Unsupported encoder type $(typeof(enc)). Please open an issue for this."))
    end
    enc(io, bytes; kwargs...)
    return nothing
end


text:
using HierarchicalTemporalMemory
using Documenter

DocMeta.setdocmeta!(HierarchicalTemporalMemory, :DocTestSetup, :(
    using HierarchicalTemporalMemory;
    const HTM=HierarchicalTemporalMemory
  ); recursive=true)

makedocs(;
    modules=[HierarchicalTemporalMemory],
    authors="Konstantinos Samaras-Tsakiris <ksamtsak@gmail.com> and contributors",
    repo="https://github.com/Oblynx/HierarchicalTemporalMemory.jl/blob/{commit}{path}#L{line}",
    sitename="HierarchicalTemporalMemory.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://oblynx.github.io/HierarchicalTemporalMemory.jl",
        assets=String[],
    ),
    pages=[
           "Home" => ["index.md","poster.md"],
        "Examples" => "examples.md",
        "Reference" => "reference.md"
    ],
)

deploydocs(;
    repo="github.com/Oblynx/HierarchicalTemporalMemory.jl",
    push_preview= true,
)


text:
using ODEInterfaceDiffEq, DiffEqProblemLibrary, DiffEqBase
using Test

jac_called = false

function Lotka(du,u,p,t)
  du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
  du[2] = -3 * u[2] + 1 * u[1] * u[2]
  nothing
end

function Lotka_jac(J,u,p,t)
  global jac_called
  jac_called = true
  J[1,1] = 1.0 - u[2]
  J[1,2] = -u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

prob = ODEProblem(ODEFunction(Lotka,jac=Lotka_jac),ones(2),(0.0,2.0))

sol =solve(prob,radau5();dt=1//2^(4))

@test jac_called == true


text:
printstyled("  ASDF\n", color=:light_green)
printstyled("    read_hdf5\n", color=:light_green)

hdf       = path*"/SampleFiles/HDF5/2days-40hz.h5"
hdf_pat   = path*"/SampleFiles/HDF5/2days-40hz.h*"
hdf_out1  = "test1.h5"
hdf_out2  = "test2.h5"
hdf_out3  = "test3.h5"
hdf_out4  = "test4.h5"
hdf_evt   = path*"/SampleFiles/HDF5/example.h5"
safe_rm(hdf_out1)
safe_rm(hdf_out2)
safe_rm(hdf_out3)
safe_rm(hdf_out4)

id  = "CI.SDD..HHZ"
idr = "C*.SDD..HH?"
ts  = "2019-07-07T23:00:00"
te  = "2019-07-08T02:00:00"

# id_to_regex
@test id_to_regex("*.*.*.*") == r".*\.*\.*\.*"
@test id_to_regex("C*.SDD..") == r"C.*\.SDD\.\."
@test id_to_regex("CI.SDD..HHZ") == r"CI\.SDD\.\.HHZ"
@test id_to_regex("C*.SDD..HH?") == r"C.*\.SDD\.\.HH."

# id_match
S2 = SeisData(SeisChannel(), SeisChannel(id="CI.SDD..HHZ"))
@test id_match(id, S2) == id_match(idr, S2) == [2]

S1 = SeisData()
read_hdf5!(S1, hdf, ts, te, id = id)
S2 = read_hdf5(hdf, ts, te, id = id)
@test S1 == S2
S2 = read_asdf(hdf, id, ts, te, true, 0)
@test S1.src[1] == abspath(hdf)
S1.src = S2.src
@test S1 == S2

# check file wildcards
S2 = read_hdf5(hdf_pat, DateTime(ts), DateTime(te), id = id)
S1.src = S2.src
@test S1 == S2

# check the default id
S2 = read_hdf5(hdf_pat, ts, te)
S1.src = S2.src
@test S1 == S2

# check that FDSN-style wildcards work
S2 = read_asdf(hdf, idr, ts, te, true, 0)
S1.src = S2.src
@test S1 == S2

# Check channel matching
S2 = SeisData(SeisChannel(id="CI.SDD..HHZ"))
read_asdf!(S2, hdf, idr, ts, te, true, 0)
S1.src = S2.src
@test S1 == S2

S2 = SeisData(SeisChannel(id="CI.SDD..HHZ",
                          fs=40.0,
                          t=[1 1562543940000000; 40 0],
                          x=randn(40)))
read_asdf!(S2, hdf, idr, ts, te, true, 0)
@test S1.x[1] == S2.x[1][41:end]
@test string(u2d(S2.t[1][1,2]*μs)) == "2019-07-07T23:59:00"

@test_throws ErrorException read_hdf5(hdf, ts, te, fmt="MatlabLol")

printstyled("    scan_hdf5\n", color=:light_green)
@test scan_hdf5(hdf) == ["CI.SDD"]
@test scan_hdf5(hdf, level="trace") == ["/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_", "/Waveforms/CI.SDD/StationXML"]
@test_throws ErrorException scan_hdf5(hdf, level="response")
@test_throws ErrorException scan_hdf5(hdf, fmt="MatlabLol")

# HDF event read
printstyled("    asdf_rqml\n", color=:light_green)
(H,R) = asdf_rqml(hdf_evt)
@test H[1].id == "20120404_0000041"
@test H[2].id == "20120404_0000038"
@test H[3].id == "20120404_0000039"
@test H[1].loc.lat ≈ 41.818
@test H[1].loc.lon ≈ 79.689
@test H[1].loc.dep ≈ 1.0
@test H[1].mag.val ≈ 4.4
@test H[1].mag.scale == "mb"
@test H[2].loc.lat ≈ 39.342
@test H[2].loc.dep ≈ 14.4
@test H[2].mag.val ≈ 4.3
@test H[2].mag.scale == "ML"

io = h5open(hdf_evt)
(H1, R1) = asdf_rqml(io)
@test H == H1
@test R == R1

# =========================================================
# HDF write
printstyled("    write_hdf5\n", color=:light_green)
@test_throws ErrorException write_hdf5(hdf_out1, S1, fmt="MatlabLol")
S2 = deepcopy(S1) + SeisChannel(id="XX.SDD..YYZ", fs=40.0)
write_hdf5(hdf_out1, S2)
safe_rm(hdf_out1)

# ASDF write test 1: can we write to a new file?
printstyled("      write to new file\n", color=:light_green)
write_hdf5( hdf_out1, S1 )
S2 = read_hdf5(hdf_out1, ts, te, id = id)
for f in datafields
  (f in (:src, :notes)) && continue
  @test isequal(getfield(S1, f), getfield(S2, f))
end

# ASDF write test 2: can we overwrite parts of an existing file?
printstyled("      add to existing file\n", color=:light_green)
ts  = "2019-07-08T00:00:00"
te  = "2019-07-08T02:00:00"
S3 = read_hdf5(hdf_out1, ts, te, id = id)
for i in 1:S3.n
  S3.x[i] .= (2.0.*S3.x[i])
end
redirect_stdout(out) do
  write_hdf5( hdf_out1, S3, ovr=true, v=3 )
  push!(S3, SeisChannel(id="YY.ZZTOP.00.LEG", fs=50.0, x=randn(1024)))

  # This should now fail since :x[2] has no :t[2]
  # @test_throws ArgumentError write_hdf5( hdf_out1, S3, v=3 )

  # This should fail
  @test_throws ErrorException write_hdf5( hdf_out1, S3, ovr=true, v=3 )

  # GphysChannel extension
  C = S3[1]
  write_hdf5( hdf_out1, C, v=3 )
  safe_rm(hdf_out1)
  deleteat!(S3, 2)

  # Force write to channel with existing net.sta
  write_hdf5(hdf_out1, S3)
  ts  = "2019-07-08T10:00:00"
  te  = "2019-07-08T12:00:00"
  S4 = read_hdf5(hdf, ts, te)
  write_hdf5(hdf_out1, S4, ovr=true, add=true)
end

# There should be 4 successful writes at this point
@test sum([occursin("wrote to file " * hdf_out1, S3.notes[1][i]) for i in 1:length(S3.notes[1])]) == 4

@test scan_hdf5(hdf_out1, level="trace") == [
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T02:00:00__hhz",
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T23:59:59.975__hhz",
  "/Waveforms/CI.SDD/StationXML"
  ]
safe_rm(hdf_out1)

# ASDF write test 3: write with gaps
printstyled("      write to new file with gaps\n", color=:light_green)
ts  = "2019-07-08T10:00:00"
te  = "2019-07-08T12:00:00"
read_hdf5!(S3, hdf, ts, te)
merge!(S3)
@test S3.n == 1 # fails if MultiStageResp merge! bug is back
write_hdf5( hdf_out1, S3 )
scan3 = scan_hdf5(hdf_out1, level="trace")
@test scan3 == [
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T02:00:00__hhz",
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T10:00:00__2019-07-08T12:00:00__hhz",
  "/Waveforms/CI.SDD/StationXML"
  ]

# ASDF write test 4: file with multiple stations and channels
if has_restricted
  S = verified_read_data("sac", path*"/SampleFiles/Restricted/20140927000000*SAC")
  redirect_stdout(out) do
    write_hdf5( hdf_out1, S, v=3 )
  end
  id = S.id[3]
  ts = "2014-09-27T09:00:00"

  # Test for intended read behavior
  printstyled("        are written traces the right length?\n", color=:light_green)
  te = "2014-09-27T10:00:00"
  C = read_hdf5(hdf_out1, ts, te, id = id)[1]
  @test length(C.x) == 360000 # Stop at last available sample

  te = "2014-09-27T09:59:59.99"
  C = read_hdf5(hdf_out1, ts, te, id = id)[1]
  @test length(C.x) == 360000 # Exact

  ts = "2014-09-27T09:00:00.01"
  te = "2014-09-27T09:59:59.90"
  C2 = read_hdf5(hdf_out1, ts, te, id = id)[1]
  @test length(C2.x) == 359990 # Exact, ten samples shorter
  @test C.x[2:5] == C2.x[1:4]

  # test 1: are the right trace names created?
  printstyled("        are the right trace names created?\n", color=:light_green)
  S1 = SeisData()
  read_hdf5!(S1, hdf, "2019-07-07T23:00:00", "2019-07-08T02:00:00", id = "CI.SDD..HHZ")

  append!(S, S1)
  safe_rm(hdf_out1)
  redirect_stdout(out) do
    write_hdf5(hdf_out1, S, add=true, ovr=true, v=3, tag="raw")
  end
  scan_expect = VERSION >= v"1.6" ? [
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-07T23:59:59.975__raw",
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T23:59:59.975__raw",
  "/Waveforms/CI.SDD/StationXML",
  "/Waveforms/JP.VONTA/JP.VONTA..E__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTA/JP.VONTA..H__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTA/JP.VONTA..N__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTA/JP.VONTA..U__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTA/StationXML",
  "/Waveforms/JP.VONTN/JP.VONTN..E__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTN/JP.VONTN..H__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTN/JP.VONTN..N__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTN/JP.VONTN..U__2014-09-27T00:00:00__2014-09-27T23:59:59.990__raw",
  "/Waveforms/JP.VONTN/StationXML"
  ] :
  [
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-07T23:59:59.975__raw",
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T23:59:59.975__raw",
  "/Waveforms/CI.SDD/StationXML",
  "/Waveforms/JP.VONTA/JP.VONTA..E__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTA/JP.VONTA..H__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTA/JP.VONTA..N__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTA/JP.VONTA..U__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTA/StationXML",
  "/Waveforms/JP.VONTN/JP.VONTN..E__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTN/JP.VONTN..H__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTN/JP.VONTN..N__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTN/JP.VONTN..U__2014-09-27T00:00:00__2014-09-27T23:59:59.99__raw",
  "/Waveforms/JP.VONTN/StationXML"
  ]
  @test scan_hdf5(hdf_out1, level="trace") == scan_expect

  # test 2: ONLY write CI.SDD..HHZ, JP.VONTA..H, JP.VONTA..N, JP.VONTA..U
  printstyled("        write with a channel sublist\n", color=:light_green)
  redirect_stdout(out) do
    write_hdf5(hdf_out2, S, chans=[2,3,4,9], add=true, ovr=true, v=3)
  end
  scan_expect = VERSION >= v"1.6" ? [
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-07T23:59:59.975__hhz",
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T23:59:59.975__hhz",
  "/Waveforms/CI.SDD/StationXML",
  "/Waveforms/JP.VONTA/JP.VONTA..H__2014-09-27T00:00:00__2014-09-27T23:59:59.990__h",
  "/Waveforms/JP.VONTA/JP.VONTA..N__2014-09-27T00:00:00__2014-09-27T23:59:59.990__n",
  "/Waveforms/JP.VONTA/JP.VONTA..U__2014-09-27T00:00:00__2014-09-27T23:59:59.990__u",
  "/Waveforms/JP.VONTA/StationXML"
  ] : [
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-07T23:59:59.975__hhz",
  "/Waveforms/CI.SDD/CI.SDD..HHZ__2019-07-08T00:00:00__2019-07-08T23:59:59.975__hhz",
  "/Waveforms/CI.SDD/StationXML",
  "/Waveforms/JP.VONTA/JP.VONTA..H__2014-09-27T00:00:00__2014-09-27T23:59:59.99__h",
  "/Waveforms/JP.VONTA/JP.VONTA..N__2014-09-27T00:00:00__2014-09-27T23:59:59.99__n",
  "/Waveforms/JP.VONTA/JP.VONTA..U__2014-09-27T00:00:00__2014-09-27T23:59:59.99__u",
  "/Waveforms/JP.VONTA/StationXML"
  ]

  @test scan_hdf5(hdf_out2, level="trace") == scan_expect

  # ensure the right trace data were written to these names
  S1 = read_hdf5(hdf_out2, "2014-09-27T09:00:00", "2014-09-27T09:59:59.99", id="*.VONTA..*", msr=false)
  S2 = S[[2,3,4]]
  for f in (:id, :loc, :fs, :gain, :resp, :t, :x)
    @test getfield(S1,f) == getfield(S2,f)
  end

  # ensure what's read back in are the correct traces with NaNs in the right places
  printstyled("        trace indexing with ovr=true\n", color=:light_green)
  S1 = read_hdf5(hdf_out2, "2014-09-27T08:00:00.00", "2014-09-27T10:00:00.00", id="*.VONTA..*")
  for i in 1:S1.n
    x = S1.x[i]
    @test length(x) == 720001
    @test isnan(last(x))
    @test eltype(x) == Float32
    for j in 1:360000
      @test isnan(x[j])
    end
  end

  # only an hour of non-NaN data; 24 hours of data created
  S1 = read_hdf5(hdf_out2, "2014-09-27T00:00:00.00", "2014-09-28T00:00:00.00", id="JP.VONTA..H")
  @test length(S1.x[1]) == 8640000
  @test length(S1.x[1])-length(findall(isnan.(S1.x[1]))) == 360000

  # data from CI.SDD should have no NaNs in this range; it's the exact range of S.x[9]
  S1 = read_hdf5(hdf_out2, "2019-07-07T23:00:00", "2019-07-08T02:00:00")
  @test S1.x[1] == S.x[9]
end

# HDF write with SeisEvent
Ev = Array{SeisEvent,1}(undef,3)
for i in 1:3
  Ev[i] = rse_wb(3)
end

printstyled("      asdf_wqml\n", color=:light_green)

printstyled("        append mode\n", color=:light_green)
SHDR = [Ev[i].hdr for i in 1:3]
SSRC = [Ev[i].source for i in 1:3]
asdf_wqml(hdf_out2, Ev[1].hdr, Ev[1].source)
asdf_wqml(hdf_out2, Ev[2])
asdf_wqml(hdf_out2, Ev)
asdf_wqml(hdf_out2, SHDR, SSRC)
(H,R) = asdf_rqml(hdf_out2)
@test length(H) == length(R) == 8

# Check that data are appended in the right order
inds = [1,2,1,2,3,1,2,3]
for i in 1:length(H)
  j = inds[i]
  compare_SeisHdr(Ev[j].hdr, H[i])
  compare_SeisSrc(Ev[j].source, R[i])
end

printstyled("        overwrite mode\n", color=:light_green)
asdf_wqml(hdf_out2, Ev[2].hdr, Ev[2].source, ovr=true, v=1)
(H,R) = asdf_rqml(hdf_out2)
@test length(H) == length(R) == 1
compare_SeisHdr(Ev[2].hdr, H[1])
compare_SeisSrc(Ev[2].source, R[1])

printstyled("        to new file\n", color=:light_green)
asdf_wqml(hdf_out4, SHDR, SSRC)
(H,R) = asdf_rqml(hdf_out4)
for i in 1:3
  compare_SeisHdr(SHDR[i], H[i])
  compare_SeisSrc(SSRC[i], R[i])
end

printstyled("        append file with incompatible QML\n", color=:light_green)

# Create incompatible QML in file
xml_evfile1 = path*"/SampleFiles/XML/fdsnws-event_2017-01-12T03-18-55Z.xml"
xml_buf = read(xml_evfile1)
io = h5open(hdf_out4, "r+")
o_delete(io, "QuakeML")
io["QuakeML"] = xml_buf
close(io)

# Now try to write to the file...does it work?
asdf_wqml(hdf_out4, SHDR, SSRC)
(H,R) = asdf_rqml(hdf_out4)
L = length(H)
H = H[L-2:L]
R = R[L-2:L]
for i in 1:3
  compare_SeisHdr(SHDR[i], H[i])
  compare_SeisSrc(SSRC[i], R[i])
end

printstyled("      write SeisEvent\n", color=:light_green)
printstyled("        to new file\n", color=:light_green)
write_hdf5(hdf_out3, Ev[1])

printstyled("        to existing file\n", color=:light_green)
@test_throws ErrorException write_hdf5(hdf_out3, Ev[2], fmt="GarfieldTheCat")
write_hdf5(hdf_out3, Ev[2], chans=[1,2])

printstyled("        to appended file\n", color=:light_green)
push!(Ev[3].data, EventChannel(id="AA.STA.00.YYY"))
write_hdf5(hdf_out3, Ev[3], chans=[3,4])

printstyled("      read_asdf_evt\n", color=:light_green)
EvCat = read_asdf_evt(hdf_out3, Ev[1].hdr.id, msr=false)

printstyled("        accuracy of SeisEvent i/o\n", color=:light_green)

printstyled("          single-event read\n", color=:light_green)
W = EvCat[1]
Ev1 = deepcopy(Ev[1])
compare_events(Ev1, W)

EvCat = read_asdf_evt(hdf_out3, Ev[2].hdr.id, msr=false)
W = EvCat[1]
Ev2 = deepcopy(Ev[2])
Ev2.data = Ev2.data[1:2]
compare_events(Ev2, W)

EvCat = read_asdf_evt(hdf_out3, Ev[3].hdr.id, msr=false)
W = EvCat[1]
Ev3 = deepcopy(Ev[3])
Ev3.data = Ev3.data[3]
compare_events(Ev3, W)

printstyled("          multi-event read\n", color=:light_green)
EvCat = read_asdf_evt(hdf_out3, msr=false)
compare_events(EvCat[1], Ev1)
compare_events(EvCat[2], Ev2)
compare_events(EvCat[3], Ev3)

if Sys.iswindows() == false
  printstyled("          multi-file read\n", color=:light_green)
  write_hdf5(hdf_out4, Ev1)

  # Check logging
  @test sum([occursin("wrote to file ", Ev1.data.notes[1][i]) for i in 1:length(Ev1.data.notes[1])]) == 2
  @test any([occursin("wrote to file " * hdf_out3, Ev1.data.notes[1][i]) for i in 1:length(Ev1.data.notes[1])])
  @test any([occursin("wrote to file " * hdf_out4, Ev1.data.notes[1][i]) for i in 1:length(Ev1.data.notes[1])])
  @test any([occursin("wrote to file " * hdf_out3, Ev1.hdr.notes[i]) for i in 1:length(Ev1.hdr.notes)])
  @test any([occursin("wrote to file " * hdf_out4, Ev1.hdr.notes[i]) for i in 1:length(Ev1.hdr.notes)])
  @test any([occursin("wrote to file " * hdf_out3, Ev1.source.notes[i]) for i in 1:length(Ev1.source.notes)])
  @test any([occursin("wrote to file " * hdf_out4, Ev1.source.notes[i]) for i in 1:length(Ev1.source.notes)])

  EC2 = read_asdf_evt("test[3-4].h5", msr=true)
  inds = [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 1]
  for i in 1:length(EC2)
    j = inds[i]
    if j > 0
      compare_SeisHdr(EC2[i].hdr, EvCat[j].hdr)
      compare_SeisSrc(EC2[i].source, EvCat[j].source)
    end
  end
end

# HDF write cleanup
safe_rm(hdf_out1)
safe_rm(hdf_out2)
safe_rm(hdf_out3)
safe_rm(hdf_out4)

# asdf_waux
# write to new file
X = randn(32768, 2)
pstr = "CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_"
asdf_waux(hdf_out1, pstr, X)

# does overwrite work?
asdf_waux(hdf_out1, pstr, X)
f = h5open(hdf_out1, "r")
aux = f["AuxiliaryData"]
@test names(aux) == [ "CI.SDD" ]
Xr = read(aux[pstr])
@test Xr == X
close(f)

# can we correct for a starting / in the string?
pstr = "//" * pstr
X = randn(32768, 2)
asdf_waux(hdf_out1, pstr, X)
f = h5open(hdf_out1, "r")
aux = f["AuxiliaryData"]
@test names(aux) == [ "CI.SDD" ]
pstr = "CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_"
@test names(aux["CI.SDD"]) == [ String(split(pstr,"/")[2]) ]
Yr = read(aux[pstr])
@test Yr == X
@test Xr != Yr
close(f)

# what about path that starts with /AuxiliaryData/?
pstr = "/AuxiliaryData/" * pstr
X = randn(32768, 2)
asdf_waux(hdf_out1, pstr, X)
f = h5open(hdf_out1, "r")
aux = f["AuxiliaryData"]
@test names(aux) == [ "CI.SDD" ]
pstr = "CI.SDD/CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_"
@test names(aux["CI.SDD"]) == [ String(split(pstr,"/")[2]) ]
Yr = read(aux[pstr])
@test Yr == X
@test Xr != Yr
close(f)


text:
function total_population(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(length(y_true))
end

function true_positive(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum((y_true .== y_pred) & (y_true .== 1)))
end

function true_negative(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum((y_true .== y_pred) & (y_true .== 0)))
end

function false_positive(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(!(y_true .== y_pred) & (y_true .== 1)))
end

function false_negative(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(!(y_true .== y_pred) & (y_true .== 0)))
end

function condition_positive(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(y_true .== 1))
end

function condition_negative(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(y_true .== 0))
end

function predicted_condition_positive(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(y_true .== 1))
end

function predicted_condition_negative(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(y_pred .== 0))
end

function accuracy_score(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tp = true_positive(y_true, y_pred)
    tn = true_negative(y_true, y_pred)
    tot_pop = total_population(y_true, y_pred)
    return((tp + tn) / tot_pop)
end

accuracy = accuracy_score

function prevalence(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    cp = condition_positive(y_true, y_pred)
    tot_pop = total_population(y_true, y_pred)
    return(cp / tot_pop)
end

function positive_predictive_value(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tp = true_positive(y_true, y_pred)
    pcp = predicted_condition_positive(y_true, y_pred)
    return(tp / pcp)
end

precision = positive_predictive_value

function false_discovery_rate(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    fp = false_positive(y_true, y_pred)
    pcp = predicted_condition_positive(y_true, y_pred)
    return(fp / pcp)
end

function negative_predictive_value(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tn = true_negative(y_true, y_pred)
    pcn = predicted_condition_negative(y_true, y_pred)
    return(tn / pcn)
end

function false_omission_rate(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    fn = false_negative(y_true, y_pred)
    pcn = predicted_condition_negative(y_true, y_pred)
    return(fn / pcn)
end

function true_positive_rate(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tp = true_positive(y_true, y_pred)
    cp = condition_positive(y_true, y_pred)
    return(tp / cp)
end

sensitivity = true_positive_rate
recall = true_positive_rate

function false_positive_rate(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    fp = false_positive(y_true, y_pred)
    cn = condition_negative(y_true, y_pred)
    return(fp / cn)
end

function false_negative_rate(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    fn = false_negative(y_true, y_pred)
    cp = condition_positive(y_true, y_pred)
    return(fn / cp)
end

function true_negative_rate(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tn = true_negative(y_true, y_pred)
    cn = condition_negative(y_true, y_pred)
    return(tn / cn)
end

specificity = true_negative_rate

function positive_likelihood_ratio(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tpr = true_positive_rate(y_true, y_pred)
    fpr = false_positive_rate(y_true, y_pred)
    return(tpr / fpr)
end

function negative_likelihood_ratio(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    fnr = false_negative_rate(y_true, y_pred)
    tnr = true_negative_rate(y_true, y_pred)
    return(fnr / tnr)
end

function diagnostic_odds_ratio(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    plr = positive_likelihood_ratio(y_true, y_pred)
    nlr = negative_likelihood_ratio(y_true, y_pred)
    return(plr / nlr)
end

function f1_score(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tp = true_positive(y_true, y_pred)
    fp = false_positive(y_true, y_pred)
    fn = false_negative(y_true, y_pred)
    return((2 * tp) / (2 * tp + fp + fn))
end

function matthews_corrcoef(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    tp = true_positive(y_true, y_pred)
    tn = true_negative(y_true, y_pred)
    fp = false_positive(y_true, y_pred)
    fn = false_negative(y_true, y_pred)
    numerator = (tp * tn) - (fp * fn)
    denominator = (tp + fp) * (tp + fn) * (tn + fp) * (tn + fn)
    return(numerator / (denominator ^ 0.5))
end


text:
"Holds the information for the various methods' tableaus."
abstract type AbstractTableau{T} end

# order(tab::T) where {T <: AbstractTableau} = tab.o
RungeKutta.nstages(tab::T) where {T <: AbstractTableau} = tab.s

# function check_symplecticity end
function symplecticity_conditions end

@define HeaderTableau begin
    name::Symbol
    o::Int
    s::Int
end


text:
module TestCore

using AbstractYieldMacros
using ContextualMacros
using Test

@testset "expandwith" begin
    ex = ContextualMacros.expandwith(
        @__MODULE__,
        quote
            @yield 1 2 3
            @yieldfrom 4 5
        end;
        yield = ctx -> QuoteNode(ctx.args),
        yieldfrom = ctx -> QuoteNode(ctx.args),
    )
    nodes = filter(x -> x isa QuoteNode, ex.args)
    @test nodes == [QuoteNode((1, 2, 3)), QuoteNode((4, 5))]
end

end  # module


text:
# # Geostrophic adjustment using Oceananigans.HydrostaticFreeSurfaceModel
#
# This example demonstrates how to simulate the one-dimensional geostrophic adjustment of a
# free surface using `Oceananigans.HydrostaticFreeSurfaceModel`. Here, we solve the hydrostatic
# Boussinesq equations beneath a free surface with a small-amplitude about rest ``z = 0``,
# with boundary conditions expanded around ``z = 0``, and free surface dynamics linearized under
# the assumption ``η / H \ll 1``, where ``η`` is the free surface displacement, and ``H`` is
# the total depth of the fluid.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## A one-dimensional domain
#
# We use a one-dimensional domain of geophysical proportions,

using Oceananigans
using Oceananigans.Units: hours, meters, kilometers

grid = RegularRectilinearGrid(size = (128, 1, 1),
                            x = (0, 1000kilometers), y = (0, 1), z = (-400meters, 0),
                            topology = (Bounded, Periodic, Bounded))

# and Coriolis parameter appropriate for the mid-latitudes on Earth,

coriolis = FPlane(f=1e-4)

# ## Building a `HydrostaticFreeSurfaceModel`
#
# We use `grid` and `coriolis` to build a simple `HydrostaticFreeSurfaceModel`,

model = HydrostaticFreeSurfaceModel(grid=grid, coriolis=coriolis)

# ## A geostrophic adjustment initial value problem
#
# We pose a geostrophic adjustment problem that consists of a partially-geostrophic
# Gaussian height field complemented by a geostrophic ``y``-velocity,

Gaussian(x, L) = exp(-x^2 / 2L^2)

U = 0.1 # geostrophic velocity
L = grid.Lx / 40 # Gaussian width
x₀ = grid.Lx / 4 # Gaussian center

vᵍ(x, y, z) = - U * (x - x₀) / L * Gaussian(x - x₀, L)

g = model.free_surface.gravitational_acceleration

η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude

ηᵍ(x, y, z) = η₀ * Gaussian(x - x₀, L)

# We use an initial height field that's twice the geostrophic solution,
# thus superimposing a geostrophic and ageostrophic component in the free
# surface displacement field:

ηⁱ(x, y, z) = 2 * ηᵍ(x, y, z)

# We set the initial condition to ``vᵍ`` and ``ηⁱ``,

set!(model, v=vᵍ, η=ηⁱ)

# ## Running a `Simulation`
#
# We pick a time-step that resolves the surface dynamics,

gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

wave_propagation_time_scale = model.grid.Δx / gravity_wave_speed

simulation = Simulation(model, Δt = 0.1 * wave_propagation_time_scale, stop_iteration = 1000)

# ## Output
#
# We output the velocity field and free surface displacement,

output_fields = merge(model.velocities, (η=model.free_surface.η,))

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = IterationInterval(10),
                                                      prefix = "geostrophic_adjustment",
                                                      force = true)

run!(simulation)

# ## Visualizing the results

using JLD2, Plots, Printf

xη = xw = xv = xnodes(model.free_surface.η)
xu = xnodes(model.velocities.u)

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)

    u = file["timeseries/u/$iter"][:, 1, 1]
    v = file["timeseries/v/$iter"][:, 1, 1]
    η = file["timeseries/η/$iter"][:, 1, 1]
    t = file["timeseries/t/$iter"]

    titlestr = @sprintf("Geostrophic adjustment at t = %.1f hours", t / hours)

    v_plot = plot(xv / kilometers, v, linewidth = 2, title = titlestr,
                  label = "", xlabel = "x (km)", ylabel = "v (m s⁻¹)", ylims = (-U, U))

    u_plot = plot(xu / kilometers, u, linewidth = 2,
                  label = "", xlabel = "x (km)", ylabel = "u (m s⁻¹)", ylims = (-2e-3, 2e-3))

    η_plot = plot(xη / kilometers, η, linewidth = 2,
                  label = "", xlabel = "x (km)", ylabel = "η (m)", ylims = (-η₀/10, 2η₀))

    plot(v_plot, u_plot, η_plot, layout = (3, 1), size = (800, 600))
end

close(file)

mp4(anim, "geostrophic_adjustment.mp4", fps = 15) # hide


text:
using ITensors, Test

@testset "Itensors.jl" begin
    @testset "$filename" for filename in (
        "test_tagset.jl",
        "test_index.jl",
        "test_indexset.jl",
        "test_itensor.jl",
        "test_contract.jl",
        "test_trg.jl",
        "test_ctmrg.jl",
        "test_siteset.jl",
        "test_mps.jl",
        "test_mpo.jl",
        "test_autompo.jl",
        "test_svd.jl",
    )
        include(filename)
    end
end


text:
const MIN_PHRED = Phred(1)
const MAX_PHRED = Phred(Int('~') - 33)

"""Convert error probability to PHRED score"""
function p_to_phred(p::Prob)
    return Phred(min(round(-10.0 * log10(p)), MAX_PHRED))
end

function p_to_phred(x::Vector{LogProb})
    return Phred[p_to_phred(p) for p in x]
end

"""Convert PHRED score to log error probability"""
@generated function phred_to_log_p(x)
    return quote
        return x / (-10.0)
    end
end

"""Convert PHRED score to error probability"""
function phred_to_p(q::Phred)
    return exp10(phred_to_log_p(q))
end

function phred_to_p(x::Vector{Phred})
    return exp10.(phred_to_log_p(x))
end

@generated function normalize(parts)
    return quote
        return parts / sum(parts)
    end
end

function cap_phreds(phreds::Vector{Phred}, max_phred::Phred)
    if max_phred < 1
        error("max phred value must be positive")
    end
    return Phred[min(p, max_phred) for p in phreds]
end


text:
@testset "check gradient for Chain" begin
    clear()

    # [1] prepare input data and its label
    TYPE = Array{Float64}
    x = randn(256, 62)
    l = rand(64, 62)
    l = l ./ sum(l, dims=1)
    x = Variable(x; type=TYPE)
    l = Variable(l; type=TYPE)

    blocks = [
        Dense(256,128,tanh; type=TYPE),
        Dense(128,128,cos;  type=TYPE),
        Dense(128,128,sin;  type=TYPE),
        Dense(128,128,cos!; type=TYPE),
        Dense(128,128,sin!; type=TYPE),
        Dense(128, 64,relu; type=TYPE),
        Maxout(64, 64; k=3, type=TYPE),
        Residual(
            Dense(64,32,sin;type=TYPE),
            Linear(32,64,   type=TYPE)
        )
    ]

    model = Chain(blocks)

    # [2] forward and backward propagation
    outs = softmax(forward(model, x); dims=1)
    COST1 = crossEntropyCost(outs, l)
    backward()

    # [3] with a samll change of a weight
    DELTA = 1e-5
    model[1].w.value[1] += DELTA
    GRAD = blocks[1].w.delta[1]

    # [4] forward and backward propagation
    outs = softmax(forward(model, x); dims=1)
    COST2 = crossEntropyCost(outs, l)
    backward()

    # [5] check if the auto-grad is true or not
    dLdW = (COST2 - COST1)/DELTA
    err  = abs((dLdW-GRAD)/(GRAD==0.0 ? 1.0 : GRAD))*100
    err  = err < 1e-3 ? 0.0 : err
    @test err < 5.0
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
module Mnlitest2

using GeneralFunctions
using OffsetArrays
using TwoFAST

export generateMnli

##Fix minus signs. Not going to repaste, but need to go over every term and look for changed signs FML

function generateMnli(Pk::Function, BiasDict22::Dict{String,Int64}, f::Float64) #where {F<:Function}
    M = OffsetArray{Array{Float64}}(undef, 0:4, 0:8, 0:8)

    b1 = BiasDict22["b1"]
    bη = BiasDict22["bη"]
    b2 = BiasDict22["b2"]
    bK² = BiasDict22["bK2"]
    bδη = BiasDict22["bδη"]
    bη² = BiasDict22["bη2"]
    bKKpara = BiasDict22["bKK∥"]
    bΠ2para = BiasDict22["bΠ2∥"]

    #@time begin
    r, xi20 = ξ(Pk,2,0)
    _, xi00 = ξ(Pk,0,0)
    _, xi1m1 = ξ(Pk,1,-1)
    _, xi3m1 = ξ(Pk,3,-1)
    _, xi2m2 = ξ(Pk,2,-2)
    _, xi4m2 = ξ(Pk,4,-2)
    _, xi10 = ξ(Pk,1,0)
    _, xi30 = ξ(Pk,3,0)
    _, xi4m1 = ξ(Pk,4,-1)
    _, xi2m1 = ξ(Pk,2,-1)
    _, xi0m1 = ξ(Pk,0,-1)
    _, xi3m2 = ξ(Pk,3,-2)
    _, xi1m2 = ξ(Pk,1,-2)
    _, xi0m2 = ξ(Pk,0,-2)
    _, xi40 = ξ(Pk,4,0)
    #end

    #@time begin
    M[0,0,0] = @. ((40*b1^2)/147. + (10*b1*b2)/21. + b2^2/4. - (20*b1^2*f)/63. - (b1*b2*f)/3. + (7*b1^2*f^2)/45. -
        (8*b1*bK²)/63. + (4*bK²^2)/45. - (32*b1*bKKpara)/189. -
        (2*b2*bKKpara)/9. + (4*b1*f*bKKpara)/15. - (8*bK²*bKKpara)/135. + (98*bKKpara^2)/675. - (20*b1*f*bδη)/63. -
        (b2*f*bδη)/3. +
        (14*b1*f^2*bδη)/45. + (4*f*bKKpara*bδη)/15. + (7*f^2*bδη^2)/45. - (16*b1*f*bη)/147. - (2*b2*f*bη)/21. +
        (16*b1*f^2*bη)/105. + (b2*f^2*bη)/9. - (2*b1*f^3*bη)/15. + (8*f*bK²*bη)/315. + (16*f^2*bK²*bη)/675. +
        (32*f*bKKpara*bη)/945. - (292*f^2*bKKpara*bη)/2025. + (4*f^2*bδη*bη)/63. - (2*f^3*bδη*bη)/15. +
        (24*f^2*bη^2)/1225. - (4*f^3*bη^2)/225. + (f^4*bη^2)/25. + (4*b1*f^2*bη²)/45. + (b2*f^2*bη²)/9. -
        (2*b1*f^3*bη²)/15. + (16*f^2*bK²*bη²)/675. - (292*f^2*bKKpara*bη²)/2025. - (2*f^3*bδη*bη²)/15. -
        (4*f^3*bη*bη²)/225. + (2*f^4*bη*bη²)/25. + (f^4*bη²^2)/25. +
        ((11700*b1 + 4725*b2 + 2730*b1*f - 10080*bK² + 10472*bKKpara + 2730*f*bδη - 5220*f*bη - 4998*f^2*bη - 4998*f^2*bη²)*
        bΠ2para)/99225. + (2717*bΠ2para^2)/33075.)*xi00^2 +
        ((-400*b1^2)/1029. - (10*b1*b2)/21. + (20*b1^2*f)/63. + (2*b1^2*f^2)/45. + (200*b1*bK²)/441. + (2*b2*bK²)/3. -
        (4*b1*f*bK²)/9. + (8*bK²^2)/63. +
        (80*b1*bKKpara)/1323. - (2*b2*bKKpara)/9. + (44*b1*f*bKKpara)/135. - (8*bK²*bKKpara)/21. + (254*bKKpara^2)/945. +
        (20*b1*f*bδη)/63. + (4*b1*f^2*bδη)/45. - (4*f*bK²*bδη)/9. + (44*f*bKKpara*bδη)/135. + (2*f^2*bδη^2)/45. +
        (160*b1*f*bη)/1029. + (2*b2*f*bη)/21. - (16*b1*f^2*bη)/147. + (4*b2*f^2*bη)/45. - (16*b1*f^3*bη)/105. -
        (40*f*bK²*bη)/441. + (172*f^2*bK²*bη)/945. - (16*f*bKKpara*bη)/1323. - (5188*f^2*bKKpara*bη)/19845. -
        (4*f^2*bδη*bη)/63. - (16*f^3*bδη*bη)/105. - (48*f^2*bη^2)/1715. + (4*f^3*bη^2)/441. + (16*f^4*bη^2)/245. -
        (20*b1*f^2*bη²)/441. + (4*b2*f^2*bη²)/45. - (16*b1*f^3*bη²)/105. + (172*f^2*bK²*bη²)/945. -
        (5188*f^2*bKKpara*bη²)/19845. - (16*f^3*bδη*bη²)/105. + (4*f^3*bη*bη²)/441. + (32*f^4*bη*bη²)/245. +
        (16*f^4*bη²^2)/245. + (4*(-10575*b1 - 13230*b2 + 12936*b1*f - 2835*bK² + 13538*bKKpara + 12936*f*bδη + 3555*f*bη -
        6762*f^2*bη - 6762*f^2*bη²)*bΠ2para)/138915. + (1286*bΠ2para^2)/46305.)*xi20^2 +
        ((8*(1125*b1^2 - 3150*b1*bK² + 2205*bK²^2 + 1050*b1*bKKpara - 1470*bK²*bKKpara + 294*bKKpara^2 - 450*b1*f*bη -
        420*b1*f^2*bη + 630*f*bK²*bη + 588*f^2*bK²*bη - 210*f*bKKpara*bη - 252*f^2*bKKpara*bη +
        81*f^2*bη^2 + 84*f^3*bη^2 + 56*f^4*bη^2 - 420*b1*f^2*bη² + 588*f^2*bK²*bη² -
        252*f^2*bKKpara*bη² + 84*f^3*bη*bη² + 112*f^4*bη*bη² + 56*f^4*bη²^2))/77175. +
        (16*(900*b1 - 1260*bK² + 469*bKKpara - 240*f*bη - 196*f^2*bη - 196*f^2*bη²)*bΠ2para)/77175. +
        (6952*bΠ2para^2)/77175.)*xi40^2

    M[2,0,0] = @. ((f^2*(15*b1^2 - 10*b1*f + 3*f^2))/270. - ((5*b1 - 3*f)*f^3*bη)/135. + (f^4*bη^2)/90.)*xi0m2*xi00 +
        ((2*b1^2)/7. + (b1*b2)/2. - (b1^2*f)/3. - (2*b1*f^2)/63. - (b1^2*f^2)/18. - (b2*f^2)/18. + (2*b1*f^3)/15. -
        f^4/50. + (4*b1*bK²)/15. - (4*f^2*bK²)/135. -
        (14*b1*bKKpara)/45. + (164*f^2*bKKpara)/2025. - (b1*f*bδη)/3. + (f^3*bδη)/15. - (16*b1*f*bη)/105. - (b2*f*bη)/6. +
        (58*b1*f^2*bη)/225. + (2*f^3*bη)/175. + (b1*f^3*bη)/15. - (2*f^4*bη)/25. - (4*f*bK²*bη)/45. +
        (14*f*bKKpara*bη)/135. + (f^2*bδη*bη)/9. + (6*f^2*bη^2)/175. - (11*f^3*bη^2)/225. - (f^4*bη^2)/50. +
        (11*b1*f^2*bη²)/75. - (f^4*bη²)/25. - (11*f^3*bη*bη²)/225. - ((165*b1 - 47*f^2 - 15*f*bη)*bΠ2para)/1575.)*-xi1m1^2 +
        ((-4*(7*b1 - 3*f)*f^3)/945. - (4*(7*b1 - 6*f)*f^3*bη)/945. + (4*f^4*bη^2)/315.)*xi2m2*xi20 +
        ((-2*(225*b1^2 - 25*b1*f^2 + 3*f^4 - 315*b1*bK² + 35*f^2*bK² + 105*b1*bKKpara - 14*f^2*bKKpara - 120*b1*f*bη -
        42*b1*f^2*bη +
        9*f^3*bη + 12*f^4*bη + 105*f*bK²*bη - 35*f*bKKpara*bη + 27*f^2*bη^2 + 14*f^3*bη^2 +
        3*f^4*bη^2 - 42*b1*f^2*bη² + 6*f^4*bη² + 14*f^3*bη*bη²))/1575. -
        (2*(180*b1 - 29*f^2 - 80*f*bη)*bΠ2para)/1575.)*xi3m1^2

    M[4,0,0] = @. ((45*b1^2 - 10*b1*f^2 + 3*f^4)/540. - (f*(5*b1 - f^2)*bη)/90. + (f^2*bη^2)/60.)*xi0m2^2 +
        ((225*b1^2 - 50*b1*f^2 + 6*f^4)/1350. - (f*(5*b1 - f^2)*bη)/45. + (f^2*bη^2)/30.)*xi2m2^2

    M[0,2,0] = @. ((4*f*bη*(-12600*b1 - 11025*b2 + 7350*b1*f + 2940*bK² + 3920*bKKpara + 7350*f*bδη + 3240*f*bη -
        2058*f^2*bη -
        2058*f^2*bη²))/231525. - (4*(-12600*b1 - 11025*b2 + 7350*b1*f + 2940*bK² + 3920*bKKpara + 7350*f*bδη + 5598*f*bη -
        2058*f^2*bη - 2058*f^2*bη²)*bΠ2para)/138915. + (1048*bΠ2para^2)/9261.)*xi00^2 +
        ((-4*f*bη*(-5040*b1 - 3087*b2 + 2058*b1*f + 2940*bK² + 392*bKKpara + 2058*f*bδη + 1296*f*bη - 294*f^2*bη -
        294*f^2*bη²))/64827. + (4*(-25200*b1 - 15435*b2 + 10290*b1*f + 14700*bK² + 1960*bKKpara + 10290*f*bδη +
        13842*f*bη -
        1470*f^2*bη - 1470*f^2*bη²)*bΠ2para)/194481. - (16360*bΠ2para^2)/64827.)*xi20^2 +
        ((-32*f*bη*(525*b1 - 735*bK² + 245*bKKpara - 135*f*bη - 98*f^2*bη - 98*f^2*bη²))/180075. +
        (32*(525*b1 - 735*bK² + 245*bKKpara - 417*f*bη - 98*f^2*bη - 98*f^2*bη²)*bΠ2para)/108045. +
        (3008*bΠ2para^2)/21609.)*xi40^2

    M[2,2,0] = @. ((f^2*(15*b1^2 - 10*b1*f + 3*f^2))/135. - (2*(5*b1 - 3*f)*f^3*bη)/135. + (f^4*bη^2)/45.)*xi0m2*xi00 +
        (-(f*(6300*b1*f + 11025*b1^2*f + 11025*b2*f - 26460*b1*f^2 + 3969*f^3 + 5880*f*bK² - 16072*f*bKKpara -
        13230*f^2*bδη +
        30240*b1*bη +
        33075*b2*bη - 22050*b1*f*bη - 3240*f^2*bη - 13230*b1*f^2*bη + 15876*f^3*bη + 17640*bK²*bη -
        20580*bKKpara*bη - 22050*f*bδη*bη - 9720*f*bη^2 + 9702*f^2*bη^2 + 3969*f^3*bη^2 +
        7938*f^3*bη² + 9702*f^2*bη*bη²))/99225. + (2*(1050*b1 + 239*f^2 - 165*f*bη)*bΠ2para)/11025.)*xi1m1^2 +
        ((-8*(7*b1 - 3*f)*f^3)/945. - (8*(7*b1 - 6*f)*f^3*bη)/945. + (8*f^4*bη^2)/315.)*xi2m2*xi20 +
        ((4*f*(175*b1*f - 21*f^3 - 245*f*bK² + 98*f*bKKpara + 840*b1*bη - 90*f^2*bη - 84*f^3*bη - 735*bK²*bη +
        245*bKKpara*bη - 270*f*bη^2 - 98*f^2*bη^2 - 21*f^3*bη^2 - 42*f^3*bη² - 98*f^2*bη*bη²))/11025. -
        (4*(525*b1 - 248*f^2 - 695*f*bη)*bΠ2para)/11025.)*xi3m1^2

    M[4,2,0] = @. (-(f^2*(7*b1 - 3*f^2))/189. - (f*(7*b1 - 2*f^2)*bη)/63. + (f^2*bη^2)/21.)*xi0m2^2 +
        ((-2*f^2*(35*b1 - 6*f^2))/945. - (2*f*(7*b1 - 2*f^2)*bη)/63. + (2*f^2*bη^2)/21.)*xi2m2^2

    M[0,4,0] = @. ((192*f^2*bη^2)/8575. - (128*f*bη*bΠ2para)/1715. + (64*bΠ2para^2)/1029.)*xi00^2 +
        ((-384*f^2*bη^2)/12005. + (256*f*bη*bΠ2para)/2401. - (640*bΠ2para^2)/7203.)*xi20^2 +
        ((576*f^2*bη^2)/60025. - (384*f*bη*bΠ2para)/12005. + (64*bΠ2para^2)/2401.)*xi40^2

    M[2,4,0] = @. ((16*f^2*bη*(f + 3*bη))/1225. - (16*f*(f + 3*bη)*bΠ2para)/735.)*xi1m1^2 +
        ((-16*f^2*bη*(f + 3*bη))/1225. + (16*f*(f + 3*bη)*bΠ2para)/735.)*xi3m1^2

    M[4,4,0] = @. ((2*f^4)/315. + (4*f^3*bη)/315. + (2*f^2*bη^2)/105.)*xi0m2^2 + ((8*f^4)/1575. + (8*f^3*bη)/315. +
        (4*f^2*bη^2)/105.)*xi2m2^2

    M[1,0,1] = @. xi10*(((-20*b1^2*f)/63. - (b1*b2*f)/3. + (4*b1*f^2)/45. + (14*b1^2*f^2)/45. + (b2*f^2)/9. -
        (2*b1*f^3)/15. + (16*f^2*bK²)/675. + (4*b1*f*bKKpara)/15. -
        (292*f^2*bKKpara)/2025. + (14*b1*f^2*bδη)/45. - (2*f^3*bδη)/15. + (64*b1*f^2*bη)/315. + (b2*f^2*bη)/9. -
        (4*f^3*bη)/125. - (4*b1*f^3*bη)/15. + (2*f^4*bη)/25. + (16*f^2*bK²*bη)/675. - (292*f^2*bKKpara*bη)/2025. -
        (2*f^3*bδη*bη)/15. - (4*f^3*bη^2)/125. + (2*f^4*bη^2)/25. - (2*b1*f^3*bη²)/15. + (2*f^4*bη²)/25. +
        (2*f^4*bη*bη²)/25. - (2*f*(15*b1 + 7*f + 7*f*bη)*bΠ2para)/525.)*xi0m1 +
        ((-8*b1^2*f)/63. - (8*b1*f^2)/315. - (4*b1^2*f^2)/45. - (4*b2*f^2)/45. + (68*b1*f^3)/525. + (8*b1*f*bK²)/45. -
        (56*f^2*bK²)/675. - (148*b1*f*bKKpara)/675. +
        (2564*f^2*bKKpara)/14175. - (4*b1*f^2*bδη)/45. + (68*f^3*bδη)/525. + (32*b1*f^2*bη)/1575. - (4*b2*f^2*bη)/45. +
        (8*f^3*bη)/875. + (124*b1*f^3*bη)/525. - (16*f^4*bη)/175. - (56*f^2*bK²*bη)/675. +
        (2564*f^2*bKKpara*bη)/14175. + (68*f^3*bδη*bη)/525. + (8*f^3*bη^2)/875. - (16*f^4*bη^2)/175. +
        (8*b1*f^3*bη²)/75. - (16*f^4*bη²)/175. - (16*f^4*bη*bη²)/175. - (124*b1*f*bΠ2para)/525. + (52*f^2*bΠ2para)/525. +
        (52*f^2*bη*bΠ2para)/525.)*xi2m1) + xi30*
        ((4*b1^2*f*xi2m1)/21. - (52*b1*f^2*xi2m1)/735. - (4*b1*f^3*xi2m1)/175. - (4*b1*f*bK²*xi2m1)/15. +
        (52*f^2*bK²*xi2m1)/525. + (8*b1*f*bKKpara*xi2m1)/75. - (296*f^2*bKKpara*xi2m1)/3675. -
        (4*f^3*bδη*xi2m1)/175. - (512*b1*f^2*bη*xi2m1)/3675. + (156*f^3*bη*xi2m1)/6125. -
        (12*b1*f^3*bη*xi2m1)/175. + (48*f^4*bη*xi2m1)/1225. + (52*f^2*bK²*bη*xi2m1)/525. -
        (296*f^2*bKKpara*bη*xi2m1)/3675. - (4*f^3*bδη*bη*xi2m1)/175. + (156*f^3*bη^2*xi2m1)/6125. +
        (48*f^4*bη^2*xi2m1)/1225. - (8*b1*f^3*bη²*xi2m1)/175. + (48*f^4*bη²*xi2m1)/1225. +
        (48*f^4*bη*bη²*xi2m1)/1225. + (116*b1*f*bΠ2para*xi2m1)/525. - (132*f^2*bΠ2para*xi2m1)/1225. -
        (132*f^2*bη*bΠ2para*xi2m1)/1225. + (32*f^2*(1 + bη)*
        (75*b1 - 105*bK² + 45*bKKpara - 27*f*bη - 20*f^2*bη - 20*f^2*bη² + 90*bΠ2para)*xi4m1)/55125.)

    M[3,0,1] = @. ((f*(25*b1^2 - 15*b1*f - 15*b1*f^2 + 9*f^3))/225. - (f^2*(10*b1 - 3*f - 3*f^2)*bη)/75. +
        (f^3*bη^2)/25.)*xi0m1*xi1m2 -
        (2*f*(125*b1^2 - 45*b1*f - 30*b1*f^2 + 18*f^3 - 120*b1*f*bη + 27*f^2*bη + 18*f^3*bη +
        27*f^2*bη^2)*xi1m2*xi2m1)/1125. -
        (4*f^2*(1 + bη)*(35*b1 - 9*f^2 - 21*f*bη)*xi2m1*xi3m2)/2625.

    M[1,2,1] = @. xi10*(((-40*b1^2*f)/63. - (2*b1*b2*f)/3. + (8*b1*f^2)/45. + (28*b1^2*f^2)/45. + (2*b2*f^2)/9. -
        (4*b1*f^3)/15. + (32*f^2*bK²)/675. + (8*b1*f*bKKpara)/15. -
        (584*f^2*bKKpara)/2025. + (28*b1*f^2*bδη)/45. - (4*f^3*bδη)/15. + (1112*b1*f^2*bη)/2205. + (2*b2*f^2*bη)/9. -
        (16*f^3*bη)/175. - (8*b1*f^3*bη)/15. + (4*f^4*bη)/25. + (32*f^2*bK²*bη)/675. - (584*f^2*bKKpara*bη)/2025. -
        (4*f^3*bδη*bη)/15. - (16*f^3*bη^2)/175. + (4*f^4*bη^2)/25. - (4*b1*f^3*bη²)/15. + (4*f^4*bη²)/25. +
        (4*f^4*bη*bη²)/25. - (4*f*(255*b1 + 7*f + 7*f*bη)*bΠ2para)/3675.)*xi0m1 +
        ((-16*b1^2*f)/63. - (16*b1*f^2)/315. - (8*b1^2*f^2)/45. - (8*b2*f^2)/45. + (136*b1*f^3)/525. + (16*b1*f*bK²)/45. -
        (112*f^2*bK²)/675. -
        (296*b1*f*bKKpara)/675. + (5128*f^2*bKKpara)/14175. - (8*b1*f^2*bδη)/45. + (136*f^3*bδη)/525. +
        (176*b1*f^2*bη)/2205. -
        (8*b2*f^2*bη)/45. + (32*f^3*bη)/1225. + (248*b1*f^3*bη)/525. - (32*f^4*bη)/175. - (112*f^2*bK²*bη)/675. +
        (5128*f^2*bKKpara*bη)/14175. + (136*f^3*bδη*bη)/525. + (32*f^3*bη^2)/1225. - (32*f^4*bη^2)/175. +
        (16*b1*f^3*bη²)/75. - (32*f^4*bη²)/175. - (32*f^4*bη*bη²)/175. - (1976*b1*f*bΠ2para)/3675. +
        (136*f^2*bΠ2para)/735. +
        (136*f^2*bη*bΠ2para)/735.)*xi2m1) + xi30*
        ((8*b1^2*f*xi2m1)/21. - (104*b1*f^2*xi2m1)/735. - (8*b1*f^3*xi2m1)/175. - (8*b1*f*bK²*xi2m1)/15. +
        (104*f^2*bK²*xi2m1)/525. + (16*b1*f*bKKpara*xi2m1)/75. - (592*f^2*bKKpara*xi2m1)/3675. -
        (8*f^3*bδη*xi2m1)/175. - (248*b1*f^2*bη*xi2m1)/735. + (624*f^3*bη*xi2m1)/8575. -
        (24*b1*f^3*bη*xi2m1)/175. + (96*f^4*bη*xi2m1)/1225. + (104*f^2*bK²*bη*xi2m1)/525. -
        (592*f^2*bKKpara*bη*xi2m1)/3675. - (8*f^3*bδη*bη*xi2m1)/175. + (624*f^3*bη^2*xi2m1)/8575. +
        (96*f^4*bη^2*xi2m1)/1225. - (16*b1*f^3*bη²*xi2m1)/175. + (96*f^4*bη²*xi2m1)/1225. +
        (96*f^4*bη*bη²*xi2m1)/1225. + (1984*b1*f*bΠ2para*xi2m1)/3675. - (432*f^2*bΠ2para*xi2m1)/1715. -
        (432*f^2*bη*bΠ2para*xi2m1)/1715. + (64*f^2*(1 + bη)*
        (105*b1 - 147*bK² + 63*bKKpara - 54*f*bη - 28*f^2*bη - 28*f^2*bη² + 153*bΠ2para)*xi4m1)/77175.)

    M[3,2,1] = @. ((2*f*(35*b1^2 - 21*b1*f - 30*b1*f^2 + 18*f^3))/315. - (2*f^2*(17*b1 - 6*f - 6*f^2)*bη)/105. +
        (4*f^3*bη^2)/35.)*xi0m1*xi1m2 -
        (4*f*(175*b1^2 - 63*b1*f - 60*b1*f^2 + 36*f^3 - 213*b1*f*bη + 54*f^2*bη + 36*f^3*bη +
        54*f^2*bη^2)*xi1m2*xi2m1)/1575. -
        (8*f^2*(1 + bη)*(49*b1 - 18*f^2 - 42*f*bη)*xi2m1*xi3m2)/3675.

    M[1,4,1] = @. xi10*(((32*f^2*bη*(25*b1 - 7*f - 7*f*bη))/6125. - (32*f*(25*b1 - 7*f - 7*f*bη)*bΠ2para)/3675.)*xi0m1 +
        (64*f*(5*b1 + f + f*bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/18375.) +
        xi30*((-32*f*(35*b1 - 13*f - 13*f*bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/42875. -
        (256*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/42875.)

    M[3,4,1] = @. ((-8*(5*b1 - 3*f)*f^3)/525. - (8*f^2*(5*b1 - 3*f - 3*f^2)*bη)/525. + (8*f^3*bη^2)/175.)*xi0m1*xi1m2 +
        (16*f^2*(10*b1*f - 6*f^2 + 25*b1*bη - 9*f*bη - 6*f^2*bη - 9*f*bη^2)*xi1m2*xi2m1)/2625. +
        (32*f^3*(3*f + 7*bη + 3*f*bη + 7*bη^2)*xi2m1*xi3m2)/6125.

    M[0,0,2] = @. ((32*f*bη*(147*b1*f + 14*bKKpara + 147*f*bδη + 6*f^2*bη + 6*f^2*bη²))/324135. -
        (32*(735*b1*f + 70*bKKpara + 735*f*bδη + 252*f*bη + 30*f^2*bη + 30*f^2*bη²)*bΠ2para)/972405. +
        (128*bΠ2para^2)/9261.)*xi20^2 +
        ((-32*f*bη*(-7*bKKpara + 4*f^2*bη + 4*f^2*bη²))/36015. -
        (32*(35*bKKpara - 21*f*bη - 20*f^2*bη - 20*f^2*bη²)*bΠ2para)/108045. - (32*bΠ2para^2)/3087.)*xi40^2 +
        xi20*(((32*f*bη*(42*b1*f + 49*bKKpara + 42*f*bδη - 24*f^2*bη - 24*f^2*bη²))/33075. -
        (32*(210*b1*f + 245*bKKpara + 210*f*bδη - 63*f*bη - 120*f^2*bη - 120*f^2*bη²)*bΠ2para)/99225. -
        (32*bΠ2para^2)/945.)*xi00 +
        ((-32*f*bη*(147*b1*f + 224*bKKpara + 147*f*bδη - 114*f^2*bη - 114*f^2*bη²))/180075. +
        (32*(735*b1*f + 1120*bKKpara + 735*f*bδη - 378*f*bη - 570*f^2*bη - 570*f^2*bη²)*bΠ2para)/540225. +
        (64*bΠ2para^2)/1715.)*xi40)

    M[2,0,2] = @. ((2*f*(10080*b1*f + 11025*b1^2*f + 11025*b2*f - 26460*b1*f^2 + 3969*f^3 + 588*f*bK² - 13426*f*bKKpara -
        13230*f^2*bδη - 8820*b1*f*bη -
        5184*f^2*bη - 13230*b1*f^2*bη + 15876*f^3*bη - 13818*bKKpara*bη - 8820*f*bδη*bη +
        7056*f^2*bη^2 + 3969*f^3*bη^2 + 7938*f^3*bη² + 7056*f^2*bη*bη²))/496125. +
        (4*f*(88*f - 441*bη)*bΠ2para)/55125.)*xi1m1^2 +
        ((-4*(7*b1 - 3*f)*f^3)/945. - (4*(7*b1 - 6*f)*f^3*bη)/945. + (4*f^4*bη^2)/315.)*xi00*xi2m2 +
        ((4*f*(420*b1*f + 630*b1*f^2 - 189*f^3 - 588*f*bK² + 756*f*bKKpara + 315*f^2*bδη + 735*b1*f*bη - 216*f^2*bη +
        315*b1*f^2*bη -
        756*f^3*bη + 1078*bKKpara*bη + 735*f*bδη*bη - 546*f^2*bη^2 - 189*f^3*bη^2 - 378*f^3*bη² -
        546*f^2*bη*bη²))/55125. + (8*f*(151*f + 98*bη)*bΠ2para)/18375.)*xi1m1*xi3m1 +
        ((-16*f*(35*b1*f - 7*f^3 - 49*f*bK² + 28*f*bKKpara - 18*f^2*bη - 28*f^3*bη + 49*bKKpara*bη - 28*f^2*bη^2 -
        7*f^3*bη^2 - 14*f^3*bη² - 28*f^2*bη*bη²))/55125. - (16*f*(58*f + 49*bη)*bΠ2para)/55125.)*xi3m1^2 +
        xi20*(((2*f^2*(15*b1^2 - 10*b1*f + 3*f^2))/675. - (4*(5*b1 - 3*f)*f^3*bη)/675. + (2*f^4*bη^2)/225.)*xi0m2 +
        ((16*(7*b1 - 3*f)*f^3)/6615. + (16*(7*b1 - 6*f)*f^3*bη)/6615. - (16*f^4*bη^2)/2205.)*xi2m2 +
        ((32*f^4)/11025. + (64*f^4*bη)/11025. + (32*f^4*bη^2)/11025.)*xi4m2)

    M[4,0,2] = @. ((-4*f^2*(7*b1 - 3*f^2))/945. + (8*f^3*bη)/315.)*xi0m2*xi2m2 + ((2*f^2*(49*b1 - 12*f^2))/6615. -
        (4*f^3*bη)/315.)*xi2m2^2

    M[0,2,2] = @. ((-80*b1^2*f)/441. - (4*b1^2*f^2)/63. + (16*b1*f*bK²)/63. - (160*b1*bKKpara)/9261. + (2*b2*bKKpara)/9. -
        (508*b1*f*bKKpara)/1323. +
        (424*bK²*bKKpara)/1323. - (2780*bKKpara^2)/9261. - (80*b1*f*bδη)/441. - (8*b1*f^2*bδη)/63. + (16*f*bK²*bδη)/63. -
        (508*f*bKKpara*bδη)/1323. - (4*f^2*bδη^2)/63. + (1072*b1*f^2*bη)/21609. - (8*b2*f^2*bη)/63. + (32*b1*f^3*bη)/147. -
        (1472*f^2*bK²*bη)/9261. + (352*f*bKKpara*bη)/64827. + (9404*f^2*bKKpara*bη)/27783. +
        (176*f^2*bδη*bη)/3087. + (32*f^3*bδη*bη)/147. + (352*f^3*bη^2)/151263. - (32*f^4*bη^2)/343. -
        (160*b1*f^2*bη²)/21609. - (8*b2*f^2*bη²)/63. + (32*b1*f^3*bη²)/147. - (1472*f^2*bK²*bη²)/9261. +
        (9404*f^2*bKKpara*bη²)/27783. + (32*f^3*bδη*bη²)/147. + (352*f^3*bη*bη²)/151263. -
        (64*f^4*bη*bη²)/343. - (32*f^4*bη²^2)/343. -
        (2*(-70560*b1 - 151263*b2 + 268422*b1*f - 102900*bK² + 274190*bKKpara + 268422*f*bδη + 22176*f*bη - 148182*f^2*bη -
        148182*f^2*bη²)*bΠ2para)/1.361367e6 - (4376*bΠ2para^2)/64827.)*xi20^2 +
        ((-16*(8085*b1*bKKpara - 11319*bK²*bKKpara + 3234*bKKpara^2 - 4620*b1*f^2*bη + 6468*f^2*bK²*bη -
        2541*f*bKKpara*bη - 3234*f^2*bKKpara*bη + 1452*f^3*bη^2 + 784*f^4*bη^2 - 4620*b1*f^2*bη² +
        6468*f^2*bK²*bη² - 3234*f^2*bKKpara*bη² + 1452*f^3*bη*bη² + 1568*f^4*bη*bη² +
        784*f^4*bη²^2))/1.663893e6 - (16*(735*b1 - 1029*bK² + 973*bKKpara - 231*f*bη - 514*f^2*bη - 514*f^2*bη²)*
        bΠ2para)/151263. - (1552*bΠ2para^2)/21609.)*xi40^2 +
        xi20*(((-32*b1^2*f)/63. - (2*b1*b2*f)/3. + (32*b1^2*f^2)/63. - (8*b1*f*bK²)/45. - (16*b1*bKKpara)/27. -
        (8*b2*bKKpara)/9. +
        (1112*b1*f*bKKpara)/945. - (16*bK²*bKKpara)/45. + (284*bKKpara^2)/405. - (32*b1*f*bδη)/63. - (2*b2*f*bδη)/3. +
        (64*b1*f^2*bδη)/63. - (8*f*bK²*bδη)/45. + (1112*f*bKKpara*bδη)/945. + (32*f^2*bδη^2)/63. +
        (992*b1*f^2*bη)/2205. + (4*b2*f^2*bη)/9. - (68*b1*f^3*bη)/105. + (176*f^2*bK²*bη)/945. +
        (176*f*bKKpara*bη)/945. - (64*f^2*bKKpara*bη)/81. + (352*f^2*bδη*bη)/2205. - (68*f^3*bδη*bη)/105. -
        (1408*f^3*bη^2)/15435. + (8*f^4*bη^2)/35. + (128*b1*f^2*bη²)/441. + (4*b2*f^2*bη²)/9. - (68*b1*f^3*bη²)/105. +
        (176*f^2*bK²*bη²)/945. - (64*f^2*bKKpara*bη²)/81. - (68*f^3*bδη*bη²)/105. -
        (1408*f^3*bη*bη²)/15435. + (16*f^4*bη*bη²)/35. + (8*f^4*bη²^2)/35. +
        (4*(-8820*b1 - 15435*b2 + 13251*b1*f - 8232*bK² + 17542*bKKpara + 13251*f*bδη + 2772*f*bη - 10596*f^2*bη -
        10596*f^2*bη²)*
        bΠ2para)/138915. + (716*bΠ2para^2)/6615.)*xi00 +
        ((16*b1^2*f)/49. - (16*b1*f*bK²)/35. + (512*b1*bKKpara)/1029. + (48*b1*f*bKKpara)/245. - (512*bK²*bKKpara)/735. +
        (528*bKKpara^2)/1715. + (16*b1*f*bδη)/49. - (16*f*bK²*bδη)/35. + (48*f*bKKpara*bδη)/245. -
        (4272*b1*f^2*bη)/12005. - (64*b1*f^3*bη)/735. + (608*f^2*bK²*bη)/1715. - (5632*f*bKKpara*bη)/36015. -
        (928*f^2*bKKpara*bη)/3087. - (176*f^2*bδη*bη)/1715. - (64*f^3*bδη*bη)/735. + (6688*f^3*bη^2)/84035. +
        (128*f^4*bη^2)/1715. - (608*b1*f^2*bη²)/2401. - (64*b1*f^3*bη²)/735. + (608*f^2*bK²*bη²)/1715. -
        (928*f^2*bKKpara*bη²)/3087. - (64*f^3*bδη*bη²)/735. + (6688*f^3*bη*bη²)/84035. +
        (256*f^4*bη*bη²)/1715. + (128*f^4*bη²^2)/1715. +
        (32*(6615*b1 + 8673*b1*f - 9261*bK² + 17626*bKKpara + 8673*f*bδη - 2079*f*bη - 8868*f^2*bη - 8868*f^2*bη²)*bΠ2para)/
        756315. + (3888*bΠ2para^2)/12005.)*xi40)

    M[2,2,2] = @. ((26460*b1^2*f + 15840*b1*f^2 + 17325*b1^2*f^2 + 17325*b2*f^2 - 41580*b1*f^3 + 6237*f^4 + 924*f^2*bK² +
        41454*b1*bKKpara - 21098*f^2*bKKpara +
        26460*b1*f*bδη - 20790*f^3*bδη - 35028*b1*f^2*bη - 7776*f^3*bη - 20790*b1*f^3*bη + 24948*f^4*bη -
        21714*f*bKKpara*bη - 13860*f^2*bδη*bη + 11088*f^3*bη^2 + 6237*f^4*bη^2 - 21168*b1*f^2*bη² +
        12474*f^4*bη² + 11088*f^3*bη*bη²)/99225. + (2*(1323*b1 + 104*f^2 - 693*f*bη)*bΠ2para)/11025.)*xi1m1^2 +
        ((-22*(7*b1 - 3*f)*f^3)/1323. - (22*(7*b1 - 6*f)*f^3*bη)/1323. + (22*f^4*bη^2)/441.)*xi00*xi2m2 +
        ((-2*(2205*b1^2*f - 660*b1*f^2 - 990*b1*f^3 + 297*f^4 + 924*f^2*bK² + 3234*b1*bKKpara - 1188*f^2*bKKpara +
        2205*b1*f*bδη -
        495*f^3*bδη -
        2793*b1*f^2*bη + 324*f^3*bη - 495*b1*f^3*bη + 1188*f^4*bη - 1694*f*bKKpara*bη -
        1155*f^2*bδη*bη + 858*f^3*bη^2 + 297*f^4*bη^2 - 1638*b1*f^2*bη² + 594*f^4*bη² +
        858*f^3*bη*bη²))/11025. - (4*(294*b1 - 233*f^2 - 154*f*bη)*bΠ2para)/3675.)*xi1m1*xi3m1 +
        ((-8*(55*b1*f^2 - 11*f^4 - 77*f^2*bK² - 147*b1*bKKpara + 44*f^2*bKKpara + 84*b1*f^2*bη - 27*f^3*bη - 44*f^4*bη +
        77*f*bKKpara*bη - 44*f^3*bη^2 - 11*f^4*bη^2 + 84*b1*f^2*bη² - 22*f^4*bη² - 44*f^3*bη*bη²))/
        11025. + (8*(147*b1 - 89*f^2 - 77*f*bη)*bΠ2para)/11025.)*xi3m1^2 +
        xi20*(((11*f^2*(15*b1^2 - 10*b1*f + 3*f^2))/945. - (22*(5*b1 - 3*f)*f^3*bη)/945. + (11*f^4*bη^2)/315.)*xi0m2 +
        ((88*(7*b1 - 3*f)*f^3)/9261. + (88*(7*b1 - 6*f)*f^3*bη)/9261. - (88*f^4*bη^2)/3087.)*xi2m2 +
        ((176*f^4)/15435. + (352*f^4*bη)/15435. + (176*f^4*bη^2)/15435.)*xi4m2)

    M[4,2,2] = @. (-(f^2*(22*b1 - 9*f^2))/189. + (2*f^3*bη)/21.)*xi0m2*xi2m2 + ((f^2*(77*b1 - 18*f^2))/1323. -
        (f^3*bη)/21.)*xi2m2^2

    M[0,4,2] = @. ((64*f*bη*(147*b1*f + 14*bKKpara + 147*f*bδη + 6*f^2*bη + 6*f^2*bη²))/252105. -
        (64*(735*b1*f + 70*bKKpara + 735*f*bδη + 252*f*bη + 30*f^2*bη + 30*f^2*bη²)*bΠ2para)/756315. +
        (256*bΠ2para^2)/7203.)*xi20^2 +
        ((-192*f*bη*(-7*bKKpara + 4*f^2*bη + 4*f^2*bη²))/84035. -
        (64*(35*bKKpara - 21*f*bη - 20*f^2*bη - 20*f^2*bη²)*bΠ2para)/84035. - (64*bΠ2para^2)/2401.)*xi40^2 +
        xi20*(((64*f*bη*(42*b1*f + 49*bKKpara + 42*f*bδη - 24*f^2*bη - 24*f^2*bη²))/25725. -
        (64*(210*b1*f + 245*bKKpara + 210*f*bδη - 63*f*bη - 120*f^2*bη - 120*f^2*bη²)*bΠ2para)/77175. -
        (64*bΠ2para^2)/735.)*xi00 +
        ((-192*f*bη*(147*b1*f + 224*bKKpara + 147*f*bδη - 114*f^2*bη - 114*f^2*bη²))/420175. +
        (64*(735*b1*f + 1120*bKKpara + 735*f*bδη - 378*f*bη - 570*f^2*bη - 570*f^2*bη²)*bΠ2para)/420175. +
        (1152*bΠ2para^2)/12005.)*xi40)

    M[2,4,2] = @. ((4*f*(15840*b1*f + 17325*b1^2*f + 17325*b2*f - 41580*b1*f^2 + 6237*f^3 + 924*f*bK² - 21098*f*bKKpara -
        20790*f^2*bδη - 13860*b1*f*bη -
        9792*f^2*bη - 20790*b1*f^2*bη + 24948*f^3*bη - 21714*bKKpara*bη - 13860*f*bδη*bη +
        11088*f^2*bη^2 + 6237*f^3*bη^2 + 12474*f^3*bη² + 11088*f^2*bη*bη²))/606375. +
        (8*f*(872*f - 2079*bη)*bΠ2para)/202125.)*xi1m1^2 +
        ((-8*(7*b1 - 3*f)*f^3)/735. - (8*(7*b1 - 6*f)*f^3*bη)/735. + (8*f^4*bη^2)/245.)*xi00*xi2m2 +
        ((8*f*(660*b1*f + 990*b1*f^2 - 297*f^3 - 924*f*bK² + 1188*f*bKKpara + 495*f^2*bδη + 1155*b1*f*bη - 408*f^2*bη +
        495*b1*f^2*bη -
        1188*f^3*bη + 1694*bKKpara*bη + 1155*f*bδη*bη - 858*f^2*bη^2 - 297*f^3*bη^2 - 594*f^3*bη² -
        858*f^2*bη*bη²))/67375. + (16*f*(769*f + 462*bη)*bΠ2para)/67375.)*xi1m1*xi3m1 +
        ((-32*f*(55*b1*f - 11*f^3 - 77*f*bK² + 44*f*bKKpara - 34*f^2*bη - 44*f^3*bη + 77*bKKpara*bη - 44*f^2*bη^2 -
        11*f^3*bη^2 - 22*f^3*bη² - 44*f^2*bη*bη²))/67375. - (32*f*(302*f + 231*bη)*bΠ2para)/202125.)*xi3m1^2 +
        xi20*(((4*f^2*(15*b1^2 - 10*b1*f + 3*f^2))/525. - (8*(5*b1 - 3*f)*f^3*bη)/525. + (4*f^4*bη^2)/175.)*xi0m2 +
        ((32*(7*b1 - 3*f)*f^3)/5145. + (32*(7*b1 - 6*f)*f^3*bη)/5145. - (32*f^4*bη^2)/1715.)*xi2m2 +
        ((64*f^4)/8575. + (128*f^4*bη)/8575. + (64*f^4*bη^2)/8575.)*xi4m2)

    M[4,4,2] = @. ((-8*f^2*(33*b1 - 17*f^2))/3465. + (272*f^3*bη)/3465.)*xi0m2*xi2m2 + ((4*f^2*(231*b1 - 68*f^2))/24255. -
        (136*f^3*bη)/3465.)*xi2m2^2

    M[2,6,2] = @. ((-256*f^3*bη)/13475. + (256*f^2*bΠ2para)/8085.)*xi1m1^2 +
        ((-192*f^3*bη)/13475. + (64*f^2*bΠ2para)/2695.)*xi1m1*xi3m1 + ((64*f^3*bη)/13475. - (64*f^2*bΠ2para)/8085.)*xi3m1^2

    M[4,6,2] = @. ((8*f^4)/693. + (16*f^3*bη)/693.)*xi0m2*xi2m2 + ((-16*f^4)/4851. - (8*f^3*bη)/693.)*xi2m2^2

    M[1,0,3] = @. xi10*((-176*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/42875. -
        (64*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/42875.) +
        xi30*(((16*f^3*bη*(1 + bη))/6125. - (16*f^2*(1 + bη)*bΠ2para)/3675.)*xi0m1 -
        (128*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/128625. + (32*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/42875.)

    M[3,0,3] = @. (-4*f^3*(5*b1 - 3*f - 2*bη - 3*f*bη - 2*bη^2)*xi1m2*xi2m1)/875. +
        ((4*f^4)/525. + (4*f^3*(1 + f)*bη)/525. + (4*f^3*bη^2)/525.)*xi0m1*xi3m2 -
        (16*f^3*(2*f + 3*bη + 2*f*bη + 3*bη^2)*xi2m1*xi3m2)/7875.

    M[1,2,3] = @. xi10*((264*b1*f^2*xi2m1)/1715. + (6*b1^2*f^2*xi2m1)/35. + (6*b2*f^2*xi2m1)/35. -
        (306*b1*f^3*xi2m1)/1225. +
        (16*f^2*bK²*xi2m1)/1225. + (44*b1*f*bKKpara*xi2m1)/175. - (1012*f^2*bKKpara*xi2m1)/3675. +
        (6*b1*f^2*bδη*xi2m1)/35. - (306*f^3*bδη*xi2m1)/1225. + (264*b1*f^2*bη*xi2m1)/1715. +
        (6*b2*f^2*bη*xi2m1)/35. - (88*f^3*bη*xi2m1)/1225. - (558*b1*f^3*bη*xi2m1)/1225. +
        (216*f^4*bη*xi2m1)/1225. + (16*f^2*bK²*bη*xi2m1)/1225. - (1012*f^2*bKKpara*bη*xi2m1)/3675. -
        (306*f^3*bδη*bη*xi2m1)/1225. - (88*f^3*bη^2*xi2m1)/1225. + (216*f^4*bη^2*xi2m1)/1225. -
        (36*b1*f^3*bη²*xi2m1)/175. + (216*f^4*bη²*xi2m1)/1225. + (216*f^4*bη*bη²*xi2m1)/1225. +
        (24*b1*f*bΠ2para*xi2m1)/175. - (4*f^2*bΠ2para*xi2m1)/147. - (4*f^2*bη*bΠ2para*xi2m1)/147. +
        (16*f^2*(1 + bη)*(270*b1 + 210*b1*f - 378*bK² + 497*bKKpara + 210*f*bδη - 126*f*bη - 252*f^2*bη - 252*f^2*bη² +
        567*bΠ2para)*xi4m1)/77175.) + xi30*
        (((2*f*(-60*b1*f + 315*b1^2*f - 105*b1*f^2 + 84*f*bK² + 420*b1*bKKpara - 196*f*bKKpara + 315*b1*f*bδη -
        105*f^2*bδη - 60*b1*f*bη +
        28*f^2*bη - 315*b1*f^2*bη + 126*f^3*bη + 84*f*bK²*bη - 196*f*bKKpara*bη -
        105*f^2*bδη*bη + 28*f^2*bη^2 + 126*f^3*bη^2 - 210*b1*f^2*bη² + 126*f^3*bη² +
        126*f^3*bη*bη²))/3675. + (4*f*(45*b1 - 37*f - 37*f*bη)*bΠ2para)/1575.)*xi0m1 + (64*b1*f^2*xi2m1)/1715. +
        (16*b1*f^3*xi2m1)/525. - (64*f^2*bK²*xi2m1)/1225. - (16*b1*f*bKKpara*xi2m1)/175. +
        (944*f^2*bKKpara*xi2m1)/11025. + (16*f^3*bδη*xi2m1)/525. + (64*b1*f^2*bη*xi2m1)/1715. -
        (64*f^3*bη*xi2m1)/3675. + (16*b1*f^3*bη*xi2m1)/175. - (64*f^4*bη*xi2m1)/1225. -
        (64*f^2*bK²*bη*xi2m1)/1225. + (944*f^2*bKKpara*bη*xi2m1)/11025. +
        (16*f^3*bδη*bη*xi2m1)/525. - (64*f^3*bη^2*xi2m1)/3675. - (64*f^4*bη^2*xi2m1)/1225. +
        (32*b1*f^3*bη²*xi2m1)/525. - (64*f^4*bη²*xi2m1)/1225. - (64*f^4*bη*bη²*xi2m1)/1225. -
        (16*b1*f*bΠ2para*xi2m1)/175. + (208*f^2*bΠ2para*xi2m1)/2205. + (208*f^2*bη*bΠ2para*xi2m1)/2205. -
        (16*f^2*(1 + bη)*(495*b1 - 693*bK² + 462*bKKpara - 231*f*bη - 252*f^2*bη - 252*f^2*bη² + 847*bΠ2para)*
        xi4m1)/282975.)

    M[3,2,3] = @. (-2*f^2*(18*b1 + 35*b1*f - 21*f^2 + 18*b1*bη - 14*f*bη - 21*f^2*bη - 14*f*bη^2)*xi1m2*xi2m1)/525. +
        ((-2*f^2*(9*b1 - 7*f^2))/315. - (2*f^2*(9*b1 - 7*f - 7*f^2)*bη)/315. + (2*f^3*bη^2)/45.)*xi0m1*xi3m2 +
        (8*f^2*(1 + bη)*(27*b1 - 14*f^2 - 21*f*bη)*xi2m1*xi3m2)/4725.

    M[1,4,3] = @. xi30*(((8*f*(-3300*b1*f + 17325*b1^2*f - 5775*b1*f^2 + 4620*f*bK² + 23100*b1*bKKpara - 10780*f*bKKpara +
        17325*b1*f*bδη -
        5775*f^2*bδη - 3300*b1*f*bη + 1512*f^2*bη - 17325*b1*f^2*bη + 6930*f^3*bη + 4620*f*bK²*bη -
        10780*f*bKKpara*bη - 5775*f^2*bδη*bη + 1512*f^2*bη^2 + 6930*f^3*bη^2 - 11550*b1*f^2*bη² +
        6930*f^3*bη² + 6930*f^3*bη*bη²))/606375. + (16*f*(11*b1 - 9*f - 9*f*bη)*bΠ2para)/1155.)*xi0m1 +
        (256*b1*f^2*xi2m1)/5145. + (64*b1*f^3*xi2m1)/1575. - (256*f^2*bK²*xi2m1)/3675. - (64*b1*f*bKKpara*xi2m1)/525. +
        (3776*f^2*bKKpara*xi2m1)/33075. + (64*f^3*bδη*xi2m1)/1575. + (256*b1*f^2*bη*xi2m1)/5145. -
        (1536*f^3*bη*xi2m1)/67375. + (64*b1*f^3*bη*xi2m1)/525. - (256*f^4*bη*xi2m1)/3675. -
        (256*f^2*bK²*bη*xi2m1)/3675. + (3776*f^2*bKKpara*bη*xi2m1)/33075. +
        (64*f^3*bδη*bη*xi2m1)/1575. - (1536*f^3*bη^2*xi2m1)/67375. - (256*f^4*bη^2*xi2m1)/3675. +
        (128*b1*f^3*bη²*xi2m1)/1575. - (256*f^4*bη²*xi2m1)/3675. - (256*f^4*bη*bη²*xi2m1)/3675. -
        (64*b1*f*bΠ2para*xi2m1)/525. + (5056*f^2*bΠ2para*xi2m1)/40425. + (5056*f^2*bη*bΠ2para*xi2m1)/40425. -
        (64*f^2*(1 + bη)*(825*b1 - 1155*bK² + 770*bKKpara - 378*f*bη - 420*f^2*bη - 420*f^2*bη² + 1400*bΠ2para)*
        xi4m1)/1.414875e6) + xi10*((352*b1*f^2*xi2m1)/1715. + (8*b1^2*f^2*xi2m1)/35. + (8*b2*f^2*xi2m1)/35. -
        (408*b1*f^3*xi2m1)/1225. + (64*f^2*bK²*xi2m1)/3675. + (176*b1*f*bKKpara*xi2m1)/525. -
        (4048*f^2*bKKpara*xi2m1)/11025. + (8*b1*f^2*bδη*xi2m1)/35. - (408*f^3*bδη*xi2m1)/1225. +
        (352*b1*f^2*bη*xi2m1)/1715. + (8*b2*f^2*bη*xi2m1)/35. - (576*f^3*bη*xi2m1)/6125. -
        (744*b1*f^3*bη*xi2m1)/1225. + (288*f^4*bη*xi2m1)/1225. + (64*f^2*bK²*bη*xi2m1)/3675. -
        (4048*f^2*bKKpara*bη*xi2m1)/11025. - (408*f^3*bδη*bη*xi2m1)/1225. - (576*f^3*bη^2*xi2m1)/6125. +
        (288*f^4*bη^2*xi2m1)/1225. - (48*b1*f^3*bη²*xi2m1)/175. + (288*f^4*bη²*xi2m1)/1225. +
        (288*f^4*bη*bη²*xi2m1)/1225. + (32*b1*f*bΠ2para*xi2m1)/175. - (48*f^2*bΠ2para*xi2m1)/1225. -
        (48*f^2*bη*bΠ2para*xi2m1)/1225. + (64*f^2*(1 + bη)*
        (14850*b1 + 11550*b1*f - 20790*bK² + 27335*bKKpara + 11550*f*bδη - 6804*f*bη - 13860*f^2*bη - 13860*f^2*bη² +
        30975*bΠ2para)*xi4m1)/1.2733875e7)

    M[3,4,3] = @. (-16*f^2*(55*b1 + 105*b1*f - 63*f^2 + 55*b1*bη - 42*f*bη - 63*f^2*bη - 42*f*bη^2)*xi1m2*xi2m1)/9625. +
        ((-8*f^2*(55*b1 - 42*f^2))/5775. - (8*f^2*(55*b1 - 42*f - 42*f^2)*bη)/5775. + (16*f^3*bη^2)/275.)*xi0m1*xi3m2 +
        (32*f^2*(1 + bη)*(55*b1 - 28*f^2 - 42*f*bη)*xi2m1*xi3m2)/28875.

    M[1,6,3] = @. xi10*((-64*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/5145. -
        (256*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/56595.) +
        xi30*(((64*f^3*bη*(1 + bη))/8085. - (64*f^2*(1 + bη)*bΠ2para)/4851.)*xi0m1 -
        (512*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/169785. +
        (128*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/56595.)

    M[3,6,3] = @. (-16*f^3*(5*b1 - 3*f - 2*bη - 3*f*bη - 2*bη^2)*xi1m2*xi2m1)/1155. +
        ((16*f^4)/693. + (16*f^3*(1 + f)*bη)/693. + (16*f^3*bη^2)/693.)*xi0m1*xi3m2 -
        (64*f^3*(2*f + 3*bη + 2*f*bη + 3*bη^2)*xi2m1*xi3m2)/10395.

    M[4,0,4] = @. (16*f^4*xi2m2^2)/11025.

    M[0,2,4] = @. ((-384*f^3*bη*(bη + bη²))/16807. + (640*f^2*(bη + bη²)*bΠ2para)/16807.)*xi20^2 +
        ((256*f^3*bη*(bη + bη²))/25725. - (256*f^2*(bη + bη²)*bΠ2para)/15435.)*xi00*xi40 +
        ((-512*f^3*bη*(bη + bη²))/50421. + (2560*f^2*(bη + bη²)*bΠ2para)/151263.)*xi20*xi40 +
        ((1152*f^3*bη*(bη + bη²))/420175. - (384*f^2*(bη + bη²)*bΠ2para)/84035.)*xi40^2

    M[2,2,4] = @. ((-32*f^2*(30*b1*f - 9*f^2 + 22*bKKpara + 15*f*bδη + 15*b1*f*bη - 36*f^2*bη - 12*f*bη^2 - 9*f^2*bη^2 -
        18*f^2*bη² - 12*f*bη*bη²))/11025. - (128*f^2*bΠ2para)/3675.)*xi1m1*xi3m1 +
        ((-16*f^2*(9*f^2 - 22*bKKpara + 36*f^2*bη + 22*f*bη^2 + 9*f^2*bη^2 + 18*f^2*bη² + 22*f*bη*bη²))/40425. +
        (32*f^2*bΠ2para)/3675.)*xi3m1^2 + ((16*f^4)/2205. + (32*f^4*bη)/2205. + (16*f^4*bη^2)/2205.)*xi00*xi4m2 +
        xi20*(((-32*(7*b1 - 3*f)*f^3)/5145. - (32*(7*b1 - 6*f)*f^3*bη)/5145. + (32*f^4*bη^2)/1715.)*xi2m2 +
        ((-128*f^4)/33957. - (256*f^4*bη)/33957. - (128*f^4*bη^2)/33957.)*xi4m2)

    M[4,2,4] = @. (272*f^4*xi2m2^2)/24255.

    M[0,4,4] = @. ((4*b1^2*f^2)/35. + (256*b1*f*bKKpara)/735. + (4168*bKKpara^2)/15435. + (8*b1*f^2*bδη)/35. +
        (256*f*bKKpara*bδη)/735. +
        (4*f^2*bδη^2)/35. + (480*b1*f^2*bη)/2401. + (8*b2*f^2*bη)/35. - (96*b1*f^3*bη)/245. + (128*f^2*bK²*bη)/5145. -
        (7136*f^2*bKKpara*bη)/15435. - (96*f^3*bδη*bη)/245. - (11232*f^3*bη^2)/184877. + (288*f^4*bη^2)/1715. +
        (480*b1*f^2*bη²)/2401. + (8*b2*f^2*bη²)/35. - (96*b1*f^3*bη²)/245. + (128*f^2*bK²*bη²)/5145. -
        (7136*f^2*bKKpara*bη²)/15435. - (96*f^3*bδη*bη²)/245. - (11232*f^3*bη*bη²)/184877. +
        (576*f^4*bη*bη²)/1715. + (288*f^4*bη²^2)/1715. +
        (48*(3773*b1*f + 5929*bKKpara + 3773*f*bδη - 2901*f^2*bη - 2901*f^2*bη²)*bΠ2para)/924385. +
        (152*bΠ2para^2)/1715.)*xi20^2 +
        ((4*b1^2*f^2)/35. + (32*b1*f*bKKpara)/105. + (16*bKKpara^2)/75. + (8*b1*f^2*bδη)/35. + (32*f*bKKpara*bδη)/105. +
        (4*f^2*bδη^2)/35. - (64*b1*f^2*bη)/735. - (16*b1*f^3*bη)/105. + (64*f^2*bK²*bη)/525. -
        (64*f^2*bKKpara*bη)/225. - (16*f^3*bδη*bη)/105. + (2496*f^3*bη^2)/94325. + (16*f^4*bη^2)/175. -
        (64*b1*f^2*bη²)/735. - (16*b1*f^3*bη²)/105. + (64*f^2*bK²*bη²)/525. - (64*f^2*bKKpara*bη²)/225. -
        (16*f^3*bδη*bη²)/105. + (2496*f^3*bη*bη²)/94325. + (32*f^4*bη*bη²)/175. + (16*f^4*bη²^2)/175. +
        (16*(8085*b1*f + 11858*bKKpara + 8085*f*bδη - 12042*f^2*bη - 12042*f^2*bη²)*bΠ2para)/848925. +
        (32*bΠ2para^2)/525.)*xi00*
        xi40 + ((16*(-11319*b1*f*bKKpara - 14014*bKKpara^2 - 11319*f*bKKpara*bδη + 9240*b1*f^2*bη + 8232*b1*f^3*bη -
        12936*f^2*bK²*bη + 22246*f^2*bKKpara*bη + 8232*f^3*bδη*bη - 2808*f^3*bη^2 - 7056*f^4*bη^2 +
        9240*b1*f^2*bη² + 8232*b1*f^3*bη² - 12936*f^2*bK²*bη² + 22246*f^2*bKKpara*bη² +
        8232*f^3*bδη*bη² - 2808*f^3*bη*bη² - 14112*f^4*bη*bη² - 7056*f^4*bη²^2))/1.663893e6 -
        (16*(11319*b1*f + 20482*bKKpara + 11319*f*bδη - 21438*f^2*bη - 21438*f^2*bη²)*bΠ2para)/1.663893e6 -
        (64*bΠ2para^2)/1029.)*
        xi20*xi40 + ((144*(7007*bKKpara^2 - 10010*b1*f^2*bη + 14014*f^2*bK²*bη - 11466*f^2*bKKpara*bη +
        3042*f^3*bη^2 + 3528*f^4*bη^2 - 10010*b1*f^2*bη² + 14014*f^2*bK²*bη² - 11466*f^2*bKKpara*bη² +
        3042*f^3*bη*bη² + 7056*f^4*bη*bη² + 3528*f^4*bη²^2))/6.0085025e7 +
        (288*(539*bKKpara - 636*f^2*bη - 636*f^2*bη²)*bΠ2para)/4.621925e6 + (144*bΠ2para^2)/8575.)*xi40^2

    M[2,4,4] = @. ((-8*f^2*(390*b1*f - 117*f^2 + 286*bKKpara + 195*f*bδη + 308*b1*bη + 195*b1*f*bη - 468*f^2*bη -
        156*f*bη^2 -
        117*f^2*bη^2 + 308*b1*bη² - 234*f^2*bη² - 156*f*bη*bη²))/13475. - (1248*f^2*bΠ2para)/13475.)*xi1m1*
        xi3m1 + ((-4*f^2*(351*f^2 - 858*bKKpara - 1694*b1*bη + 1404*f^2*bη + 858*f*bη^2 + 351*f^2*bη^2 - 1694*b1*bη² +
        702*f^2*bη² + 858*f*bη*bη²))/148225. + (312*f^2*bΠ2para)/13475.)*xi3m1^2 +
        ((52*f^4)/2695. + (104*f^4*bη)/2695. + (52*f^4*bη^2)/2695.)*xi00*xi4m2 +
        xi20*(((-312*(7*b1 - 3*f)*f^3)/18865. - (312*(7*b1 - 6*f)*f^3*bη)/18865. + (936*f^4*bη^2)/18865.)*xi2m2 +
        ((-416*f^4)/41503. - (832*f^4*bη)/41503. - (416*f^4*bη^2)/41503.)*xi4m2)

    M[4,4,4] = @. (3858*f^4*xi2m2^2)/175175.

    M[0,6,4] = @. ((-960*f^3*bη*(bη + bη²))/26411. + (1600*f^2*(bη + bη²)*bΠ2para)/26411.)*xi20^2 +
        ((128*f^3*bη*(bη + bη²))/8085. - (128*f^2*(bη + bη²)*bΠ2para)/4851.)*xi00*xi40 +
        ((-1280*f^3*bη*(bη + bη²))/79233. + (6400*f^2*(bη + bη²)*bΠ2para)/237699.)*xi20*xi40 +
        ((576*f^3*bη*(bη + bη²))/132055. - (192*f^2*(bη + bη²)*bΠ2para)/26411.)*xi40^2

    M[2,6,4] = @. ((-16*f^2*(30*b1*f - 9*f^2 + 22*bKKpara + 15*f*bδη + 15*b1*f*bη - 36*f^2*bη - 12*f*bη^2 - 9*f^2*bη^2 -
        18*f^2*bη² - 12*f*bη*bη²))/3465. - (64*f^2*bΠ2para)/1155.)*xi1m1*xi3m1 +
        ((-8*f^2*(9*f^2 - 22*bKKpara + 36*f^2*bη + 22*f*bη^2 + 9*f^2*bη^2 + 18*f^2*bη² + 22*f*bη*bη²))/12705. +
        (16*f^2*bΠ2para)/1155.)*xi3m1^2 + ((8*f^4)/693. + (16*f^4*bη)/693. + (8*f^4*bη^2)/693.)*xi00*xi4m2 +
        xi20*(((-16*(7*b1 - 3*f)*f^3)/1617. - (16*(7*b1 - 6*f)*f^3*bη)/1617. + (16*f^4*bη^2)/539.)*xi2m2 +
        ((-320*f^4)/53361. - (640*f^4*bη)/53361. - (320*f^4*bη^2)/53361.)*xi4m2)

    M[4,6,4] = @. (424*f^4*xi2m2^2)/24255.

    M[4,8,4] = @. (32*f^4*xi2m2^2)/6435.

    M[3,2,5] = @. (160*f^4*(1 + bη)*xi2m1*xi3m2)/14553.

    M[1,4,5] = @. (-8*f^2*(1 + bη)*(15*b1*f + 22*bKKpara + 15*f*bδη - 18*f^2*bη - 18*f^2*bη² + 12*bΠ2para)*xi10*xi4m1)/
        2079. + xi30*((-40*f^2*(21*b1*f + 32*bKKpara + 21*f*bδη + 63*b1*f*bη - 36*f^2*bη + 32*bKKpara*bη +
        21*f*bδη*bη - 36*f^2*bη^2 + 42*b1*f*bη² - 36*f^2*bη² - 36*f^2*bη*bη² + 18*bΠ2para +
        18*bη*bΠ2para)*xi2m1)/14553. - (32*f^2*(1 + bη)*(-13*bKKpara + 12*f^2*bη + 12*f^2*bη² - 13*bΠ2para)*
        xi4m1)/21021.)

    M[3,4,5] = @. (920*f^4*(1 + bη)*xi2m1*xi3m2)/21021.

    M[1,6,5] = @.  (-16*f^2*(1 + bη)*(15*b1*f + 22*bKKpara + 15*f*bδη - 18*f^2*bη - 18*f^2*bη² + 12*bΠ2para)*xi10*xi4m1)/
        3465. + xi30*((-16*f^2*(21*b1*f + 32*bKKpara + 21*f*bδη + 63*b1*f*bη - 36*f^2*bη + 32*bKKpara*bη +
        21*f*bδη*bη - 36*f^2*bη^2 + 42*b1*f*bη² - 36*f^2*bη² - 36*f^2*bη*bη² + 18*bΠ2para +
        18*bη*bΠ2para)*xi2m1)/4851. - (64*f^2*(1 + bη)*(-13*bKKpara + 12*f^2*bη + 12*f^2*bη² - 13*bΠ2para)*
        xi4m1)/35035.)

    M[3,6,5] = @. (544*f^4*(1 + bη)*xi2m1*xi3m2)/10395.

    M[3,8,5] = @. (128*f^4*(1 + bη)*xi2m1*xi3m2)/6435.

    M[2,4,6] = @. ((80*f^4*(1 + 4*bη + bη^2))/11011. + (160*f^4*bη²)/11011.)*xi3m1^2 +
        ((80*f^4)/11011. + (160*f^4*bη)/11011. + (80*f^4*bη^2)/11011.)*xi20*xi4m2

    M[0,6,6] = @. ((-32*f^2*(bη + bη²)*(21*b1*f + 32*bKKpara + 21*f*bδη - 18*f^2*bη - 18*f^2*bη²))/4851. -
        (64*f^2*(bη + bη²)*bΠ2para)/539.)*xi20*xi40 +
        ((-32*f^2*(bη + bη²)*(-15*bKKpara + 8*f^2*bη + 8*f^2*bη²))/24255. + (32*f^2*(bη + bη²)*bΠ2para)/1617.)*xi40^2

    M[2,6,6] = @. ((664*f^4*(1 + 4*bη + bη^2))/38115. + (1328*f^4*bη²)/38115.)*xi3m1^2 +
        ((664*f^4)/38115. + (1328*f^4*bη)/38115. + (664*f^4*bη^2)/38115.)*xi20*xi4m2

    M[2,8,6] = @. ((64*f^4*(1 + 4*bη + bη^2))/6435. + (128*f^4*bη²)/6435.)*xi3m1^2 +
        ((64*f^4)/6435. + (128*f^4*bη)/6435. + (64*f^4*bη^2)/6435.)*xi20*xi4m2

    M[1,6,7] = @. (224*f^4*(1 + bη)*(bη + bη²)*xi30*xi4m1)/6435.

    M[1,8,7] = @. (256*f^4*(1 + bη)*(bη + bη²)*xi30*xi4m1)/6435.

    M[0,8,8] = @. ((128*f^4*bη^2)/6435. + (256*f^4*bη*bη²)/6435. + (128*f^4*bη²^2)/6435.)*xi40^2
    #end

    return r, M
end
#generateMnli(x->x, f=1.0, BiasDict22=Dict([("b1",0), ("bη", 0), ("b2", 0), ("bK2", 0), ("bδη", 0),
#                           ("bη2", 0), ("bKK∥", 0), ("bΠ2∥", 0)]));
precompile(ξ, (Function, Int64, Int64, Int64))
precompile(ξ, (Function, Int64, Int64))
precompile(generateMnli, (Function, Dict{String,Int64}, Float64))
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
function wrap_api_call(spec::SpecFunc, args; with_func_ptr = false)
    ex = :($(spec.name)($((with_func_ptr ? [args; first(func_ptrs(spec))] : args)...)))
    ex = if with_func_ptr
        ex
    else
        dispatch_call(spec, ex)
    end
    wrap_return(
        ex,
        spec.return_type,
        idiomatic_julia_type(spec.return_type),
    )
end

function dispatch_handle(spec::SpecFunc)
    maybe_handle = !isempty(children(spec)) ? first(children(spec)).type : :nothing
    if maybe_handle in spec_handles.name
        handle = handle_by_name(maybe_handle)
        handle_id = wrap_identifier(handle)
        hierarchy = parent_hierarchy(handle)
        if handle.name == :VkDevice || handle.name == :VkInstance
            # to avoid name conflicts
            handle_id
        elseif :VkDevice in hierarchy
            :(device($handle_id))
        elseif :VkInstance in hierarchy
            :(instance($handle_id))
        end
    else
        :nothing
    end
end

dispatch_call(spec::SpecFunc, ex) = spec.name in (:vkGetInstanceProcAddr, :vkGetDeviceProcAddr) ? ex : :(@dispatch $(dispatch_handle(spec)) $ex)

function wrap_enumeration_api_call(spec::SpecFunc, exs::Expr...; free = [])
    if must_repeat_while_incomplete(spec)
        if !isempty(free)
            free_block = quote
                if _return_code == VK_INCOMPLETE
                    $(map(x -> :(Libc.free($(x.name))), free)...)
                end
            end
            [:(@repeat_while_incomplete $(Expr(:block, exs..., free_block)))]
        else
            [:(@repeat_while_incomplete $(Expr(:block, exs...)))]
        end
    else
        exs
    end
end

function APIFunction(spec::SpecFunc, with_func_ptr)
    p = Dict(
        :category => :function,
        :name => function_name(spec),
        :relax_signature => true,
    )

    count_ptr_index = findfirst(x -> (is_length(x) || is_size(x)) && x.requirement == POINTER_REQUIRED, children(spec))
    queried_params = getindex(children(spec), findall(is_implicit_return, children(spec)))
    if !isnothing(count_ptr_index)
        count_ptr = children(spec)[count_ptr_index]
        queried_params =
            getindex(children(spec), findall(x -> x.len == count_ptr.name && !x.is_constant, children(spec)))

        first_call_args = map(@λ(begin
            &count_ptr => count_ptr.name
            GuardBy(in(queried_params)) => :C_NULL
            x => vk_call(x)
        end), children(spec))

        i = 0
        second_call_args = map(@λ(begin
            :C_NULL && Do(i += 1) => queried_params[i].name
            x => x
        end), first_call_args)

        p[:body] = concat_exs(
            initialize_ptr(count_ptr),
            wrap_enumeration_api_call(
                spec,
                wrap_api_call(spec, first_call_args; with_func_ptr),
                (is_length(count_ptr) ? initialize_array : initialize_ptr).(queried_params, count_ptr)...,
                wrap_api_call(spec, second_call_args; with_func_ptr),
                ;
                free = filter(is_data, queried_params),
            )...,
            wrap_implicit_return(spec, queried_params; with_func_ptr),
        )

        args = filter(!in(vcat(queried_params, count_ptr)), children(spec))

        ret_type = @match length(queried_params) begin
            if any(is_data_with_retrievable_size, queried_params)
            end => Expr(:curly, :Tuple, idiomatic_julia_type.([unique(len.(queried_params)); queried_params])...)
            1 => idiomatic_julia_type(first(queried_params))
            _ => Expr(
                :curly,
                :Tuple,
                (idiomatic_julia_type(param) for param ∈ queried_params)...,
            )
        end
    elseif !isempty(queried_params)
        call_args = map(@λ(begin
            x && GuardBy(in(queried_params)) => x.name
            x => vk_call(x)
        end), children(spec))

        args = filter(!in(filter(x -> x.requirement ≠ POINTER_REQUIRED, queried_params)), children(spec))

        ret_type = @match length(queried_params) begin
            1 => idiomatic_julia_type(first(queried_params))
            _ => Expr(:curly, :Tuple, (idiomatic_julia_type(param) for param ∈ queried_params)...)
        end

        if spec.type == FTYPE_QUERY && length(queried_params) == 1 && begin
                t = ptr_type(only(queried_params).type)
                is_vulkan_type(t) && !in(t, spec_handles.name) && any(Base.Fix1(in, t), spec_structs.extends)
            end

            param = only(queried_params)
            t = ptr_type(param.type)
            intermediate_t = struct_name(t)
            var = wrap_identifier(param.name)
            p[:body] = quote
                $var = initialize($intermediate_t, next_types...)
                $(param.name) = Ref(Base.unsafe_convert($t, $var))
                GC.@preserve $var begin
                    $(wrap_api_call(spec, call_args; with_func_ptr))
                    $intermediate_t($(param.name)[], Any[$var])
                end
            end

            add_func_args!(p, spec, args; with_func_ptr)
            push!(p[:args], :(next_types::Type...))
            p[:return_type] = wrap_return_type(spec, ret_type)
            return APIFunction(spec, with_func_ptr, p)
        else
            p[:body] = concat_exs(
                map(initialize_ptr, queried_params)...,
                wrap_api_call(spec, call_args; with_func_ptr),
                wrap_implicit_return(spec, queried_params; with_func_ptr),
            )
        end
    else
        p[:short] = true
        p[:body] = :($(wrap_api_call(spec, map(vk_call, children(spec)); with_func_ptr)))

        args = children(spec)

        ret_type = idiomatic_julia_type(spec.return_type)
    end

    add_func_args!(p, spec, args; with_func_ptr)
    p[:return_type] = wrap_return_type(spec, ret_type)
    APIFunction(spec, with_func_ptr, p)
end

"""
Extend functions that create (or allocate) one or several handles,
by exposing the parameters of the associated CreateInfo structures.
`spec` must have one or several CreateInfo arguments.
"""
function APIFunction(spec::CreateFunc, with_func_ptr)
    @assert !isnothing(spec.create_info_param) "Cannot extend handle constructor with no create info parameter."
    def = APIFunction(spec.func, false)
    p_func = def.p
    p_info = Constructor(StructDefinition{false}(spec.create_info_struct)).p

    args = [p_func[:args]; p_info[:args]]
    kwargs = [p_func[:kwargs]; p_info[:kwargs]]

    info_expr = reconstruct_call(p_info; is_decl = false)
    info_index = findfirst(==(spec.create_info_param), filter(!is_optional, children(spec.func)))
    deleteat!(args, info_index)

    func_call_args::Vector{ExprLike} = name.(p_func[:args])
    func_call_args[info_index] = info_expr

    if with_func_ptr
        append!(args, func_ptr_args(spec.func))
        append!(func_call_args, func_ptrs(spec.func))
    end

    body = reconstruct_call(Dict(:name => name(def), :args => func_call_args, :kwargs => name.(p_func[:kwargs])))

    p = Dict(
        :category => :function,
        :name => p_func[:name],
        :args => args,
        :kwargs => kwargs,
        :short => true,
        :body => body,
        :relax_signature => true,
        # Do not include the return type in generated code, but keep the return type information for the promotion to a high-level function.
        :_return_type => p_func[:return_type],
    )
    APIFunction(spec, with_func_ptr, p)
end

function contains_api_structs(def::Union{APIFunction,Constructor})
    any(x -> x ≠ promote_hl(x), def.p[:args])
end

is_promoted(ex) = ex == promote_hl(ex)

function promote_hl(def::APIFunction)
    promoted = APIFunction(def, def.with_func_ptr, promote_hl(def.p))
    type = def.p[def.spec isa CreateFunc ? :_return_type : :return_type]
    wrap_body = :(next_types::Type...) in promoted.p[:args] ? promote_return_hl_next_types : promote_return_hl
    merge!(promoted.p,
        Dict(
            :short => false,
            :body => wrap_body(type, promoted.p[:body])
        )
    )
    promote_return_type_hl!(promoted, type)
    promoted
end

function promote_return_type_hl!(promoted, type)
    T = hl_type(type)
    include_rtype = false
    rtype = @match type begin
        :(ResultTypes.Result{$S,$E}) => (include_rtype = true; :(ResultTypes.Result{$(promote_return_type_hl!(nothing, S)), $E}))
        :(Vector{$_T}) => :(Vector{$(promote_return_type_hl!(nothing, _T))})
        _ => T
    end
    include_rtype && (promoted.p[:return_type] = rtype)
    rtype
end

function promote_hl(def::Constructor)
    Constructor(promote_hl(def.p), def.to, def.from)
end

function promote_hl(arg::ExprLike)
    id, type = @match arg begin
        :($id::$t) => (id, t)
        _ => return arg
    end
    type = postwalk(type) do ex
        if ex isa Symbol && startswith(string(ex), '_')
            Symbol(string(ex)[2:end]) # remove underscore prefix
        else
            ex
        end
    end
    :($id::$type)
end

function promote_hl(p::Dict)
    args = promote_hl.(p[:args])
    modified_args = [arg for (arg, new_arg) in zip(p[:args], args) if arg ≠ new_arg]
    call_args = map(p[:args]) do arg
        id, type = @match arg begin
            :($id::$t) => (id, t)
            :($id::$t...) => (:($id...), nothing)
            id => (id, nothing)
        end
        if arg in modified_args
            T = @match type begin
                :(AbstractArray{<:$t}) => :(Vector{$t})
                t => t
            end
            :(convert($T, $id))
        else
            id
        end
    end
    p = Dict(
        :category => :function,
        :name => promote_name_hl(p[:name]),
        :args => args,
        :kwargs => p[:kwargs],
        :body => reconstruct_call(Dict(:name => p[:name], :args => call_args, :kwargs => name.(p[:kwargs]))),
        :short => true,
        :relax_signature => true,
    )
end

function function_name(sym::Symbol, is_high_level = false)
    sym = nc_convert(SnakeCaseLower, remove_vk_prefix(sym))
    is_high_level ? sym : Symbol('_', sym)
end
function_name(spec::Spec, is_high_level = false) = function_name(spec.name, is_high_level)

function promote_name_hl(name)
    str = String(name)
    startswith(str, '_') ? Symbol(str[2:end]) : name
end

function promote_return_hl_next_types(type, ex)
    rex = promote_return_hl(type, ex)
    call = isexpr(rex, :block) ? last(rex.args) : rex
    push!(call.args, :(next_types_hl...))
    quote
        next_types_hl = next_types
        next_types = intermediate_type.(next_types)
        $((isexpr(rex, :block) ? rex.args : [rex])...)
    end
end

function promote_return_hl(type, ex)
    T = hl_type(type)
    @match type begin
        :Cvoid => ex
        :(ResultTypes.Result{$S,$E}) => Expr(:block, :(val = @propagate_errors $ex), promote_return_hl(S, :val))
        :(Vector{$T}) => begin
            rex = promote_return_hl(T, ex)
            rex isa Symbol ? rex : broadcast_ex(rex)
        end
        GuardBy(is_intermediate) => :($T($ex))
        _ => ex
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
struct SMPSStage{T <: AbstractFloat, M <: AbstractMatrix}
    id::Int
    model::LPData{T,M}
    uncertain::LPData{T,M}
    technology::UnitRange{Int}
    recourse::UnitRange{Int}

    function SMPSStage(model::LPData{T,M}) where {T <: AbstractFloat, M <: AbstractMatrix}
        return new{T, M}(1, model, LPData(T,M,0,0,0), 0:0, 0:0)
    end

    function SMPSStage(id::Integer, model::LPData{T,M}, uncertain::LPData{T,M}, technology::UnitRange{Int}, recourse::UnitRange{Int}) where {T <: AbstractFloat, M <: AbstractMatrix}
        return new{T, M}(id, model, uncertain, technology, recourse)
    end
end