text:
module TestAlign

using Base.Test

using Bio
using Bio.Seq
using Bio.Align
using BGZFStreams
using TestFunctions
using YAML

# Generate a random valid alignment of a sequence of length n against a sequence
# of length m. If `glob` is true, generate a global alignment, if false, a local
# alignment.
function random_alignment(m, n, glob=true)
    match_ops = [OP_MATCH, OP_SEQ_MATCH, OP_SEQ_MISMATCH]
    insert_ops = [OP_INSERT, OP_SOFT_CLIP, OP_HARD_CLIP]
    delete_ops = [OP_DELETE, OP_SKIP]
    ops = vcat(match_ops, insert_ops, delete_ops)

    # This is just a random walk on a m-by-n matrix, where steps are either
    # (+1,0), (0,+1), (+1,+1). To make somewhat more realistic alignments, it's
    # biased towards going in the same direction. Local alignments have a random
    # start and end time, global alignments always start at (0,0) and end at
    # (m,n).

    # probability of choosing the same direction as the last step
    straight_pr = 0.9

    op = OP_MATCH
    if glob
        i = 0
        j = 0
        i_end = m
        j_end = n
    else
        i = rand(1:m-1)
        j = rand(1:n-1)
        i_end = rand(i+1:m)
        j_end = rand(j+1:n)
    end

    path = AlignmentAnchor[AlignmentAnchor(i, j, OP_START)]
    while (glob && i < i_end && j < j_end) || (!glob && (i < i_end || j < j_end))
        straight = rand() < straight_pr

        if i == i_end
            if !straight
                op = rand(delete_ops)
            end
            j += 1
        elseif j == j_end
            if !straight
                op = rand(inset_ops)
            end
            i += 1
        else
            if !straight
                op = rand(ops)
            end

            if isdeleteop(op)
                j += 1
            elseif isinsertop(op)
                i += 1
            else
                i += 1
                j += 1
            end
        end
        push!(path, AlignmentAnchor(i, j, op))
    end

    return path
end

# Make an Alignment from a path returned by random_alignment. Converting from
# path to Alignment is just done by removing redundant nodes from the path.
function anchors_from_path(path)
    anchors = AlignmentAnchor[]
    for k in 1:length(path)
        if k == length(path) || path[k].op != path[k+1].op
            push!(anchors, path[k])
        end
    end
    return anchors
end

# Generate a random range within `range`.
function randrange(range)
    x = rand(range)
    y = rand(range)
    if x < y
        return x:y
    else
        return y:x
    end
end

@testset "Alignments" begin
    @testset "Operations" begin
        @testset "Constructors and Conversions" begin
            ops = Set(Align.char_to_op)
            for op in ops
                if op != Align.OP_INVALID
                    @test Operation(Char(op)) == op
                end
            end
            @test_throws Exception Char(Align.OP_INVALID)
            @test_throws Exception Operation('m')
            @test_throws Exception Operation('7')
            @test_throws Exception Operation('A')
            @test_throws Exception Operation('\n')
        end
    end

    @testset "AlignmentAnchor" begin
        anchor = AlignmentAnchor(1, 2, OP_MATCH)
        @test string(anchor) == "AlignmentAnchor(1, 2, 'M')"
    end

    @testset "Alignment" begin
        # alignments with nonsense operations
        @test_throws Exception Alignment(AlignmentAnchor[
            Operation(0, 0, OP_START),
            Operation(100, 100, convert(Operation, 0xfa))])

        # test bad alignment anchors by swapping nodes in paths
        for _ in 1:100
            path = random_alignment(rand(1000:10000), rand(1000:10000))
            anchors = anchors_from_path(path)
            n = length(anchors)
            n < 3 && continue
            i = rand(2:n-1)
            j = rand(i+1:n)
            anchors[i], anchors[j] = anchors[j], anchors[i]
            @test_throws Exception Alignment(anchors)
        end

        # test bad alignment anchors by swapping operations
        for _ in 1:100
            path = random_alignment(rand(1000:10000), rand(1000:10000))
            anchors = anchors_from_path(path)
            n = length(anchors)
            n < 3 && continue
            i = rand(2:n-1)
            j = rand(i+1:n)
            u = anchors[i]
            v = anchors[j]
            if (ismatchop(u.op) && ismatchop(v.op)) ||
               (isinsertop(u.op) && isinsertop(v.op)) ||
               (isdeleteop(u.op) && isdeleteop(v.op))
                continue
            end
            anchors[i] = AlignmentAnchor(u.seqpos, u.refpos, v.op)
            anchors[j] = AlignmentAnchor(v.seqpos, v.refpos, u.op)
            @test_throws Exception Alignment(anchors)
        end

        # cigar string round-trip
        for _ in 1:100
            path = random_alignment(rand(1000:10000), rand(1000:10000))
            anchors = anchors_from_path(path)
            aln = Alignment(anchors)
            cig = cigar(aln)
            @test Alignment(cig, aln.anchors[1].seqpos + 1,
                            aln.anchors[1].refpos + 1) == aln
        end
    end

    @testset "AlignedSequence" begin
        #               0   4        9  12 15     19
        #               |   |        |  |  |      |
        #     query:     TGGC----ATCATTTAACG---CAAG
        # reference: AGGGTGGCATTTATCAG---ACGTTTCGAGAC
        #               |   |   |    |     |  |   |
        #               4   8   12   17    20 23  27
        anchors = [
            AlignmentAnchor( 0,  4, OP_START),
            AlignmentAnchor( 4,  8, OP_MATCH),
            AlignmentAnchor( 4, 12, OP_DELETE),
            AlignmentAnchor( 9, 17, OP_MATCH),
            AlignmentAnchor(12, 17, OP_INSERT),
            AlignmentAnchor(15, 20, OP_MATCH),
            AlignmentAnchor(15, 23, OP_DELETE),
            AlignmentAnchor(19, 27, OP_MATCH)
        ]
        query = "TGGCATCATTTAACGCAAG"
        alnseq = AlignedSequence(query, anchors)
        @test Bio.Align.first(alnseq) ==  5
        @test Bio.Align.last(alnseq)  == 27
        # OP_MATCH
        for (seqpos, refpos) in [(1, 5), (2, 6), (4, 8), (13, 18), (19, 27)]
            @test seq2ref(alnseq, seqpos) == (refpos, OP_MATCH)
            @test ref2seq(alnseq, refpos) == (seqpos, OP_MATCH)
        end
        # OP_INSERT
        @test seq2ref(alnseq, 10) == (17, OP_INSERT)
        @test seq2ref(alnseq, 11) == (17, OP_INSERT)
        # OP_DELETE
        @test ref2seq(alnseq,  9) == ( 4, OP_DELETE)
        @test ref2seq(alnseq, 10) == ( 4, OP_DELETE)
        @test ref2seq(alnseq, 23) == (15, OP_DELETE)

        seq = dna"ACGG--TGAAAGGT"
        ref = dna"-CGGGGA----TTT"
        alnseq = AlignedSequence(seq, ref)
        @test Bio.Align.first(alnseq) == 1
        @test Bio.Align.last(alnseq)  == 9
        @test alnseq.aln.anchors == [
             AlignmentAnchor( 0, 0, '0')
             AlignmentAnchor( 1, 0, 'I')
             AlignmentAnchor( 4, 3, '=')
             AlignmentAnchor( 4, 5, 'D')
             AlignmentAnchor( 5, 6, 'X')
             AlignmentAnchor( 9, 6, 'I')
             AlignmentAnchor(11, 8, 'X')
             AlignmentAnchor(12, 9, '=')
        ]
    end
end

# generate test cases from two aligned sequences
function alnscore{S,T}(::Type{S}, affinegap::AffineGapScoreModel{T}, alnstr::AbstractString, clip::Bool)
    gap_open = affinegap.gap_open
    gap_extend = affinegap.gap_extend
    lines = split(chomp(alnstr), '\n')
    a, b = lines[1:2]
    m = length(a)
    @assert m == length(b)

    if length(lines) == 2
        start = 1
        while start ≤ m && a[start] == ' ' || b[start] == ' '
            start += 1
        end
        stop = start
        while stop + 1 ≤ m && !(a[stop+1] == ' ' || b[stop+1] == ' ')
            stop += 1
        end
    elseif length(lines) == 3
        start = search(lines[3], '^')
        stop = rsearch(lines[3], '^')
    else
        error("invalid alignment string")
    end

    score = T(0)
    gap_extending_a = false
    gap_extending_b = false
    for i in start:stop
        if a[i] == '-'
            score += gap_extending_a ? gap_extend : (gap_open + gap_extend)
            gap_extending_a = true
        elseif b[i] == '-'
            score += gap_extending_b ? gap_extend : (gap_open + gap_extend)
            gap_extending_b = true
        else
            score += affinegap.submat[a[i],b[i]]
            gap_extending_a = false
            gap_extending_b = false
        end
    end
    sa = S(replace(a, r"\s|-", ""))
    sb = S(replace(b, r"\s|-", ""))
    return sa, sb, score, clip ? string(a[start:stop], '\n', b[start:stop]) : string(a, '\n', b)
end

function alnscore{T}(affinegap::AffineGapScoreModel{T}, alnstr::AbstractString; clip=true)
    return alnscore(AbstractString, affinegap, alnstr, clip)
end

function alndistance{S,T}(::Type{S}, cost::CostModel{T}, alnstr::AbstractString)
    lines = split(chomp(alnstr), '\n')
    @assert length(lines) == 2
    a, b = lines
    m = length(a)
    @assert length(b) == m
    dist = T(0)
    for i in 1:m
        if a[i] == '-'
            dist += cost.deletion
        elseif b[i] == '-'
            dist += cost.insertion
        else
            dist += cost.submat[a[i],b[i]]
        end
    end
    return S(replace(a, r"\s|-", "")), S(replace(b, r"\s|-", "")), dist
end

function alndistance{T}(cost::CostModel{T}, alnstr::AbstractString)
    return alndistance(AbstractString, cost, alnstr)
end

function alignedpair(alnres)
    aln = alignment(alnres)
    a = aln.a
    b = aln.b
    anchors = a.aln.anchors
    buf = IOBuffer()
    print_seq(buf, a, anchors)
    println(buf)
    print_ref(buf, b, anchors)
    return String(buf)
end

function print_seq(io, seq, anchors)
    for i in 2:length(anchors)
        if ismatchop(anchors[i].op) || isinsertop(anchors[i].op)
            for j in anchors[i-1].seqpos+1:anchors[i].seqpos
                print(io, seq.seq[j])
            end
        elseif isdeleteop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                print(io, '-')
            end
        end
    end
end

function print_ref(io, ref, anchors)
    for i in 2:length(anchors)
        if ismatchop(anchors[i].op) || isdeleteop(anchors[i].op)
            for j in anchors[i-1].refpos+1:anchors[i].refpos
                print(io, ref[j])
            end
        elseif isinsertop(anchors[i].op)
            for _ in anchors[i-1].seqpos+1:anchors[i].seqpos
                print(io, '-')
            end
        end
    end
end

@testset "PairwiseAlignment" begin
    @testset "SubstitutionMatrix" begin
        # DNA
        @test EDNAFULL[DNA_A,DNA_A] ===  5
        @test EDNAFULL[DNA_G,DNA_G] ===  5
        @test EDNAFULL[DNA_A,DNA_G] === -4
        @test EDNAFULL[DNA_G,DNA_A] === -4
        @test EDNAFULL[DNA_M,DNA_T] === -4
        @test EDNAFULL[DNA_M,DNA_C] ===  1

        # amino acid
        @test BLOSUM62[AA_A,AA_R] === -1
        @test BLOSUM62[AA_R,AA_A] === -1
        @test BLOSUM62[AA_R,AA_R] ===  5
        @test BLOSUM62[AA_O,AA_R] ===  0  # default
        @test BLOSUM62[AA_R,AA_O] ===  0  # default

        # update
        myblosum = copy(BLOSUM62)
        @test myblosum[AA_A,AA_R] === -1
        myblosum[AA_A,AA_R] = 10
        @test myblosum[AA_A,AA_R] === 10

        @test BLOSUM62[AA_O,AA_R] ===  0  # default
        myblosum[AA_O,AA_R] = -3
        @test myblosum[AA_O,AA_R] === -3

        submat = SubstitutionMatrix(DNANucleotide, rand(Float64, 15, 15))
        @test isa(submat, SubstitutionMatrix{DNANucleotide,Float64})

        submat = DichotomousSubstitutionMatrix(5, -4)
        @test isa(submat, DichotomousSubstitutionMatrix{Int})
        submat = convert(SubstitutionMatrix{DNANucleotide,Int}, submat)
        @test submat[DNA_A,DNA_A] ===  5
        @test submat[DNA_C,DNA_C] ===  5
        @test submat[DNA_A,DNA_C] === -4
        @test submat[DNA_C,DNA_A] === -4

        try
            print(IOBuffer(), EDNAFULL)
            print(IOBuffer(), BLOSUM62)
            # no error
            @test true
        catch
            @test false
        end
    end

    @testset "AffineGapScoreModel" begin
        # predefined substitution matrix
        for affinegap in [AffineGapScoreModel(BLOSUM62, -10, -1),
                          AffineGapScoreModel(BLOSUM62, gap_open=-10, gap_extend=-1),
                          AffineGapScoreModel(BLOSUM62, gap_open_penalty=10, gap_extend_penalty=1)]
            @test affinegap.gap_open == -10
            @test affinegap.gap_extend == -1
            @test typeof(affinegap) == AffineGapScoreModel{Int}
        end

        # matrix
        submat = SubstitutionMatrix(DNANucleotide, rand(Float64, 15, 15))
        for affinegap in [AffineGapScoreModel(submat, -3, -1),
                          AffineGapScoreModel(submat, gap_open=-3, gap_extend=-1),
                          AffineGapScoreModel(submat, gap_open_penalty=3, gap_extend_penalty=1)]
            @test affinegap.gap_open == -3
            @test affinegap.gap_extend == -1
            @test typeof(affinegap) == AffineGapScoreModel{Float64}
        end

        affinegap = AffineGapScoreModel(match=3, mismatch=-3, gap_open=-5, gap_extend=-2)
        @test affinegap.gap_open == -5
        @test affinegap.gap_extend == -2
        @test typeof(affinegap) == AffineGapScoreModel{Int}
    end

    @testset "CostModel" begin
        submat = SubstitutionMatrix(DNANucleotide, rand(Int, 15, 15))
        for cost in [CostModel(submat, 5, 6),
                     CostModel(submat, insertion=5, deletion=6)]
            @test cost.insertion == 5
            @test cost.deletion == 6
            @test typeof(cost) == CostModel{Int}
        end

        cost = CostModel(match=0, mismatch=3, insertion=5, deletion=6)
        @test cost.insertion == 5
        @test cost.deletion == 6
        @test typeof(cost) == CostModel{Int}
    end

    @testset "Alignment" begin
        anchors = [
            AlignmentAnchor(0, 0, OP_START),
            AlignmentAnchor(3, 3, OP_SEQ_MATCH)
        ]
        seq = AlignedSequence("ACG", anchors)
        ref = "ACG"
        aln = PairwiseAlignment(seq, ref)
        @test collect(aln) == [('A', 'A'), ('C', 'C'), ('G', 'G')]
        result = PairwiseAlignmentResult(3, true, seq, ref)
        @test isa(result, PairwiseAlignmentResult) == true
        @test isa(alignment(result), PairwiseAlignment) == true
        @test score(result) == 3
        @test hasalignment(result) == true
    end

    @testset "count_<ops>" begin
        # anchors are derived from an alignment:
        #   seq: ACG---TGCAGAATTT
        #        |     || || ||
        #   ref: AAAATTTGAAGTAT--
        a = dna"ACGTGCAGAATTT"
        b = dna"AAAATTTGAAGTAT"
        anchors = [
            AlignmentAnchor( 0,  0, '0'),
            AlignmentAnchor( 1,  1, '='),
            AlignmentAnchor( 3,  3, 'X'),
            AlignmentAnchor( 3,  6, 'D'),
            AlignmentAnchor( 5,  8, '='),
            AlignmentAnchor( 6,  9, 'X'),
            AlignmentAnchor( 8, 11, '='),
            AlignmentAnchor( 9, 12, 'X'),
            AlignmentAnchor(11, 14, '='),
            AlignmentAnchor(13, 14, 'I')
        ]
        aln = PairwiseAlignment(AlignedSequence(a, anchors), b)
        @test count_matches(aln) == 7
        @test count_mismatches(aln) == 4
        @test count_insertions(aln) == 2
        @test count_deletions(aln) == 3
        @test count_aligned(aln) == 16
    end

    @testset "Interfaces" begin
        seq = dna"ACGTATAGT"
        ref = dna"ATCGTATTGGT"
        # seq:  1 A-CGTATAG-T  9
        #         | ||||| | |
        # ref:  1 ATCGTATTGGT 11
        model = AffineGapScoreModel(EDNAFULL, gap_open=-4, gap_extend=-1)
        result = pairalign(GlobalAlignment(), seq, ref, model)
        @test isa(result, PairwiseAlignmentResult)
        aln = alignment(result)
        @test isa(aln, PairwiseAlignment)
        @test seq2ref(aln, 1) == (1, OP_SEQ_MATCH)
        @test seq2ref(aln, 2) == (3, OP_SEQ_MATCH)
        @test seq2ref(aln, 3) == (4, OP_SEQ_MATCH)
        @test ref2seq(aln, 1) == (1, OP_SEQ_MATCH)
        @test ref2seq(aln, 2) == (1, OP_DELETE)
        @test ref2seq(aln, 3) == (2, OP_SEQ_MATCH)
    end

    @testset "GlobalAlignment" begin
        affinegap = AffineGapScoreModel(
            match=0,
            mismatch=-6,
            gap_open=-5,
            gap_extend=-3
        )

        function testaln(alnstr)
            a, b, s, alnpair = alnscore(affinegap, alnstr)
            aln = pairalign(GlobalAlignment(), a, b, affinegap)
            @test score(aln) == s
            @test alignedpair(aln) == alnpair
            aln = pairalign(GlobalAlignment(), a, b, affinegap, score_only=true)
            @test score(aln) == s
        end

        @testset "empty sequences" begin
            aln = pairalign(GlobalAlignment(), "", "", affinegap)
            @test score(aln) == 0
        end

        @testset "complete match" begin
            testaln("""
            ACGT
            ACGT
            """)
        end

        @testset "mismatch" begin
            testaln("""
            ACGT
            AGGT
            """)

            testaln("""
            ACGT
            AGGA
            """)
        end

        @testset "insertion" begin
            testaln("""
            ACGTT
            ACGT-
            """)

            testaln("""
            ACGTTT
            ACGT--
            """)

            testaln("""
            ACCGT
            AC-GT
            """)

            testaln("""
            ACCCGT
            AC--GT
            """)

            testaln("""
            AACGT
            A-CGT
            """)

            testaln("""
            AAACGT
            A--CGT
            """)
        end

        @testset "deletion" begin
            testaln("""
            ACGT-
            ACGTT
            """)

            testaln("""
            ACGT-
            ACGTT
            """)

            testaln("""
            ACGT--
            ACGTTT
            """)

            testaln("""
            AC-GT
            ACCGT
            """)

            testaln("""
            AC--GT
            ACCCGT
            """)

            testaln("""
            A-CGT
            AACGT
            """)

            testaln("""
            A--CGT
            AAACGT
            """)
        end

        @testset "banded" begin
            a, b, s, alnpair = alnscore(affinegap, """
            ACGT
            ACGT
            """)
            aln = pairalign(GlobalAlignment(), a, b, affinegap, banded=true)
            @test score(aln) == s
            @test alignedpair(aln) == alnpair

            a, b, s, alnpair = alnscore(affinegap, """
            ACGT
            AGGT
            """)
            aln = pairalign(GlobalAlignment(), a, b, affinegap, banded=true)
            @test score(aln) == s
            @test alignedpair(aln) == alnpair

            a, b, s, alnpair = alnscore(affinegap, """
            ACG--T
            ACGAAT
            """)
            aln = pairalign(GlobalAlignment(), a, b, affinegap, banded=true, lower_offset=0, upper_offset=0)
            @test score(aln) == s
            @test alignedpair(aln) == alnpair
        end
    end

    @testset "SemiGlobalAlignment" begin
        affinegap = AffineGapScoreModel(
            match=0,
            mismatch=-6,
            gap_open=-5,
            gap_extend=-3
        )

        function testaln(alnstr)
            a, b, s, alnpair = alnscore(affinegap, alnstr, clip=false)
            aln = pairalign(SemiGlobalAlignment(), a, b, affinegap)
            @test score(aln) == s
            @test alignedpair(aln) == alnpair
            aln = pairalign(SemiGlobalAlignment(), a, b, affinegap, score_only=true)
            @test score(aln) == s
        end

        @testset "complete match" begin
            testaln("""
            ACGT
            ACGT
            """)
        end

        @testset "partial match" begin
            testaln("""
            --ACTT---
            TTACGTAGT
              ^^^^
            """)

            testaln("""
            --AC-TTG-
            TTACGTTGT
              ^^^^^^
            """)

            testaln("""
            --ACTAGT---
            TTAC--GTTGT
              ^^^^^^
            """)
        end
    end

    @testset "OverlapAlignment" begin
        affinegap = AffineGapScoreModel(
            match=3,
            mismatch=-6,
            gap_open=-5,
            gap_extend=-3
        )

        function testaln(alnstr)
            a, b, s, alnpair = alnscore(affinegap, alnstr, clip=false)
            aln = pairalign(OverlapAlignment(), a, b, affinegap)
            @test score(aln) == s
            @test alignedpair(aln) == alnpair
            aln = pairalign(OverlapAlignment(), a, b, affinegap, score_only=true)
            @test score(aln) == s
        end

        @testset "complete match" begin
            testaln("""
            ACGT
            ACGT
            """)
        end

        @testset "partial match" begin
            testaln("""
            ---ACGGTGATTAT
            GATACGGTGA----
               ^^^^^^^
            """)

            testaln("""
            ---AACGT-GATTAT
            GATAACGGAGA----
               ^^^^^^^^
            """)

            testaln("""
            GATACGGTGA----
            ---ACGGTGATTAT
               ^^^^^^^
            """)

            testaln("""
            GATAACGGAGA----
            ---AACGT-GATTAT
               ^^^^^^^^
            """)
        end
    end

    @testset "LocalAlignment" begin
        @testset "zero matching score" begin
            affinegap = AffineGapScoreModel(
                match=0,
                mismatch=-6,
                gap_open=-5,
                gap_extend=-3
            )

            function testaln(alnstr)
                a, b, s, alnpair = alnscore(affinegap, alnstr)
                aln = pairalign(LocalAlignment(), a, b, affinegap)
                @test score(aln) == s
                @test alignedpair(aln) == alnpair
                aln = pairalign(LocalAlignment(), a, b, affinegap, score_only=true)
                @test score(aln) == s
            end

            @testset "empty sequences" begin
                aln = pairalign(LocalAlignment(), "", "", affinegap)
                @test score(aln) == 0
            end

            @testset "complete match" begin
                testaln("""
                ACGT
                ACGT
                """)
            end

            @testset "partial match" begin
                testaln("""
                ACGT
                AGGT
                  ^^
                """)

                testaln("""
                   ACGT
                AACGTTT
                      ^
                """)
            end

            @testset "no match" begin
                a = "AA"
                b = "TTTT"
                aln = pairalign(LocalAlignment(), a, b, affinegap)
                @test score(aln) == 0
            end
        end

        @testset "positive matching score" begin
            affinegap = AffineGapScoreModel(
                match=5,
                mismatch=-6,
                gap_open=-5,
                gap_extend=-3
            )

            function testaln(alnstr)
                a, b, s, alnpair = alnscore(affinegap, alnstr)
                aln = pairalign(LocalAlignment(), a, b, affinegap)
                @test score(aln) == s
                @test alignedpair(aln) == alnpair
                aln = pairalign(LocalAlignment(), a, b, affinegap, score_only=true)
                @test score(aln) == s
            end

            @testset "complete match" begin
                testaln("""
                ACGT
                ACGT
                ^^^^
                """)
            end

            @testset "partial match" begin
                testaln("""
                ACGT
                AGGT
                  ^^
                """)

                testaln("""
                 ACGT
                AACGTTT
                 ^^^^
                """)

                testaln("""
                  AC-GT
                AAACTGTTT
                """)
            end

            @testset "no match" begin
                a = "AA"
                b = "TTTT"
                aln = pairalign(LocalAlignment(), a, b, affinegap)
                @test score(aln) == 0
            end
        end
    end

    @testset "EditDistance" begin
        mismatch = 1
        submat = DichotomousSubstitutionMatrix(0, mismatch)
        insertion = 1
        deletion = 2
        cost = CostModel(submat, insertion, deletion)

        function testaln(alnstr)
            a, b, dist = alndistance(cost, alnstr)
            aln = pairalign(EditDistance(), a, b, cost)
            @test distance(aln) == dist
            @test alignedpair(aln) == chomp(alnstr)
            aln = pairalign(EditDistance(), a, b, cost, distance_only=true)
            @test distance(aln) == dist
        end

        @testset "empty sequences" begin
            aln = pairalign(EditDistance(), "", "", cost)
            @test distance(aln) == 0
        end

        @testset "complete match" begin
            testaln("""
            ACGT
            ACGT
            """)
        end

        @testset "mismatch" begin
            testaln("""
            AGGT
            ACGT
            """)

            testaln("""
            AGGT
            ACGT
            """)
        end

        @testset "insertion" begin
            testaln("""
            ACGTT
            ACG-T
            """)

            testaln("""
            ACGTT
            -CG-T
            """)
        end

        @testset "deletion" begin
            testaln("""
            AC-T
            ACGT
            """)

            testaln("""
            -C-T
            ACGT
            """)
        end
    end

    @testset "LevenshteinDistance" begin
        @testset "empty sequences" begin
            aln = pairalign(LevenshteinDistance(), "", "")
            @test distance(aln) == 0
        end

        @testset "complete match" begin
            a = "ACGT"
            b = "ACGT"
            aln = pairalign(LevenshteinDistance(), a, b)
            @test distance(aln) == 0
        end
    end

    @testset "HammingDistance" begin
        function testaln(alnstr)
            a, b = split(chomp(alnstr), '\n')
            dist = sum([x != y for (x, y) in zip(a, b)])
            aln = pairalign(HammingDistance(), a, b)
            @test distance(aln) == dist
            @test alignedpair(aln) == chomp(alnstr)
            aln = pairalign(HammingDistance(), a, b, distance_only=true)
            @test distance(aln) == dist
        end

        @testset "empty sequences" begin
            aln = pairalign(HammingDistance(), "", "")
            @test distance(aln) == 0
        end

        @testset "complete match" begin
            testaln("""
            ACGT
            ACGT
            """)
        end

        @testset "mismatch" begin
            testaln("""
            ACGT
            AGGT
            """)

            testaln("""
            ACGT
            AGGA
            """)
        end

        @testset "indel" begin
            @test_throws Exception pairalign(HammingDistance(), "ACGT", "ACG")
            @test_throws Exception pairalign(HammingDistance(), "ACG", "ACGT")
        end
    end
end

@testset "High-throughput Sequencing" begin
    get_bio_fmt_specimens()

    @testset "AuxDataDict" begin
        dict = AuxDataDict()
        @test length(dict) == 0
        @test isempty(dict)
        @test_throws KeyError dict["NM"]

        dict = AuxDataDict(
            "X1" => UInt8(1),
            "X2" => UInt16(2),
            "X3" => UInt32(3),
            "X4" => Int8(4),
            "X5" => Int16(5),
            "X6" => Int32(6),
            "X7" => Float32(7),
            "X8" => "eight",
            "X9" => Int32[9])
        @test length(dict) == 9
        @test !isempty(dict)
        @test dict["X1"] === UInt8(1)
        @test dict["X2"] === UInt16(2)
        @test dict["X3"] === UInt32(3)
        @test dict["X4"] === Int8(4)
        @test dict["X5"] === Int16(5)
        @test dict["X6"] === Int32(6)
        @test dict["X7"] === Float32(7)
        @test dict["X8"] == "eight"
        @test typeof(dict["X8"]) == String
        @test dict["X9"] == Int32[9]
        @test typeof(dict["X9"]) == Vector{Int32}

        dict = AuxDataDict("NM" => 0x01, "XY" => Int32(100), "XZ" => [0x11, 0x23])
        @test length(dict) == 3
        @test dict["NM"] === 0x01
        @test dict["XY"] === Int32(100)
        @test dict["XZ"] == [0x11, 0x23]
        @test eltype(dict["XZ"]) == UInt8

        dict = AuxDataDict("NM" => 0x01, "MD" => "8T1T39")
        @test length(dict) == 2
        @test dict["NM"] === 0x01
        @test dict["MD"] == "8T1T39"
        dict["NM"] = 0x00
        @test dict["NM"] === 0x00
        dict["MD"] = "50"
        @test dict["MD"] == "50"
        @test collect(dict) == ["NM" => 0x00, "MD" => "50"]
        dict["XY"] = "foobar"
        @test dict["XY"] == "foobar"
        @test collect(dict) == ["NM" => 0x00, "MD" => "50", "XY" => "foobar"]
        delete!(dict, "NM")
        @test length(dict) == 2
        @test collect(dict) == ["MD" => "50", "XY" => "foobar"]
    end

    @testset "SAM" begin
        samdir = Pkg.dir("Bio", "test", "BioFmtSpecimens", "SAM")

        @testset "SAMHeader" begin
            h = SAMHeader()
            @test isa(h, Associative)
            @test isempty(h)
            h["HD"] = Dict("VN" => "100.100", "SO" => "unknown")
            @test length(h) == 1
            @test h["HD"]["VN"] == "100.100"
            h["CO"] = ["comment1", "comment2"]
            @test length(h) == 2
            @test h["CO"] == ["comment1", "comment2"]
            delete!(h, "CO")
            @test length(h) == 1
        end

        @testset "Record" begin
            rec = SAMRecord()
            @test !ismapped(rec)

            # default values
            @test seqname(rec) == "*"
            @test flag(rec) == SAM_FLAG_UNMAP
            @test refname(rec) == "*"
            @test leftposition(rec) == 0
            @test rightposition(rec) == -1
            @test nextrefname(rec) == "*"
            @test nextleftposition(rec) == 0
            @test templatelength(rec) == 0
            @test cigar(rec) == "*"
            @test alignment(rec) == Alignment(AlignmentAnchor[])
            @test sequence(rec) == "*"
            @test_throws ArgumentError seqlength(rec)
            @test qualities(rec) == "*"
            @test Align.alignment_length(rec) === 0

            buf = IOBuffer()
            show(buf, rec)
            @test startswith(takebuf_string(buf), "Bio.Align.SAMRecord:")

            buf = IOBuffer()
            showcompact(buf, rec)
            @test takebuf_string(buf) == "*\tunmapped"

            # set & delete tags
            rec = SAMRecord()
            @test !haskey(rec, "MN")
            rec["MN"] = 0x01
            @test rec["MN"] === 0x01
            @test haskey(rec, "MN")
            @test !haskey(rec, "XY")
            rec["XY"] = "foobar"
            @test rec["XY"] == "foobar"
            @test haskey(rec, "XY")
            delete!(rec, "MN")
            @test !haskey(rec, "MN")
        end

        @testset "Reader" begin
            reader = open(SAMReader, joinpath(samdir, "ce#1.sam"))
            @test isa(reader, SAMReader)
            @test eltype(reader) === SAMRecord
            @test startswith(repr(reader), "Bio.Align.SAMReader:")

            # header
            h = header(reader)
            @test h["SQ"] == [Dict("SN" => "CHROMOSOME_I", "LN" => "1009800")]

            # first record
            rec = SAMRecord()
            read!(reader, rec)
            @test ismapped(rec)
            @test refname(rec) == "CHROMOSOME_I"
            @test leftposition(rec) == 2
            @test rightposition(rec) == 102
            @test seqname(rec) == "SRR065390.14978392"
            @test sequence(rec)  == "CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA"
            @test seqlength(rec) == 100
            @test qualities(rec) == "#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"
            @test flag(rec) == 16
            @test cigar(rec) == "27M1D73M"
            @test alignment(rec) == Alignment([
                AlignmentAnchor(  0,   1, OP_START),
                AlignmentAnchor( 27,  28, OP_MATCH),
                AlignmentAnchor( 27,  29, OP_DELETE),
                AlignmentAnchor(100, 102, OP_MATCH)])
            @test rec["XG"] == 1
            @test rec["XM"] == 5
            @test rec["XN"] == 0
            @test rec["XO"] == 1
            @test rec["AS"] == -18
            @test rec["XS"] == -18
            @test rec["YT"] == "UU"
            @test eof(reader)
            close(reader)

            # iterator
            @test length(collect(open(SAMReader, joinpath(samdir, "ce#1.sam")))) == 1
            @test length(collect(open(SAMReader, joinpath(samdir, "ce#2.sam")))) == 2

            # IOStream
            @test length(collect(SAMReader(open(joinpath(samdir, "ce#1.sam"))))) == 1
            @test length(collect(SAMReader(open(joinpath(samdir, "ce#2.sam"))))) == 2
        end

        @testset "Round trip" begin
            for specimen in YAML.load_file(joinpath(samdir, "index.yml"))
                filepath = joinpath(samdir, specimen["filename"])
                mktemp() do path, io
                    # copy
                    reader = open(SAMReader, filepath)
                    writer = SAMWriter(io, header(reader))
                    records = SAMRecord[]
                    for rec in reader
                        push!(records, rec)
                        write(writer, rec)
                    end
                    close(reader)
                    close(writer)

                    # read again
                    reader = open(SAMReader, path)
                    @test collect(reader) == records
                    close(reader)
                end
            end
        end
    end

    @testset "BAM" begin
        bamdir = Pkg.dir("Bio", "test", "BioFmtSpecimens", "BAM")

        @testset "Record" begin
            rec = BAMRecord()
            @test !ismapped(rec)

            # default values
            @test refname(rec) == "*"
            @test refindex(rec) == 0
            @test leftposition(rec) == 0
            @test rightposition(rec) == -1
            @test mappingquality(rec) == 0
            @test flag(rec) == SAM_FLAG_UNMAP
            @test nextrefname(rec) == "*"
            @test nextrefindex(rec) == 0
            @test nextleftposition(rec) == 0
            @test templatelength(rec) == 0
            @test seqname(rec) == ""
            @test cigar(rec) == ""
            @test sequence(rec) == dna""
            @test seqlength(rec) == 0
            @test alignment(rec) == Alignment(AlignmentAnchor[])
            @test qualities(rec) == UInt8[]
            @test Align.alignment_length(rec) === 0

            buf = IOBuffer()
            show(buf, rec)
            @test startswith(takebuf_string(buf), "Bio.Align.BAMRecord:")

            buf = IOBuffer()
            showcompact(buf, rec)
            @test takebuf_string(buf) == "(empty name)\tunmapped"

            # set & delete tags
            rec = BAMRecord()
            @test !haskey(rec, "MN")
            rec["MN"] = 0x01
            @test rec["MN"] === 0x01
            @test haskey(rec, "MN")
            @test !haskey(rec, "XY")
            rec["XY"] = "foobar"
            @test rec["XY"] == "foobar"
            @test haskey(rec, "XY")
            delete!(rec, "MN")
            @test !haskey(rec, "MN")
        end

        @testset "Reader" begin
            reader = open(BAMReader, joinpath(bamdir, "ce#1.bam"))
            @test isa(reader, BAMReader)
            @test eltype(reader) === BAMRecord
            @test startswith(repr(reader), "Bio.Align.BAMReader{IOStream}:")

            # header
            h = header(reader)
            @test h["SQ"] == [Dict("SN" => "CHROMOSOME_I", "LN" => "1009800")]

            # first record
            rec = BAMRecord()
            read!(reader, rec)
            @test ismapped(rec)
            @test refname(rec) == "CHROMOSOME_I"
            @test refindex(rec) == 1
            @test leftposition(rec) == 2
            @test rightposition(rec) == 102
            @test seqname(rec) == "SRR065390.14978392"
            @test sequence(rec) == dna"""
            CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCT
            AAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
            """
            @test seqlength(rec) == 100
            @test eltype(qualities(rec)) == Int8
            @test qualities(rec) == [Int(x) - 33 for x in "#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"]
            @test flag(rec) == 16
            @test cigar(rec) == "27M1D73M"
            @test alignment(rec) == Alignment([
                AlignmentAnchor(  0,   1, OP_START),
                AlignmentAnchor( 27,  28, OP_MATCH),
                AlignmentAnchor( 27,  29, OP_DELETE),
                AlignmentAnchor(100, 102, OP_MATCH)])
            @test rec["XG"] == 1
            @test rec["XM"] == 5
            @test rec["XN"] == 0
            @test rec["XO"] == 1
            @test rec["AS"] == -18
            @test rec["XS"] == -18
            @test rec["YT"] == "UU"
            @test eof(reader)
            close(reader)

            # iterator
            @test length(collect(open(BAMReader, joinpath(bamdir, "ce#1.bam")))) == 1
            @test length(collect(open(BAMReader, joinpath(bamdir, "ce#2.bam")))) == 2

            # IOStream
            @test length(collect(BAMReader(open(joinpath(bamdir, "ce#1.bam"))))) == 1
            @test length(collect(BAMReader(open(joinpath(bamdir, "ce#2.bam"))))) == 2
        end

        @testset "Round trip" begin
            for specimen in YAML.load_file(joinpath(bamdir, "index.yml"))
                filepath = joinpath(bamdir, specimen["filename"])
                mktemp() do path, _
                    # copy
                    reader = open(BAMReader, filepath)
                    writer = BAMWriter(
                        BGZFStream(path, "w"),
                        header(reader, true))
                    records = BAMRecord[]
                    for rec in reader
                        push!(records, rec)
                        write(writer, rec)
                    end
                    close(reader)
                    close(writer)

                    # read again
                    reader = open(BAMReader, path)
                    @test collect(reader) == records
                    close(reader)
                end
            end
        end

        @testset "Random access" begin
            filepath = joinpath(bamdir, "GSE25840_GSM424320_GM06985_gencode_spliced.head.bam")
            reader = open(BAMReader, filepath, index=filepath * ".bai")

            # expected values are counted using samtools
            for (refname, interval, expected) in [
                    ("chr1", 1000:10000,      21),
                    ("chr1", 8000:10000,      20),
                    ("chr1", 766_000:800_000, 142),
                    ("chr1", 786_000:800_000, 1),
                    ("chr1", 796_000:800_000, 0)]
                intsect = intersect(reader, refname, interval)
                @test eltype(intsect) == BAMRecord
                @test count(_ -> true, intsect) == expected
                # check that the intersection iterator is stateless
                @test count(_ -> true, intsect) == expected
            end

            # randomized tests
            for n in 1:50
                refindex = 1
                refname = "chr1"
                range = randrange(1:1_000_000)
                seekstart(reader)
                # linear scan
                expected = collect(filter(reader) do rec
                    isoverlapping(rec, refindex, range)
                end)
                # indexed scan
                actual = collect(intersect(reader, refname, range))
                @test actual == expected
            end

            close(reader)
        end
    end
end

end # TestAlign


text:
@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$.
Use the same tableaus for $\tilde{a}^{1}$ and $\tilde{a}^{3}$, so that $\tilde{s} = s$, as well as

```math
\begin{aligned}
\begin{array}{c|cc}
& \tfrac{1}{2} b^{1} \\
& \vdots \\
& \tfrac{1}{2} b^{1} \\
\hline
a^{2} & \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} b^{3} \\
& \vdots \\
& \tfrac{1}{2} b^{3} \\
\hline
a^{4} & \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} b^{1} \\
c^{1} & \vdots \\
& \tfrac{1}{2} b^{1} \\
\hline
\tilde{a}^{2} & \tfrac{1}{2} (1 + R(\infty)) \, b^{1} \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} b^{3} \\
c^{3} & \vdots \\
& \tfrac{1}{2} b^{3} \\
\hline
\tilde{a}^{4} & \tfrac{1}{2} (1 + R(\infty)) \, b^{3} \\
\end{array}
\end{aligned}
```

Set $\omega = [0, ..., 0, 1]$ and
```math
\delta_{ij} = \begin{cases}
+1 & j = i , \\
-1 & j = \tilde{s} , \\
 0 & \text{else} ,
\end{cases}
```
so that $\Lambda_{1} = \Lambda_{2} = ... = \Lambda_{\tilde{s}}$.

This methods is constructed to satisfy the constraint on the projective stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$.
Note, however, that it violates the symplecticity conditions $b^{1}_{i} b^{4}_{j} = b^{1}_{i} a^{4}_{ij} + b^{4}_{j} \tilde{a}^{1}_{ji}$ and $b^{2}_{i} b^{3}_{j} = b^{2}_{i} \tilde{a}^{3}_{ij} + b^{3}_{j} a^{2}_{ji}$.

"""
function TableauVSPARKInternalProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    s = q.s
    o = min(q.o, p.o)

    α_q = zeros(T, s, s)
    α_p = zeros(T, s, s)

    for i in 1:s
        α_q[i,:] .= q.b ./ 2
        α_p[i,:] .= p.b ./ 2
    end

    β_q = q.b .* (1 + R∞) / 2
    β_p = p.b .* (1 + R∞) / 2

    c_λ = q.c
    d_λ = q.b
    ω_λ = zeros(T, 1, s+1)
    δ_λ = zeros(T, s-1, s)

    ω_λ[1,s+1] = 1

    for i in 1:s-1
        δ_λ[i,i] = +1
        δ_λ[i,s] = -1
    end

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            q.a, p.a, α_q, α_p,
                            q.a, p.a, α_q, α_p,
                            q.b, p.b, β_q, β_p,
                            q.c, p.c, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            q.a, p.a, α_q, α_p,
                            q.a, p.a, α_q, α_p,
                            q.b, p.b, β_q, β_p,
                            q.c, p.c, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauVSPARKGLRKpInternal(s)
    glrk = TableauGauss(s)
    TableauVSPARKInternalProjection(Symbol("GLRK($s)pInternal"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$, and set

```math
\begin{aligned}
\begin{array}{c|cc}
& \tfrac{1}{2} b^{1} \\
& \vdots \\
& \tfrac{1}{2} b^{1} \\
\hline
a^{2} & \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} b^{4} \\
& \vdots \\
& \tfrac{1}{2} b^{4} \\
\hline
a^{4} & \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} b^{1} \\
c^{1} & \vdots \\
& \tfrac{1}{2} b^{1} \\
\hline
\tilde{a}^{2} & \tfrac{1}{2} (1 + R(\infty)) \, b^{1} \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} b^{3} \\
c^{3} & \vdots \\
& \tfrac{1}{2} b^{3} \\
\hline
\tilde{a}^{4} & \tfrac{1}{2} (1 + R(\infty)) \, b^{3} \\
\end{array}
\end{aligned}
```

Note that by this definition $\tilde{s} = s$.
The coefficients $\tilde{a}^{1}$ and $\tilde{a}^{3}$ are determined by the (modified) symplecticity conditions, specifically $a^{4}_{ij} = b^{3}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{1}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$, where $b^{2}$ has been replaced with $b^{1}$ and $b^{4}$ with $b^{3}$, respectively.
Set $\omega = [0, ..., 0, 1]$ and
```math
\delta_{ij} = \begin{cases}
+1 & j = i , \\
-1 & j = \tilde{s} , \\
 0 & \text{else} ,
\end{cases}
```
so that $\Lambda_{1} = \Lambda_{2} = ... = \Lambda_{\tilde{s}}$.

Note that this method satisfies the symplecticity conditions $b^{1}_{i} b^{4}_{j} = b^{1}_{i} a^{4}_{ij} + b^{4}_{j} \tilde{a}^{1}_{ji}$ and $b^{2}_{i} b^{3}_{j} = b^{2}_{i} \tilde{a}^{3}_{ij} + b^{3}_{j} a^{2}_{ji}$ only if $R(\infty) = 1$ due to the definitions of $b^{2}$ and $b^{4}$.
Moreover, it does usually not satisfy the constraint on the projective stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$, exactly, but only approximately with bounded error, thus implying a residual in the symplecticity equation even if $R(\infty) = 1$.

"""
function TableauVSPARKModifiedInternalProjection(name, q, p, d=[]; R∞=1)
    T = eltype(q.a)

    @assert q.s == p.s

    s = q.s
    o = min(q.o, p.o)

    α_q = zeros(T, s, s)
    α_p = zeros(T, s, s)

    for i in 1:s
        α_q[i,:] .= q.b ./ 2
        α_p[i,:] .= p.b ./ 2
    end

    β_q = q.b .* (1 + R∞) / 2
    β_p = p.b .* (1 + R∞) / 2

#    a_q̃, a_p̃ = get_ã_vspark_primary(α_q, β_q, q.b, α_p, β_p, p.b)
    a_q̃, a_p̃ = get_ã_vspark_primary(α_q, q.b, q.b, α_p, p.b, p.b)

    c_λ = q.c
    d_λ = q.b
    ω_λ = zeros(T, 1, s+1)
    δ_λ = zeros(T, s-1, s)

    ω_λ[1,s+1] = 1

    for i in 1:s-1
        δ_λ[i,i] = +1
        δ_λ[i,s] = -1
    end

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            q.a, p.a, α_q, α_p,
                            a_q̃, a_p̃, α_q, α_p,
                            q.b, p.b, β_q, β_p,
                            q.c, p.c, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            q.a, p.a, α_q, α_p,
                            a_q̃, a_p̃, α_q, α_p,
                            q.b, p.b, β_q, β_p,
                            q.c, p.c, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauVSPARKGLRKpModifiedInternal(s)
    glrk = TableauGauss(s)
    TableauVSPARKModifiedInternalProjection(Symbol("GLRK($s)pModifiedInternal"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$.
For the projection, choose the tableau with $\tilde{s} = 1$ and $\rho = 0$, such that $\tilde{Q}_{n,1} = \tfrac{1}{2} ( q_{n} + q_{n+1})$, $\tilde{P}_{n,1} = \tfrac{1}{2} ( p_{n} + p_{n+1})$, i.e.,

```math
\begin{aligned}
\begin{array}{c|cc}
& \tfrac{1}{2} \\
& \vdots \\
& \tfrac{1}{2} \\
\hline
a^{2} & \\
\end{array}
&&
\begin{array}{c|cc}
& \tfrac{1}{2} \\
& \vdots \\
& \tfrac{1}{2} \\
\hline
a^{4} & \\
\end{array}
&&
\begin{array}{c|cc}
\tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{2} & \tfrac{1}{2} (1 + R(\infty))\\
\end{array}
&&
\begin{array}{c|cc}
\tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{4} & \tfrac{1}{2} (1 + R(\infty))\\
\end{array}
\end{aligned}
```

The coefficients $\tilde{a}^{1}$ and $\tilde{a}^{3}$ are determined by the symplecticity conditions, specifically $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$, and $\omega = [0, 1]$.

"""
function TableauVSPARKMidpointProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}
    @assert q.s == p.s
    s = q.s

    g = TableauGauss(1)
    α̃ = g.a
    β = g.b
    γ = g.c

    α = 0.5 * ones(T, s, 1)

    q_ã = zeros(T, 1, s)
    for i in 1:s
        q_ã[1,i] = q.b[i] / β[1] * ( β[1] - α[i,1] )
    end

    p_ã = zeros(T, 1, s)
    for i in 1:s
        p_ã[1,i] = p.b[i] / β[1] * ( β[1] - α[i,1] )
    end

    β .= (1 + R∞) / 2
    dλ = ones(T, 1)
    ω  = reshape(T[0  1], (1,2))
    δ  = zeros(T, 0, 1)

    if length(d) == 0
        return TableauVSPARKprimary(name, min(q.o, p.o),
                            q.a, p.a, α, α,
                            q_ã, p_ã, α̃, α̃,
                            q.b, p.b, β, β,
                            q.c, p.c, γ, dλ,
                            ω, δ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, min(q.o, p.o),
                            q.a, p.a, α, α,
                            q_ã, p_ã, α̃, α̃,
                            q.b, p.b, β, β,
                            q.c, p.c, γ, dλ,
                            ω, δ, d)
    end
end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and midpoint projection."
function TableauVSPARKLobattoIIIAIIIBpMidpoint(s)
    TableauVSPARKMidpointProjection(Symbol("LobattoIIIAIIIB($s)pMidpoint"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and midpoint projection."
function TableauVSPARKLobattoIIIBIIIApMidpoint(s)
    TableauVSPARKMidpointProjection(Symbol("LobattoIIIBIIIA($s)pMidpoint"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and midpoint projection."
function TableauVSPARKGLRKpMidpoint(s)
    glrk = TableauGauss(s)
    TableauVSPARKMidpointProjection(Symbol("GLRK($s)pMidpoint"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$.
For the projection, choose the tableau with $\tilde{s} = 1$ and $\rho = 0$, such that $\tilde{Q}_{n,1} = \tfrac{1}{2} ( q_{n} + q_{n+1})$, $\tilde{P}_{n,1} = \tfrac{1}{2} ( p_{n} + p_{n+1})$, i.e.,

```math
\begin{aligned}
\begin{array}{c|cc}
\tfrac{1}{2} & \tfrac{1}{2} b^{1} \\
\hline
\tilde{a}^{1} & \\
\end{array}
&&
\begin{array}{c|cc}
\tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{2} & \tfrac{1}{2} ( 1 + R (\infty) ) \\
\end{array}
&&
\begin{array}{c|cc}
\tfrac{1}{2} & \tfrac{1}{2} b^{3} \\
\hline
\tilde{a}^{3} & \\
\end{array}
&&
\begin{array}{c|cc}
\tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{4} & \tfrac{1}{2} ( 1 + R (\infty) ) \\
\end{array}
\end{aligned}
```

The coefficients $a^{2}$ and $a^{4}$ are determined by the symplecticity conditions, specifically $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$, and $\omega = [0, 1]$.

"""
function TableauVSPARKModifiedMidpointProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}
    @assert q.s == p.s
    s = q.s

    g = TableauGauss(1)
    α̃ = g.a
    β = g.b
    γ = g.c

    q_ã = reshape(q.b ./ 2, 1, s)
    p_ã = reshape(p.b ./ 2, 1, s)

    α_q = zeros(T, s, 1)
    α_p = zeros(T, s, 1)

    for i in 1:s
        α_q[i,1] = β[1] / p.b[i] * ( p.b[i] - p_ã[1,i] )
        α_p[i,1] = β[1] / q.b[i] * ( q.b[i] - q_ã[1,i] )
    end

    β .= (1 + R∞) / 2
    dλ = ones(T, 1)
    ω  = reshape(T[0  1], (1,2))
    δ  = zeros(T, 0, 1)

    if length(d) == 0
        return TableauVSPARKprimary(name, min(q.o, p.o),
                            q.a, p.a, α_q, α_p,
                            q_ã, p_ã, α̃, α̃,
                            q.b, p.b, β, β,
                            q.c, p.c, γ, dλ,
                            ω, δ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, min(q.o, p.o),
                            q.a, p.a, α_q, α_p,
                            q_ã, p_ã, α̃, α̃,
                            q.b, p.b, β, β,
                            q.c, p.c, γ, dλ,
                            ω, δ, d)
    end
end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and midpoint projection."
function TableauVSPARKLobattoIIIAIIIBpModifiedMidpoint(s)
    TableauVSPARKModifiedMidpointProjection(Symbol("LobattoIIIAIIIB($s)pModifiedMidpoint"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and midpoint projection."
function TableauVSPARKLobattoIIIBIIIApModifiedMidpoint(s)
    TableauVSPARKModifiedMidpointProjection(Symbol("LobattoIIIBIIIA($s)pModifiedMidpoint"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and midpoint projection."
function TableauVSPARKGLRKpModifiedMidpoint(s)
    glrk = TableauGauss(s)
    TableauVSPARKModifiedMidpointProjection(Symbol("GLRK($s)pModifiedMidpoint"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$.
For the projection, choose the tableau with $\tilde{s} = 2$ and $\rho = 1$, such that $\tilde{Q}_{n,1} = q_{n}$, $\tilde{Q}_{n,2} = q_{n+1}$, $\tilde{P}_{n,1} = p_{n}$, $\tilde{P}_{n,2} = p_{n+1}$, i.e.,

```math
\begin{aligned}
\begin{array}{c|cc}
0 & 0 \\
1 & b^{1} \\
\hline
\tilde{a}^{1} & \\
\end{array}
&&
\begin{array}{c|cc}
& 0 & 0 \\
& \tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{2} & \tfrac{1}{2} & \tfrac{1}{2} \\
\end{array}
&&
\begin{array}{c|cc}
0 & 0 \\
1 & b^{3} \\
\hline
\tilde{a}^{3} & \\
\end{array}
&&
\begin{array}{c|cc}
& 0 & 0 \\
& \tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{4} & \tfrac{1}{2} & \tfrac{1}{2} \\
\end{array}
\end{aligned}
```

The coefficients $a^{2}$ and $a^{4}$ are determined by the symplecticity conditions, specifically $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$.
Further choose $\omega = [1, 1, 0]$ and $\delta = [-1, R_{\infty}]$, so that $\tilde{\Lambda}_{n,1} = R_{\infty} \tilde{\Lambda}_{n,2}$ and
```math
\tilde{P}_{n,1} - \vartheta (\tilde{Q}_{n,1}) + R_{\infty} ( \tilde{P}_{n,2} - \vartheta (\tilde{Q}_{n,2}) ) = 0 .
```
Due to the particular choice of projective stages, this is equivalent to
```math
p_{n} - \vartheta (q_{n}) + R_{\infty} ( p_{n+1} - \vartheta (q_{n+1}) ) = 0 ,
```
so that the constraint $\phi(q_{n+1}, p_{n+1}) = 0$ is satisfied if $\phi(q_{n}, p_{n}) = 0$.
Note that the choice of $\tilde{a}^{2}$ and $\tilde{a}^{4}$ violates the symplecticity condition $b^{2}_{i} b^{4}_{j} = b^{2}_{i} \tilde{a}^{4}_{ij} + b^{4}_{j} \tilde{a}^{2}_{ji}$.

"""
function TableauVSPARKSymmetricProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    a_q = q.a
    a_p = p.a

    α_q = zeros(T, q.s, 2)
    α_q[:,1] .= 0.5

    α_p = zeros(T, p.s, 2)
    α_p[:,1] .= 0.5

    a_q̃ = vcat(zero(q.b)', q.b')
    a_p̃ = vcat(zero(p.b)', p.b')

    α_q̃ = [[0.0  0.0]
           [0.5  0.5]]
    α_p̃ = [[0.0  0.0]
           [0.5  0.5]]

    b_q = q.b
    b_p = p.b
    β_q = [0.5, 0.5]
    β_p = [0.5, 0.5]

    c_q = q.c
    c_p = p.c
    c_λ = [ 0.0, 1.0]
    d_λ = [ 0.5, 0.5*R∞]

    ω_λ = reshape(T[1  1  0], (1,3))
    δ_λ = reshape(T[-1   R∞], (1,2))

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end
end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and symmetric projection."
function TableauVSPARKLobattoIIIAIIIBpSymmetric(s)
    TableauVSPARKSymmetricProjection(Symbol("LobattoIIIAIIIB($s)pSymmetric"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and symmetric projection."
function TableauVSPARKLobattoIIIBIIIApSymmetric(s)
    TableauVSPARKSymmetricProjection(Symbol("LobattoIIIBIIIA($s)pSymmetric"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauVSPARKGLRKpSymmetric(s)
    glrk = TableauGauss(s)
    TableauVSPARKSymmetricProjection(Symbol("GLRK($s)pSymmetric"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$.
For the projection, choose the Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{4}, b^{4})$ and $(\tilde{a}^{2}, b^{2})$, respectively, and choose $\tilde{a}^{1}$ and $\tilde{a}^{3}$ such that the projective stages correspond to the initial condition and the solution, i.e.,

```math
\begin{aligned}
\begin{array}{c|cc}
0 & 0 \\
1 & b^{1} \\
\hline
\tilde{a}^{1} & \\
\end{array}
&&
\begin{array}{c|cc}
0 & \tfrac{1}{2} & 0 \\
1 & \tfrac{1}{2} & 0 \\
\hline
\tilde{a}^{2} & \\
\end{array}
&&
\begin{array}{c|cc}
0 & 0 \\
1 & b^{3} \\
\hline
\tilde{a}^{3} & \\
\end{array}
&&
\begin{array}{c|cc}
0 & 0 & 0 \\
1 & \tfrac{1}{2} & \tfrac{1}{2} \\
\hline
\tilde{a}^{4} & \\
\end{array}
\end{aligned}
```

and compute $a^{2}$ and $a^{4}$ by the symplecticity conditions, that is $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$ and $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$.
Finally choose $\omega = [0, 0, 1]$ and $\delta = [-1, R_{\infty}]$, implying that $\rho = 1$.
By construction, this method satisfies all symplecticity conditions, but the constraint on the projection stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$, is not satisfied exactly, but only approximately, although with bounded error.

"""
function TableauVSPARKLobattoIIIAIIIBProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    loba = TableauLobattoIIIA(2)
    lobb = TableauLobattoIIIB(2)

    a_q = q.a
    a_p = p.a
    b_q = q.b
    b_p = p.b

    a_q̃ = vcat(zero(q.b)', q.b')
    a_p̃ = vcat(zero(p.b)', p.b')

    α_q̃ = loba.a
    α_p̃ = lobb.a

    β_q = loba.b
    β_p = lobb.b

    α_q, α_p = get_α_vspark_primary(a_q̃, b_q, β_q, a_p̃, b_p, β_p)

    c_q = q.c
    c_p = p.c
    c_λ = [ 0.0, 1.0]
    d_λ = [ 0.5, 0.5]

    ω_λ = reshape(T[0  0  1], (1,3))
    δ_λ = reshape(T[-1   R∞], (1,2))

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIA-IIIB projection."
function TableauVSPARKLobattoIIIAIIIBpLobattoIIIAIIIB(s)
    TableauVSPARKLobattoIIIAIIIBProjection(Symbol("LobattoIIIAIIIB($s)pLobattoIIIAIIIB(2)"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIA-IIIB projection."
function TableauVSPARKLobattoIIIBIIIApLobattoIIIAIIIB(s)
    TableauVSPARKLobattoIIIAIIIBProjection(Symbol("LobattoIIIBIIIA($s)pLobattoIIIAIIIB(2)"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and Lobatto-IIIA-IIIB projection."
function TableauVSPARKGLRKpLobattoIIIAIIIB(s)
    glrk = TableauGauss(s)
    TableauVSPARKLobattoIIIAIIIBProjection(Symbol("GLRK($s)pLobattoIIIAIIIB(2)"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

This methods is the same as `TableauVSPARKLobattoIIIAIIIBProjection`, except for using Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{2}, b^{2})$, and $(\tilde{a}^{4}, b^{4})$ respectively, instead of the other way around.

"""
function TableauVSPARKLobattoIIIBIIIAProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    loba = TableauLobattoIIIA(2)
    lobb = TableauLobattoIIIB(2)

    a_q = q.a
    a_p = p.a
    b_q = q.b
    b_p = p.b

    a_q̃ = vcat(zero(q.b)', q.b')
    a_p̃ = vcat(zero(p.b)', p.b')

    α_q̃ = lobb.a
    α_p̃ = loba.a

    β_q = lobb.b
    β_p = loba.b

    α_q, α_p = get_α_vspark_primary(a_q̃, b_q, β_q, a_p̃, b_p, β_p)

    c_q = q.c
    c_p = p.c
    c_λ = [ 0.0, 1.0]
    d_λ = [ 0.5, 0.5]

    ω_λ = reshape(T[0  0  1], (1,3))
    δ_λ = reshape(T[-1   R∞], (1,2))

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIB-IIIA projection."
function TableauVSPARKLobattoIIIAIIIBpLobattoIIIBIIIA(s)
    TableauVSPARKLobattoIIIBIIIAProjection(Symbol("LobattoIIIAIIIB($s)pLobattoIIIBIIIA(2)"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIB-IIIA projection."
function TableauVSPARKLobattoIIIBIIIApLobattoIIIBIIIA(s)
    TableauVSPARKLobattoIIIBIIIAProjection(Symbol("LobattoIIIBIIIA($s)pLobattoIIIBIIIA(2)"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and Lobatto-IIIB-IIIA projection."
function TableauVSPARKGLRKpLobattoIIIBIIIA(s)
    glrk = TableauGauss(s)
    TableauVSPARKLobattoIIIBIIIAProjection(Symbol("GLRK($s)pLobattoIIIBIIIA(2)"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$.
For the projection, choose the Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{4}, b^{4})$ and $(\tilde{a}^{2}, b^{2})$, respectively.

The coefficients $\tilde{a}^{1}$ and $\tilde{a}^{3}$ are determined by the relations
```math
\begin{aligned}
\sum \limits_{j=1}^{s} \tilde{a}^{1}_{ij} (c_{j}^{1})^{k-1} &= \frac{(c_{i}^{2})^k}{k} , \qquad &
\sum \limits_{j=1}^{s} \tilde{a}^{3}_{ij} (c_{j}^{3})^{k-1} &= \frac{(c_{i}^{4})^k}{k} , \qquad &
i &= 1 , \, ... , \, \tilde{s} , \qquad &
k &= 1 , \, ... , \, s .
\end{aligned}
```

The coefficients $a^{2}$ and $a^{4}$ by the symplecticity conditions, that is $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$ and $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$.
Finally choose $\omega = [0, 0, 1]$ and $\delta = [-1, R_{\infty}]$, implying that $\rho = 1$.
By construction, this method satisfies all symplecticity conditions, but the constraint on the projection stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$, is not satisfied exactly, but only approximately, although with bounded error.

"""
function TableauVSPARKModifiedLobattoIIIAIIIBProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    s = q.s
    o = min(q.o, p.o)

    loba = TableauLobattoIIIA(2)
    lobb = TableauLobattoIIIB(2)

    a_q = q.a
    a_p = p.a
    b_q = q.b
    b_p = p.b

    β_q = loba.b
    β_p = lobb.b

    a_q̃ = get_lobatto_glrk_coefficients(s, 2, T).a
    a_p̃ = get_lobatto_glrk_coefficients(s, 2, T).a

    α_q̃ = loba.a
    α_p̃ = lobb.a

    α_q, α_p = get_α_vspark_primary(a_q̃, b_q, β_q, a_p̃, b_p, β_p)

    c_q = q.c
    c_p = p.c
    c_λ = loba.c
    d_λ = loba.b

    ω_λ = reshape(T[0  0  1], (1,3))
    δ_λ = reshape(T[-1   R∞], (1,2))

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIA-IIIB projection."
function TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIAIIIB(s)
    TableauVSPARKModifiedLobattoIIIAIIIBProjection(Symbol("LobattoIIIAIIIB($s)pModifiedLobattoIIIAIIIB(2)"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIA-IIIB projection."
function TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIAIIIB(s)
    TableauVSPARKModifiedLobattoIIIAIIIBProjection(Symbol("LobattoIIIBIIIA($s)pModifiedLobattoIIIAIIIB(2)"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauVSPARKGLRKpModifiedLobattoIIIAIIIB(s)
    glrk = TableauGauss(s)
    TableauVSPARKModifiedLobattoIIIAIIIBProjection(Symbol("GLRK($s)pModifiedLobattoIIIAIIIB(2)"), glrk, glrk; R∞=(-1)^s)
end

@doc raw"""

This methods is the same as `TableauVSPARKModifiedLobattoIIIAIIIBProjection`, except for using Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{2}, b^{2})$, and $(\tilde{a}^{4}, b^{4})$ respectively, instead of the other way around.

"""
function TableauVSPARKModifiedLobattoIIIBIIIAProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    s = q.s
    o = min(q.o, p.o)

    loba = TableauLobattoIIIA(2)
    lobb = TableauLobattoIIIB(2)

    a_q = q.a
    a_p = p.a
    b_q = q.b
    b_p = p.b

    β_q = loba.b
    β_p = lobb.b

    a_q̃ = get_lobatto_glrk_coefficients(s, 2, T).a
    a_p̃ = get_lobatto_glrk_coefficients(s, 2, T).a

    α_q̃ = lobb.a
    α_p̃ = loba.a

    α_q, α_p = get_α_vspark_primary(a_q̃, b_q, β_q, a_p̃, b_p, β_p)

    c_q = q.c
    c_p = p.c
    c_λ = loba.c
    d_λ = loba.b

    ω_λ = reshape(T[0  0  1], (1,3))
    δ_λ = reshape(T[-1   R∞], (1,2))

    if length(d) == 0
        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIB-IIIA projection."
function TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIBIIIA(s)
    TableauVSPARKModifiedLobattoIIIBIIIAProjection(Symbol("LobattoIIIAIIIB($s)pModifiedLobattoIIIBIIIA(2)"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIB-IIIA projection."
function TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIBIIIA(s)
    TableauVSPARKModifiedLobattoIIIBIIIAProjection(Symbol("LobattoIIIBIIIA($s)pModifiedLobattoIIIBIIIA(2)"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauVSPARKGLRKpModifiedLobattoIIIBIIIA(s)
    glrk = TableauGauss(s)
    TableauVSPARKModifiedLobattoIIIBIIIAProjection(Symbol("GLRK($s)pModifiedLobattoIIIBIIIA(2)"), glrk, glrk; R∞=(-1)^s)
end

function TableauVSPARKSymplecticProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=+1) where {T}

    la = TableauLobattoIIIA(2)
    lb = TableauLobattoIIIB(2)

    @assert q.s == p.s
    @assert q.b == p.b
    @assert q.c == p.c

    @assert la.s == lb.s
    @assert la.b == lb.b

    s = q.s
    σ = la.s
    ρ = 0

    a_q = q.a
    a_p = p.a
    b_q = q.b
    b_p = p.b

    β_q = [0.5, R∞*0.5]
    β_p = [0.5, R∞*0.5]

    α_q = zeros(T, s, 2)
    α_q[:,1] .= 0.5

    α_p = zeros(T, s, 2)
    α_p[:,1] .= 0.5

    a_q̃ = zeros(T, σ, s)
    a_p̃ = zeros(T, σ, s)
    for i in 1:σ
        for j in 1:s
            a_q̃[i,j] = b_q[j] / β_p[i] * (β_p[i] - α_p[j,i])
            a_p̃[i,j] = b_p[j] / β_q[i] * (β_q[i] - α_q[j,i])
        end
    end

    α_q̃ = la.a
    α_p̃ = lb.a

    c_q = q.c
    c_p = p.c
    c_λ = la.c
    d_λ = [0.5, 0.5]

    ω_λ = [0.5 0.5 0.0
           0.0 0.0 1.0]

    δ_λ = zeros(T, ρ, σ)

    if length(d) == 0
        return TableauVSPARKprimary(name, min(q.o, p.o),
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == s

        return TableauVSPARKprimary(name, min(q.o, p.o),
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

function TableauVSPARKGLRKpSymplectic(s)
    glrk = TableauGauss(s)
    TableauVSPARKSymplecticProjection(Symbol("GLRK($s)pSymplectic"), glrk, glrk; R∞=(-1)^s)
end

function TableauVSPARKLobABCCD(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)

    a_q = lobd.a
    a_p = lobd.a

    α_q = lobb.a
    α_p = lobb.a

    a_q̃ = loba.a
    a_p̃ = loba.a

    α_q̃ = lob.a
    α_p̃ = lobc.a

    b_q = lobd.b
    b_p = lobd.b

    β_q̃ = lobd.b
    β_p̃ = lobd.b

    c_q = lobd.c
    c_p = lobd.c

    γ = lobd.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobABCCD($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobABCCE(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobe.a
    a_p = lobe.a

    α_q = lobb.a
    α_p = lobb.a

    a_q̃ = loba.a
    a_p̃ = loba.a

    α_q̃ = lob.a
    α_p̃ = lobc.a

    b_q = lobd.b
    b_p = lobd.b

    β_q̃ = lobd.b
    β_p̃ = lobd.b

    c_q = lobd.c
    c_p = lobd.c

    γ = lobd.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobABCCE($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobABED(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobe.a
    a_p = lobe.a

    α_q = lobb.a
    α_p = lobb.a

    a_q̃ = loba.a
    a_p̃ = loba.a

    α_q̃ = lobd.a
    α_p̃ = lobd.a

    b_q = lob.b
    b_p = lob.b

    β_q̃ = lob.b
    β_p̃ = lob.b

    c_q = lob.c
    c_p = lob.c

    γ   = lob.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobABED($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobABDE(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobd.a
    a_p = lobd.a

    α_q = lobb.a
    α_p = lobb.a

    a_q̃ = loba.a
    a_p̃ = loba.a

    α_q̃ = lobe.a
    α_p̃ = lobe.a

    b_q = lob.b
    b_p = lob.b

    β_q̃ = lob.b
    β_p̃ = lob.b

    c_q = lob.c
    c_p = lob.c

    γ   = lob.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobABDE($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobABD(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobd.a
    a_p = lobd.a

    α_q = lobb.a
    α_p = lobb.a

    a_q̃ = loba.a
    a_p̃ = loba.a

    α_q̃ = loba.a
    α_p̃ = lobb.a

    b_q = lob.b
    b_p = lob.b

    β_q̃ = lob.b
    β_p̃ = lob.b

    c_q = lob.c
    c_p = lob.c

    γ   = lob.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobABD($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobABE(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobe.a
    a_p = lobe.a

    α_q = lobb.a
    α_p = lobb.a

    a_q̃ = loba.a
    a_p̃ = loba.a

    α_q̃ = loba.a
    α_p̃ = lobb.a

    b_q = lob.b
    b_p = lob.b

    β_q̃ = lob.b
    β_p̃ = lob.b

    c_q = lob.c
    c_p = lob.c

    γ   = lob.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobABE($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobDE(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobd.a
    a_p = lobd.a

    α_q = lobe.a
    α_p = lobe.a

    a_q̃ = lobe.a
    a_p̃ = lobe.a

    α_q̃ = lobd.a
    α_p̃ = lobd.a

    b_q = lob.b
    b_p = lob.b

    β_q̃ = lob.b
    β_p̃ = lob.b

    c_q = lob.c
    c_p = lob.c

    γ   = lob.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobDE($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

function TableauVSPARKLobED(s=2, T=Float64)
    lob  = TableauLobattoIIIC̄(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    lobd = TableauLobattoIIID(s)
    lobe = TableauLobattoIIIE(s)

    a_q = lobe.a
    a_p = lobe.a

    α_q = lobd.a
    α_p = lobd.a

    a_q̃ = lobd.a
    a_p̃ = lobd.a

    α_q̃ = lobe.a
    α_p̃ = lobe.a

    b_q = lob.b
    b_p = lob.b

    β_q̃ = lob.b
    β_p̃ = lob.b

    c_q = lob.c
    c_p = lob.c

    γ   = lob.c

    d = ones(T, s)
    ω = get_lobatto_ω_matrix(s)
    δ = zeros(T, 0, s)

    return TableauVSPARKprimary(Symbol("LobED($s)"), 2*s-2,
                        a_q, a_p, α_q, α_p,
                        a_q̃, a_p̃, α_q̃, α_p̃,
                        b_q, b_p, β_q̃, β_p̃,
                        c_q, c_p, γ, d,
                        ω, δ)
end

ext: jl


text:
"""
    mul!!(C, A, B, [α, β]) -> C′
"""
mul!!(C, A, B) = may(mul!, C, A, B)
mul!!(C, A, B, α, β) = may(mul!, C, A, B, α, β)

pure(::typeof(mul!)) = NoBang.mul
_asbb(::typeof(mul!)) = mul!!
possible(::typeof(mul!), C, A, B) =
    ismutable(C) && _matmuleltype(A, B) <: eltype(C)
possible(::typeof(mul!), C, A, B, α, β) =
    ismutable(C) && _matmuleltype(C, A, B, α, β) <: eltype(C)

# Estimate `eltype` of `C`.  This is how it's done in LinearAlgebra.jl
# but maybe it's better to use the approach of
# https://github.com/tpapp/AlgebraResultTypes.jl ?
_matprod(x, y) = x * y + x * y
_matprod(c, a, b, α, β) = a * b * α + a * b * α + c * β
_matmuleltype(A, B) = Base.promote_op(_matprod, eltype(A), eltype(B))
_matmuleltype(C, A, B, α, β) =
    Base.promote_op(_matprod, eltype(C), eltype(A), eltype(B), typeof(α), typeof(β))

"""
    lmul!!(A, B) -> B′
"""
lmul!!(A, B) = may(lmul!, A, B)

pure(::typeof(lmul!)) = *
_asbb(::typeof(lmul!)) = lmul!!
possible(::typeof(lmul!), A, B) =
    ismutable(B) && _matmuleltype(A, B) <: eltype(B)

"""
    rmul!!(A, B) -> A′
"""
rmul!!(A, B) = may(rmul!, A, B)

pure(::typeof(rmul!)) = *
_asbb(::typeof(rmul!)) = rmul!!
possible(::typeof(rmul!), A, B) =
    ismutable(A) && _matmuleltype(A, B) <: eltype(A)


text:
"""
    true_positives(targets, outputs, [encoding]) -> Union{Int, Dict}

Count how many positive predicted outcomes in `outputs` are also
marked as positive outcomes in `targets`. Which value denotes
"positive" depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_positive`](@ref),
[`condition_positive`](@ref),
[`true_positive_rate`](@ref)

# Examples

```jldoctest
julia> true_positives(1, 1, LabelEnc.ZeroOne()) # single observation
1

julia> true_positives([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
2

julia> true_positives([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 1
  :b => 2
  :c => 0
```
"""
true_positives(target, output, encoding::BinaryLabelEncoding) =
    Int(isposlabel(target, encoding) & isposlabel(output, encoding))

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

"""
    true_negatives(targets, outputs, [encoding]) -> Union{Int, Dict}

Count how many negative predicted outcomes in `outputs` are also
marked as negative outcomes in `targets`. Which value(s) denote
"negative" depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_negative`](@ref),
[`condition_negative`](@ref),
[`true_negative_rate`](@ref)

# Examples

```jldoctest
julia> true_negatives(0, 0, LabelEnc.ZeroOne()) # single observation
1

julia> true_negatives([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
1

julia> true_negatives([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 2
  :b => 3
  :c => 4
```
"""
true_negatives(target, output, encoding::BinaryLabelEncoding) =
    Int(isneglabel(target, encoding) & isneglabel(output, encoding))

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

"""
    false_positives(targets, outputs, [encoding]) -> Union{Int, Dict}

Count how many positive predicted outcomes in `outputs` are
actually marked as negative outcomes in `targets`. These
occurrences are also known as `type_1_errors`. Which value
denotes "positive" depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_positive`](@ref),
[`condition_negative`](@ref),
[`false_positive_rate`](@ref)

# Examples

```jldoctest
julia> false_positives(0, 1, LabelEnc.ZeroOne()) # single observation
1

julia> false_positives([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
1

julia> false_positives([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 2
  :b => 1
  :c => 0
```
"""
false_positives(target, output, encoding::BinaryLabelEncoding) =
    Int(isneglabel(target, encoding) & isposlabel(output, encoding))

const type_1_errors = false_positives

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

"""
    false_negatives(targets, outputs, [encoding]) -> Union{Int, Dict}

Count how many negative predicted outcomes in `outputs` are
actually marked as positive outcomes in `targets`. These
occurrences are also known as `type_2_errors`. Which value
denotes "positive" depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_negative`](@ref),
[`condition_positive`](@ref),
[`false_negative_rate`](@ref)

# Examples

```jldoctest
julia> false_negatives(1, 0, LabelEnc.ZeroOne()) # single observation
1

julia> false_negatives([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
1

julia> false_negatives([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 1
  :b => 0
  :c => 2
```
"""
false_negatives(target, output, encoding::BinaryLabelEncoding) =
    Int(isposlabel(target, encoding) & isneglabel(output, encoding))

const type_2_errors = false_negatives

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

"""
    condition_positive(targets, outputs, [encoding]) -> Union{Int, Dict}

Count the number of positive outcomes in `targets`. Which value
denotes "positive" depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_positive`](@ref),
[`condition_negative`](@ref),
[`prevalence`](@ref)

# Examples

```jldoctest
julia> condition_positive(1, 0, LabelEnc.ZeroOne()) # single observation
1

julia> condition_positive([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
3

julia> condition_positive([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 2
  :b => 2
  :c => 2
```
"""
condition_positive(target, output, encoding::BinaryLabelEncoding) =
    Int(isposlabel(target, encoding))

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

"""
    condition_negative(targets, outputs, [encoding]) -> Union{Int, Dict}

Count the number of negative outcomes in `targets`. Which values
denote "negative" depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_negative`](@ref),
[`condition_positive`](@ref)

# Examples

```jldoctest
julia> condition_negative(0, 1, LabelEnc.ZeroOne()) # single observation
1

julia> condition_negative([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
2

julia> condition_negative([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 4
  :b => 4
  :c => 4
```
"""
condition_negative(target, output, encoding::BinaryLabelEncoding) =
    Int(isneglabel(target, encoding))

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

"""
    predicted_positive(targets, outputs, [encoding]) -> Union{Int, Dict}

Count the number of positive predicted outcomes in `outputs`.
Which value denotes "positive" depends on the given (or inferred)
`encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_negative`](@ref),
[`condition_positive`](@ref)

# Examples

```jldoctest
julia> predicted_positive(0, 1, LabelEnc.ZeroOne()) # single observation
1

julia> predicted_positive([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
3

julia> predicted_positive([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 3
  :b => 3
  :c => 0
```
"""
predicted_positive(target, output, encoding::BinaryLabelEncoding) =
    Int(isposlabel(output, encoding))

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

"""
    predicted_negative(targets, outputs, [encoding]) -> Union{Int, Dict}

Count the number of negative predicted outcomes in `outputs`.
Which values denote "negative" depends on the given (or inferred)
`encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`predicted_positive`](@ref),
[`condition_negative`](@ref)

# Examples

```jldoctest
julia> predicted_negative(1, 0, LabelEnc.ZeroOne()) # single observation
1

julia> predicted_negative([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
2

julia> predicted_negative([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 3
  :b => 3
  :c => 6
```
"""
predicted_negative(target, output, encoding::BinaryLabelEncoding) =
    Int(isneglabel(output, encoding))

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

"""
    correctly_classified(targets, outputs, [encoding]) -> Union{Int, Dict}

Count the number of predicted outcomes in `outputs` that agree
with the expected outcomes in `targets` under a two-class
interpretation. Which value(s) denote "positive" or "negative"
depends on the given (or inferred) `encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`incorrectly_classified`](@ref),
[`accuracy`](@ref)

# Examples

```jldoctest
julia> correctly_classified(0, 0, LabelEnc.ZeroOne()) # single observation
1

julia> correctly_classified([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
3

julia> correctly_classified([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 3
  :b => 5
  :c => 4
```
"""
correctly_classified(target, output, encoding::BinaryLabelEncoding) =
    Int(isposlabel(target, encoding) == isposlabel(output, encoding))

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

"""
    incorrectly_classified(targets, outputs, [encoding]) -> Union{Int, Dict}

Count the number of predicted outcomes in `outputs` that are
misclassified according to the expected outcomes in `targets`
under a two-class interpretation. Which value(s) denote
"positive" or "negative" depends on the given (or inferred)
`encoding`.
$SCALAR_DESC

$COUNT_ENCODING_DESCR

# Arguments

$COUNT_ARGS

# See also

[`correctly_classified`](@ref),
[`accuracy`](@ref)

# Examples

```jldoctest
julia> incorrectly_classified(0, 1, LabelEnc.ZeroOne()) # single observation
1

julia> incorrectly_classified([1,0,1,1,0], [1,1,1,0,0]) # multiple observations
2

julia> incorrectly_classified([:a,:a,:b,:b,:c,:c], [:a,:b,:b,:b,:a,:a]) # multi-class
Dict{Symbol,$Int} with 3 entries:
  :a => 3
  :b => 1
  :c => 2
```
"""
incorrectly_classified(target, output, encoding::BinaryLabelEncoding) =
    Int(xor(isposlabel(target, encoding), isposlabel(output, encoding)))

const misclassified = incorrectly_classified

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

# only for internal use to define `prevalence`
_length_targets(target, output, encoding) = 1
_length_targets(target::AbstractArray, output, encoding) = length(targets)

# --------------------------------------------------------------------
# Generate common fallback functions
for fun in (:true_positives,  :true_negatives,
            :false_positives, :false_negatives,
            :condition_positive, :condition_negative,
            :predicted_positive, :predicted_negative,
            :correctly_classified, :incorrectly_classified)

    # Convenience syntax for using native labels
    @eval function ($fun)(targets, outputs, encoding::AbstractVector)
        ($fun)(targets, outputs, LabelEnc.NativeLabels(encoding))
    end

    # Generic fallback. Tries to infer label encoding
    @eval function ($fun)(targets, outputs)
        encoding = _labelenc(targets, outputs)
        ($fun)(targets, outputs, encoding)
    end

    # BinaryLabelEncoding: Generate shared accumulator
    @eval function ($fun)(
            targets::AbstractArray,
            outputs::AbstractArray,
            encoding::BinaryLabelEncoding)
        @_dimcheck length(targets) == length(outputs)
        result::Int = 0
        @inbounds for I in eachindex(targets, outputs)
            result += ($fun)(targets[I], outputs[I], encoding)
        end
        result
    end

    # Multiclass LabelEncoding: Implementation for scalar
    @eval function ($fun)(target, output, encoding::LabelEncoding)
        labels = label(encoding)
        result = Dict{eltype(labels),Int}()
        for l in labels
            result[l] = ($fun)(target, output, LabelEnc.OneVsRest(l))
        end
        result
    end

    # Multiclass LabelEncoding: Generate shared accumulator
    @eval function ($fun)(
            targets::AbstractArray,
            outputs::AbstractArray,
            encoding::LabelEncoding)
        @_dimcheck length(targets) == length(outputs)
        labels = label(encoding)
        results = zeros(Int, length(labels))
        ovr = [LabelEnc.OneVsRest(l) for l in labels]
        @inbounds for I in eachindex(targets, outputs)
            target = targets[I]
            output = outputs[I]
            for j in eachindex(results, ovr)
                results[j] += ($fun)(target, output, ovr[j])
            end
        end
        Dict(Pair.(labels, results))
    end

    # Disable support for OneOfK (otherwise behaviour is misleading)
    @eval function ($fun)(target, output, encoding::LabelEnc.OneOfK)
        throw(ArgumentError("encoding LabelEnc.OneOfK not yet supported"))
    end

    @eval function ($fun)(targets::AbstractArray, outputs::AbstractArray, encoding::LabelEnc.OneOfK)
        throw(ArgumentError("encoding LabelEnc.OneOfK not yet supported"))
    end
end


text: # utility functions for interaction with GA API endpoints


text:
include("../../global_curved.jl")
include("odefun.jl")
using Plots

function main()
  sim_years = 1000.

  Vp = 1e-9 # plate rate
  ρ = 2.670
  cs = 3.464
  σn = 50
  RSamin = 0.01
  RSamax = 0.025
  RSb = 0.015
  RSDc = 0.016
  RSf0 = 0.6
  RSV0 = 1e-6
  RSVinit = 1e-9
  RSH1 = 15
  RSH2 = 18

  μshear = cs^2 * ρ
  η = μshear / (2 * cs)

  # N + 1 is the number of grid points in each dimension
  N = 200
  δNp = N + 1

  # SBP interior order
  SBPp   = 2

  # mesh file side set type to actually boundary condition type
  bc_map = [BC_DIRICHLET, BC_DIRICHLET, BC_NEUMANN, BC_NEUMANN,
            BC_JUMP_INTERFACE]

  (verts, EToV, EToF, FToB, EToDomain) = read_inp_2d("meshes/1_1_block.inp";
            bc_map=bc_map)

   # number of elements and faces
   (nelems, nfaces) = (size(EToV, 2), size(FToB, 1))
   @show (nelems, nfaces)

   # EToN0 is the base mesh size (e.g., before refinement)
   EToN0 = zeros(Int64, 2, nelems)
   EToN0[1, :] .= N
   EToN0[2, :] .= N

   # Exact solution
   Lx = 80
   Ly = 80

     # Set up the local grid dimensions
     Nr = EToN0[1, :]
     Ns = EToN0[2, :]

     # Dictionary to store the operators
     OPTYPE = typeof(locoperator(2, 8, 8))
     lop = Dict{Int64,OPTYPE}() # Be extra careful about the () here

     el_x = 10e12 # set it to be infinity to have even spread
     el_y = 10e12
     xt = (r,s) -> (el_x .* tan.(atan((Lx )/el_x).* (0.5*r .+ 0.5))  , el_x .* sec.(atan((Lx )/el_x).* (0.5*r .+ 0.5)).^2 * atan((Lx)/el_x) * 0.5 ,zeros(size(s)))
     yt = (r,s) -> (el_y .* tan.(atan((Ly )/el_y).* (0.5*s .+ 0.5))  , zeros(size(r)), el_y .* sec.(atan((Ly )/el_y).*(0.5*s .+ 0.5)) .^2 * atan((Ly )/el_y) * 0.5 )

     # create metrics
     metrics = create_metrics(SBPp, Nr[1], Ns[1], xt, yt) # not quite sure about this part

     # create local operator
     LFtoB = [BC_DIRICHLET,BC_DIRICHLET,BC_NEUMANN,BC_NEUMANN]
     lop = Dict{Int64,OPTYPE}() # This step to create a dict is essential
     lop[1] = locoperator(SBPp, Nr[1], Ns[1], metrics, LFtoB) # this function might not be correct

     # obtain M
     factorization = (x) -> cholesky(Symmetric(x))
     M = SBPLocalOperator1(lop, Nr[1], Ns[1], factorization)

     # obtain ge that stores boundary data
     ge = zeros((Nr[1]+1) * (Ns[1]+1))
     e = 1

     # boundary conditions
     bc_Dirichlet = (lf, x, y, e) -> zeros(size(x))
     bc_Neumann   = (lf, x, y, nx, ny, e) -> zeros(size(x))
     # set right-hand side and solve
     locbcarray_mod!(ge, lop[e], LFtoB, bc_Dirichlet, bc_Neumann,(e))
     u = M.F[e] \ ge

  Δτ = zeros(N+1)

  # Assemble fault variables/data
  RSa = zeros(N+1)

  xf = lop[1].facecoord[1][1]
  yf = lop[1].facecoord[2][1]
  for n = 1:N+1
      RSa[n] = RSamin - (RSamin - RSamax) *
        min(1, max(0, (RSH1 - yf[n])/(RSH1 - RSH2)))
  end

  τz0 = σn * RSamax * asinh(RSVinit / (2 * RSV0) *
                                 exp.((RSf0 + RSb * log.(RSV0 / RSVinit)) /
                                      RSamax)) + η * RSVinit

  θ = (RSDc ./ RSV0) .* exp.((RSa ./ RSb) .* log.((2 .* RSV0 ./ RSVinit) .*
      sinh.((τz0 .- η .* RSVinit) ./ (RSa .* σn))) .- RSf0 ./ RSb)

  ψ0 = RSf0 .+ RSb .* log.(RSV0 .* θ ./ RSDc)

  ψδ = zeros(2δNp)
  ψδ[1:δNp] .= ψ0

  #δ = zeros(δNp)
  #pyplot()
  #display(plot(yf, δ))
  #sleep(1)

  # set up parameters sent to right hand side
  odeparam = (reject_step = [false],
              Vp=Vp,
              lop=lop,
              F = M.F[e],
              u=u,
              Δτ = Δτ,
              ge = ge,
              μshear=μshear,
              RSa=RSa,
              RSb=RSb,
              σn=σn,
              η=η,
              RSV0=RSV0,
              τz0=τz0,
              RSDc=RSDc,
              RSf0=RSf0,
              LFtoB = LFtoB
             )

  dψV = zeros(2δNp)
  tspan = (0, sim_years * year_seconds)
  prob = ODEProblem(odefun, ψδ, tspan, odeparam)
  function stepcheck(_, p, _)
    if p.reject_step[1]
      p.reject_step[1] = false
      println("reject")
      return true
    end
    return false
  end

  sol = solve(prob, Tsit5(); isoutofdomain=stepcheck, dt=year_seconds,
              atol = 1e-5, rtol = 1e-3, save_everystep=true,
              internalnorm=(x, _)->norm(x, Inf))

  return (sol, yf)
end

function plot_slip(S, yf, stride_time)

  m = length(yf)
  no_time_steps = size(S.t)
  slip_final = S.u[end][end]

  for i = 1:stride_time:no_time_steps[1]

    slip_t = S.u[i][m+1:end] # slip at time t
    #pyplot()
    display(plot(slip_t, -yf, xtickfont=font(18),
    ytickfont=font(18),
    guidefont=font(18),
    legendfont=font(18), ylabel = "Depth (km)", xlabel = "Slip (m)", xlims = (0, slip_final)))
    sleep(0.1)
  end

  #nothing
end

(S, yf) = main()
plot_slip(S, yf, 10)

max_stars_repo_name: brittany-erickson/HybridSBP

max_issues_repo_name: brittany-erickson/HybridSBP

max_forks_repo_name: brittany-erickson/HybridSBP


text:
try
    eval(Expr(:import,:Lists))
catch err
    Pkg.clone("https://github.com/Samayel/Lists.jl")
    Pkg.build("Lists")
end

try
    eval(Expr(:import,:Multicombinations))
catch err
    Pkg.clone("https://github.com/jlep/Multicombinations.jl")
    Pkg.build("Multicombinations")
end

file = joinpath(Pkg.dir("Nemo"), "src", "flint", "fmpq_rel_series.jl")
content = readstring(file)
content = replace(content, "fmpq_poly_scalar_divexact_", "fmpq_poly_scalar_div_")
open(file, "w") do f
    write(f, content)
end