text:
# Constructors
# ============
#
# Constructor methods for Biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Construct a blank sequence of length `len`.
function BioSequence{A}(len::Integer) where {A<:Alphabet}
    return BioSequence{A}(Vector{UInt64}(undef, seq_data_len(A, len)), 1:len, false)
end

# Construct a void sequence.
function BioSequence()
    return BioSequence{VoidAlphabet}(Vector{UInt64}(), 0:-1, false)
end

BioSequence(::Type{DNA}) = DNASequence()
BioSequence(::Type{RNA}) = RNASequence()
BioSequence(::Type{AminoAcid}) = AminoAcidSequence()
BioSequence(::Type{Char}) = CharSequence()

# Create a sequence from a string or a vector.
function BioSequence{A}(
        src::Union{AbstractString, AbstractVector},
        startpos::Integer = 1,
        stoppos::Integer = length(src)) where {A<:Alphabet}
    len = stoppos - startpos + 1
    seq = BioSequence{A}(len)
    return encode_copy!(seq, 1, src, startpos, len)
end

# Create a sequence that is a subsequence of another sequence.
function BioSequence{A}(
        other::BioSequence{A},
        part::UnitRange{<:Integer}) where {A<:Alphabet}
    checkbounds(other, part)
    start = other.part.start + part.start - 1
    stop = start + length(part) - 1
    subseq = BioSequence{A}(other.data, start:stop, true)
    other.shared = true
    return subseq
end
function BioSequence(other::BioSequence{A}, part::UnitRange) where {A<:Alphabet}
    return BioSequence{A}(other, part)
end

# Concatenate multiple sequences.
function BioSequence{A}(chunks::BioSequence{A}...) where {A<:Alphabet}
    len = 0
    for chunk in chunks
        len += length(chunk)
    end
    seq = BioSequence{A}(len)
    offset = 1
    for chunk in chunks
        copyto!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

# Create a 4 bit DNA/RNA sequences from a 2 bit DNA/RNA sequences, and vice-versa.
for (alpha, alphb) in [(DNAAlphabet{4}, DNAAlphabet{2}), # DNA to DNA
                       (DNAAlphabet{2}, DNAAlphabet{4}),
                       (RNAAlphabet{4}, RNAAlphabet{2}), # RNA to RNA
                       (RNAAlphabet{2}, RNAAlphabet{4}),
                       (DNAAlphabet{2}, RNAAlphabet{4}), # DNA to RNA
                       (DNAAlphabet{4}, RNAAlphabet{2}),
                       (RNAAlphabet{4}, DNAAlphabet{2}), # RNA to DNA
                       (RNAAlphabet{2}, DNAAlphabet{4})]

    @eval function (::Type{BioSequence{$alpha}})(seq::BioSequence{$alphb})
        newseq = BioSequence{$alpha}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end
end

for (alpha, alphb) in [(DNAAlphabet{2}, RNAAlphabet{2}),
                       (RNAAlphabet{2}, DNAAlphabet{2}),
                       (DNAAlphabet{4}, RNAAlphabet{4}),
                       (RNAAlphabet{4}, DNAAlphabet{4})]

    @eval function (::Type{BioSequence{$alpha}})(seq::BioSequence{$alphb})
        newseq = BioSequence{$alpha}(seq.data, seq.part, true)
        seq.shared = true
        return newseq
    end
end

function Base.repeat(chunk::BioSequence{A}, n::Integer) where {A<:Alphabet}
    seq = BioSequence{A}(length(chunk) * n)
    offset = 1
    for _ in 1:n
        copyto!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

# Concatenation and Base.repeat operators.
Base.:*(chunk::BioSequence{A}, chunks::BioSequence{A}...) where {A} =
    BioSequence{A}(chunk, chunks...)
Base.:^(chunk::BioSequence, n::Integer) = repeat(chunk, n)

function Base.similar(seq::BioSequence{A}, len::Integer = length(seq)) where {A<:Alphabet}
    return BioSequence{A}(len)
end


text:
using CoffeeBeetles
using Test

@testset "CoffeeBeetles.jl" begin
    main()
    files = ["figure$i.pdf" for i in 4:7]
    append!(files, ["table1.txt", "table2.csv"])
    @test all(isfile(file) for file in files)
end

ext: jl


text:
# Autogenerated wrapper script for Ccluster_jll for x86_64-apple-darwin
export libccluster

using Arb_jll
using FLINT_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Ccluster")
JLLWrappers.@declare_library_product(libccluster, "@rpath/libccluster-1.1.7.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Arb_jll, FLINT_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libccluster,
        "lib/libccluster.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ReachabilityAnalysis: logarithmic_norm

@testset "Logarithmic norm" begin

    # Example 1
    A = [1 2; 3 -4]

    @test logarithmic_norm(A, 1) == 4
    @test logarithmic_norm(A, Inf) == 3

    # Example 2 (see Example 2 in "Explicit Error Bounds for Carleman Linearization.
    # Marcelo Forets and Amaury Pouly (2018)").
    A = [0 1; -1 -2]
    Aop2 = kron(A, A)

    @test logarithmic_norm(A, Inf) == 1
    @test logarithmic_norm(Aop2, Inf) == 9;

    @test logarithmic_norm(A, 1) == 1
    @test logarithmic_norm(Aop2, 1) == 9

    @test logarithmic_norm(A, 2) == 0
    @test logarithmic_norm(Aop2, 2) ≈ 4.23606797749979
end


text:
"""
    DepthwiseConvDims

Concrete subclass of `ConvDims` for a depthwise convolution.  Differs primarily due to
characterization by C_in, C_mult, rather than C_in, C_out.  Useful to be separate from
DenseConvDims primarily for channel calculation differences.
"""
struct DepthwiseConvDims{N, K, S, P, D} <: ConvDims{N}
    input_size::NTuple{N, Int}

    kernel_size::NTuple{K, Int}
    channels_in::Int
    channels_multiplier::Int

    stride::NTuple{S, Int}
    padding::NTuple{P, Int}
    dilation::NTuple{D, Int}
    flipkernel::Bool
end

function DepthwiseConvDims(
    x_size::NTuple{M}, w_size::NTuple{M};
    stride = 1, padding = 0, dilation = 1, flipkernel::Bool = false,
) where M
    sstride, ppadding, ddilation = check_spdf(
        x_size, w_size, stride, padding, dilation)

    # Ensure channels are equal
    if x_size[end-1] != w_size[end]
        xs = x_size[end-1]
        ws = w_size[end]
        throw(DimensionMismatch("Input channels must match! ($xs vs. $ws)"))
    end

    DepthwiseConvDims(
        x_size[1:(end - 2)],
        w_size[1:(end - 2)], x_size[end - 1], w_size[end - 1],
        sstride, ppadding, ddilation, flipkernel)
end

function DepthwiseConvDims(x::AbstractArray, w::AbstractArray; kwargs...)
    if ndims(x) != ndims(w)
        throw(DimensionMismatch("Rank of x and w must match! ($(ndims(x)) vs. $(ndims(w)))"))
    end
    return DepthwiseConvDims(size(x), size(w); kwargs...)
end

# Useful for constructing a new DepthwiseConvDims that has only a few elements different
# from the original progenitor object.
@inline DepthwiseConvDims(
    c::DepthwiseConvDims; I=input_size(c), K=kernel_size(c),
    C_in=channels_in(c), C_m=channel_multiplier(c), S=stride(c),
    P=padding(c), D=dilation(c), F=flipkernel(c),
) = DepthwiseConvDims(
    I,
    K, C_in, C_m,
    S, P, D, F)

@inline channels_in(c::DepthwiseConvDims) = c.channels_in
@inline channels_out(c::DepthwiseConvDims) = c.channels_in * c.channels_multiplier
@inline channel_multiplier(c::DepthwiseConvDims) = c.channels_multiplier

@inline input_size(c::DepthwiseConvDims) = c.input_size
@inline kernel_size(c::DepthwiseConvDims) = c.kernel_size

@inline stride(c::DepthwiseConvDims) = c.stride
@inline padding(c::DepthwiseConvDims) = c.padding
@inline dilation(c::DepthwiseConvDims) = c.dilation
@inline flipkernel(c::DepthwiseConvDims) = c.flipkernel

# This one is basically the same as for DenseConvDims, we only change a few lines for kernel channel count
function check_dims(x::NTuple{M}, w::NTuple{M}, y::NTuple{M}, cdims::DepthwiseConvDims) where {M}
    # First, check that channel counts are all correct:
    @assert x[M-1] == channels_in(cdims) DimensionMismatch("Data input channel count ($(x[M-1]) vs. $(channels_in(cdims)))")
    @assert y[M-1] == channels_out(cdims) DimensionMismatch("Data output channel count ($(y[M-1]) vs. $(channels_out(cdims)))")
    @assert w[M-1] == channel_multiplier(cdims) DimensionMismatch("Kernel multiplier channel count ($(w[M-1]) vs. $(channel_multiplier(cdims))")
    @assert w[M] == channels_in(cdims) DimensionMismatch("Kernel input channel count ($(w[M]) vs. $(channels_in(cdims)))")

    # Next, check that the spatial dimensions match up
    @assert x[1:M-2] == input_size(cdims) DimensionMismatch("Data input spatial size ($(x[1:M-2]) vs. $(input_size(cdims)))")
    @assert y[1:M-2] == output_size(cdims) DimensionMismatch("Data output spatial size ($(y[1:M-2]) vs. $(output_size(cdims)))")
    @assert w[1:M-2] == kernel_size(cdims) DimensionMismatch("Kernel spatial size ($(w[1:M-2]) vs. $(kernel_size(cdims)))")

    # Finally, check that the batch size matches
    @assert x[M] == y[M] DimensionMismatch("Batch size ($(x[M]) vs. $(y[M]))")
end


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module HandTests

dir = @__DIR__
hand_test_implementations = Tuple{String, Float64}[
    ("$dir/../../../../src/julia/modules/Zygote/ZygoteHAND.jl", 1e-8)
]

map!(tup -> (abspath(tup[1]), tup[2]), hand_test_implementations, hand_test_implementations)

using Test
include("../../../../src/julia/shared/load.jl")
include("../../../../src/julia/runner/load.jl")
include("load.jl")
import ADPerfTest
import TestLoader
using HandData
using TestUtils

@testset "Hand Module Tests" begin

@testset "Hand Module Test ($(basename(module_path)))" for (module_path, tolerance) in hand_test_implementations
    module_dir, module_filename = splitdir(module_path)
    module_name, module_ext = splitext(module_filename)
    # Assert that the path to the module is correct
    ext_r = @test module_ext == ".jl"
    ispath_r = @test (ispath(module_path) && !isdir(module_path))
    if isa(ext_r, Test.Pass) && isa(ispath_r, Test.Pass)
        need_modify_load_path = !(module_dir ∈ LOAD_PATH)
        if need_modify_load_path
            push!(LOAD_PATH, module_dir)
        end
        try
            @testset "Hand Simple Tests" begin
                test = TestLoader.get_hand_test(module_name)
                # Module loads
                @test isa(test, ADPerfTest.Test{HandInput, HandOutput})
                input = load_hand_input("$dir/../../../../data/hand/test.txt", false)
                output = empty_hand_output()
                test.prepare!(test.context, input)
                correct_objective = [ 1.65193147941611551e-01, -1.74542769272742593e-01, 1.54751161622253441e-01, -1.25651749731793605e-01, -4.25102935355075040e-02, -1.30665781132340175e-01 ]
                correct_jacobian = [
                    -1.95512959341282745e-02 1.18421319736437044e-02 -1.08140043450366714e-02 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -2.21537299179855318e-02 2.98023294354714807e-02 -8.89718384720838799e-03 -2.24887887916698764e-03;
                    -6.81027576369377941e-03 1.74108797580776152e-03 -6.54731655106268462e-02 0.00000000000000000e+00 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 1.42669175415592216e-02 2.58519935565392558e-02 2.30412622537484605e-03 3.09895230315512067e-03;
                    -4.09359630708754085e-02 2.57836165362687400e-02 2.95529253068245066e-03 0.00000000000000000e+00 0.00000000000000000e+00 -1.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -3.12350905279757118e-02 -9.40429974684049881e-03 -1.84753145105870202e-02 -3.09211880497457705e-04;
                    -7.10032003143440671e-02 5.14892578985870278e-02 -6.33408854067304267e-02 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -3.45326893768281365e-02 -1.95437595587887036e-02 -1.29102628630682510e-02 1.99738182488361877e-03 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00;
                    -2.81346929141111299e-02 -4.50591977492771356e-02 -8.98189200086859235e-02 0.00000000000000000e+00 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 2.42024596962444292e-02 -1.84392814127634017e-02 9.76606245804552756e-03 -1.42570914120403679e-03 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00;
                    -2.54484136783155838e-02 6.23981738234076627e-02 2.07440902671436785e-02 0.00000000000000000e+00 0.00000000000000000e+00 -1.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -3.10675692993888850e-03 6.52740542794425036e-03 8.49553578389919579e-03 -4.35490319929073976e-04 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00
                ]
                test.calculate_objective!(test.context, 1)
                test.output!(output, test.context)
                # Objective is calculated correctly
                @test correct_objective ≈ output.objective atol=tolerance
                test.calculate_objective!(test.context, 3)
                test.output!(output, test.context)
                # Objective is calculated correctly for times = 3
                @test correct_objective ≈ output.objective atol=tolerance
                test.calculate_jacobian!(test.context, 1)
                test.output!(output, test.context)
                # Jacobian is calculated correctly
                @test correct_jacobian ≈ output.jacobian atol=tolerance
                test.calculate_jacobian!(test.context, 3)
                test.output!(output, test.context)
                # Jacobian is calculated correctly for times = 3
                @test correct_jacobian ≈ output.jacobian atol=tolerance
                # Objective runs multiple times
                @test can_objective_run_multiple_times!(test.context, test.calculate_objective!)
                # Jacobian runs multiple times
                @test can_objective_run_multiple_times!(test.context, test.calculate_jacobian!)
            end

            @testset "Hand Complicated Tests" begin
                test = TestLoader.get_hand_test(module_name)
                # Module loads
                @test isa(test, ADPerfTest.Test{HandInput, HandOutput})
                input = load_hand_input("$dir/../../../../data/hand/hand_complicated.txt", true)
                output = empty_hand_output()
                test.prepare!(test.context, input)
                correct_objective = [ 1.56187661696463698e-01, -1.49300526003322220e-01, 1.72238089826454832e-01, -9.88770451849596554e-02, -1.61238035462101248e-02, -1.97586768465579654e-01 ]
                correct_jacobian = [
                    3.19994705650605837e-03 5.94121036324768426e-03 -3.16780414139280486e-02 2.51273846313546725e-02 -3.43751692490045363e-02 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -2.67300459304900973e-03 -2.05375964421033422e-04 -0.00000000000000000e+00 -0.00000000000000000e+00;
                    -5.83664316352128232e-03 -2.25087272196717869e-02 -1.17603133455238852e-02 -6.73431005210687033e-03 -7.93824304538508912e-02 0.00000000000000000e+00 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 2.51903206295353913e-03 6.35746365607068199e-03 0.00000000000000000e+00 0.00000000000000000e+00;
                    5.26539698350125818e-03 -7.98030902470658887e-04 -4.00357012848860105e-02 4.48594779455911905e-02 1.08995867664168408e-02 0.00000000000000000e+00 0.00000000000000000e+00 -1.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -1.65741767881614870e-03 2.11520284492772506e-03 -0.00000000000000000e+00 -0.00000000000000000e+00;
                    6.09513107611836524e-03 4.72340259356063275e-03 -2.21424420388885296e-02 3.97751453573731983e-02 -8.70476400071599210e-02 -1.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 -3.33849263402929107e-03 1.28209220660817119e-03 -0.00000000000000000e+00 -0.00000000000000000e+00 1.37205123955743335e-03 -5.45177232375653083e-03 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00;
                    -1.30273083127474543e-02 5.73648537699011918e-03 -9.03798668988376082e-03 -1.01785098875940283e-02 -4.37081649532052238e-02 0.00000000000000000e+00 -1.00000000000000000e+00 0.00000000000000000e+00 9.53458231126915585e-04 4.01972269593312698e-03 0.00000000000000000e+00 0.00000000000000000e+00 -9.33621518070085026e-04 -1.38222424349224984e-03 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00;
                    -1.97054905573765815e-02 8.52371287096032049e-03 1.63019777955282150e-03 7.14646422504415374e-02 2.94403102374449152e-02 0.00000000000000000e+00 0.00000000000000000e+00 -1.00000000000000000e+00 3.24947244952083786e-03 8.49942589390221431e-04 -0.00000000000000000e+00 -0.00000000000000000e+00 -5.87537490240527642e-03 2.95494237613404438e-03 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00 -0.00000000000000000e+00
                ]
                test.calculate_objective!(test.context, 1)
                test.output!(output, test.context)
                # Objective is calculated correctly
                @test correct_objective ≈ output.objective atol=tolerance
                test.calculate_objective!(test.context, 3)
                test.output!(output, test.context)
                # Objective is calculated correctly for times = 3
                @test correct_objective ≈ output.objective atol=tolerance
                test.calculate_jacobian!(test.context, 1)
                test.output!(output, test.context)
                # Jacobian is calculated correctly
                @test correct_jacobian ≈ output.jacobian atol=tolerance
                test.calculate_jacobian!(test.context, 3)
                test.output!(output, test.context)
                # Jacobian is calculated correctly for times = 3
                @test correct_jacobian ≈ output.jacobian atol=tolerance
                # Objective runs multiple times
                @test can_objective_run_multiple_times!(test.context, test.calculate_objective!)
                # Jacobian runs multiple times
                @test can_objective_run_multiple_times!(test.context, test.calculate_jacobian!)
            end
        finally
            pop!(LOAD_PATH)
        end
    end
end

end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
"""
    Ensembles

This module provides the main functions [`run_ensemble`](@ref Ensembles.run_ensemble).
This serves to run multiple trajectories for a given simulation type, sampling
from an initial distribution.
"""
module Ensembles

using SciMLBase: remake, EnsembleThreads, EnsembleProblem, solve
using RecursiveArrayTools: ArrayPartition
using TypedTables: TypedTables

using NQCBase: NQCBase
using NQCDynamics:
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    DynamicsUtils,
    DynamicsMethods,
    RingPolymers,
    NonadiabaticDistributions
using NQCDynamics.NonadiabaticDistributions:
    NuclearDistribution,
    CombinedDistribution

export run_ensemble

function sample_distribution(sim::AbstractSimulation, distribution::NuclearDistribution, i)
    u = getindex(distribution, i)
    DynamicsMethods.DynamicsVariables(sim, u.v, u.r)
end

function sample_distribution(sim::RingPolymerSimulation{<:DynamicsMethods.ClassicalMethods.ThermalLangevin}, distribution::NuclearDistribution, i)
    u = getindex(distribution, i)
    DynamicsMethods.DynamicsVariables(sim, RingPolymers.RingPolymerArray(u.v), RingPolymers.RingPolymerArray(u.r))
end

function sample_distribution(sim::AbstractSimulation, distribution::CombinedDistribution, i)
    u = getindex(distribution.nuclear, i)
    DynamicsMethods.DynamicsVariables(sim, u.v, u.r, distribution.electronic)
end

include("selections.jl")
include("outputs.jl")
include("reductions.jl")

"""
    run_ensemble(sim::AbstractSimulation, tspan, distribution;
        selection=nothing,
        output=(sol,i)->(sol,false),
        reduction=(u,data,I)->(append!(u,data),false),
        ensemble_algorithm=EnsembleThreads(),
        algorithm=DynamicsMethods.select_algorithm(sim),
        kwargs...
        )

Run multiple trajectories for timespan `tspan` sampling from `distribution`.
The DifferentialEquations ensemble interface is used which allows us to specify
functions to modify the output and how it is reduced across trajectories.

# Keywords

* `selection` should be an `AbstractVector` containing the indices to sample from the `distribution`
By default, `nothing` leads to random sampling.
* `output` should be a function as described in the DiffEq docs that takes the solution and returns the output.
* `algorithm` is the algorithm used to integrate the equations of motion.
* `ensemble_algorithm` is the algorithm from DifferentialEquations which determines which
form of parallelism is used.

This function wraps the EnsembleProblem from DifferentialEquations and passes the `kwargs`
to the `solve` function.
"""
function run_ensemble(
    sim::AbstractSimulation, tspan, distribution;
    selection::Union{Nothing,AbstractVector}=nothing,
    output=(sol,i)->(sol,false),
    reduction::Symbol=:append,
    ensemble_algorithm=EnsembleThreads(),
    algorithm=DynamicsMethods.select_algorithm(sim),
    saveat=[],
    trajectories=1,
    kwargs...
)

    if output isa Symbol
        output = (output,)
    end
    trajectories = convert(Int, trajectories)

    output_func = get_output_func(output)

    kwargs = NQCBase.austrip_kwargs(;kwargs...)
    saveat = austrip.(saveat)

    u0 = sample_distribution(sim, distribution, 1)
    problem = DynamicsMethods.create_problem(u0, austrip.(tspan), sim)

    reduction = select_reduction(reduction)
    u_init = get_u_init(reduction, saveat, kwargs, tspan, u0, output_func)
    prob_func = choose_selection(distribution, selection, sim, output, saveat, trajectories)

    ensemble_problem = EnsembleProblem(problem; prob_func, output_func, reduction, u_init)
    sol = solve(ensemble_problem, algorithm, ensemble_algorithm; u_init, saveat, trajectories, kwargs...)

    if output isa Tuple
        return [TypedTables.Table(
                    t=vals.t,
                    [(;zip(output, val)...) for val in vals.saveval]
                ) for vals in prob_func.values
            ]
    else
        return sol.u
    end
end

function choose_selection(distribution, selection, sim, output::Tuple, saveat, trajectories)
    selection = Selection(distribution, selection)
    return SelectWithCallbacks(selection, DynamicsMethods.get_callbacks(sim), output, trajectories, saveat=saveat)
end

function choose_selection(distribution, selection, sim, output, saveat, trajectories)
    return Selection(distribution, selection)
end

get_output_func(::Tuple) = (sol,i)->(sol,false)
get_output_func(output) = output

end # module


text:
function wrap_macro2enum!(ctx::AbstractContext,
                          cursors::Vector{CLMacroDefinition},
                          constant_type::String)

    enum_sym = symbol_safe(constant_type)
    enum_type = Cint
    name2value = Tuple{Symbol,enum_type}[]

    for cursor in cursors
        if constant_type == "DAQmxAttribute"
            item_name = spelling(cursor)[7:end]
        elseif constant_type == "DAQmxConstant"
            item_name = spelling(cursor)[length("DAQmx_Val_")+1:end]
        else
            item_name = prod(split(name(cursor), "_")[3:end])
        end
        if startswith(item_name, r"[0-9]")
            item_name = "_" * item_name
        end
        item_sym = symbol_safe(item_name)
        tokens = tokenize(cursor)

        # Macros definitions include bitshift operators, which we have to check for...
        for i in 1:length(tokens)
            if typeof(tokens[i]) == Literal && (i == length(tokens) || tokens[i+1].text ∉ ["<<",">>"])
                push!(name2value, (item_sym, parse(enum_type, tokens[i].text)))
                break
            elseif typeof(tokens[i]) == Literal && i < length(tokens) && tokens[i+1].text ∈ ["<<",">>"]
                bitshift = tokens[i].text * tokens[i + 1].text * tokens[i + 2].text
                push!(name2value, (item_sym, eval(Meta.parse(bitshift))))
                break
            end
        end
    end

    expr = Expr(:macrocall, Symbol("@cenum"), nothing, Expr(:(::), enum_sym, enum_type))
    enum_pairs = Expr(:block)
    ctx.common_buffer[enum_sym] = ExprUnit(expr)

    for (name,value) in name2value
        push!(enum_pairs.args, :($name = $value))
    end

    push!(expr.args, enum_pairs)

    return ctx
end

function wrap_dealias!(ctx::AbstractContext, cursor::CLFunctionDecl)
    func_type = type(cursor)
    func_name = isempty(ctx.force_name) ? Symbol(spelling(cursor)) : ctx.force_name
    ret_type = :Cint
    args = function_args(cursor)
    arg_types = [argtype(func_type, i) for i in 0:length(args)-1]
    arg_reps = clang2julia.(arg_types)

    for (i, arg) in enumerate(arg_reps)
        if isa(arg, Symbol)
            if String(arg) ∈ keys(typedef_map)
                arg_reps[i] = Symbol(typedef_map[String(arg)])
            end
        end
        if Meta.isexpr(arg, :curly) && first(arg.args) == :NTuple
            if String(last(arg.args)) ∈ keys(typedef_map)
                arg_reps[i] = Expr(:curly, :Ref, Symbol(typedef_map[String(last(arg.args))]))
            else
                arg_reps[i] = Expr(:curly, :Ref, last(arg.args))
            end
        end
        if Meta.isexpr(arg, :curly) && first(arg.args) == :Ptr
            if String(last(arg.args)) ∈ keys(typedef_map)
                arg_reps[i] = Expr(:curly, :Ref, Symbol(typedef_map[String(last(arg.args))]))
            else
                arg_reps[i] = Expr(:curly, :Ref, last(arg.args))
            end
        end
    end

    arg_count = 0
    arg_names = map(args) do x
        n = name_safe(name(x))
        s = !isempty(n) ? n : "arg"*string(arg_count+=1)
        Symbol(s)
    end

    signature = efunsig(Symbol(spelling(cursor)[6:end]), arg_names, arg_reps)
    body = eccall(func_name, Symbol(ctx.libname), ret_type, arg_names, arg_reps)
    push!(ctx.api_buffer, Expr(:function, signature, Expr(:block, body)))

    return ctx
end

function eccall(func_name::Symbol, libname::Symbol, rtype, args, types)
    :(ccall(($(QuoteNode(func_name)), $(QuoteNode(libname))),
        $rtype,
        $(Expr(:tuple, types...)),
        $(args...))
     )
end

function is_ptr_type_expr(@nospecialize t)
    (t === :Cstring || t === :Cwstring) && return true
    isa(t, Expr) || return false
    t = t::Expr
    t.head === :curly && t.args[1] === :Ptr
end

function efunsig(name::Symbol, args::Vector{Symbol}, types)
    x = [is_ptr_type_expr(t) ? a : Expr(:(::), a, t) for (a,t) in zip(args,types)]
    Expr(:call, name, x...)
end

function get_comment(cursor::CLCursor)
    stringptr = clang_Cursor_getRawCommentText(cursor.cursor)
    if stringptr.data == C_NULL
        return nothing
    else
        cstring = clang_getCString(stringptr)
        return unsafe_string(cstring)
    end
end

function sort_cursors(cursors::Vector{CLMacroDefinition})
    cursor_names = spelling.(cursors)
    cursor_dict = Dict(zip(cursor_names, cursors))
    sort!(unique!(cursor_names))
    sorted_cursors = [cursor_dict[cursor_names[j]] for j in 1:length(cursor_names)]
    return sorted_cursors
end


text:
"""
An iterative exchange algorithm for penalized least squares regression.
The exchange algorithm operates similar to cyclic coordinate descent,
except that it swaps (exchanges) predictors to improve the residual sum of squares.
"""
module ExchangeLeastsq

#using GLMNet
using RegressionTools
using PLINK
using Distances: euclidean, chebyshev
using StatsBase: sample
using OpenCL
using DataFrames
using Gadfly

export exlstsq
export cv_exlstsq

#"An alias for the `OpenCL` module name"
#const cl = OpenCL

# typealias for Float32/Float64
typealias Float Union{Float32, Float64}

include("common.jl")
include("gpu.jl")
include("gwas.jl")
include("cv.jl")
include("exlstsq.jl")
#include("test.jl")

end # end module ExchangeLeastsq


text:
module BigIntegers
    mutable struct BigInteger
        str::AbstractString
        positive::Bool

        BigInteger(str::AbstractString, positive::Bool) = new(str == "" ? "0" : str, positive)
        BigInteger(x::Int) = new(string(abs(x)), x ≥ 0)
        BigInteger(x::BigInteger) = copy(x)
        function BigInteger(x::AbstractString)
            if x == ""
                new("0", true)
            else
                if x[1] == '-'
                    new(x[2:end], false)
                else
                    new(x, true)
                end
            end
        end
    end

    pad(s::AbstractString, n::Int, right::Bool = true)::String = right ? string(s, repeat("0", n)) : string(repeat("0", n), s)

    function align_strings(a::AbstractString, b::AbstractString, right::Bool = true)::Tuple{String, String}
        if length(a) > length(b)
            b = pad(b, length(a) - length(b), right)
        elseif length(b) > length(a)
            a = pad(a, length(b) - length(a), right)
        end
        return a, b
    end

    function mul_s(a::BigInteger, b::BigInteger)::BigInteger
        if length(a) + length(b) < 19
            return BigInteger(string(parse(Int, a.str) * parse(Int, b.str)), a.positive == b.positive)
        end
        g, h = reverse(a.str), reverse(b.str)
        result = BigInteger(0)
        for i ∈ 1:length(g)
            carry = 0
            sub_result = repeat('0', i - 1)
            for j ∈ 1:length(h)
                carry += parse(Int, g[i]) * parse(Int, h[j])
                sub_result = string(carry % 10, sub_result)
                carry ÷= 10
            end
            result += BigInteger(lstrip(string(carry, sub_result), '0'), true)
        end
        result.positive = a.positive == b.positive
        return result
    end

    Base.show(io::IO, x::BigInteger)::Nothing = show(io, string(x.positive ? "" : "-", x.str))

    Base.copy(x::BigInteger)::BigInteger = BigInteger(x.str, x.positive)

    Base.length(x::BigInteger)::Int = length(x.str)

    Base.getindex(x::BigInteger, i::Int)::Char = getindex(x.str, i)

    Base.getindex(x::BigInteger, r::UnitRange{Int})::SubString = getindex(x.str, r)

    Base.lastindex(x::BigInteger)::Int = lastindex(x.str)

    Base.:-(x::BigInteger)::BigInteger = BigInteger(x.str, !x.positive)

    Base.abs(x::BigInteger)::BigInteger = BigInteger(x.str, true)

    Base.isless(x::BigInteger, y::BigInteger)::Bool = x < y

    Base.hash(x::BigInteger)::UInt = hash(x.str)

    function Base.:<(x::BigInteger, y::BigInteger)::Bool
        if x.positive == y.positive
            if length(x) == length(y)
                if x.str == y.str
                    return false
                else
                    for i ∈ 1:length(x)
                        if x[i] == y[i]
                            continue
                        else
                            return !((x[i] < y[i]) ⊻ x.positive)
                        end
                    end
                end
            else
                return !((length(x) < length(y)) ⊻ x.positive)
            end
        else
            return !x.positive
        end
    end

    Base.:<(x::BigInteger, y::Int)::Bool = x < BigInteger(y)
    Base.:<(x::Int, y::BigInteger)::Bool = BigInteger(x) < y

    Base.:(==)(x::BigInteger, y::BigInteger)::Bool = x.str == y.str && x.positive == y.positive
    Base.:(==)(x::BigInteger, y::Int)::Bool = x == BigInteger(y)

    function Base.:+(x::BigInteger, y::BigInteger)::BigInteger
        if x.positive == y.positive
            a, b = align_strings(x.str, y.str, false)
            result = ""
            carry = 0
            for i ∈ 1:Int(ceil(length(a) / 18))
                carry += parse(Int, a[max(end - 18 * i + 1, 1):end - 18 * (i - 1)]) + parse(Int, b[max(end - 18 * i + 1, 1):end - 18 * (i - 1)])
                result = string(repeat('0', max(18 - length(string(carry)), 0)), string(carry)[max(end - 17, 1):end], result)
                carry ÷= 10 ^ 18
            end
            if carry == 0
                result = lstrip(result, '0')
            else
                result = string(carry, result)
            end
            return BigInteger(result, x.positive)
        else
            if x.positive
                return x - -y
            else
                return y - -x
            end
        end
    end

    Base.:+(x::BigInteger, y::Int)::BigInteger = x + BigInteger(y)
    Base.:+(x::Int, y::BigInteger)::BigInteger = BigInteger(x) + y

    function Base.:-(x::BigInteger, y::BigInteger)::BigInteger
        if x.positive == y.positive
            if x == y
                return BigInteger(0)
            else
                if x.positive
                    if x > y
                        positive = true
                        a, b = x.str, y.str
                    else
                        positive = false
                        a, b = y.str, x.str
                    end
                else
                    if x > y
                        positive = true
                        a, b = y.str, x.str
                    else
                        positive = false
                        a, b = x.str, y.str
                    end
                end
                a, b = align_strings(a, b, false)
                result = ""
                carry = 0
                for i ∈ 1:Int(ceil(length(a) / 18))
                    carry += parse(Int, a[max(end - 18 * i + 1, 1):end - 18 * (i - 1)]) - parse(Int, b[max(end - 18 * i + 1, 1):end - 18 * (i - 1)])
                    if carry < 0
                        result = string(repeat('0', 18 - length(string(10 ^ 18 + carry))), 10 ^ 18 + carry, result)
                        carry = -1
                    else
                        result = string(repeat('0', 18 - length(string(carry))), carry, result)
                        carry = 0
                    end
                end
                result = lstrip(result, '0')
                return BigInteger(result, positive)
            end
        else
            if x.positive
                return x + -y
            else
                return -(-x + y)
            end
        end
    end

    Base.:-(x::BigInteger, y::Int)::BigInteger = x - BigInteger(y)
    Base.:-(x::Int, y::BigInteger)::BigInteger = BigInteger(x) - y

    function Base.:*(x::BigInteger, y::BigInteger)::BigInteger
        x_len, y_len = length(x), length(y)
        if x_len + y_len < 19 || x_len == 1 || y_len == 1 return mul_s(x, y) end
        n = x_len > y_len ? x_len : y_len
        m = n ÷ 2
        a, b = BigInteger(x[1:end - m], true), BigInteger(x[max(end - m + 1, 1):end], true)
        c, d = BigInteger(y[1:end - m], true), BigInteger(y[max(end - m + 1, 1):end], true)
        ac, bd, ab_cd = a * c, b * d, (a + b) * (c + d)
        r = BigInteger(pad(ac.str, 2 * m), true)
        s = BigInteger(pad((ab_cd - ac - bd).str, m), true)
        res = r + s + bd
        res.positive = x.positive == y.positive
        return res
    end

    Base.:*(x::BigInteger, y::Int)::BigInteger = x * BigInteger(y)
    Base.:*(x::Int, y::BigInteger)::BigInteger = BigInteger(x) * y

    function Base.:÷(x::BigInteger, y::BigInteger)::BigInteger
        if abs(x) < abs(y)
            result = BigInteger(0)
        else
            value = copy(y)
            multipliers = Dict(BigInteger(1) => copy(value))
            i = BigInteger(1)
            while abs(x) > abs(value)
                value *= 2
                i *= 2
                multipliers[i] = copy(value)
            end
            value = copy(multipliers[BigInteger(1)])
            result = BigInteger(1)
            for i ∈ reverse(sort(collect(keys(multipliers))))
                if value + multipliers[i] ≤ abs(x)
                    result += i
                    value += multipliers[i]
                end
            end
        end
        return x.positive == y.positive ? result : -(result + 1)
    end

    Base.:÷(x::BigInteger, y::Int)::BigInteger = x ÷ BigInteger(y)

    function Base.:%(x::BigInteger, y::BigInteger)::BigInteger
        if x.str == y.str
            return BigInteger(0)
        elseif abs(x) > abs(y)
            return x - x ÷ y * y
        elseif x.positive == y.positive
            return copy(x)
        else
            return x + y
        end
    end

    Base.:%(x::BigInteger, y::Int)::BigInteger = x % BigInteger(y)

    function Base.:^(x::BigInteger, p::BigInteger)
        value = copy(x)
        result = BigInteger(1)
        while p > 0
            if p % 2 == 1
                result *= value
            end
            value *= value
            p ÷= 2
        end
        return result
    end

    Base.:^(x::BigInteger, p::Int) = x ^ BigInteger(p)
    Base.:^(x::Int, p::BigInteger) = BigInteger(x) ^ p

    function Base.powermod(x::BigInteger, p::BigInteger, m::BigInteger)
        value = copy(x)
        result = BigInteger(1)
        while p > 0
            if p % 2 == 1
                result *= value
                result %= m
            end
            value *= value
            value %= m
            p ÷= 2
        end
        return result
    end

    Base.powermod(x::BigInteger, p::Int, m::BigInteger) = powermod(x, BigInteger(p), m)
    Base.powermod(x::BigInteger, p::BigInteger, m::Int) = powermod(x, p, BigInteger(m))
    Base.powermod(x::BigInteger, p::Int, m::Int) = powermod(x, BigInteger(p), BigInteger(m))
    Base.powermod(x::Int, p::BigInteger, m::BigInteger) = powermod(BigInteger(x), p, m)
    Base.powermod(x::Int, p::Int, m::BigInteger) = powermod(BigInteger(x), BigInteger(p), m)
    Base.powermod(x::Int, p::BigInteger, m::Int) = powermod(BigInteger(x), p, BigInteger(m))
end

module Calculus
    """
        sum_arithmetic_series(a::Real, d::Real, n::Int)

    Get sum of first `n` terms of arithmetic series

    # Arguments
    - `n::Real`: First number in series
    - `d::Real`: Common difference of series
    - `n::Int`: Number of terms in series

    # Examples
    ```
    julia> sum_arithmetic_series(1, 2, 3)
    9.0
    julia> sum_arithmetic_series(0.5, 0.75, 3)
    3.75
    ```
    """
    sum_arithmetic_series(a::Real, d::Real, n::Int)::Float64 = 0.5n * (2a + d * (n - 1))

    """
        sum_geometric_series(a::Real, r::Real, n::Int)

    Get sum of `n` first terms of geometric series

    # Arguments

    - `a::Real`: First number in series
    - `r::Real`: Common ratio of series
    - `n::Int`: Number of terms in series

    # Examples
    ```
    julia> sum_geometric_series(1, 2, 3)
    7.0
    julia> sum_geometric_series(0.5, 0.75, 3)
    1.15625
    ```
    """
    sum_geometric_series(a::Real, r::Real, n::Int)::Float64 = a * (r ^ n - 1) / (r - 1)

    """
        fibonacci_numbers(n::Int, T::Type=Int)

    Generate `n` first Fibonacci numbers [OEIS A000045](https://oeis.org/A000045)

    # Arguments

    - `n::Int`: Number of Fibonacci numbers
    - `T::Type`: Type of numbers

    # Examples
    ```
     julia> fibonacci_numbers(10)
     10-element Vector{Int64}:
      1
      1
      2
      3
      5
      8
     13
     21
     34
     55
    julia> fibonacci_numbers(10, BigInteger)
    11-element Vector{BigInteger}:
     "1"
     "1"
     "2"
     "3"
     "5"
     "8"
     "13"
     "21"
     "34"
     "55"
    ```
    """
    function fibonacci_numbers(n::Int, T::Type=Int)::Vector{T}
        f_prev, f_curr = T(1), T(1)
        result = [f_prev, f_curr]
        for i ∈ 1:n - 2
            f_prev, f_curr = f_curr, f_prev + f_curr
            push!(result, f_curr)
        end
        return result
    end

    """
        get_fibonacci(n::Int)

    Get `n`-th Fibonacci number

    # Arguments

    - `n::Int`: Index of Fibonacci number

    # Examples
    ```
    julia> get_fibonacci(10)
    55
    ```
    """
    get_fibonacci(n::Int)::Int = trunc(Int, ((1 + √5) / 2) ^ n / √5 + 0.5)

    """
        fibonacci_index(n::Int)

    Get index of `n` in Fibonacci sequence

    # Examples
    ```
    julia> fibonacci_index(55)
    10
    julia> fibonacci_index(54)
    9
    ```
    """
    fibonacci_index(n::Int)::Int = trunc(Int, log((1 + √5) / 2, √5n + 0.5))

    """
        len_collatz_chain(n::Int, hashmap::Dict{Int, Int})

    Get length of Collatz chain

    # Arguments
    - `n::Int`: Starting number
    - `hashmap::Dict{Int, Int}`: Collatz chains lengths

    # Examples
    ```
    julia> len_collatz_chain(13, Dict(1 => 1))
    10
    ```
    """
    function len_collatz_chain(n::Int, hashmap::Dict{Int, Int})::Int
        if !haskey(hashmap, n)
            if n % 2 == 0
                hashmap[n] = 1 + len_collatz_chain(n ÷ 2, hashmap)
            else
                hashmap[n] = 2 + len_collatz_chain((3n + 1) ÷ 2, hashmap)
            end
        end
        return hashmap[n]
    end

    """
        triangular_numbers(n::Int)

    Generate `n` first numbers in triangular numbers sequence [OEIS A000217](<https://oeis.org/A000217>)

    # Arguments

    - `n::Int`: Number of triangular numbers

    # Examples
    ```
    julia> triangular_numbers(10)
    10-element Vector{Int64}:
      1
      3
      6
     10
     15
     21
     28
     36
     45
     55
    ```
    """
    function triangular_numbers(n::Int)::Vector{Int}
        result = [1]
        for i ∈ 2:n
            push!(result, result[end] + i)
        end
        return result
    end

    """
        get_triangular(n::Int)

    Get `n`-th triangular number

    # Arguments

    - `n::Int`: Index of triangular number

    # Examples
    ```
    julia> get_triangular(10)
    55
    ```
    """
    get_triangular(n::Int)::Int = n * (n + 1) ÷ 2

    """
        triangular_index(n::Int)

    Get index of `n` in triangular sequence

    # Examples
    ```
    julia> triangular_index(55)
    10
    julia> triangular_index(54)
    9
    ```
    """
    triangular_index(n::Int)::Int = (√(8n + 1) - 1) ÷ 2

    """
        is_triangular(n::Int)

    Check if number is triangular

    # Examples
    ```
    julia> is_triangular(55)
    true
    julia> is_triangular(54)
    false
    ```
    """
    is_triangular(n::Int)::Bool = √(8n + 1) % 2 == 1

    """
        square_numbers(n::Int)

    Generate `n` first numbers in square numbers sequence [OEIS A000290](<https://oeis.org/A000290>)

    # Arguments

    - `n::Int`: Number of square numbers

    # Examples
    ```
    julia> square_numbers(10)
    10-element Vector{Int64}:
       1
       4
       9
      16
      25
      36
      49
      64
      81
     100
    ```
    """
    square_numbers(n::Int)::Vector{Int} = (1:n) .^ 2

    """
        get_square(n::Int)

    Get `n`-th square number

    # Arguments

    - `n::Int`: Index of square number

    # Examples
    ```
    julia> get_square(10)
    100
    ```
    """
    get_square(n::Int)::Int = n * n

    """
        square_index(n::Int)

    Get index of `n` in square sequence

    # Examples
    ```
    julia> square_index(100)
    10
    julia> square_index(99)
    9
    ```
    """
    square_index(n::Int)::Int = trunc(Int, √n)

    """
        is_square(n::Int)

    Check if number is square

    # Examples
    ```
    julia> is_square(100)
    true
    julia> is_square(99)
    false
    ```
    """
    is_square(n::Int)::Bool = √n % 1 == 0

        """
        pentagonal_numbers(n::Int)

    Generate `n` first numbers in pentagonal numbers sequence [OEIS A000326](<https://oeis.org/A000326>)

    # Arguments

    - `n::Int`: Number of pentagonal numbers

    # Examples
    ```
    julia> pentagonal_numbers(10)
    10-element Vector{Int64}:
       1
       5
      12
      22
      35
      51
      70
      92
     117
     145
    ```
    """
    function pentagonal_numbers(n::Int)::Vector{Int}
        result = [1]
        for i ∈ 4:3:3n
            push!(result, result[end] + i)
        end
        return result
    end

    """
        get_pentagonal(n::Int)

    Get `n`-th pentagonal number

    # Arguments

    - `n::Int`: Index of pentagonal number

    # Examples
    ```
    julia> get_pentagonal(10)
    145
    ```
    """
    get_pentagonal(n::Int)::Int = n * (3 * n - 1) ÷ 2

    """
        pentagonal_index(n::Int)

    Get index of `n` in pentagonal sequence

    # Examples
    ```
    julia> pentagonal_index(145)
    10
    julia> pentagonal_index(144)
    9
    ```
    """
    pentagonal_index(n::Int)::Int = (√(24n + 1) + 1) ÷ 6

    """
        is_pentagonal(n::Int)

    Check if number is pentagonal

    # Examples
    ```
    julia> is_pentagonal(145)
    true
    julia> is_pentagonal(144)
    false
    ```
    """
    is_pentagonal(n::Int)::Bool = √(24n + 1) % 6 == 5

        """
        hexagonal_numbers(n::Int)

    Generate `n` first numbers in hexagonal numbers sequence [OEIS A000384](<https://oeis.org/A000384>)

    # Arguments

    - `n::Int`: Number of hexagonal numbers

    # Examples
    ```
    julia> hexagonal_numbers(10)
    10-element Vector{Int64}:
       1
       6
      15
      28
      45
      66
      91
     120
     153
     190
    ```
    """
    function hexagonal_numbers(n::Int)::Vector{Int}
        result = [1]
        for i ∈ 5:4:4n
            push!(result, result[end] + i)
        end
        return result
    end

    """
        get_hexagonal(n::Int)

    Get `n`-th hexagonal number

    # Arguments

    - `n::Int`: Index of hexagonal number

    # Examples
    ```
    julia> get_hexagonal(10)
    190
    ```
    """
    get_hexagonal(n::Int)::Int = n * (2n - 1)

    """
        hexagonal_index(n::Int)

    Get index of `n` in hexagonal sequence

    # Examples
    ```
    julia> hexagonal_index(190)
    10
    julia> hexagonal_index(189)
    9
    ```
    """
    hexagonal_index(n::Int)::Int = (√(8n + 1) + 1) ÷ 4

    """
        is_hexagonal(n::Int)

    Check if number is hexagonal

    # Examples
    ```
    julia> is_hexagonal(190)
    true
    julia> is_hexagonal(189)
    false
    ```
    """
    is_hexagonal(n::Int)::Bool = √(8n + 1) % 4 == 3

    """
        champernowne_digits(n::Int)

    Get number of digits formed by `n`-digit numbers in Champernowne constant

    # Examples
    julia> champernowne_digits(3)
    2889
    """
    champernowne_digits(n::Int)::Int = (10 ^ n * (9n - 1) + 1) ÷ 9
end

module Math
    using ..BigIntegers: BigInteger

    """
        fast_max(a::Int, b::Int)

    Get max between two numbers

    # Examples:
    julia> fast_max(100, 50)
    100
    """
    fast_max(a::Int, b::Int)::Int = a > b ? a : b

    """
        fast_max(a::BigInteger, b::BigInteger)

    Get max between two numbers

    # Examples:
    julia> fast_max(BigInteger(48), BigInteger(49))
    "49"
    """
    fast_max(a::BigInteger, b::BigInteger)::BigInteger = a > b ? a : b

    """
        fast_min(a::Int, b::Int)

    Get min between two numbers

    # Examples:
    julia> fast_min(100, 50)
    50
    """
    fast_min(a::Int, b::Int)::Int = a < b ? a : b

    """
        fast_min(a::BigInteger, b::BigInteger)

    Get min between two numbers

    # Examples:
    julia> fast_min(BigInteger(48), BigInteger(49))
    "48"
    """
    fast_min(a::BigInteger, b::BigInteger)::BigInteger = a < b ? a : b
end

module Grid
    using ..Math: fast_max

    """
        adjacent_product(grid::Vector{Vector{Int}}, n::Int, coordinate::Tuple{Int, Int},
                         direction::Tuple{Int, Int})

    Get product of `n` adjacent numbers in `grid`

    # Arguments:
    - `grid::Vector{Vector{Int}}`: Matrix-like array
    - `n::Int`: Number of adjacent numbers in product
    - `coordinate::Tuple{Int, Int}`: Start position in product
    - `direction::Tuple{Int, Int}`: Direction for getting adjacent numbers (vertical, horizontal, diagonal)

    # Examples:
    julia> grid = [[8, 2, 22, 97], [49, 49, 99, 40], [81, 49, 31, 73], [52, 70, 95, 23]]
    4-element Vector{Vector{Int64}}:
     [8, 2, 22, 97]
     [49, 49, 99, 40]
     [81, 49, 31, 73]
     [52, 70, 95, 23]
    Horizontal (8 * 2 * 22 * 97):
    julia> adjacent_product(grid, 4, (1, 1), (1, 0))
    34144
    Vertical (97 * 40 * 73):
    julia> adjacent_product(grid, 3, (4, 1), (0, 1))
    283240
    Up-left to down-right diagonal (49 * 31 * 23)
    julia> adjacent_product(grid, 3, (2, 2), (1, 1))
    34937
    Up-right to down-left diagonal (97 * 99)
    julia> adjacent_product(grid, 2, (4, 1), (-1, 1))
    9603
    """
    function adjacent_product(grid::Vector{Vector{Int}}, n::Int, coordinate::Tuple{Int, Int},
                              direction::Tuple{Int, Int})::Int
        product = 1
        for i ∈ 0:n - 1
            number = grid[coordinate[2] + direction[2] * i][coordinate[1] + direction[1] * i]
            if number == 0
                return 0
            end
            product *= number
        end
        return product
    end

    """
        grid::Vector{Vector{Int}}

    Get triangle maximum path sum

    # Arguments:
    - `grid::Vector{Vector{Int}}`: Triangle-like list

    # Examples:
    julia> max_path_sum([[3], [7, 4], [2, 4, 6], [8, 5, 9, 3]])
    23
    """
    function max_path_sum(grid::Vector{Vector{Int}})::Int
        while length(grid) > 1
            for i ∈ 1:length(grid) - 1
                grid[end - 1][i] += fast_max(grid[end][i:i + 1]...)
            end
            grid = grid[1:end - 1]
        end
        return grid[1][1]
    end
end

module Numbers
    using ..BigIntegers: BigInteger
    using ..Calculus: sum_geometric_series
    """
        is_palindrome(n::Int)

    Check if Integer number is palindrome

    # Arguments
    - `n::Int`: Number

    # Examples
    ```
    julia> is_palindrome(123321)
    true
    julia> is_palindrome(123456)
    false
    ```
    """
    function is_palindrome(n::Int)::Bool
        n = string(n)
        return n == reverse(n)
    end

    """
        is_palindrome(n::BigInteger)

    Check if BigInteger number is palindrome

    # Arguments
    - `n::BigInteger`: Number

    # Examples
    ```
    julia> is_palindrome(BigInteger(123321))
    true
    julia> is_palindrome(BigInteger(123456))
    false
    ```
    """
    is_palindrome(n::BigInteger)::Bool = n.str == reverse(n.str)

    """
        count_divisors(n::Int)

    Get number of divisors of `n`

    # Arguments
    - `n::Int`: Number

    # Examples
    ```
    julia> count_divisors(23)
    2
    julia> count_divisors(284)
    6
    ```
    """
    function count_divisors(n::Int)::Int
        divisors = 1
        i = 2
        while i * i ≤ n
            f = 1
            while n % i == 0
                n /= i
                f += 1
            end
            i += 1
            divisors *= f
        end
        if n > 1; divisors *= 2 end
        return divisors
    end

    """
        sum_proper_factors(n::Int, primes::Vector{Int} = prime_numbers(trunc(Int, √n)))

    Get sum of proper divisors of `n`

    # Arguments
    - `n::Int`: Number
    - `primes::Vector{Int} = prime_numbers(trunc(Int, √n))`: Vector of prime numbers below at least √n

    # Examples
    ```
    julia> sum_proper_factors(23)
    2
    julia> sum_proper_factors(284, prime_numbers(trunc(Int, √284)))
    220
    ```
    """
    function sum_proper_divisors(n::Int, primes::Vector{Int} = prime_numbers(trunc(Int, √n)))::Int
        result = 1
        number = n
        for prime ∈ primes
            j = 1
            while number % prime == 0
                number ÷= prime
                j += 1
            end
            result *= sum_geometric_series(1, prime, j)
        end
        if number > 1; result *= number + 1 end
        return result - n
    end

    """
        get_digits(n::Int)

    Get array of digits in number

    # Arguments:
    - `n::Int`: Number

    # Examples
    ```
    julia> get_digits(123)
    3-element Vector{Int64}:
     1
     2
     3
    ```
    """
    get_digits(n::Int)::Vector{Int} = parse.(Int, collect(string(n)))

    """
        get_digits(n::BigInteger)

    Get array of digits in number

    # Arguments:
    - `n::BigInt`: Number

    # Examples
    ```
    julia> get_digits(BigInteger(123456))
    6-element Vector{Int64}:
     1
     2
     3
     4
     5
     6
    ```
    """
    get_digits(n::BigInteger)::Vector{Int} = parse.(Int, collect(n.str))

    """
        sum_digits(n::Int)

    Get sum of digits in number

    # Arguments

    - `n::Int`: Number

    # Examples
    ```
    julia> sum_digits(123)
    6
    ```
    """
    digits_sum(n::Int)::Int = sum(get_digits(n))

    """
        sum_digits(n::BigInteger)

    Get sum of digits in number

    # Arguments

    - `n::BigInt`: Number

    # Examples
    ```
    julia> sum_digits(BigInteger(123456))
    21
    ```
    """
    digits_sum(n::BigInteger)::Int = sum(get_digits(n))

    """
        multiplicative_order(a::Int, n::Int)

    Get multiplicative order of `n` modulo `m`

    # Arguments:
    - `n::Int`: Number
    - `m::Int`: Number

    # Examples
    ```
    multiplicative_order(16, 7)
    3
    ```
    """
    function multiplicative_order(n::Int, m::Int)::Int
        i, k = 1, n
        while k ≠ 1 && i < m
            k = k * n % m
            i += 1
        end
        return i
    end

    """
        is_s_number(a::Int, b::Int)

    Check if number is S-number, that is a perfect square and its square root can be obtained by
    splitting the decimal representation of `n` into 2 or more numbers then adding the numbers.

    # Arguments
    - `a::Int`: Number
    - `b::Int`: Square of a

    # Examples
    julia> is_s_number(99, 9801)
    true
    julia> is_s_number(5, 25)
    false
    """
    function is_s_number(a::Int, b::Int)::Bool
        if a > b return false elseif a == b return true end
        m = 10
        while m < b
            quotient, remainder = b ÷ m, b % m
            if remainder < a && is_s_number(a - remainder, quotient) return true end
            m *= 10
        end
        return false
    end
end

module Primes
    """
        prime_sieve(n::Int)

    Get the sieve of primes below `n`

    # Arguments:
    - `n::Int`: Index of last prime in list

    # Examples:
    ```
    julia> prime_sieve(10)
    10-element Vector{Bool}:
     0
     1
     1
     0
     1
     0
     1
     0
     0
     0
    ```
    """
    function prime_sieve(n::Int)::Vector{Bool}
        sieve = repeat([true], n)
        sieve[1] = false
        sieve[4:2:n] .= false
        for i ∈ 3:2:trunc(Int, √n)
            if sieve[i]; sieve[i*i:2i:n] .= false end
        end
        return sieve
    end

    """
        prime_numbers(n::Int)

    Generate prime numbers below `n`

    # Arguments
    - `n::Int`: Number below of which prime numbers will generate

    # Examples
    ```
    julia> prime_numbers(20)
    8-element Vector{Int64}:
      2
      3
      5
      7
     11
     13
     17
     19
    ```
    """
    prime_numbers(n::Int)::Vector{Int} = findall(prime_sieve(n))

    """
        is_prime(n::Int)

    Check if number is prime

    # Arguments
    - `n::Int`: Number

    # Examples
    ```
    julia> is_prime(23)
    true
    julia> is_prime(22)
    false
    ```
    """
    function is_prime(n::Int)::Bool
        if n < 2
            return false
        elseif n < 4
            return true
        elseif n % 2 == 0 || n % 3 == 0
            return false
        end
        for i ∈ 6:6:trunc(Int, √n + 1)
            if n % (i - 1) == 0 || n % (i + 1) == 0; return false end
        end
        return true
    end

    """
        smallest_prime_factor(n::Int, primes:: Vector{Int} = prime_numbers(trunc(Int, √n)))

    Get smallest prime factor of number

    # Arguments
    - `n::Int`: Number
    - `primes:: Vector{Int} = prime_numbers(trunc(Int, √n))`: Vector of primes below at least √n

    # Examples
    ```
    julia> smallest_prime_factor(23)
    23
    julia> n = 22
    julia> primes = prime_numbers(trunc(Int, √n))
    julia> smallest_prime_factor(n, primes)
    2
    ```
    """
    function smallest_prime_factor(n::Int, primes:: Vector{Int} = prime_numbers(trunc(Int, √n)))::Int
        for prime ∈ primes
            if n % prime == 0; return prime end
        end
        return n
    end
end

module Geometry
    """
        triangle_area(a::Vector{Int}, b::Vector{Int}, c::Vector{Int})

    Get area of triangle using coordinates of points

    # Examples:
    julia> triangle_area([0, 0], [0, 3], [4, 0])
    6.0
    julia> triangle_area([-23, 11], [19, 95], [-12, 114])
    1701.0
    """
    triangle_area(a::Vector{Int}, b::Vector{Int}, c::Vector{Int})::Float64 = abs((a[1] - c[1]) * (b[2] - a[2]) - (a[1] - b[1]) * (c[2] - a[2])) / 2
end

module Game
    """
        dice_combinations(p::Int, n::Int, s::Int)

    Get number of combinations of totals `p` for `n` `s`-sided dices

    # Arguments:
    - `p::Int`: Total sum of dices
    - `n::Int`: Number of dices
    - `s::Int`: Number of sides for each dice

    # Examples:
    julia> dice_combinations(6, 6, 6)
    1
    julia> dice_combinations(16, 8, 4)
    3823
    """
    dice_combinations(p::Int, n::Int, s::Int)::Int = sum([(-1) ^ k * binomial(n, k) * binomial(p - 1 - s * k, n - 1) for k ∈ 0:(p - n) ÷ s])
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
@testset "Multivariate Basic Structural Model" begin
    # Multivariate
    front_rear_df = CSV.File(StateSpaceModels.FRONT_REAR_SEAT_KSI) |> DataFrame
    front_rear = [front_rear_df.front_seat_KSI front_rear_df.rear_seats_KSI]

    model = MultivariateBasicStructural(front_rear, 12)
    fit!(model)
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)
    # simualting
    scenarios = simulate_scenarios(model, 10, 10_000)
    test_scenarios_adequacy_with_forecast(forec, scenarios)
end


text:
import Mosek

function run_rmp_mosek(frmp)
    task = Mosek.maketask()
    Mosek.putstreamfunc(task, Mosek.MSK_STREAM_LOG, print)

    Mosek.readdata(task, frmp)

    # Set parameters
    Mosek.putintparam(task, Mosek.MSK_IPAR_NUM_THREADS, 1)
    Mosek.putintparam(task, Mosek.MSK_IPAR_PRESOLVE_USE, 0)
    Mosek.putintparam(task, Mosek.MSK_IPAR_INTPNT_BASIS, 0)

    # Solve
    t = @elapsed Mosek.optimize(task)

    println("Solve time: $t")
    println("Barrier iterations: ", Mosek.getintinf(task, Mosek.MSK_IINF_INTPNT_ITER))

    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_rmp_mosek(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
"""
    photosynthesis(object, meteo, constants = Constants())
    photosynthesis!(object, meteo, constants = Constants())

Generic photosynthesis model for photosynthetic organs. Computes the assimilation and
stomatal conductance according to the models set in (or each component in) `object`.

The models used are defined by the types of the `photosynthesis` and `stomatal_conductance`
fields of `leaf`. For exemple to use the implementation of the Farquhar–von Caemmerer–Berry
(FvCB) model (see [`photosynthesis`](@ref)), the `leaf.photosynthesis` field should be of type
[`Fvcb`](@ref).

# Arguments

- `object`: a [`Component`](@ref) struct ([`AbstractComponentModel`](@ref)), or a Dict/Array of.
- `meteo::Union{AbstractAtmosphere,Weather}`: meteorology structure, see [`Atmosphere`](@ref) or
[`Weather`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

# Using Fvcb model:
leaf = LeafModels(photosynthesis = Fvcb(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Tₗ = 25.0, PPFD = 1000.0, Cₛ = 400.0, Dₗ = meteo.VPD)

photosynthesis(leaf, meteo)

# ---Using several components---

leaf2 = copy(leaf)
leaf2.status.PPFD = 800.0

photosynthesis([leaf,leaf2],meteo)

# ---Using several meteo time-steps---

w = Weather([Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65),
             Atmosphere(T = 25.0, Wind = 1.5, P = 101.3, Rh = 0.55)], (site = "Test site,))

photosynthesis(leaf, w)

# ---Using several meteo time-steps and several components---

photosynthesis(Dict(:leaf1 => leaf, :leaf2 => leaf2), w)

# Using a model file:

model = read_model("a-model-file.yml")

# Initialising the mandatory variables:
init_status!(model, Tₗ = 25.0, PPFD = 1000.0, Cₛ = 400.0, Dₗ = meteo.VPD)

# Running a simulation for all component types in the same scene:
photosynthesis!(model, meteo)
model["Leaf"].status.A

```
"""
function photosynthesis(leaf::AbstractComponentModel, meteo::AbstractAtmosphere, constants = Constants())
    leaf_tmp = copy(leaf)
    photosynthesis!(leaf_tmp, meteo, constants)
    leaf_tmp.status
end

# photosynthesis over several objects (e.g. all leaves of a plant) in an Array
function photosynthesis!(object::O, meteo::AbstractAtmosphere, constants = Constants()) where O <: AbstractArray{<:AbstractComponentModel}

    for i in values(object)
        photosynthesis!(i, meteo, constants)
    end

return nothing
end

# photosynthesis over several objects (e.g. all leaves of a plant) in a kind of Dict.
function photosynthesis!(object::O, meteo::AbstractAtmosphere, constants = Constants()) where {O <: AbstractDict{N,<:AbstractComponentModel} where N}
    for (k, v) in object
        photosynthesis!(v, meteo, constants)
    end

return nothing
end

# same as the above but non-mutating
function photosynthesis(
    object::O,
    meteo::AbstractAtmosphere,
    constants = Constants()
    ) where O <: Union{AbstractArray{<:AbstractComponentModel},AbstractDict{N,<:AbstractComponentModel} where N}

    # Copy the objects only once before the computation for performance reasons:
    object_tmp = copy(object)

    # Computation:
    photosynthesis!(object_tmp, meteo, constants)

    return DataFrame(object_tmp)
end

# photosynthesis over several meteo time steps (called Weather) and possibly several components:
function photosynthesis!(
    object::T,
    meteo::Weather,
    constants = Constants()
    ) where T <: Union{AbstractArray{<:AbstractComponentModel},AbstractDict{N,<:AbstractComponentModel} where N}

    # Check if the meteo data and the status have the same length (or length 1)
    check_status_wheather(object, meteo)

    # Computing for each time-steps:
    for (i, meteo_i) in enumerate(meteo.data)
        # Each object in a time-step:
        for obj in object
            photosynthesis!(obj[i], meteo_i, constants)
        end
    end
end

# If we call weather with one component only, put it in an Array and call the function above
function photosynthesis!(object::AbstractComponentModel, meteo::Weather, constants = Constants())
    photosynthesis!([object], meteo, constants)
end

# photosynthesis over several meteo time steps (same as above) but non-mutating
function photosynthesis(
    object::T,
    meteo::Weather,
    constants = Constants()
    ) where T <: Union{AbstractComponentModel,AbstractDict{N,<:AbstractComponentModel} where N}

    object_tmp = copy(object)
    photosynthesis!(object_tmp, meteo, constants)

    return object_tmp
end

# The function that finally calls assimilation!:
function photosynthesis!(leaf::AbstractComponentModel, meteo::AbstractAtmosphere, constants = Constants())
    is_init = is_initialised(leaf, leaf.photosynthesis, leaf.stomatal_conductance)
    !is_init && error("Some variables must be initialized before simulation (see info message for more details)")
    return assimilation!(leaf, meteo, constants)
end

# The same function (calls assimilation!) but for LeafModels without meteo, so we have to check
# first if the leaf as several time-steps or not. NB: we could do it each time in the function
# just above but I prefer make a new method for performance reasons.
function photosynthesis!(leaf::AbstractComponentModel, meteo::Nothing = nothing, constants = Constants()) where {I,E,A,Gs}
    is_init = is_initialised(leaf, leaf.photosynthesis, leaf.stomatal_conductance)
    !is_init && error("Some variables must be initialized before simulation (see info message for more details)")

    if typeof(leaf.status) == MutableNamedTuples.MutableNamedTuple
        assimilation!(leaf, meteo, constants)
    else
        # We have several time-steps here
        for i in 1:length(leaf.status)
            assimilation!(leaf[i], meteo, constants)
        end

    end

end

"""
    check_status_meteo(component,weather)

Checks if a component status and the wheather have the same length, or
if they can be recycled (length 1)
"""
function check_status_wheather(
        component::LeafModels{I,E,A,Gs,<:Vector{MutableNamedTuples.MutableNamedTuple}},
        weather::Weather
    ) where {I,E,A,Gs}
   length(component.status) != length(weather.data) &&
    @error "component status should have the same number of time-steps than weather"
end

# here we only have one time-step in the status, so we use it for all time-steps
function check_status_wheather(
        component::LeafModels{I,E,A,Gs,<:MutableNamedTuples.MutableNamedTuple},
        weather::Weather
    )  where {I,E,A,Gs}
    nothing
end

# for several components as an array
function check_status_wheather(component::T, weather::Weather) where T <: AbstractArray{<:AbstractComponentModel}
    for i in component
        check_status_wheather(i, weather)
    end
end

# for several components as a Dict
function check_status_wheather(component::T, weather::Weather) where T <: AbstractDict{N,<:AbstractComponentModel} where N
    for (key, val) in component
        check_status_wheather(val, weather)
    end
end


text:
using Flux

function create_model(feature_height, feature_width)

	return Chain(
		input -> reshape(input, size(input)[1:2]..., 1, size(input)[3]),
		# Conv feature_height×4×256 (ReLU, Output: 1×589×256)
		Conv((feature_height, 4), 1 => 256, relu, pad=(0, 5)),
		# Max Pooling 2 (Output: 294×256)
		MaxPool((1, 2)),
		# Conv 4×256 (ReLU, Output: 291×256)
		Conv((1, 4), 256 => 256, relu),
		# Batch Normalization
		BatchNorm(256),
		# Max Pooling 2 (Output: 145×256)
		MaxPool((1, 2)),
		# Conv 4×384 (ReLU, Output: 142×384)
		Conv((1, 4), 256 => 384, relu),
		# Max Pooling 2 (Output 71×384)
		MaxPool((1, 2)),
		# Batch Normalization
		BatchNorm(384),
		# Conv 71×50 (Sigmoid, Output: 50×1)
		Conv((1, 71), 384 => 50, sigmoid),
		input -> reshape(input, getindex(size(input), [3, 4])...),
	)

end

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute()::Int
    for prime ∈ reverse(prime_numbers(7_654_321))
        str_prime = string(prime)
        if sort(parse.(Int, ∪(str_prime))) == collect(1:length(str_prime))
            return prime
        end
    end
end

compute()

@benchmark compute()


text:
"""
    ht_p2_p1_veclap

The manufactured-solution colliding flow example from Elman et al 2014. The
Hood-Taylor formulation with quadratic triangles for the velocity and continuous
pressure on linear triangles.

The formulation is the one derived in Elman, et al., Finite elements and fast
iterative solvers, p. 132. In other words, it is the vector Laplacian version.
"""
module ht_p2_p1_veclap

using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities, ir_identity, attribute, VecAttrib
using MeshSteward: T6block, T6toT3
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary, transform
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs, eldofvals
using Elfel.FEIterators: jacjac, location
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(2 * A, 2 * A, N, N), "velocity")
    ir = baseincrel(vmesh)
    transform(ir, x -> x .- A)
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(Uxh, Uyh, Ph, tndof, mu)
    function integrate!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(pel, pqp)
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * dot(gradNux[i], gradNux[j])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * dot(gradNuy[i], gradNuy[j])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(Uxh), FEIterator(Uyh), FEIterator(Ph))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(Uxh, qargs), QPIterator(Uyh, qargs), QPIterator(Ph, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    @time integrate!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(Ph)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(Ph)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(Ph, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(Uxh, Uyh)
    function integrate!(elits, qpits, trueux, trueuy)
        uxnedof, uynedof = ndofsperel.(elits)
        E = 0.0
        for el in zip(elits...)
            uxel, uyel = el
            uxdofvals = eldofvals(uxel)
            uydofvals = eldofvals(uyel)
            for qp in zip(qpits...)
                uxqp, uyqp = qp
                Jac, J = jacjac(uxel, uxqp)
                JxW = J * weight(uxqp)
                Np = bfun(uxqp)
                uxt = trueux(location(uxel, uxqp)...)
                uyt = trueuy(location(uyel, uyqp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:uxnedof
                    uxa += (uxdofvals[j] * Np[j])
                    uya += (uydofvals[j] * Np[j])
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elits = (FEIterator(Uxh), FEIterator(Uyh),)
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(Uxh, qargs), QPIterator(Uyh, qargs),)
    return integrate!(elits, qpits, trueux, trueuy)
end

function run(N)
    vmesh, pmesh = genmesh(N)
    # Velocity spaces
    Uxh = FESpace(Float64, vmesh, FEH1_T6(), 1)
    Uyh = FESpace(Float64, vmesh, FEH1_T6(), 1)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(Uxh, 0, i, 1, trueux(locs[i]...))
            setebc!(Uyh, 0, i, 1, trueuy(locs[i]...))
        end
    end
    # Pressure space
    Ph = FESpace(Float64, pmesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(pmesh); nearestto = [0.0, 0.0])
    setebc!(Ph, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([Uxh, Uyh, Ph])
    tndof = ndofs(Uxh) + ndofs(Uyh) + ndofs(Ph)
    tnunk = nunknowns(Uxh) + nunknowns(Uyh) + nunknowns(Ph)
    # Assemble the coefficient matrix
    K = assembleK(Uxh, Uyh, Ph, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [Uxh, Uyh, Ph])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([Uxh, Uyh, Ph], U)
    # Postprocessing
    makeattribute(Ph, "p", 1)
    makeattribute(Uxh, "ux", 1)
    makeattribute(Uyh, "uy", 1)
    ep = evaluate_pressure_error(Ph)
    ev = evaluate_velocity_error(Uxh, Uyh)
    vtkwrite("ht_p2_p1_veclap-p", baseincrel(pmesh), [(name = "p",), ])
    vtkwrite("ht_p2_p1_veclap-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    geom = geometry(Ph.mesh)
    ir = baseincrel(Ph.mesh)
    pt = VecAttrib([truep(geom[i]...) for i in 1:length(geom)])
    ir.right.attributes["pt"] = pt
    vtkwrite("ht_p2_p1_veclap-pt", baseincrel(pmesh), [(name = "pt",), ])
    return (ep, ev)
end

end

using .ht_p2_p1_veclap
let
    N = 2
    for loop in 1:5
        ep, ev = ht_p2_p1_veclap.run(N)
        @show N, ep, ev
        N = N * 2
    end
end


text:
function plot(h::Hamiltonian{<:Lattice}; resolution = (1000, 1000), kw...)
    scene = hamiltonianplot(h; resolution = resolution, scale_plot = false, kw...)
    plot = scene[end]
    plot[:tooltips][] && addtooltips!(scene, h)
    scale!(scene)
    return scene
end

@recipe(HamiltonianPlot, hamiltonian) do scene
    Theme(
        allintra = false, allcells = true, showsites = true, showlinks = true,
        shadedsites = false, shadedlinks = false, dimming = 0.75,
        siteradius = 0.2, siteborder = 3, siteborderdarken = 1.0, linkdarken = 0.0,
        linkthickness = 6, linkoffset = missing, linkradius = 0.1,
        tooltips = true, digits = 3,
        _tooltips_rowcolhar = Vector{Tuple{Int,Int,Int}}[],
        colors = map(t -> RGBAf0(t...),
            ((0.960,0.600,.327), (0.410,0.067,0.031),(0.940,0.780,0.000),
            (0.640,0.760,0.900),(0.310,0.370,0.650),(0.600,0.550,0.810),
            (0.150,0.051,0.100),(0.870,0.530,0.640),(0.720,0.130,0.250))),
        light = Vec{3,Float32}[[0, 0, 10], [0, 10, 0], [10, 0, 0], [10, 10, 10], [-10, -10, -10]]
    )
end

function plot!(plot::HamiltonianPlot)
    h = to_value(plot[1])
    lat = h.lattice
    colors = cycle(plot[:colors][])
    sublats = Elsa.sublats(lat)

    mdist = meandist(h)
    mdist > 0 || (mdist = 1)
    plot[:siteradius][] *= mdist/2
    plot[:linkradius][] *= mdist/2

    lo = plot[:linkoffset][]
    plot[:linkoffset][] =
        lo === missing ? √max(1-(plot[:linkradius][]/plot[:siteradius][])^2, 0) : lo * mdist

    # plot sites
    plot[:showsites][] &&
        for (n, har) in enumerate(h.harmonics), (ssrc, csrc) in zip(sublats, colors)
            iszero(har.dn) || plot[:allcells][] || break
            csrc´ = iszero(har.dn) ? csrc : transparent(csrc, 1 - plot[:dimming][])
            itr = siterange(lat, ssrc)
            plotsites!(plot, lat, itr, har.dn, n, csrc´)
        end

    # plot links
    plot[:showlinks][] &&
        for (n, har) in enumerate(h.harmonics)
            iszero(har.dn) || plot[:allcells][] || break
            for (ssrc, csrc) in zip(sublats, colors)
                csrc´ = iszero(har.dn) ? csrc : transparent(csrc, 1 - plot[:dimming][])
                csrc´ = darken(csrc´, plot[:linkdarken][])
                for (sdst, cdst) in zip(sublats, colors)
                    itr = Elsa.eachindex_nz(har, siterange(lat, sdst), siterange(lat, ssrc))
                    plotlinks!(plot, lat, itr, har.dn, n, csrc´)
                end
            end
        end

    return plot
end

function plotsites!(plot, lat, srange, dn, n, color)
    allsites = Elsa.sites(lat)
    br = lat.bravais.matrix
    sites = [padright(allsites[i] + br * dn, Val(3)) for i in srange]
    plot[:tooltips][] && (tt = [(site, 0, n) for site in srange])
    if !isempty(sites)
        plot[:shadedsites][] ? plotsites_hi!(plot, sites, color) : plotsites_lo!(plot, sites, color)
        plot[:tooltips][] && push!(plot[:_tooltips_rowcolhar][], tt)
    end
    return plot
end

function plotsites_lo!(plot, sites, color)
    scatter!(plot, sites;
        color = color,
        markersize = 2 * plot[:siteradius][],
        strokewidth = plot[:siteborder][],
        strokecolor = darken(color, plot[:siteborderdarken][]))
    return nothing
end

function plotsites_hi!(plot, sites, color)
    meshscatter!(plot, sites;
        color = color,
        markersize = plot[:siteradius][], light = plot[:light][])
    return nothing
end

function plotlinks!(plot, lat, itr, dn, n, color)
    links = Pair{SVector{3,Float32},SVector{3,Float32}}[]
    plot[:tooltips][] && (tt = Tuple{Int,Int,Int}[])
    sites = Elsa.sites(lat)
    br = lat.bravais.matrix
    for (row, col) in itr
        iszero(dn) && row == col && continue
        rdst = padright(sites[row] + br * dn, Val(3))
        rsrc = padright(sites[col], Val(3))
        rdst = iszero(dn) ? (rdst + rsrc) / 2 : rdst
        rsrc = rsrc + plot[:siteradius][] * plot[:linkoffset][] * normalize(rdst - rsrc)
        push!(links, rsrc => rdst)
        plot[:tooltips][] && push!(tt, (row, col, n))
    end
    if !isempty(links)
        plot[:shadedlinks][] ? plotlinks_hi!(plot, links, color) : plotlinks_lo!(plot, links, color)
        plot[:tooltips][] && push!(plot[:_tooltips_rowcolhar][], tt)
    end
    return plot
end

function plotlinks_lo!(plot, links, color)
    # linesegments!(plot, links; color = darken(color, plot[:siteborderdarken][]), linewidth = 2+plot[:linkthickness][])
    linesegments!(plot, links; color = color, linewidth = plot[:linkthickness][])
    return nothing
end

function plotlinks_hi!(plot, links, color)
    positions = [(r1 + r2) / 2 for (r1, r2) in links]
    rotvectors = [r2 - r1 for (r1, r2) in links]
    radius = plot[:linkradius][]
    scales = [Vec3f0(radius, radius, norm(r2 - r1)/2) for (r1, r2) in links]
    cylinder = GeometryTypes.Cylinder(Point3f0(0., 0., -1.0), Point3f0(0., 0, 1.0), Float32(1))
    meshscatter!(plot, positions;
        color = color, marker = cylinder, markersize = scales, rotations = rotvectors,
        light = plot[:light][])
    return nothing
end

function addtooltips!(scene, h)
    sceneplot = scene[end]
    visible = Node(false)
    N = Elsa.blockdim(h)
    poprect = lift(scene.events.mouseposition) do mp
        FRect((mp .+ 5), 1,1)
    end
    textpos = lift(scene.events.mouseposition) do mp
        Vec3f0((mp .+ 5 .+ (0, -50))..., 0)
    end
    popup = poly!(campixel(scene), poprect, raw = true, color = RGBAf0(1,1,1,0), visible = visible)
    rect = popup[end]
    translate!(popup, Vec3f0(0, 0, 10000))
    text!(popup, " ", textsize = 30, position = textpos, align = (:center, :center),
        color = :black, strokewidth = 4, strokecolor = :white, raw = true, visible = visible)
    text_field = popup[end]
    on(scene.events.mouseposition) do event
        subplot, idx = mouse_selection(scene)
        layer = findfirst(isequal(subplot), sceneplot.plots)
        if layer !== nothing && idx > 0
            idx´ = fix_linesegments_bug(idx, subplot)
            txt = popuptext(sceneplot, layer, idx´, h)
            text_field[1] = txt
            visible[] = true
        else
            visible[] = false
        end
        return
    end
    return scene
end

# idx returned by mouse_selection seems wrong by a factor 2 in LineSegments subplot
fix_linesegments_bug(idx, subplot::LineSegments) = Int(idx/2) # catches odd idx
fix_linesegments_bug(idx, subplot) = idx

function popuptext(sceneplot, layer, idx, h)
    cache = sceneplot[:_tooltips_rowcolhar][]
    checkbounds(Bool, cache, layer) || return string("Bug in layer: ", layer, " of ", length(cache))
    checkbounds(Bool, cache[layer], idx) || return string("Bug in idx: ", idx, " of ", length(cache[layer]))
    (row, col_or_zero, haridx) = cache[layer][idx]
    if col_or_zero == 0
        col = iszero(col_or_zero) ? row : col_or_zero
        har = h.harmonics[1]
    else
        col = col_or_zero
        har = h.harmonics[haridx]
    end
    element = round.(har.h[row, col], digits = sceneplot[:digits][])
    isreal = all(o -> imag(o) ≈ 0, element)
    txt = isreal ? matrixstring(real.(element)) : matrixstring(element)
    if col_or_zero == 0
        txt´ = string("Onsite[$col] : ", txt)
    else
        txt´ = string("Hopping[$row, $col] : ", txt)
    end
    return txt´
end

matrixstring(x::Number) = string(x)
function matrixstring(s::SMatrix)
    ss = repr("text/plain", s)
    pos = findfirst(isequal('\n'), ss)
    return pos === nothing ? ss : ss[pos:end]
end


text:
using GeometricAlgebra
const ga = GeometricAlgebra
using Test

@basis "---+"
include("definitions.jl")
include("../operators.jl")


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Rotations from an Earth-Centered, Earth-Fixed (ECEF) reference frame to
#   another ECEF reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export rECEFtoECEF

"""
    rECEFtoECEF([T,] ECEFo, ECEFf, JD_UTC::Number, eop_data)

Compute the rotation from an Earth-Centered, Earth-Fixed (`ECEF`) reference
frame to another ECEF reference frame at the Julian Day [UTC] `JD_UTC`. The
rotation description that will be used is given by `T`, which can be `DCM` or
`Quaternion`. The origin ECEF frame is selected by the input `ECEFo` and the
destination ECEF frame is selected by the input `ECEFf`. The model used to
compute the rotation is specified by the selection of the origin and destination
frames. Currently, there are two models supported: IAU-76/FK5 and IAU-2006 with
2010 conventions (CIO approach only).

# Rotation description

The rotations that aligns the origin ECEF frame with the destination ECEF frame
can be described by Direction Cosine Matrices or Quaternions. This is selected
by the parameter `T`.

The possible values are:

* `DCM`: The rotation will be described by a Direction Cosine Matrix.
* `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported yet.**

# ECEF Frame

The supported ECEF frames for both origin `ECEFo` and destination `ECEFf` are:

* `ITRF()`: ECEF will be selected as the International Terrestrial Reference
            Frame (ITRF).
* `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
           frame.
* `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
            System (TIRS).

# EOP Data

The conversion between the supported ECEF frames **always** depends on EOP Data
(see `get_iers_eop` and `read_iers_eop`). If IAU-76/FK5 model is used, then the
type of `eop_data` must be `EOPData_IAU1980`. Otherwise, if IAU-2006/2010 model
is used, then the type of `eop_data` must be `EOPData_IAU2000A`.

# Returns

The rotation description represented by `T` that rotates the ECEF reference
frame into alignment with the ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(:IAU1980);

julia> rECEFtoECEF(PEF(), ITRF(), DatetoJD(1986,6,19,21,35,0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  1.0          0.0         4.35684e-7
  0.0          1.0         1.44762e-6
 -4.35684e-7  -1.44762e-6  1.0

julia> rECEFtoECEF(Quaternion, PEF(), ITRF(), DatetoJD(1986,6,19,21,35,0), eop_IAU1980)
Quaternion{Float64}:
  + 0.9999999999997147 - 7.236343481310813e-7.i + 2.1765518308012794e-7.j + 0.0.k

julia> eop_IAU2000A = get_iers_eop(:IAU2000A);

julia> rECEFtoECEF(TIRS(), ITRF(), DatetoJD(1986,6,19,21,35,0), eop_IAU2000A)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  1.0          3.08408e-11  -4.3531e-7
 -3.14708e-11  1.0          -1.44727e-6
  4.3531e-7    1.44727e-6    1.0

julia> rECEFtoECEF(Quaternion, TIRS(), ITRF(), DatetoJD(1986,6,19,21,35,0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.9999999999997146 - 7.236343481345639e-7.i + 2.176551830689726e-7.j + 1.5577911634233308e-11.k
```
"""
@inline rECEFtoECEF(T_ECEFo::T_ECEFs,
                    T_ECEFf::T_ECEFs,
                    JD_UTC::Number,
                    eop_data::EOPData_IAU1980) =
    rECEFtoECEF(DCM, T_ECEFo, T_ECEFf, JD_UTC, eop_data)

@inline rECEFtoECEF(T_ECEFo::T_ECEFs_IAU_2006, T_ECEFf::T_ECEFs_IAU_2006,
                    JD_UTC::Number, eop_data::EOPData_IAU2000A) =
    rECEFtoECEF(DCM, T_ECEFo, T_ECEFf, JD_UTC, eop_data)

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                 ITRF <=> PEF
# ==============================================================================

function rECEFtoECEF(T::T_ROT, ::Val{:ITRF}, ::Val{:PEF}, JD_UTC::Number,
                     eop_data::EOPData_IAU1980)
    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p = eop_data.x(JD_UTC)*π/648000
    y_p = eop_data.y(JD_UTC)*π/648000

    # Return the rotation.
    return rITRFtoPEF_fk5(T, x_p, y_p)
end

@inline rECEFtoECEF(T::T_ROT, T_ECEFo::Val{:PEF}, T_ECEFf::Val{:ITRF},
                    JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECEFtoECEF(T, T_ECEFf, T_ECEFo, JD_UTC, eop_data))

################################################################################
#                                IAU-2006/2010
################################################################################

#                                ITRF <=> TIRS
# ==============================================================================

function rECEFtoECEF(T::T_ROT, ::Val{:ITRF}, ::Val{:TIRS}, JD_UTC::Number,
                     eop_data::EOPData_IAU2000A)

    arcsec2rad = π/648000

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC)*arcsec2rad
    y_p = eop_data.y(JD_UTC)*arcsec2rad

    # Return the rotation.
    rITRFtoTIRS_iau2006(T, JD_UTC, x_p, y_p)
end

@inline rECEFtoECEF(T::T_ROT, T_ECEFo::Val{:TIRS}, T_ECEFf::Val{:ITRF},
                    JD_UTC::Number, eop_data::EOPData_IAU2000A) =
    inv_rotation(rECEFtoECEF(T, T_ECEFf, T_ECEFo, JD_UTC, eop_data))


text:
using Test
using CompScienceMeshes
using BEAST

# G0 is the interior of the metal (no simulation required here)
# G1 is the unbounded exterior
# G2 is the bounded interior

fn = joinpath(@__DIR__, "assets", "thick_cladding.msh")
G01 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma01")
G02 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma02")
G12 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma12")

G = CompScienceMeshes.read_gmsh_mesh(fn)

@test numcells(G01) + numcells(G02) + numcells(G12) == numcells(G)

# Build region boundaries with normals pointing inward:
G1 = weld(-G01, -G12)
@test CompScienceMeshes.isoriented(G1)

G2 = weld(-G02, G12)
@test CompScienceMeshes.isoriented(G2)

isclosed(m) = (numcells(boundary(m)) == 0)

@test isclosed(G1)
@test isclosed(G2)

# Test the validity of the embedding matrices
E = skeleton(G,1)
E1 = skeleton(G1,1)
A1 = CompScienceMeshes.embedding(E1,E)

nd = BEAST.nedelec(G,E)
@test numfunctions(nd) == numcells(E)
nd1 = BEAST.nedelec(G1,E1)
@test numfunctions(nd1) == numcells(E1)
Z_direct = assemble(BEAST.Identity(), nd1, nd)
Z_expand = assemble(BEAST.Identity(), nd1, nd1) * A1
@test norm(Z_direct - Z_expand, Inf) ≤ √(eps(1.0))