text:
module TestAlign

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Bio
using Bio.Seq
using Bio.Align
using TestFunctions

# 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

@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(seqpos, alnseq) == (refpos, OP_MATCH)
            @test ref2seq(refpos, alnseq) == (seqpos, OP_MATCH)
        end
        # OP_INSERT
        @test seq2ref(10, alnseq) == (17, OP_INSERT)
        @test seq2ref(11, alnseq) == (17, OP_INSERT)
        # OP_DELETE
        @test ref2seq( 9, alnseq) == ( 4, OP_DELETE)
        @test ref2seq(10, alnseq) == ( 4, OP_DELETE)
        @test ref2seq(23, alnseq) == (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 bytestring(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 "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

end # TestAlign


text:
Test.@test_throws ErrorException PredictMD.generate_examples(PredictMD.maketempdir())
Test.@test_throws ErrorException PredictMD.generate_examples(PredictMD.maketempdir(); scripts = true)

ext: jl


text:
module VortexHelperBowlPuffer
using ..Ahorn, Maple

@mapdef Entity "VortexHelper/BowlPuffer" BowlPuffer(x::Integer, y::Integer, noRespawn::Bool = false, explodeTimer::Number = 1.0)

const placements = Ahorn.PlacementDict(
    "Pufferfish Bowl (Vortex Helper)" => Ahorn.EntityPlacement(
        BowlPuffer,
		"point"
    ),
	"Pufferfish Bowl (No Respawn) (Vortex Helper)" => Ahorn.EntityPlacement(
        BowlPuffer,
		"point",
		Dict{String, Any}(
			"noRespawn" => true
		)
    )
)

pufferBowlSprite = "objects/VortexHelper/pufferBowl/idle00"

function Ahorn.selection(entity::BowlPuffer)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 11, y - 11, 21, 19)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BowlPuffer, room::Maple.Room)
    Ahorn.drawSprite(ctx, pufferBowlSprite, 0, -3)
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
"""
    struct LagrangianRefFE{D} <: ReferenceFE{D}
      data::GenericRefFE{D}
      face_own_nodeids::Vector{Vector{Int}}
      own_nodes_permutations::Vector{Vector{Int}},
    end

Type representing a Lagrangian finite element. In addition to all the information
provided by a `ReferenceFE`, this type also provides "node-based" information in the following
fields:

- `face_own_nodeids::Vector{Vector{Int}}`: nodes owned by each face
- `own_nodes_permutations::Vector{Vector{Int}}`: permutations of the nodes when the vertices are permuted

For this type

-  `get_dofs(reffe)` returns a `LagrangianDofBasis`
-  `get_prebasis(reffe)` returns a `MonomialBasis`
-  `ReferenceFE{N}(reffe,faceid) where N` returns a `LagrangianRefFE{N}`
"""
struct LagrangianRefFE{D} <: ReferenceFE{D}
  data::GenericRefFE{D}
  face_own_nodeids::Vector{Vector{Int}}
  own_nodes_permutations::Vector{Vector{Int}}
  @doc """
      LagrangianRefFE(
        polytope::Polytope{D},
        prebasis::MonomialBasis,
        dofs::LagrangianDofBasis,
        face_own_nodeids::Vector{Vector{Int}},
        own_nodes_permutations::Vector{Vector{Int}},
        reffaces...) where D

  Low level (inner) constructor of `LagrangianRefFE`.
  """
  function LagrangianRefFE(
    polytope::Polytope{D},
    prebasis::MonomialBasis,
    dofs::LagrangianDofBasis,
    face_own_nodeids::Vector{Vector{Int}},
    own_nodes_permutations::Vector{Vector{Int}},
    reffaces...) where D

    face_own_dofids = _generate_nfacedofs(face_own_nodeids,dofs.node_and_comp_to_dof)
    own_dofs_permutations = _find_own_dof_permutaions(own_nodes_permutations,dofs.node_and_comp_to_dof,face_own_nodeids,face_own_dofids)

    data = GenericRefFE(
      polytope,prebasis,dofs,face_own_dofids;
      own_dofs_permutations = own_dofs_permutations,
      reffaces = reffaces)

    new{D}(data,face_own_nodeids,own_nodes_permutations)
  end
end

num_dofs(reffe::LagrangianRefFE) = reffe.data.ndofs

get_polytope(reffe::LagrangianRefFE) = reffe.data.polytope

get_prebasis(reffe::LagrangianRefFE) = reffe.data.prebasis

get_dofs(reffe::LagrangianRefFE) = reffe.data.dofs

get_face_own_dofids(reffe::LagrangianRefFE) = reffe.data.face_own_dofids

get_own_dofs_permutations(reffe::LagrangianRefFE) = reffe.data.own_dofs_permutations

get_shapefuns(reffe::LagrangianRefFE) = reffe.data.shapefuns

"""
    get_node_coordinates(reffe::LagrangianRefFE) -> Vector{Point{D,Float64}}

Returns the nodal coordinates of the underlying `LagrangianDofBasis`.
"""
get_node_coordinates(reffe::LagrangianRefFE) = reffe.data.dofs.nodes

"""
    get_dof_to_node(reffe::LagrangianRefFE) -> Vector{Int}

Returns the field `dof_to_node` of the underlying `LagrangianDofBasis`.
"""
get_dof_to_node(reffe::LagrangianRefFE) = reffe.data.dofs.dof_to_node

"""
    get_dof_to_comp(reffe::LagrangianRefFE) -> Vector{Int}

Returns the field `dof_to_comp` of the underlying `LagrangianDofBasis`.
"""
get_dof_to_comp(reffe::LagrangianRefFE) = reffe.data.dofs.dof_to_comp

"""
    get_node_and_comp_to_dof(reffe::LagrangianRefFE) -> Vector

Returns the field `node_and_comp_to_dof` of the underlying `LagrangianDofBasis`.
"""
get_node_and_comp_to_dof(reffe::LagrangianRefFE) = reffe.data.dofs.node_and_comp_to_dof

function ReferenceFE{N}(reffe::LagrangianRefFE,iface::Integer) where N
  ReferenceFE{N}(reffe.data,iface)
end

function ReferenceFE{D}(reffe::LagrangianRefFE{D},iface::Integer) where D
  @assert iface==1 "Only one D-face"
  reffe
end

# Helpers for LagrangianRefFE

function _generate_nfacedofs(nfacenodes,node_and_comp_to_dof)
  faces = 1:length(nfacenodes)
  T = eltype(node_and_comp_to_dof)
  comps = 1:n_components(T)
  nfacedofs = [Int[] for i in faces]
  for face in faces
    nodes = nfacenodes[face]
    # Node major
    for comp in comps
      for node in nodes
        comp_to_dofs = node_and_comp_to_dof[node]
        dof = comp_to_dofs[comp]
        push!(nfacedofs[face],dof)
      end
    end
  end
  nfacedofs
end

function _find_own_dof_permutaions(node_perms,node_and_comp_to_dof,nfacenodeids,nfacedofsids)
  dof_perms = Vector{Int}[]
  T = eltype(node_and_comp_to_dof)
  ncomps = n_components(T)
  idof_to_dof = nfacedofsids[end]
  inode_to_node = nfacenodeids[end]
  for inode_to_pinode in node_perms
    ninodes = length(inode_to_pinode)
    nidofs = ncomps*ninodes
    idof_to_pidof = fill(INVALID_PERM,nidofs)
    for (inode,ipnode) in enumerate(inode_to_pinode)
      if ipnode == INVALID_PERM
        continue
      end
      node = inode_to_node[inode]
      pnode = inode_to_node[ipnode]
      comp_to_pdof = node_and_comp_to_dof[pnode]
      comp_to_dof = node_and_comp_to_dof[node]
      for comp in 1:ncomps
        dof = comp_to_dof[comp]
        pdof = comp_to_pdof[comp]
        idof = findfirst(i->i==dof,idof_to_dof)
        ipdof = findfirst(i->i==pdof,idof_to_dof)
        idof_to_pidof[idof] = ipdof
      end
    end
    push!(dof_perms,idof_to_pidof)
  end
  dof_perms
end

# Construction of LagrangianRefFE from Polytopes

"""
    LagrangianRefFE(::Type{T},p::Polytope,orders) where T
    LagrangianRefFE(::Type{T},p::Polytope,order::Int) where T

Builds a `LagrangianRefFE` object on top of the given polytope. `T` is the type of
the value of the approximation space (e.g., `T=Float64` for scalar-valued problems,
`T=VectorValue{N,Float64}` for vector-valued problems with `N` components). The arguments `order` or `orders`
are for the polynomial order of the resulting space, which allows isotropic or anisotropic orders respectively
(provided that the cell topology allows the given anisotropic order). The argument `orders` should be an
indexable collection of `D` integers (e.g., a tuple or a vector), being `D` the number of space dimensions.

In order to be able to use this function, the type of the provided polytope `p` has to implement the
following additional methods. They have been implemented for `ExtrusionPolytope` in the library. They
need to be implemented for new polytope types in order to build Lagangian reference elements on top of them.

- [`compute_monomial_basis(::Type{T},p::Polytope,orders) where T`](@ref)
- [`compute_own_nodes(p::Polytope,orders)`](@ref)
- [`compute_face_orders(p::Polytope,face::Polytope,iface::Int,orders)`](@ref)

The following methods are also used in the construction of the `LagrangianRefFE` object. A default implementation
of them is available in terms of the three previous methods. However, the user can also implement them for
new polytope types increasing customization possibilities.

- [`compute_nodes(p::Polytope,orders)`](@ref)
- [`compute_own_nodes_permutations(p::Polytope, interior_nodes)`](@ref)
- [`compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T`](@ref)

"""
function LagrangianRefFE(::Type{T},p::Polytope{D},orders) where {T,D}
  prebasis = compute_monomial_basis(T,p,orders)
  nodes, face_own_nodeids = compute_nodes(p,orders)
  dofs = LagrangianDofBasis(T,nodes)
  interior_nodes = dofs.nodes[face_own_nodeids[end]]
  own_nodes_permutations = compute_own_nodes_permutations(p, interior_nodes)
  reffaces = compute_lagrangian_reffaces(T,p,orders)
  LagrangianRefFE(p,prebasis,dofs,face_own_nodeids,own_nodes_permutations,reffaces...)
end

function MonomialBasis(::Type{T},p::Polytope,orders) where T
  compute_monomial_basis(T,p,orders)
end

function LagrangianDofBasis(::Type{T},p::Polytope,orders) where T
  nodes, _ = compute_nodes(p,orders)
  LagrangianDofBasis(T,nodes)
end

# Constructors taking Int

function LagrangianRefFE(::Type{T},p::Polytope{D},order::Int) where {T,D}
  orders = tfill(order,Val{D}())
  LagrangianRefFE(T,p,orders)
end

function MonomialBasis(::Type{T},p::Polytope{D},order::Int) where {D,T}
  orders = tfill(order,Val{D}())
  MonomialBasis(T,p,orders)
end

function LagrangianDofBasis(::Type{T},p::Polytope{D},order::Int) where {T,D}
  orders = tfill(order,Val{D}())
  LagrangianDofBasis(T,p,orders)
end

# Queries needed to be implemented for polytopes in order to use them
# for building LagrangianRefFEs in a seamless way

"""
    compute_monomial_basis(::Type{T},p::Polytope,orders) where T -> MonomialBasis

Returns the monomial basis of value type `T` and order per direction described by `orders`
on top of the polytope `p`.
"""
function compute_monomial_basis(::Type{T},p::Polytope,orders) where T
  @abstractmethod
end

"""
    compute_own_nodes(p::Polytope{D},orders) where D -> Vector{Point{D,Float64}}

Returns the coordinates of the nodes owned by the interior of the polytope
associated with a Lagrangian space with the order per direction described by `orders`.
"""
function compute_own_nodes(p::Polytope,orders)
  @abstractmethod
end

"""
    compute_face_orders(p::Polytope,face::Polytope,iface::Int,orders)

Returns a vector or a tuple with the order per direction at the face `face`
of the polytope `p` when restricting the order per direction `orders` to this face.
`iface` is the face id of `face` in the numeration restricted to the face dimension.
"""
function compute_face_orders(p::Polytope,face::Polytope,iface::Int,orders)
  @abstractmethod
end

"""
    compute_nodes(p::Polytope,orders)

When called

    node_coords, face_own_nodeids = compute_nodes(p,orders)

Returns `node_coords`, the nodal coordinates of all the Lagrangian nodes associated with the order per direction
`orders`, and `face_own_nodeids`, being a vector of vectors indicating which nodes are owned by each of
the faces of the polytope `p`.
"""
function compute_nodes(p::Polytope,orders)
  _compute_nodes(p,orders)
end

"""
    compute_own_nodes_permutations(
      p::Polytope, own_nodes_coordinates) -> Vector{Vector{Int}}

Returns a vector of vectors with the permutations of the nodes owned by the interior of the
polytope.
"""
function compute_own_nodes_permutations(p::Polytope, interior_nodes)
  _compute_node_permutations(p, interior_nodes)
end

"""
    compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T

Returns a tuple of length `D` being the number of space dimensions.
The entry `d+1` of this tuple contains a vector of `LagrangianRefFE`
one for each face of dimension `d` on the boundary of the polytope.
"""
function compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T
  _compute_lagrangian_reffaces(T,p,orders)
end

# Default implementations

function _compute_nodes(p,orders)
  if any( orders .== 0)
    _compute_constant_nodes(p,orders)
  elseif all(orders .== 1)
    _compute_linear_nodes(p)
  else
    _compute_high_order_nodes(p,orders)
  end
end

function _compute_constant_nodes(p,orders)
  @assert all( orders .== 0) "If an order is 0 in some direction, it should be 0 also in the others"
  x = compute_own_nodes(p,orders)
  facenodes = [Int[] for i in 1:num_faces(p)]
  push!(facenodes[end],1)
  x, facenodes
end

function _compute_linear_nodes(p)
  x = get_vertex_coordinates(p)
  facenodes = [Int[] for i in 1:num_faces(p)]
  for i in 1:num_vertices(p)
    push!(facenodes[i],i)
  end
  x, facenodes
end

function _compute_high_order_nodes(p::Polytope{D},orders) where D
  nodes = Point{D,Float64}[]
  facenodes = [Int[] for i in 1:num_faces(p)]
  _compute_high_order_nodes_dim_0!(nodes,facenodes,p)
  for d in 1:(num_dims(p)-1)
    _compute_high_order_nodes_dim_d!(nodes,facenodes,p,orders,Val{d}())
  end
  _compute_high_order_nodes_dim_D!(nodes,facenodes,p,orders)
  (nodes, facenodes)
end

function _compute_high_order_nodes_dim_0!(nodes,facenodes,p)
  x = get_vertex_coordinates(p)
  k = 1
  for vertex in 1:num_vertices(p)
    push!(nodes,x[vertex])
    push!(facenodes[vertex],k)
    k += 1
  end
end

@noinline function _compute_high_order_nodes_dim_d!(nodes,facenodes,p,orders,::Val{d}) where d
  x = get_vertex_coordinates(p)
  offset = get_offset(p,d)
  k = length(nodes)+1
  for iface in 1:num_faces(p,d)
    face = Polytope{d}(p,iface)
    face_ref_x = get_vertex_coordinates(face)
    face_prebasis = MonomialBasis(Float64,face,1)
    change = inv(evaluate(face_prebasis,face_ref_x))
    face_shapefuns = change_basis(face_prebasis,change)
    face_vertex_ids = get_faces(p,d,0)[iface]
    face_x = x[face_vertex_ids]
    face_orders = compute_face_orders(p,face,iface,orders)
    face_interior_nodes = compute_own_nodes(face,face_orders)
    face_high_x = evaluate(face_shapefuns,face_interior_nodes)*face_x
    for xi in 1:length(face_high_x)
      push!(nodes,face_high_x[xi])
      push!(facenodes[iface+offset],k)
      k += 1
    end
  end
end

function _compute_high_order_nodes_dim_D!(nodes,facenodes,p,orders)
  k = length(nodes)+1
  p_high_x = compute_own_nodes(p,orders)
  for xi in 1:length(p_high_x)
    push!(nodes,p_high_x[xi])
    push!(facenodes[end],k)
    k += 1
  end
end

_compute_node_permutations(::Polytope{0}, interior_nodes) = [[1]]

function _compute_node_permutations(p, interior_nodes)
  vertex_to_coord = get_vertex_coordinates(p)
  lbasis = MonomialBasis(Float64,p,1)
  change = inv(evaluate(lbasis,vertex_to_coord))
  lshapefuns = change_basis(lbasis,change)
  perms = get_vertex_permutations(p)
  map = evaluate(lshapefuns,interior_nodes)
  pvertex_to_coord = similar(vertex_to_coord)
  node_perms = Vector{Int}[]
  tol = 1.0e-10
  for vertex_to_pvertex in perms
    node_to_pnode = fill(INVALID_PERM,length(interior_nodes))
    pvertex_to_coord[vertex_to_pvertex] = vertex_to_coord
    pinterior_nodes = map*pvertex_to_coord
    for node in 1:length(interior_nodes)
      x = interior_nodes[node]
      pnode = findfirst(i->norm(i-x)<tol,pinterior_nodes)
      if pnode != nothing
         node_to_pnode[node] = pnode
      end
    end
    push!(node_perms,node_to_pnode)
  end
  node_perms
end

_compute_lagrangian_reffaces(::Type{T},p::Polytope{0},orders) where T = ()

function _compute_lagrangian_reffaces(::Type{T},p::Polytope{D},orders) where {T,D}
  reffaces = [ LagrangianRefFE{d}[]  for d in 0:D ]
  p0 = Polytope{0}(p,1)
  reffe0 = LagrangianRefFE(T,p0,())
  for vertex in 1:num_vertices(p)
    push!(reffaces[0+1],reffe0)
  end
  offsets = get_offsets(p)
  for d in 1:(num_dims(p)-1)
    offset = offsets[d+1]
    for iface in 1:num_faces(p,d)
      face = Polytope{d}(p,iface)
      face_orders = compute_face_orders(p,face,iface,orders)
      refface = LagrangianRefFE(T,face,face_orders)
      push!(reffaces[d+1],refface)
    end
  end
  tuple(reffaces...)
end

# Particular implementation for ExtrusionPolytope

function compute_monomial_basis(::Type{T},p::ExtrusionPolytope{D},orders) where {D,T}
  extrusion = Tuple(p.extrusion.array)
  terms = _monomial_terms(extrusion,orders)
  MonomialBasis{D}(T,orders,terms)
end

function compute_own_nodes(p::ExtrusionPolytope{D},orders) where D
  extrusion = Tuple(p.extrusion.array)
  if all(orders .== 0)
    _interior_nodes_order_0(p)
  else
    _interior_nodes(extrusion,orders)
  end
end

function _interior_nodes_order_0(p)
  x = get_vertex_coordinates(p)
  x0 = sum(x) / length(x)
  [x0,]
end

function compute_face_orders(p::ExtrusionPolytope,face::ExtrusionPolytope{D},iface::Int,orders) where D
  d = num_dims(face)
  offset = get_offset(p,d)
  nface = p.nfaces[iface+offset]
  face_orders = _extract_nonzeros(nface.extrusion,orders)
  face_orders
end

function compute_nodes(p::ExtrusionPolytope{D},orders) where D
  _nodes, facenodes = _compute_nodes(p,orders)
  if any( orders .== 0)
    return (_nodes, facenodes)
  end
  terms = _coords_to_terms(_nodes,orders)
  nodes = _terms_to_coords(terms,orders)
  (nodes, facenodes)
end

# Helpers for the ExtrusionPolytope-related implementation

function _monomial_terms(extrusion::NTuple{D,Int},orders) where D
  terms = CartesianIndex{D}[]
  if D == 0
    push!(terms,CartesianIndex(()))
    return terms
  end
  _check_orders(extrusion,orders)
  M = mutable(VectorValue{D,Int})
  term = zero(M)
  _orders = M(orders)
  k = 0
  _add_terms!(terms,term,extrusion,_orders,D,k)
  terms
end

function _interior_nodes(extrusion::NTuple{D,Int},orders) where D
  _check_orders(extrusion,orders)
  terms = CartesianIndex{D}[]
  M = mutable(VectorValue{D,Int})
  term = zero(M)
  _orders = M(orders)
  k = 1
  _add_terms!(terms,term,extrusion,_orders,D,k)
  _terms_to_coords(terms,orders)
end

function _check_orders(extrusion,orders)
  D = length(extrusion)
  @assert length(orders) == D "container of orders not long enough"
  _orders = collect(orders)
  if extrusion[D] == HEX_AXIS
    _orders[D] = 0
  end
  for d in (D-1):-1:1
    if (extrusion[d] == HEX_AXIS || d == 1) && _orders[d+1] == 0
      _orders[d] = 0
    end
  end
  nz = _orders[_orders .!= 0]
  if length(nz) > 1
    @assert all(nz .== nz[1]) "The provided anisotropic order is not compatible with polytope topology"
  end
  nothing
end

function _add_terms!(terms,term,extrusion,orders,dim,k)
  _term = copy(term)
  _orders = copy(orders)
  indexbase = 1
  for i in k:(_orders[dim]-k)
    _term[dim] = i + indexbase
    if dim > 1
      if (extrusion[dim] == TET_AXIS) && i != 0
        _orders .-= 1
      end
      _add_terms!(terms,_term,extrusion,_orders,dim-1,k)
    else
      push!(terms,CartesianIndex(Tuple(_term)))
    end
  end
end

function _coords_to_terms(coords::Vector{<:Point{D}},orders) where D
  indexbase = 1
  terms = CartesianIndex{D}[]
  P = Point{D,Int}
  t = zero(mutable(P))
  for x in coords
    for d in 1:D
      t[d] = round(x[d]*orders[d]) + indexbase
    end
    term = CartesianIndex(Tuple(t))
    push!(terms,term)
  end
  terms
end

function  _terms_to_coords(terms::Vector{CartesianIndex{D}},orders) where D
  P = Point{D,Float64}
  indexbase = 1
  nodes = P[]
  x = zero(mutable(P))
  for t in terms
    for d in 1:D
      x[d] = (t[d] - indexbase) / orders[d]
    end
    node = P(x)
    push!(nodes,node)
  end
  nodes
end

function _extract_nonzeros(mask,values)
  b = Int[]
  for (m,n) in zip(mask,values)
    if (m != 0)
      push!(b, n)
    end
  end
  return Tuple(b)
end


text:
using WebSockets

function todict(rc::Tuple{Request, WebSocket})
  req, client = rc
  req′ = todict(req)
  req′[:socket] = client
  return req′
end

function wcatch(app, req)
  try
    app(req)
  catch e
    println(STDERR, "Error handling websocket connection:")
    showerror(STDERR, e, catch_backtrace())
  end
end

function wclose(req)
  close(req[:socket])
end

function echo(req)
  sock = req[:socket]
  while isopen(sock)
    try
      write(sock, read(sock))
    end
  end
end


text:
#
# This module contains predefined colors for elements in a world like energy and
# nothing. We need this for solving an issue with invisible energy. This is
# because GR library has only 1256 colors and not 24bit :( Only 56 slots are
# available for color redefining. Other 1200 colors are for organisms. It also
# contains dots related constants. Maximum supported amouny of colors is 4096.
#
# @author DeadbraiN
#
module Dots
  export INDEX_MAX_COLOR
  export INDEX_MAX_ORG_COLOR

  export COLOR_EMPTY
  export INDEX_EMPTY
  export COLOR_ENERGY
  export INDEX_ENERGY
  export COLOR_TEXT
  export INDEX_TEXT
  export INDEX_FIRST_COLOR
  export COLORS
  #
  # Dots directions
  #
  export DIRECTION_NO
  export DIRECTION_LEFT
  export DIRECTION_UP
  export DIRECTION_RIGHT
  export DIRECTION_DOWN
  #
  # Evergy move possible directions
  #
  const DIRECTION_NO        = 0
  const DIRECTION_LEFT      = 1
  const DIRECTION_UP        = 2
  const DIRECTION_RIGHT     = 3
  const DIRECTION_DOWN      = 4
  #
  # Describes one color and it's index
  #
  type Color
    index::UInt16
    rgb::Array{Float64, 1}
  end
  #
  # Maximum index of color we may use
  #
  const INDEX_MAX_ORG_COLOR = UInt16(1200)
  const INDEX_MAX_COLOR     = UInt16(1256)
  #
  # Color and index of empty space
  #
  const COLOR_EMPTY   = Float64[0.0, 0.0, 0.0]
  const INDEX_EMPTY   = UInt16(0)
  #
  # Color and index of energy block
  #
  const COLOR_ENERGY  = Float64[0.0, 1.0, 0.0]
  const INDEX_ENERGY  = UInt16(1201)
  #
  # Text coloring
  #
  const COLOR_TEXT    = Float64[1.0, 1.0, 1.0]
  const INDEX_TEXT    = UInt16(1202)
  #
  # Index of first color
  #
  const INDEX_FIRST_COLOR = UInt16(1)
  #
  # All available colors from above
  #
  const COLORS = Color[
    Color(INDEX_ENERGY, COLOR_ENERGY),
    Color(INDEX_EMPTY,  COLOR_EMPTY),
    Color(INDEX_TEXT,   COLOR_TEXT)
  ]
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
struct Add{T} <: PushBinaryOp{T} end
execute_op(i::Add{T}, a::T, b::T) where {T} = b+a
oldopname(::Type{<:Add}) = "+"

struct Sub{T} <: PushBinaryOp{T} end
execute_op(i::Sub{T}, a::T, b::T) where {T} = b-a # Push order is that element at top of stack (a) is subtracted from the next element (b)
oldopname(::Type{<:Sub}) = "-"

struct Mul{T} <: PushBinaryOp{T} end
execute_op(i::Mul{T}, a::T, b::T) where {T} = b*a
oldopname(::Type{<:Mul}) = "*"

struct Div{T} <: PushBinaryOp{T} end
@inline function execute!(i::Div{T}, interp::AbstractPushInterpreter) where {T}
    stackhas(interp, T, 2) || return(nothing)
    last(stack(interp, T)) == 0 && return(nothing)
    a, b = pop2!(interp, T)
    res = b/a
    push!(interp, res)
end
oldopname(::Type{<:Div}) = "/"

struct Mod{T} <: PushBinaryOp{T} end
@inline function execute!(i::Mod{T}, interp::AbstractPushInterpreter) where {T}
    stackhas(interp, T, 2) || return(nothing)
    last(stack(interp, T)) == 0 && return(nothing)
    a, b = pop2!(interp, T)
    res = div(b, a)
    push!(interp, res)
end
oldopname(::Type{<:Mod}) = "%"

struct Max{T} <: PushBinaryOp{T} end
execute_op(i::Max{T}, a::T, b::T) where {T} = max(a, b)

struct Min{T} <: PushBinaryOp{T} end
execute_op(i::Min{T}, a::T, b::T) where {T} = min(a, b)


text:
struct UnionFind
    par::Vector{Int}
    size::Vector{Int}
    UnionFind(n) = new(collect(1:n), ones(Int, n))
end

root!(uf, x) = uf.par[x] == x ? x : uf.par[x] = root!(uf, uf.par[x])

same(uf, x, y) = root!(uf, x) == root!(uf, y)

function unite!(uf, x, y)
    x = root!(uf, x)
    y = root!(uf, y)
    x == y && return false

    if uf.size[x] > uf.size[y]
        uf.par[y] = x
        uf.size[x] += uf.size[y]
    else
        uf.par[x] = y
        uf.size[y] += uf.size[x]
    end
    return true
end
size!(uf, x) = uf.size[root!(uf, x)]

function main()
    N, K, L = parse.(Int, split(readline()))

    uf₁ = UnionFind(N)
    uf₂ = UnionFind(N)

    for _ in 1:K
        p, q = parse.(Int, split(readline()))
        unite!(uf₁, p, q)
    end

    for _ in 1:L
        r, s = parse.(Int, split(readline()))
        unite!(uf₂, r, s)
    end

    d = Dict{Tuple{Int, Int}, Int}()
    for i in 1:N
        key = (root!(uf₁, i), root!(uf₂, i))
        if haskey(d, key)
            d[key] += 1
        else
            d[key] = 1
        end
    end

    for i in 1:N-1
        print(d[(root!(uf₁, i), root!(uf₂, i))], " ")
    end

    println(d[(root!(uf₁, N), root!(uf₂, N))])
end

main()


text:
using QuantEcon

include(string(pwd(),"/hmc/diffs_estrut.jl"))
include(string(pwd(),"/hmc/diff_kalman.jl"))
include(string(pwd(),"/src/simulation.jl"))

aa = [1/3,0.99,6,2/3,1,1,1,1.5,0.5/4,0.5]

teste_diff_mod = diff_mod(aa,1)

yy,shocks = simulate_dsge(GAMMA_0_foo(aa),GAMMA_1_foo(aa),PSI_foo(aa),PI_foo(aa),500)

data = yy[:,1]

sol = gensys(GAMMA_0_foo(aa),GAMMA_1_foo(aa),PSI_foo(aa),PI_foo(aa); verbose = false)

G = zeros(1,4)
G[1,2] = 1

A = sol.Theta1
R = [0] .+ 1e-8
Q = aa[6]^2*sol.Theta2*sol.Theta2'

dA = teste_diff_mod[:,1:16]
dG = teste_diff_mod[:,17:20]
dQ = teste_diff_mod[:,21:30]
dR = zeros(size(aa,1),1)

dx0 = zeros(4)
S0 = zeros(4,4)

g_x_hat_l = zeros(size(aa,1),4)
d_Sl = zeros(size(aa,1),size(S0,1)^2)

kalman_res = Kalman(A,G,Q,R) #create a Kalman filter instance

y_mean = mean(data,dims=1)
y_var = var(data,dims=1)
y_var = reshape(y_var,size(y_var,2))
#y_var = repeat([1],p)

x_hat = repeat(y_mean,4) #initial mean of the state
x_var = diagm(repeat(y_var,4))#variance initial of state
#x_var = x_var*x_var'
set_state!(kalman_res,x_hat,x_var)

med = kalman_res.cur_x_hat
varian = kalman_res.cur_sigma
#println(det(varian))
eta = data[2,:] - kalman_res.G*med #mean loglike
P = kalman_res.G*varian*kalman_res.G' + kalman_res.R
teste_cond = 1/cond(P)

gain = (A*varian*G')*inv(P)

d_p = diff_P(G,varian,dG,d_Sl,dR)
d_gain = diff_K(A,P,G,varian,dA,dG,d_p,d_Sl)
g_eta = grad_y(med,g_x_hat_l,G,dG)
g_x_hat = grad_x_hat(med,dA,A,g_x_hat_l,eta,g_eta, gain,d_gain)

p_inv = inv(P)

-1/2*vec(p_inv)'*d_p' - eta'*p_inv*g_eta' -1/2*(kron(eta'*p_inv, eta*p_inv)*d_p') #updating the kalman estimates

m = size(P,1)
n = size(A,1)
Knn = commutation_matrix(n,n)

b1 = (kron(A*varian,I(n)) + kron(I(n),A*varian)*Knn)*dA'
b2 = (kron(gain*P,I(n)) + kron(I(n),gain*P)*Kmn)*d_gain'
b1 + kron(A,A)*d_Sl' -b2 - kron(gain,gain)*d_p' + duplication_matrix(4)*dQ'

d_Sl = diff_S(A,P,gain,varian,dA,d_p,d_gain,d_Sl,dQ)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using ReservoirComputing

#model parameters
res_size = 30
radius = 1.2
activation = tanh
degree = 6
sparsity = 0.5
sigma = 0.1
beta = 0.0
alpha = 1.0
nla_type = NLADefault()
in_size = 3
extended_states = false
h_steps = 2

W = init_reservoir_givendeg(res_size, radius, degree)

train_len = 50
predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#test input layers functions
input_layer = [init_input_layer, init_dense_input_layer, init_sparse_input_layer, min_complex_input, irrational_sign_input]

for t in input_layer

    if t == init_sparse_input_layer
        W_in = t(res_size, in_size, sigma, sparsity)
    else
        W_in = t(res_size, in_size, sigma)
    end
    esn = ESN(W,
        train,
        W_in,
        activation = activation,
        alpha = alpha,
        nla_type = nla_type,
        extended_states = extended_states)

    @test size(W_in, 1) == res_size
    @test size(W_in, 2) == in_size
end


text:
module NumericalFluxes

export Rusanov, CentralGradPenalty, CentralNumericalFluxDiffusive,
       CentralNumericalFluxNonDiffusive

using StaticArrays
using GPUifyLoops: @unroll
import ..DGmethods: BalanceLaw, Grad, Vars, vars_state, vars_diffusive,
                    vars_aux, vars_gradient, boundary_state!, wavespeed,
                    flux_nondiffusive!, flux_diffusive!, diffusive!, num_state,
                    num_gradient, gradvariables!

"""
    GradNumericalPenalty

Any `P <: GradNumericalPenalty` should define methods for:

   diffusive_penalty!(gnf::P, bl::BalanceLaw, diffF, nM, QM, QdiffM, QauxM, QP,
                      QdiffP, QauxP, t)
   diffusive_boundary_penalty!(gnf::P, bl::BalanceLaw, l_Qvisc, nM, l_GM, l_QM,
                               l_auxM, l_GP, l_QP, l_auxP, bctype, t)

"""
abstract type GradNumericalPenalty end

function diffusive_penalty! end
function diffusive_boundary_penalty! end

"""
    CentralGradPenalty <: GradNumericalPenalty

"""
struct CentralGradPenalty <: GradNumericalPenalty end

function diffusive_penalty!(::CentralGradPenalty, bl::BalanceLaw,
                            VF, nM, diffM, QM, aM, diffP, QP, aP, t)
  FT = eltype(QM)

  @inbounds begin
    ndim = 3
    ngradstate = num_gradient(bl,FT)
    n_Δdiff = similar(VF, Size(ndim, ngradstate))
    @unroll for j = 1:ngradstate
      @unroll for i = 1:ndim
        n_Δdiff[i, j] = nM[i] * (diffP[j] - diffM[j]) / 2
      end
    end
    diffusive!(bl, Vars{vars_diffusive(bl,FT)}(VF),
               Grad{vars_gradient(bl,FT)}(n_Δdiff),
               Vars{vars_state(bl,FT)}(QM), Vars{vars_aux(bl,FT)}(aM),
               t)
  end
end

function diffusive_boundary_penalty!(nf::CentralGradPenalty, bl::BalanceLaw,
                                     VF, nM, diffM, QM, aM, diffP, QP, aP,
                                     bctype, t, Q1, aux1)
  FT = eltype(diffP)
  boundary_state!(nf, bl, Vars{vars_state(bl,FT)}(QP),
                  Vars{vars_aux(bl,FT)}(aP), nM,
                  Vars{vars_state(bl,FT)}(QM),
                  Vars{vars_aux(bl,FT)}(aM), bctype, t,
                  Vars{vars_state(bl,FT)}(Q1),
                  Vars{vars_aux(bl,FT)}(aux1))

  gradvariables!(bl, Vars{vars_gradient(bl,FT)}(diffP),
                 Vars{vars_state(bl,FT)}(QP),
                 Vars{vars_aux(bl,FT)}(aP), t)

  diffusive_penalty!(nf, bl, VF, nM, diffM, QM, aM, diffP, QP, aP, t)
end

"""
    NumericalFluxNonDiffusive

Any `N <: NumericalFluxNonDiffusive` should define the a method for

    numerical_flux_nondiffusive!(nf::N, bl::BalanceLaw, F, nM, QM, QauxM, QP,
                                 QauxP, t)

where
- `F` is the numerical flux array
- `nM` is the unit normal
- `QM`/`QP` are the minus/positive state arrays
- `t` is the time

An optional method can also be defined for

    numerical_boundary_flux_nondiffusive!(nf::N, bl::BalanceLaw, F, nM, QM,
                                          QauxM, QP, QauxP, bctype, t)

"""
abstract type NumericalFluxNonDiffusive end

function numerical_flux_nondiffusive! end

function numerical_boundary_flux_nondiffusive!(nf::NumericalFluxNonDiffusive,
                                               bl::BalanceLaw,
                                               F::MArray{Tuple{nstate}}, nM, QM,
                                               auxM, QP, auxP, bctype, t,
                                               Q1, aux1) where {nstate}
  FT = eltype(F)
  boundary_state!(nf, bl, Vars{vars_state(bl,FT)}(QP),
                  Vars{vars_aux(bl,FT)}(auxP), nM,
                  Vars{vars_state(bl,FT)}(QM),
                  Vars{vars_aux(bl,FT)}(auxM), bctype, t,
                  Vars{vars_state(bl,FT)}(Q1),
                  Vars{vars_aux(bl,FT)}(aux1))
  numerical_flux_nondiffusive!(nf, bl, F, nM, QM, auxM, QP, auxP, t)
end

"""
    Rusanov <: NumericalFluxNonDiffusive

The Rusanov (aka local Lax-Friedrichs) numerical flux.

# Usage

    Rusanov()

Requires a `flux_nondiffusive!` and `wavespeed` method for the balance law.
"""
struct Rusanov <: NumericalFluxNonDiffusive end

function numerical_flux_nondiffusive!(::Rusanov, bl::BalanceLaw, F::MArray, nM,
                                      QM, auxM, QP, auxP, t)
  FT = eltype(F)
  nstate = num_state(bl,FT)

  λM = wavespeed(bl, nM, Vars{vars_state(bl,FT)}(QM),
                 Vars{vars_aux(bl,FT)}(auxM), t)
  FM = similar(F, Size(3, nstate))
  fill!(FM, -zero(eltype(FM)))
  flux_nondiffusive!(bl, Grad{vars_state(bl,FT)}(FM),
                     Vars{vars_state(bl,FT)}(QM),
                     Vars{vars_aux(bl,FT)}(auxM), t)

  λP = wavespeed(bl, nM, Vars{vars_state(bl,FT)}(QP),
                 Vars{vars_aux(bl,FT)}(auxP), t)
  FP = similar(F, Size(3, nstate))
  fill!(FP, -zero(eltype(FP)))
  flux_nondiffusive!(bl, Grad{vars_state(bl,FT)}(FP),
                     Vars{vars_state(bl,FT)}(QP),
                     Vars{vars_aux(bl,FT)}(auxP), t)

  λ  =  max(λM, λP)

  @unroll for s = 1:nstate
    @inbounds F[s] += (nM[1] * (FM[1, s] + FP[1, s]) + nM[2] * (FM[2, s] + FP[2, s]) +
                       nM[3] * (FM[3, s] + FP[3, s]) + λ * (QM[s] - QP[s])) / 2
  end
end

"""
    CentralNumericalFluxNonDiffusive() <: NumericalFluxNonDiffusive

The central numerical flux for nondiffusive terms

# Usage

    CentralNumericalFluxNonDiffusive()

Requires a `flux_nondiffusive!` method for the balance law.
"""
struct CentralNumericalFluxNonDiffusive <: NumericalFluxNonDiffusive end

function numerical_flux_nondiffusive!(::CentralNumericalFluxNonDiffusive,
                                      bl::BalanceLaw, F::MArray, nM,
                                      QM, auxM, QP, auxP, t)
  FT = eltype(F)
  nstate = num_state(bl,FT)

  FM = similar(F, Size(3, nstate))
  fill!(FM, -zero(eltype(FM)))
  flux_nondiffusive!(bl, Grad{vars_state(bl,FT)}(FM),
                     Vars{vars_state(bl,FT)}(QM),
                     Vars{vars_aux(bl,FT)}(auxM), t)

  FP = similar(F, Size(3, nstate))
  fill!(FP, -zero(eltype(FP)))
  flux_nondiffusive!(bl, Grad{vars_state(bl,FT)}(FP),
                     Vars{vars_state(bl,FT)}(QP),
                     Vars{vars_aux(bl,FT)}(auxP), t)

  @unroll for s = 1:nstate
    @inbounds F[s] += (nM[1] * (FM[1, s] + FP[1, s]) + nM[2] * (FM[2, s] + FP[2, s]) +
                       nM[3] * (FM[3, s] + FP[3, s])) / 2
  end
end

"""
    NumericalFluxDiffusive

Any `N <: NumericalFluxDiffusive` should define the a method for

    numerical_flux_diffusive!(nf::N, bl::BalanceLaw, F, nM, QM, QdiffM, QauxM, QP,
                              QdiffP, QauxP, t)

where
- `F` is the numerical flux array
- `nM` is the unit normal
- `QM`/`QP` are the minus/positive state arrays
- `QdiffM`/`QdiffP` are the minus/positive diffusive state arrays
- `QdiffM`/`QdiffP` are the minus/positive auxiliary state arrays
- `t` is the time

An optional method can also be defined for

    numerical_boundary_flux_diffusive!(nf::N, bl::BalanceLaw, F, nM, QM, QdiffM,
                                       QauxM, QP, QdiffP, QauxP, bctype, t)

"""
abstract type NumericalFluxDiffusive end

function numerical_flux_diffusive! end

function numerical_boundary_flux_diffusive!(nf::NumericalFluxDiffusive,
                                            bl::BalanceLaw,
                                            F::MArray{Tuple{nstate}}, nM,
                                            QM, QVM, auxM, QP, QVP, auxP,
                                            bctype, t, Q1, QV1,
                                            aux1) where {nstate}
  FT = eltype(F)
  boundary_state!(nf, bl, Vars{vars_state(bl,FT)}(QP),
                  Vars{vars_diffusive(bl,FT)}(QVP),
                  Vars{vars_aux(bl,FT)}(auxP), nM,
                  Vars{vars_state(bl,FT)}(QM),
                  Vars{vars_diffusive(bl,FT)}(QVM),
                  Vars{vars_aux(bl,FT)}(auxM), bctype, t,
                  Vars{vars_state(bl,FT)}(Q1),
                  Vars{vars_diffusive(bl,FT)}(QV1),
                  Vars{vars_aux(bl,FT)}(aux1))
  numerical_flux_diffusive!(nf, bl, F, nM, QM, QVM, auxM, QP, QVP, auxP, t)
end

"""
    CentralNumericalFluxDiffusive <: NumericalFluxDiffusive

The central numerical flux for diffusive terms

# Usage

    CentralNumericalFluxDiffusive()

Requires a `flux_diffusive!` for the balance law.
"""
struct CentralNumericalFluxDiffusive <: NumericalFluxDiffusive end

function numerical_flux_diffusive!(::CentralNumericalFluxDiffusive,
                                   bl::BalanceLaw, F::MArray, nM,
                                   QM, QVM, auxM,
                                   QP, QVP, auxP,
                                   t)
  FT = eltype(F)
  nstate = num_state(bl,FT)

  FM = similar(F, Size(3, nstate))
  fill!(FM, -zero(eltype(FM)))
  flux_diffusive!(bl, Grad{vars_state(bl,FT)}(FM),
                  Vars{vars_state(bl,FT)}(QM),
                  Vars{vars_diffusive(bl,FT)}(QVM),
                  Vars{vars_aux(bl,FT)}(auxM), t)

  FP = similar(F, Size(3, nstate))
  fill!(FP, -zero(eltype(FP)))
  flux_diffusive!(bl, Grad{vars_state(bl,FT)}(FP),
                  Vars{vars_state(bl,FT)}(QP),
                  Vars{vars_diffusive(bl,FT)}(QVP),
                  Vars{vars_aux(bl,FT)}(auxP), t)

  @unroll for s = 1:nstate
    @inbounds F[s] += (nM[1] * (FM[1, s] + FP[1, s]) + nM[2] * (FM[2, s] + FP[2, s]) +
                       nM[3] * (FM[3, s] + FP[3, s])) / 2
  end
end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using GeckoScreen
using Test
using DataFrames
using CSV
using Query

testcsv(path::AbstractString; kwargs...) = CSV.read(path, DataFrame; kwargs...)

@testset "GeckoScreen.jl" begin
    # Write your tests here.

    gmap = readGeCKOlibrary()
    first(gmap,6)
    #@test size(readGeCKOlibrary, 1) == 34

end


text:
using Documenter, Literate, Aeroelasticity

# Pre-install matplotlib
import Plots; Plots.pyplot()

# Generate examples
include("generate.jl")

# Build Documentation
makedocs(;
    modules = [Aeroelasticity],
    pages = [
        "Home" => "index.md",
        "Getting Started" => joinpath("examples", "guide.md"),
        "Examples" => [
            joinpath("examples", "section-stability.md"),
            joinpath("examples", "section-simulation.md"),
        ],
        "Model Documentation" => [
            "Aerodynamic Models" => [
                "Steady Thin Airfoil Theory" => joinpath("models", "aerodynamics", "steady.md"),
                "Quasi-Steady Thin Airfoil Theory" => joinpath("models", "aerodynamics", "quasisteady.md"),
                "Wagner's Function" => joinpath("models", "aerodynamics", "wagner.md"),
                "Peters' Finite State" => joinpath("models", "aerodynamics", "peters.md"),
            ],
            "Structural Models" => [
                "Typical Section" => joinpath("models", "structures", "section.md"),
                "Geometrically Exact Beam Theory" => joinpath("models", "structures", "gxbeam.md"),
            ],
            "Coupled Models" => [
                "Steady Thin Airfoil Theory + Typical Section" => joinpath("models", "couplings", "steady-section.md"),
                "Quasi-Steady Thin Airfoil Theory + Typical Section" => joinpath("models", "couplings", "quasisteady-section.md"),
                "Wagner's Function + Typical Section" => joinpath("models", "couplings", "wagner-section.md"),
                "Peters' Finite State + Typical Section" => joinpath("models", "couplings", "peters-section.md"),
            ],
        ],
        "Developer Guide" => "developer.md",
        "API Reference" => [
            "Public" => joinpath("reference", "public.md"),
            "Internals" => joinpath("reference", "internals.md"),
        ],
    ],
    sitename = "Aeroelasticity.jl",
    authors = "Taylor McDonnell <taylormcd@byu.edu>",
)

# Deploy Documentation
deploydocs(
    repo = "github.com/byuflowlab/Aeroelasticity.jl.git",
    devbranch = "main",
)


text:
# Copies figure 18 to figures directory:

run(`cp ../static/figure_water_period_TTV_trim.png  ..`)


text:
function nz_stream_foodweb()
  data_path = joinpath(@__DIR__, "../..", "data", "nz_stream")
  files = readdir(data_path)
  documents = filter(f -> endswith(f, ".csv"), files)
  Ns = UnipartiteNetwork{String}[]
  for doc in documents
    content = readdlm(joinpath(@__DIR__, "../..", "data", "nz_stream", doc), ',')
    species_names = convert(Array{String}, content[:,1][2:end])
    interaction_matrix = map(Int64, content[2:end,2:end]).>0
    push!(Ns, UnipartiteNetwork(interaction_matrix', species_names))
  end
  return Ns
end

function web_of_life(name)
  file_pathname = name * ".csv"
  data_path = joinpath(@__DIR__, "../..", "data", "weboflife")
  files = readdir(data_path)
  @assert file_pathname in files
  content = readdlm(joinpath(@__DIR__, "../..", "data", "weboflife", file_pathname), ',')
  bottom_species_names = convert(Array{String}, vec(content[:,1][2:end]))
  top_species_names = convert(Array{String}, vec(content[1,:][2:end]))
  int_mat = content[2:end,2:end]
  repl = findall(int_mat .== "#")
  int_mat[repl] .= 0.0
  interaction_matrix = map(Int64, round.(permutedims(int_mat)))
  ntype = BipartiteQuantitativeNetwork
  if maximum(interaction_matrix) == 1
    interaction_matrix = interaction_matrix .> 0
    ntype = BipartiteNetwork
  end
  this_net = ntype(interaction_matrix, top_species_names, bottom_species_names)
  forbidden_names = ["Numbers of flowers", "Num. of hosts sampled", "#", "Frequency of occurrences", "Number of flowers", "Number of droppings analysed"]
  sp_1 = filter(x -> !(x ∈ forbidden_names), species(this_net; dims=1))
  sp_2 = filter(x -> !(x ∈ forbidden_names), species(this_net; dims=2))
  return this_net[sp_1, sp_2]
end

function web_of_life()
  wol_ref = joinpath(@__DIR__, "../..", "data", "weboflife", "references.csv")
  wol_infos = readdlm(wol_ref, ',')
  names = [Symbol(replace(n, " " => "_")) for n in wol_infos[1,:]]
  infos = [NamedTuple{tuple(names...)}(tuple(wol_infos[i,:]...)) for i in 2:size(wol_infos,1)]
  return infos
end


text:
using CUDA
using Adapt
using CUDA.CUSPARSE
using SparseArrays

if CUSPARSE.version() >= v"11.4.1"
    @testset "generic mv!" for T in [Float32, Float64]
        A = sprand(T, 10, 10, 0.1)
        x = rand(Complex{T}, 10)
        y = similar(x)
        dx = adapt(CuArray, x)
        dy = adapt(CuArray, y)

        dA = adapt(CuArray, A)
        mv!('N', T(1.0), dA, dx, T(0.0), dy, 'O')
        @test Array(dy) ≈ A * x

        dA = CuSparseMatrixCSR(dA)
        mv!('N', T(1.0), dA, dx, T(0.0), dy, 'O')
        @test Array(dy) ≈ A * x
    end
end

if CUSPARSE.version() >= v"11.4.1" # lower CUDA version doesn't support these algorithms

    @testset "mm algo=$algo" for algo in [
        CUSPARSE.CUSPARSE_SPMM_ALG_DEFAULT,
        CUSPARSE.CUSPARSE_SPMM_CSR_ALG1,
        CUSPARSE.CUSPARSE_SPMM_CSR_ALG2,
        CUSPARSE.CUSPARSE_SPMM_CSR_ALG3,
    ]
        @testset "mm $T" for T in [Float32, Float64, ComplexF32, ComplexF64]
            A = sprand(T, 10, 10, 0.1)
            B = rand(T, 10, 2)
            C = rand(T, 10, 2)
            dA = CuSparseMatrixCSR(A)
            dB = CuArray(B)
            dC = CuArray(C)

            alpha = 1.2
            beta = 1.3
            mm!('N', 'N', alpha, dA, dB, beta, dC, 'O', algo)
            @test alpha * A * B + beta * C ≈ collect(dC)
        end
    end

    @testset "mv algo=$algo" for algo in [
        CUSPARSE.CUSPARSE_SPMV_ALG_DEFAULT,
        CUSPARSE.CUSPARSE_SPMV_CSR_ALG1,
        CUSPARSE.CUSPARSE_SPMV_CSR_ALG2,
    ]
        @testset "mv $T" for T in [Float32, Float64, ComplexF32, ComplexF64]
            A = sprand(T, 10, 10, 0.1)
            B = rand(T, 10)
            C = rand(T, 10)
            dA = CuSparseMatrixCSR(A)
            dB = CuArray(B)
            dC = CuArray(C)

            alpha = 1.2
            beta = 1.3
            mv!('N', alpha, dA, dB, beta, dC, 'O', algo)
            @test alpha * A * B + beta * C ≈ collect(dC)
        end
    end
end # version >= 11.4.1


text:
using Decimals
using Base.Test

# Equality
@test isequal(Decimal(0, 2, -3), Decimal(0, 2, -3))
@test ~isequal(Decimal(0, 2, -3), Decimal(0, 2, 3))
@test isequal(Decimal(0, 2, -3), 0.002)
@test isequal(Decimal(1, 2, 0), -2)
@test ~isequal(Decimal(1, 2, 0), 2)

@test isequal(0.1, 0.1)
@test isequal(2, 2)

@test Decimal(0, 2, -3) == Decimal(0, 2, -3)
@test Decimal(0, 2, -3) != Decimal(0, 2, 3)
@test Decimal(0, 2, -3) == 0.002

@test -2 == Decimal(1, 2, 0)
@test 2 != Decimal(1, 2, 0)

@test Decimal(1, 2, 0) == -2
@test Decimal(1, 2, 0) != 2

@test 0.1 == 0.1
@test 2 == 2

# Triple equals
# Is there a way to redefine === in Julia 4.x?
# @test Decimal(0, 2, -3) === Decimal(0, 2, -3)
# @test decimal("0.2") === 0.2


text:
using Distributions, Turing, Random, Parameters
ProjDir = @__DIR__
cd(ProjDir)
include("LNR.jl")

@model model(data,Nr) = begin
    μ = Array{Real,1}(undef,Nr)
    μ ~ [Normal(0,3)]
    σ ~ Truncated(Cauchy(0, 1), 0, Inf)
    for i in 1:length(data)
        data[i] ~ LNR(μ, σ, 0.0)
    end
end

Random.seed!(343)
Nreps = 30
error_count = fill(0.0,Nreps)
Nr = 3

for i in 1:Nreps
    μ = rand(Normal(0 ,3), Nr)
    σ = rand(Uniform(.5, 2))
    dist = LNR(μ=μ, σ=σ, ϕ=0.0)
    data = rand(dist, 30)
    chain = sample(model(data, Nr), NUTS(1000, .8), 2000, discard_adapt = false, progress = false)
    error_count[i] = sum(get(chain,:numerical_error)[1])
end

describe(error_count)


text:
export Units

"""
    Units(unitname::String)

`Units` objects `show` a FastHenry .units command.
"""
struct Units <: Element
  unitname :: String
  Units() = new("") # flag for internal use only
  function Units(x)
    if ~issubset(Set([x]),Set(["km", "m", "cm", "mm", "um", "in", "mils"]))
      throw(ArgumentError("valid units are km, m, cm, mm, um, in, mils"))
    end
    new(x)
  end
end

Base.:(==)(x::Units, y::Units) = x.unitname == y.unitname
# mm are default if no unit line is present.


text:
# DESCRIPTION
# Simulates open-loop dynamics using a dynamic inverse controller.
#
# INPUT
# - acfun: name of system dynamics function
# - finp: name of control input function
# - Tsim: total simulation time
# - dt: time step
# - x0: initial state vector
# - u0: initial control vector
#
# OUTPUT
# - state: state vector history
# - time: time vector history
#
# ------------------------------------------------------------------------------

function simulate(acfun,finp,Tsim,dt,x0,u0)
    # set up simulation
    # vector of times
    time=collect(0:dt:Tsim)
    # initialize state time history
    state=zeros(NSTATES,length(time))
    # set initial condition
    state[:,1]=x0
    # simulate open-loop system
    for i=1:length(time)-1
        # integrate
        sol=rk4(acfun,finp,time[i],dt,state[:,i],u0)
        state[:,i+1]=sol
    end
    #
    return state, time
end

max_stars_repo_name: umbertosaetti/J-SimpleHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
# Translating
# https://proplot.readthedocs.io/en/latest/basics.html#Figures-and-subplots
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

# start of example
data = cumsum(2(rand(100, 5) .- 0.5), dims=1)

# Simple plot
f, axs = proplot.subplots(ncols=2)
axs[1].plot(data, lw=2)
axs[1].format(xticks=20, xtickminor=false)
[ax.format(
    suptitle="Simple subplot grid", title="Title",
    xlabel="x axis", ylabel="y axis"
) for ax in axs]
f.savefig("figs/Figures-and-subplots-1.svg", transparent=false)

# Complex grid
array = [  # the "picture"; 1 == subplot a, 2 == subplot b, etc.
    [1, 1, 2, 2],
    [0, 3, 3, 0],
]
f, axs = proplot.subplots(array, axwidth=1.8)
[ax.format(
    abc=true, abcloc="ul", suptitle="Complex subplot grid",
    xlabel="xlabel", ylabel="ylabel"
) for ax in axs]
axs[3].plot(data, lw=2)
f.savefig("figs/Figures-and-subplots-2.svg", transparent=false)

# Really complex grid
array = [  # the "picture"
    [1, 1, 2],
    [1, 1, 6],
    [3, 4, 4],
    [3, 5, 5],
]
f, axs = proplot.subplots(array, width=5, span=false)
[ax.format(
    suptitle="Really complex subplot grid",
    xlabel="xlabel", ylabel="ylabel", abc=true
) for ax in axs]
axs[1].plot(data, lw=2)
f.savefig("figs/Figures-and-subplots-3.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
using Test
using OpenQuantumSystems
using LinearAlgebra
# using Traceur

@testset "molecules" begin

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))

    fc1 = franckCondonFactors(3, 0.0)
    fc2 = [
        1.0+0.0im 0.0+0.0im 0.0+0.0im
        0.0+0.0im 1.0+0.0im 0.0+0.0im
        0.0+0.0im 0.0+0.0im 1.0+0.0im
    ]
    @test 1e-12 > D(fc1, fc2)

    fc1 = franckCondonFactors(3, 1.0)
    fc2 = [
        0.7788007830714049+0.0im -0.5506953149031837+0.0im 0.2753476574515917+0.0im
        0.5506953149031838+0.0im 0.3894003915357024+0.0im -0.5841005873035535+0.0im
        0.27534765745159184+0.0im 0.5841005873035536+0.0im 0.0973500978839258+0.0im
    ]
    @test 1e-12 > D(fc1, fc2)

    fc1 = franckCondonFactors(3, 2.0im)
    fc2 = franckCondonFactors(3, -2.0im)
    fc2 .= conj(fc2)
    @test 1e-12 > D(fc1, fc2)

    @test vibrationalIndices([2]) == [[1], [2]]

    inds = [[1, 1], [1, 2], [2, 1], [2, 2]]
    @test vibrationalIndices([2, 2]) == inds

    inds = [
        [1, 1],
        [1, 2],
        [1, 3],
        [2, 1],
        [2, 2],
        [2, 3],
        [3, 1],
        [3, 2],
        [3, 3],
        [4, 1],
        [4, 2],
        [4, 3],
    ]
    @test vibrationalIndices([4, 3]) == inds
    # @trace vibrationalIndices([4, 3])
    # println("dispatch 2")

    @test electronicIndices(3) == [[1, 1, 1], [2, 1, 1], [1, 2, 1], [1, 1, 2]]

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.4, 1.0)
    Energy = [0.0, 200.0]
    fcFactors = [franckCondonFactors(3, mode1.shift), franckCondonFactors(3, mode2.shift)]
    mol = Molecule([mode1, mode2], 3, Energy)
    @test mol.modes == [mode1, mode2]
    @test mol.Nvib == 3
    @test 1e-12 > D(mol.fcFactors, fcFactors)
    @test mol.E == Energy

    @test getMolStateEnergy(mol, 1, [1, 1]) ≈ 0.3
    @test getMolStateEnergy(mol, 1, [3, 1]) ≈ 0.7
    @test getMolStateEnergy(mol, 2, [2, 2]) ≈ 200.9

    @test getMolStateFC(mol, 1, [1, 1], 1, [1, 1]) == 1.0
    @test getMolStateFC(mol, 1, [2, 1], 1, [1, 2]) == 0.0
    @test getMolStateFC(mol, 1, [1, 1], 2, [1, 1]) ≈ 0.6065306597126334
    @test getMolStateFC(mol, 2, [3, 1], 2, [1, 2]) == 0.0

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.4, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1, mode2], 3, Energy)
    mol1.E = [0.0, 400.0]
    mol2 = updateMolecule(mol1)
    @test mol2.E == [0.0, 400.0]
    mol2.E = [0.0, 200.0]
    updateMolecule!(mol2)
    mol1.E = [0.0, 200.0]
    mol1 = updateMolecule(mol1)
    @test mol1.E == mol2.E

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.4, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1, mode2], 3, Energy)

    @test getMolShifts(mol1) == [1.0, 2.0]
    @test getMolFrequencies(mol1) == [0.2, 0.4]

end


text:
struct TanhModulation <: InfectionModulation
  weight_detected::Float64
  weight_deaths::Float64
  loc::Float64
  scale::Float64
  limit_value::Float64
  TanhModulation(;weight_detected::Real, weight_deaths::Real, loc::Real, scale::Real, limit_value::Real) = new(weight_detected, weight_deaths, loc, scale, limit_value)
end

function (f::TanhModulation)(state::SimState, params::SimParams, event::Event)
  @assert kind(event) == TransmissionEvent

  ck = contactkind(event)
  if ConstantKernelContact !== ck && SporadicContact !== ck && FriendshipContact !== ck
    return true # do not affect other types of contact than "outer" ones
  end

  num_detected = numdetected(state.stats)
  num_deaths = numdead(state.stats)
  fear = num_detected * f.weight_detected + num_deaths * f.weight_deaths
  x = (fear - f.loc) / f.scale
  scaling = ((1 - f.limit_value) / 2)
  base = (1 - (1 - f.limit_value) / 2)
  rand(state.rng) < -tanh(x) * scaling + base
end

# This all to avoid using @eval and others
const modulations = Dict{String, Type{T} where T}(
    "TanhModulation" => TanhModulation
)

make_infection_modulation(name::AbstractString; args...) = modulations[name](;args...)


text:
@inline function add_fpfp_db(x::T, y::T) where {T<:AbstractFloat}
    return DoubleFloat(add_fpfp_dd(x, y))
end

@inline function sub_fpfp_db(x::T, y::T) where {T<:AbstractFloat}
    return DoubleFloat(sub_fpfp_dd(x, y))
end

@inline function mul_fpfp_db(x::T, y::T) where {T<:AbstractFloat}
    return DoubleFloat(mul_fpfp_dd(x, y))
end

@inline function dvi_fpfp_db(x::T, y::T) where {T<:AbstractFloat}
    return DoubleFloat(dvi_fpfp_dd(x, y))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, FileWatching

# This script does the following
# Sets up n unix pipes
# For the odd pipes, a byte is written to the write end at intervals specified in intvls
# Nothing is written into the even numbered pipes
# Odd numbered pipes are tested for reads
# Even numbered pipes are tested for timeouts
# Writable ends are always tested for writability before a write

n = 20
intvls = [2, .2, .1, .005]

pipe_fds = Vector{Any}(uninitialized, n)
for i in 1:n
    @static if Sys.iswindows()
        pipe_fds[i] = Vector{Libc.WindowsRawSocket}(uninitialized, 2)
        0 == ccall(:wsasocketpair, Cint, (Cint, Cuint, Cint, Ptr{Libc.WindowsRawSocket}), 1, 1, 6, pipe_fds[i]) || error(Libc.FormatMessage())
    else
        pipe_fds[i] = Array{RawFD}(uninitialized, 2)
        @test 0 == ccall(:pipe, Cint, (Ptr{RawFD},), pipe_fds[i])
    end
end

function pfd_tst_reads(idx, intvl)
    global ready += 1
    wait(ready_c)
    t_elapsed = @elapsed begin
        start_evt2 = Condition()
        evt2 = @async (notify(start_evt2); poll_fd(pipe_fds[idx][1], intvl; readable=true, writable=false))
        wait(start_evt2); yield() # make sure the async poll_fd is pumping events
        evt = poll_fd(pipe_fds[idx][1], intvl; readable=true, writable=false)
    end
    @test !evt.timedout
    @test evt.readable
    @test !evt.writable
    @test evt === wait(evt2)

    # println("Expected ", intvl, ", actual ", t_elapsed, ", diff ", t_elapsed - intvl)
    # Disabled since this assertion fails randomly, notably on build VMs (issue #12824)
    # @test t_elapsed <= (intvl + 1)

    dout = Vector{UInt8}(uninitialized, 1)
    @static if Sys.iswindows()
        1 == ccall(:recv, stdcall, Cint, (Ptr{Void}, Ptr{UInt8}, Cint, Cint), pipe_fds[idx][1], dout, 1, 0) || error(Libc.FormatMessage())
    else
        @test 1 == ccall(:read, Csize_t, (Cint, Ptr{UInt8}, Csize_t), pipe_fds[idx][1], dout, 1)
    end
    @test dout[1] == Int8('A')
end

function pfd_tst_timeout(idx, intvl)
    global ready += 1
    wait(ready_c)
    t_elapsed = @elapsed begin
        start_evt2 = Condition()
        evt2 = @async (notify(start_evt2); poll_fd(pipe_fds[idx][1], intvl; readable=true, writable=false))
        wait(start_evt2); yield() # make sure the async poll_fd is pumping events
        evt = poll_fd(pipe_fds[idx][1], intvl; readable=true, writable=false)
        @test evt.timedout
        @test !evt.readable
        @test !evt.writable
        @test evt === wait(evt2)
    end

    # Disabled since these assertions fail randomly, notably on build VMs (issue #12824)
    # @test intvl <= t_elapsed
    # @test t_elapsed <= (intvl + 1)
end

# Odd numbers trigger reads, even numbers timeout
for (i, intvl) in enumerate(intvls)
    @sync begin
        global ready = 0
        global ready_c = Condition()
        t = Vector{Task}(uninitialized, n)
        for idx in 1:n
            if isodd(idx)
                t[idx] = @async pfd_tst_reads(idx, intvl)
            else
                t[idx] = @async pfd_tst_timeout(idx, intvl)
            end
        end

        while ready < n
            sleep(0.1)
        end
        ready = 0
        # tickle only the odd ones, but test for writablity for everyone
        for idx in 1:n
            event = poll_fd(pipe_fds[idx][2], 0.001; readable=true, writable=true)
            @test !event.timedout
            @test !event.readable
            @test event.writable

            if isodd(idx)
                @static if Sys.iswindows()
                    1 == ccall(:send, stdcall, Cint, (Ptr{Void}, Ptr{UInt8}, Cint, Cint), pipe_fds[idx][2], "A", 1, 0) || error(Libc.FormatMessage())
                else
                    @test 1 == ccall(:write, Csize_t, (Cint, Ptr{UInt8}, Csize_t), pipe_fds[idx][2], "A", 1)
                end
            end
        end
        notify(ready_c, all=true)
        for idx in 1:n
            wait(t[idx])
        end
    end
end

for i in 1:n
    for j = 1:2
        @static if Sys.iswindows()
            0 == ccall(:closesocket, stdcall, Cint, (Ptr{Void},), pipe_fds[i][j]) || error(Libc.FormatMessage())
        else
            @test 0 == ccall(:close, Cint, (Cint,), pipe_fds[i][j])
        end
    end
end

# issue #12473
# make sure 1-shot timers work
let a = []
    Timer(t -> push!(a, 1), 0.01, 0)
    sleep(0.2)
    @test a == [1]
end

# make sure repeating timers work
@noinline function make_unrooted_timer(a)
    t = Timer(0.0, 0.1)
    finalizer(t -> a[] += 1, t)
    wait(t)
    e = @elapsed for i = 1:5
        wait(t)
    end
    @test 1.5 > e >= 0.4
    @test a[] == 0
    nothing
end
let a = Ref(0)
    make_unrooted_timer(a)
    gc()
    @test a[] == 1
end

for f in (watch_file, poll_file)
    local f
    @test_throws ArgumentError f("adir\0bad")
end

#issue #12992
function test_12992()
    pfw = PollingFileWatcher(@__FILE__, 0.01)
    close(pfw)
    pfw = PollingFileWatcher(@__FILE__, 0.01)
    close(pfw)
    pfw = PollingFileWatcher(@__FILE__, 0.01)
    close(pfw)
    gc()
    gc()
end

# Make sure multiple close is fine
function test2_12992()
    pfw = PollingFileWatcher(@__FILE__, 0.01)
    close(pfw)
    close(pfw)
    pfw = PollingFileWatcher(@__FILE__, 0.01)
    close(pfw)
    close(pfw)
    pfw = PollingFileWatcher(@__FILE__, 0.01)
    close(pfw)
    close(pfw)
    gc()
    gc()
end

test_12992()
test_12992()
test_12992()

test2_12992()
test2_12992()
test2_12992()

#######################################################################
# This section tests file watchers.                                   #
#######################################################################
dir = mktempdir()
file = joinpath(dir, "afile.txt")
# like touch, but lets the operating system update the timestamp
# for greater precision on some platforms (windows)
@test close(open(file,"w")) === nothing

function test_file_poll(channel,interval,timeout_s)
    rc = poll_file(file, interval, timeout_s)
    put!(channel,rc)
end

function test_timeout(tval)
    t_elapsed = @elapsed begin
        channel = Channel(1)
        @async test_file_poll(channel, 10, tval)
        tr = take!(channel)
    end
    @test tr[1] === Base.Filesystem.StatStruct() && tr[2] === EOFError()
    @test tval <= t_elapsed
end

function test_touch(slval)
    tval = slval*1.1
    channel = Channel(1)
    @async test_file_poll(channel, tval/3, tval)
    sleep(tval/3)  # one poll period
    f = open(file,"a")
    write(f,"Hello World\n")
    close(f)
    tr = take!(channel)
    @test ispath(tr[1]) && ispath(tr[2])
end

function test_watch_file_timeout(tval)
    watch = @async watch_file(file, tval)
    @test wait(watch) == FileWatching.FileEvent(false, false, true)
end

function test_watch_file_change(tval)
    watch = @async watch_file(file, tval)
    sleep(tval/3)
    open(file, "a") do f
        write(f, "small change\n")
    end
    @test wait(watch) == FileWatching.FileEvent(false, true, false)
end

function test_monitor_wait(tval)
    fm = FileMonitor(file)
    @async begin
        sleep(tval)
        f = open(file,"a")
        write(f,"Hello World\n")
        close(f)
    end
    fname, events = wait(fm)
    close(fm)
    if Sys.islinux() || Sys.iswindows() || Sys.isapple()
        @test fname == basename(file)
    else
        @test fname == ""  # platforms where F_GETPATH is not available
    end
    @test events.changed
end

function test_monitor_wait_poll()
    pfw = PollingFileWatcher(file, 5.007)
    @async begin
        sleep(2.5)
        f = open(file,"a")
        write(f,"Hello World\n")
        close(f)
    end
    (old, new) = wait(pfw)
    close(pfw)
    @test new.mtime - old.mtime > 2.5 - 1.5 # mtime may only have second-level accuracy (plus add some hysteresis)
end

test_timeout(0.1)
test_timeout(1)
test_touch(6)
test_monitor_wait(0.1)
test_monitor_wait(0.1)
test_monitor_wait_poll()
test_monitor_wait_poll()
test_watch_file_timeout(0.1)
test_watch_file_change(6)

@test_throws Base.UVError watch_file("____nonexistent_file", 10)
@test(@elapsed(
    @test(poll_file("____nonexistent_file", 1, 3.1) ===
          (Base.Filesystem.StatStruct(), EOFError()))) > 3)

rm(file)
rm(dir)

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
# energy (epsxc) only
function XC_c_pbe_E( rho, grho )

    ga = 0.0310906908696548950
    be = 0.06672455060314922
    third = 1.0/3.0
    pi34 = 0.6203504908994
    xkf = 1.919158292677513
    xks = 1.128379167095513

    rs = pi34/rho^third
    ec = XC_c_pw_E( rho )

    kf = xkf/rs
    ks = xks * sqrt(kf)
    t = sqrt(grho) / (2.0 * ks * rho)

    expe = exp(-ec/ga)
    af = be / ga * (1.0 / (expe - 1.0) )

    y = af * t * t

    xy = (1.0 + y) / (1.0 + y + y * y)

    qy = y * y * (2.0 + y) / (1.0 + y + y * y)^2

    s1 = 1.0 + be / ga * t * t * xy
    h0 = ga * log(s1)

    sc = rho * h0

    return sc

end


text:
#
#
# Hello World in Many Programming Languages!!!
#
# Authors:
# @author Ruben Andre Barreiro
#
#

# Prints "Hello World!!!" to the terminal window
println("Hello World!!!")

max_stars_repo_path: programming, markup & markdown languages/Julia/HelloWorld.jl


text:
#####
##### Builder
#####

f(x) = begin
    y = 10 + x
    z = 20 + y
    q = 20 + z
    return q + 50
end

function g(x)
    try
        if x > 1
            x + g(x - 1)
        else
            return 1
        end
        while true
            println("Nice!")
        end
        return
    catch e
        return 0
    end
end

function fn(x, y)
    z = 10
    if z > 10
        n = 10
        return x + y
    else
        return x + y + z
    end
end

@testset "Builder -- iterate" begin
    ir = code_info(f, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    local c = 1
    for (v, st) in p
        @test v == var(c)
        println(getindex(p, v))
        c += 1
    end
    ir = code_info(g, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    c = 1
    for (v, st) in p
        @test v == var(c)
        c += 1
    end
end

@testset "Builder -- setindex!" begin
    ir = code_info(f, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    len = length(p.from.code)
    for (v, st) in p
        p[v] = st
    end
    ir = finish(p)
    @test length(p.to.code) == len
    for (v, st) in CodeInfoTools.Builder(ir)
        st isa Expr && st.head == :call || continue
        @test st.args[1] == Base.:(+)
    end
    ir = code_info(g, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    len = length(p.from.code)
    for (v, st) in p
        p[v] = st
    end
    ir = finish(p)
    @test length(p.to.code) == len
    for (v, st) in CodeInfoTools.Builder(ir)
        st isa Expr && st.head == :call || continue
        @test (st.args[1] in (g, Base.:(+), Base.:(-), Base.:(>), Base.println))
    end
end

@testset "Builder -- push!" begin
    ir = code_info(f, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    c = deepcopy(p.from.code)
    for (v, st) in p
    end
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == length(c) + 1
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == length(c) + 2
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == length(c) + 3
    @test unwrap(p[end - 3]) isa Core.ReturnNode
    ir = code_info(g, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    c = deepcopy(p.from.code)
    for (v, st) in p
    end
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == length(c) + 1
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == length(c) + 2
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == length(c) + 3
    @test unwrap(p[end - 3]) isa Core.ReturnNode
end

@testset "Builder -- pushfirst!" begin
    ir = code_info(f, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == 1
    pushfirst!(p, Expr(:call, rand))
    @test length(p.to.code) == 2
    pushfirst!(p, Expr(:call, rand))
    @test length(p.to.code) == 3
    ir = code_info(g, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    push!(p, Expr(:call, rand))
    @test length(p.to.code) == 1
    pushfirst!(p, Expr(:call, rand))
    @test length(p.to.code) == 2
    pushfirst!(p, Expr(:call, rand))
    @test length(p.to.code) == 3
end

@testset "Builder -- delete!" begin
    ir = code_info(f, Tuple{Int})
    p = CodeInfoTools.Builder(ir)
    for (v, st) in p
    end
    display(p)
    println()
    delete!(p, var(1))
    display(p)
    println()
    insert!(p, var(2), Expr(:call, rand))
    display(p)
    println()
    display(renumber(p.to))
    println()
    display(finish(p))
    println()
end


text:
function init_pot_Hcoul( LF, center )
    Npoints = LF.Nx * LF.Ny * LF.Nz
    Vpot = zeros( Float64,Npoints )
    for ip = 1:Npoints
        r = norm( LF.lingrid[:,ip] - center[:] )
        Vpot[ip] = -1.0/r
    end
    return Vpot
end