text:
@testset "Checked arithmetic" begin

    using Base: checked_abs, checked_neg, checked_add, checked_sub, checked_mul,
                checked_cld, checked_div, checked_fld, checked_mod, checked_rem
    using Base.Checked: add_with_overflow, sub_with_overflow, mul_with_overflow
    using HalfIntegers: checked_twice

    @testset "Base.checked_abs" begin
        for T in halfinttypes
            @eval @test checked_abs($T(5//2)) ==ₜ $T(5//2)
            @eval @test checked_abs($T(-3//2)) ==ₜ $T(3//2)
            @eval @test_throws OverflowError checked_abs(typemin($T))
            @eval @test abs(typemin($T)) ==ₜ typemin($T) # unchecked `abs` overflows
        end
        for T in halfuinttypes
            @eval @test checked_abs($T(5//2)) ==ₜ $T(5//2)
        end
        @test checked_abs(BigHalfInt(-3//2)) ==ₜ BigHalfInt(3//2)
        @test checked_abs(BigHalfInt(5//2)) ==ₜ BigHalfInt(5//2)
    end

    @testset "Base.checked_neg" begin
        for T in halfinttypes
            @eval @test checked_neg(typemax($T)) ==ₜ typemin($T) + onehalf($T)
            @eval @test_throws OverflowError checked_neg(typemin($T))
            @eval @test -typemin($T) ==ₜ typemin($T) # unchecked `-` overflows
        end
        for T in halfuinttypes
            @eval @test checked_neg(zero($T)) ==ₜ zero($T)
            @eval @test_throws OverflowError checked_neg(onehalf($T))
            @eval @test -onehalf($T) ==ₜ typemax($T) # unchecked `-` overflows
        end
        @test checked_neg(BigHalfInt(-3//2)) ==ₜ BigHalfInt(3//2)
        @test checked_neg(BigHalfInt(5//2)) ==ₜ BigHalfInt(-5//2)
    end

    @testset "HalfIntegers.checked_twice" begin
        for T in (inttypes..., uinttypes...)
            @eval @test checked_twice($T(5)) ==ₜ $T(10)
            @eval @test_throws OverflowError checked_twice(typemax($T))
            @eval @test checked_twice(typemax(Half{$T})) ==ₜ typemax($T)
        end
        @test twice(typemax(Int)) ==ₜ -2 # unchecked `twice` overflows
        @test twice(typemax(UInt)) ==ₜ typemax(UInt) - one(UInt) # unchecked `twice` overflows
        @test checked_twice(BigInt(100)) ==ₜ BigInt(200)
        @test checked_twice(BigHalfInt(100)) ==ₜ BigInt(200)
        @test checked_twice(0.375) ==ₜ 0.75
        @test checked_twice(3//8) ==ₜ 3//4
    end

    @testset "Base.checked_add" begin
        for T in (halfinttypes..., halfuinttypes...)
            @eval @test checked_add($T(5//2), $T(3//2)) ==ₜ $T(4)
            @eval @test_throws OverflowError checked_add(typemax($T), one($T))
            @eval @test typemax($T) + onehalf($T) == typemin($T) # unchecked `+` overflows
        end
        @test checked_add(BigHalfInt(5//2), BigHalfInt(-3//2)) ==ₜ BigHalfInt(1)
        # Mixed HalfInteger types
        @test_throws InexactError checked_add(HalfInt(-1), HalfUInt(3/2))
        @test_throws InexactError checked_add(HalfInt(-2), HalfUInt(3/2))
        @test_broken HalfInt(-1) + HalfUInt(3/2) ==ₜ HalfUInt(1/2) # This should work, since it matches the Integer behavior
        @test_broken HalfInt(-2) + HalfUInt(3/2) ==ₜ typemax(HalfUInt) # This should work, since it matches the Integer behavior
        # Mixed Integer/HalfInteger
        @test checked_add(HalfInt8(50), Int8(10)) ==ₜ HalfInt8(60)
        @test_throws OverflowError checked_add(HalfInt8(50), Int8(100))
        @test_throws OverflowError checked_add(HalfInt8(-50), Int8(100))
        @test HalfInt8(50) + Int8(100) ==ₜ HalfInt8(22) # (100 + 200) % Int8 == 44
        @test HalfInt8(-50) + Int8(100) ==ₜ HalfInt8(50) # (-100 + 200) % Int8 == 100
    end

    @testset "Base.checked_sub" begin
        for T in (halfinttypes..., halfuinttypes...)
            @eval @test checked_sub($T(5//2), $T(3//2)) ==ₜ $T(1)
            @eval @test_throws OverflowError checked_sub(typemin($T), one($T))
            @eval @test typemin($T) - onehalf($T) == typemax($T) # unchecked `-` overflows
        end
        @test checked_sub(BigHalfInt(5//2), BigHalfInt(-3//2)) ==ₜ BigHalfInt(4)
        # Mixed HalfInteger types
        @test_throws InexactError checked_sub(HalfUInt(3/2), HalfInt(-1))
        @test_throws InexactError checked_sub(HalfInt(-1/2), HalfUInt(0))
        @test_broken HalfUInt(3/2) - HalfInt(-1) ==ₜ HalfUInt(5/2) # This should work, since it matches the Integer behavior
        @test_broken HalfInt(-1/2) - HalfUInt(0) ==ₜ typemax(HalfUInt) # This should work, since it matches the Integer behavior
        # Mixed Integer/HalfInteger
        @test checked_sub(HalfInt8(50), Int8(10)) ==ₜ HalfInt8(40)
        @test_throws OverflowError checked_sub(HalfInt8(50), Int8(100))
        @test_throws OverflowError checked_sub(HalfInt8(-50), Int8(100))
        @test HalfInt8(50) - Int8(100) ==ₜ HalfInt8(-50) # (100 - 200) % Int8 == -100
        @test HalfInt8(-50) - Int8(100) ==ₜ HalfInt8(-22) # (-100 - 200) % Int8 == -44
    end

    @testset "Base.checked_mul" begin
        for T in (halfinttypes..., halfuinttypes...)
            @eval @test typemax($T) * one($T) ==ₜ float(typemax($T)) # `*(::HalfInteger, ::HalfInteger)` uses float arithmetic, so no overflow
            @eval @test checked_mul(typemax($T), one($T)) ==ₜ float(typemax($T))
        end
        @test checked_mul(BigHalfInt(7/2), BigHalfInt(3/2)) ==ₜ BigFloat(5.25)
        # Mixed HalfInteger types
        @test checked_mul(HalfInt8(7/2), HalfUInt16(3/2)) ==ₜ 5.25
        # Mixed HalfInteger/Integer
        @test Int8(100) * HalfInt8(1) == HalfInt8(-28) # (100*2) % Int8 == -56
        @test HalfInt8(50) * Int8(2) == HalfInt8(-28) # (100*2) % Int8 == -56
        @test_throws OverflowError checked_mul(Int8(100), HalfInt8(1))
        @test_throws OverflowError checked_mul(HalfInt8(50), Int8(2))
    end

    @testset "Base.$checked_op" for (op, checked_op) = ((:cld,:checked_cld), (:div,:checked_div), (:fld,:checked_fld))
        for T in inttypes
            @eval @test_throws DivideError $op(typemin(Half{$T}), Half{$T}(-1/2)) # op is already checked
            @eval @test_throws DivideError $checked_op(typemin(Half{$T}), Half{$T}(-1/2))
            @eval @test $checked_op(Half{$T}(6), Half{$T}(5/2)) ==ₜ $op(Half{$T}(6), Half{$T}(5/2))
            @eval @test $checked_op(Half{$T}(-6), Half{$T}(5/2)) ==ₜ $op(Half{$T}(-6), Half{$T}(5/2))
        end
        for T in uinttypes
            @eval @test $checked_op(Half{$T}(6), Half{$T}(5/2)) ==ₜ $op(Half{$T}(6), Half{$T}(5/2))
        end
        @eval @test $checked_op(BigHalfInt(6), BigHalfInt(5/2)) ==ₜ $op(BigHalfInt(6), BigHalfInt(5/2))
        @eval @test $checked_op(-BigHalfInt(6), BigHalfInt(5/2)) ==ₜ $op(-BigHalfInt(6), BigHalfInt(5/2))
        # Mixed HalfInteger/Integer
        for T in (inttypes..., uinttypes...)
            @eval @test $checked_op($T(6), Half{$T}(5/2)) ==ₜ $op($T(6), Half{$T}(5/2))
            @eval @test $checked_op(Half{$T}(15/2), $T(5)) ==ₜ $op(Half{$T}(15/2), $T(5))
            @eval @test_throws OverflowError $checked_op(typemax($T), Half{$T}(5/2))
            @eval @test_throws OverflowError $checked_op(Half{$T}(15/2), typemax($T))
            @eval @test_broken try; $op(typemax($T), Half{$T}(5/2)); catch e; e; end isa OverflowError
            @eval @test_broken try; $op(Half{$T}(15/2), typemax($T)); catch e; e; end isa OverflowError
        end
        @eval @test $checked_op(BigInt(6), BigHalfInt(5/2)) ==ₜ $op(BigInt(6), BigHalfInt(5/2))
        @eval @test $checked_op(-BigInt(6), BigHalfInt(5/2)) ==ₜ $op(-BigInt(6), BigHalfInt(5/2))
        @eval @test $checked_op(BigHalfInt(15/2), BigInt(5)) ==ₜ $op(BigHalfInt(15/2), BigInt(5))
        @eval @test $checked_op(-BigHalfInt(15/2), BigInt(5)) ==ₜ $op(-BigHalfInt(15/2), BigInt(5))
    end

    @testset "Base.$checked_op" for (op, checked_op) = ((:mod,:checked_mod), (:rem,:checked_rem))
        for T in (inttypes..., :BigInt)
            @eval @test $checked_op(Half{$T}(6), Half{$T}(5/2)) ==ₜ $op(Half{$T}(6), Half{$T}(5/2))
            @eval @test $checked_op(Half{$T}(-6), Half{$T}(5/2)) ==ₜ $op(Half{$T}(-6), Half{$T}(5/2))
        end
        for T in uinttypes
            @eval @test $checked_op(Half{$T}(6), Half{$T}(5/2)) ==ₜ $op(Half{$T}(6), Half{$T}(5/2))
        end
        # Mixed HalfInteger/Integer
        for T in (inttypes..., uinttypes...)
            @eval @test $checked_op($T(6), Half{$T}(5/2)) ==ₜ $op($T(6), Half{$T}(5/2))
            @eval @test $checked_op(Half{$T}(15/2), $T(5)) ==ₜ $op(Half{$T}(15/2), $T(5))
            @eval @test_throws OverflowError $checked_op(typemax($T), Half{$T}(5/2))
            @eval @test_throws OverflowError $checked_op(Half{$T}(15/2), typemax($T))
            @eval @test_broken try; $op(typemax($T), Half{$T}(5/2)); catch e; e; end isa OverflowError
            @eval @test_broken try; $op(Half{$T}(15/2), typemax($T)); catch e; e; end isa OverflowError
        end
        @eval @test $checked_op(BigInt(6), BigHalfInt(5/2)) ==ₜ $op(BigInt(6), BigHalfInt(5/2))
        @eval @test $checked_op(-BigInt(6), BigHalfInt(5/2)) ==ₜ $op(-BigInt(6), BigHalfInt(5/2))
        @eval @test $checked_op(BigHalfInt(15/2), BigInt(5)) ==ₜ $op(BigHalfInt(15/2), BigInt(5))
        @eval @test $checked_op(-BigHalfInt(15/2), BigInt(5)) ==ₜ $op(-BigHalfInt(15/2), BigInt(5))
    end

    @testset "Base.Checked.add_with_overflow" begin
        for T in (halfinttypes..., halfuinttypes...)
            @eval @test add_with_overflow($T(5//2), $T(3//2)) ==ₜ ($T(4), false)
            @eval @test add_with_overflow(typemax($T), onehalf($T)) ==ₜ (typemin($T), true)
        end
        @test add_with_overflow(BigHalfInt(5//2), BigHalfInt(-3//2)) ==ₜ (BigHalfInt(1), false)
    end

    @testset "Base.Checked.sub_with_overflow" begin
        for T in (halfinttypes..., halfuinttypes...)
            @eval @test sub_with_overflow($T(5//2), $T(3//2)) ==ₜ ($T(1), false)
            @eval @test sub_with_overflow(typemin($T), onehalf($T)) ==ₜ (typemax($T), true)
        end
        @test sub_with_overflow(BigHalfInt(5//2), BigHalfInt(-3//2)) ==ₜ (BigHalfInt(4), false)
    end

    @testset "Base.Checked.mul_with_overflow" begin
        for T in (halfinttypes..., halfuinttypes...)
            @eval @test mul_with_overflow(typemax($T), one($T)) ==ₜ (float(typemax($T)), false)
        end
        @test mul_with_overflow(BigHalfInt(7/2), BigHalfInt(3/2)) ==ₜ (BigFloat(5.25), false)
        # Mixed HalfInteger types
        @test mul_with_overflow(HalfInt8(7/2), HalfUInt16(3/2)) ==ₜ (5.25, false)
        # Mixed HalfInteger/Integer
        @test mul_with_overflow(Int(100), HalfInt(3/2)) ==ₜ (HalfInt(150), false)
        @test mul_with_overflow(Int8(100), HalfInt8(3/2)) ==ₜ (HalfInt8(22), true)
        @test mul_with_overflow(HalfInt(101/2), Int(2)) ==ₜ (HalfInt(101), false)
        @test mul_with_overflow(HalfInt8(101/2), Int8(2)) ==ₜ (HalfInt8(-27), true)
    end

    @static if VERSION ≥ v"1.1"
        # In Julia < 1.1, zero(Half{T}) and one(Half{T}) do not work for
        # T ∈ [SafeInt8, SafeInt16, SafeUInt8, SafeUInt16] due to a bug in Julia, cf.
        # https://github.com/JuliaLang/julia/issues/32203

        using SaferIntegers

        safeinttypes  = (:SafeInt, :SafeInt8, :SafeInt16, :SafeInt32, :SafeInt64, :SafeInt128)
        safeuinttypes = (:SafeUInt, :SafeUInt8, :SafeUInt16, :SafeUInt32, :SafeUInt64, :SafeUInt128)

        @testset "SaferIntegers" begin
            @test_throws OverflowError Half{SafeInt8}(typemax(Int8))
            @test_throws OverflowError typemax(Half{SafeInt8}) + one(HalfInt8)
            @test_throws OverflowError typemax(Half{SafeInt8}) + one(Int8)
            @test_throws OverflowError typemin(Half{SafeInt8}) - one(HalfInt8)
            @test_throws OverflowError typemin(Half{SafeInt8}) - one(Int8)
            for T in safeinttypes
                @eval @test HalfInteger(one($T)) === one(Half{$T})
                @eval @test_throws OverflowError Half{$T}(typemax($T))
                @eval @test_throws OverflowError Half{$T}(typemin($T))
                @eval @test_throws OverflowError Half{$T}(typemax($T)//one($T))
                @eval @test_throws OverflowError Half{$T}(typemin($T)//one($T))
                @eval @test_throws OverflowError Half{$T}(complex(typemax($T)))
                @eval @test_throws OverflowError Half{$T}(complex(typemin($T)))
                @eval @test_throws OverflowError Half{$T}(complex(typemax($T)//one($T)))
                @eval @test_throws OverflowError Half{$T}(complex(typemin($T)//one($T)))
                @eval @test_throws OverflowError -typemin(Half{$T})
                @eval @test_throws OverflowError abs(typemin(Half{$T}))
                @eval @test_throws OverflowError flipsign(typemin(Half{$T}), -1)
                @eval @test_throws OverflowError typemax(Half{$T}) + onehalf(Half{$T})
                @eval @test_throws OverflowError typemin(Half{$T}) - onehalf(Half{$T})
                @eval @test_throws OverflowError twice(typemax($T))
                @eval @test_throws OverflowError twice(typemin($T))
                @eval @test_throws OverflowError round(typemax(Half{$T}))
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundNearest)
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundNearestTiesAway)
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundNearestTiesUp)
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundUp)
            end
            for T in safeuinttypes
                @eval @test HalfInteger(one($T)) === one(Half{$T})
                @eval @test_throws OverflowError Half{$T}(typemax($T))
                @eval @test_throws OverflowError Half{$T}(typemax($T)//one($T))
                @eval @test_throws OverflowError Half{$T}(complex(typemax($T)))
                @eval @test_throws OverflowError Half{$T}(complex(typemax($T)//one($T)))
                @eval @test_throws OverflowError -onehalf(Half{$T})
                @eval @test_throws OverflowError flipsign(onehalf(Half{$T}), -1)
                @eval @test_throws OverflowError typemax(Half{$T}) + onehalf(Half{$T})
                @eval @test_throws OverflowError typemin(Half{$T}) - onehalf(Half{$T})
                @eval @test_throws OverflowError twice(typemax($T))
                @eval @test_throws OverflowError round(typemax(Half{$T}))
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundNearest)
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundNearestTiesAway)
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundNearestTiesUp)
                @eval @test_throws OverflowError round(typemax(Half{$T}), RoundUp)
            end
        end
    end
end


text:
function test1(name, n)
    while n > 0
        println("$name $n")
        n -= 1
        yield()
    end
    false
end

t1 = @task test1("foo", 5)
t2 = @task test1("bar", 6)
t3 = @task test1("baz", 4)
schedule(t1)
schedule(t2)
schedule(t3)
wait(t2)

ext: jl


text: const numTests = Ref(20)


text:
using HTSLib
using HTSLib: htslib  # low-level API module
using Test
using FormatSpecimens
import HttpIO

# Low-level API
# -------------

@testset "Low-level API" begin
    # check struct size (these sizes are derived from C's sizeof operator)
    if Sys.isunix() && Sys.ARCH == :x86_64  # LP64
        @test sizeof(htslib.kstring_t) == 24
        @test sizeof(htslib.htsFormat) == 32
        @test sizeof(htslib.htsFile) == 128
        @test sizeof(htslib.bam1_core_t) == 48
        @test sizeof(htslib.bam1_t) == 80
    end

    @testset "cigar" begin
        bamfile = joinpath(path_of_format("BAM"), "ce#1.bam")
        reader = open(HTSReadWriter, bamfile)
        record = read(reader)
        GC.@preserve record begin
            record_ptr = pointer(record)
            n_cigar = record.n_cigar
            @test n_cigar == 3
            cigar1 = unsafe_load(htslib.bam_get_cigar(record_ptr), 1)
            cigar2 = unsafe_load(htslib.bam_get_cigar(record_ptr), 2)
            cigar3 = unsafe_load(htslib.bam_get_cigar(record_ptr), 3)
            @test Char(htslib.bam_cigar_opchar(htslib.bam_cigar_op(cigar1))) == 'M'
            @test Char(htslib.bam_cigar_opchar(htslib.bam_cigar_op(cigar2))) == 'D'
            @test Char(htslib.bam_cigar_opchar(htslib.bam_cigar_op(cigar3))) == 'M'
            @test htslib.bam_cigar_oplen(cigar1) == 27
            @test htslib.bam_cigar_oplen(cigar2) == 1
            @test htslib.bam_cigar_oplen(cigar3) == 73
        end
    end
end

# High-level API
# --------------

@testset "HTSLib version" begin
    @test HTSLib.HTSLIB_VERSION == v"1.10.2"
end

@testset "SAM header" begin
    bamfile = joinpath(path_of_format("BAM"), "ce#1.bam")
    reader = open(HTSReadWriter, bamfile)
    header_str = string(header(reader))

    @test isa(header_str, String)
    @test string(SamHeader(header_str)) == header_str
end

@testset "BAM reading" begin
    @testset "Read local bam" begin
        bamfile = joinpath(path_of_format("BAM"), "ce#1.bam")
        reader = open(HTSReadWriter, bamfile)

        @test eltype(reader) == HTSRecord
        @test isa(header(reader), SamHeader)

        record = HTSRecord()
        read!(reader, record)

        @test refid(record) == 1
        @test leftposition(record) == 2
        @test rightposition(record) == 102
        @test sequence(record, 3) == 'T'
        @test seqlength(record) == 100
        @test sequence(record) == collect("CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA")
        @test quality(record) == [Int(x) - 33 for x in "#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"]
        @test eltype(quality(record)) == UInt8
    end

    @testset "iterator" begin
        bamfile1 = joinpath(path_of_format("BAM"), "ce#1.bam")
        bamfile2 = joinpath(path_of_format("BAM"), "ce#2.bam")

        @test length(collect(HTSReadWriter(bamfile1))) == 1
        @test length(collect(HTSReadWriter(bamfile2))) == 2
    end
end

@testset "BAM reading from url" begin
    ## TODO: test invalid url
    bamfile1 = "https://raw.githubusercontent.com/BioJulia/FormatSpecimens.jl/master/BAM/ce#1.bam"
    bamfile2 = "https://raw.githubusercontent.com/BioJulia/FormatSpecimens.jl/master/BAM/ce#2.bam"

    @test length(collect(HTSReadWriter(bamfile1))) == 1
    @test length(collect(HTSReadWriter(bamfile2))) == 2
end

@testset "BAM/CRAM reading from gsurl" begin
    # skip test if gcloud is not installed
    if isnothing(Sys.which("gcloud"))
        return
    end

    gcsclient = HttpIO.GCSClient()

    cram = "gs://gatk-test-data/wgs_cram/NA12878_20k_hg38/NA12878.cram"
    crai = "gs://gatk-test-data/wgs_cram/NA12878_20k_hg38/NA12878.crai"

    bam = "gs://gatk-test-data/wgs_bam/NA12878_20k_hg38/NA12878.bam"
    bai = "gs://gatk-test-data/wgs_bam/NA12878_20k_hg38/NA12878.bai"

    @testset "Simple BAM read" begin
        bamio = HTSReadWriter(bam; index=bai, gcsclient=gcsclient)
        record = HTSRecord()
        read!(bamio, record)
        @test leftposition(record) == 10001
        @test rightposition(record) == 10110
    end

    @testset "Simple CRAM read" begin
        cramio = HTSReadWriter(cram; gcsclient=gcsclient) # TODO: hts_idx_load2 does not recognize cram index file
        record = HTSRecord()
        read!(cramio, record)
        @test leftposition(record) == 10001
        @test rightposition(record) == 10110
    end

    @testset "Bam random access" begin
        bamio = HTSReadWriter(bam; index=bai, gcsclient=gcsclient)

        @testset "single region" begin
            bamitr1 = HTSRegionsIterator(bamio, "chr1:1000000-2000000")
            bamitr2 = HTSRegionsIterator(bamio, "chr1", 1000000, 2000000)

            @test length(collect(bamitr1)) == 25
            @test length(collect(bamitr2)) == 25

            @test !isnothing(tryread!(bamio, HTSRecord()))
        end

        @testset "unmapped reads" begin
            bamitr = HTSRegionsIterator(bamio, "*")
            @test length(collect(bamitr)) == 442

            ## we have iterated unmapped reads at end of the file
            @test isnothing(tryread!(bamio, HTSRecord()))
        end

        @testset "GC" begin
            itr = open(HTSReadWriter, bam; index=bai, gcsclient=gcsclient) do reader
                itr = HTSRegionsIterator(reader, "chr1", 1000000, 2000000)
                @test !isnothing(read(itr))
                itr
            end
            # reader has been closed, itr is invalid now, but it should be safe to finalize it.
            finalize(itr)
        end
    end

    @testset "CRAM random access" begin
        # TODO: add index
        #cramio = HTSReadWriter(cram; gcsclient=gcsclient)

        # @testset "from file start" begin
        #     itr = open(HTSReadWriter, cram; gcsclient=gcsclient) do reader
        #         itr = HTSRegionsIterator(reader, ".")
        #         @test !isnothing(read(itr))
        #         itr
        #     end
        #     finalize(itr)
        # end
    end
end

@testset "Write hts file" begin
    @testset "bam" begin
        bamfile = joinpath(path_of_format("BAM"), "ce#1.bam")
        bamin = open(HTSReadWriter, bamfile)
        (path, io) = mktemp()
        bamout = HTSReadWriter(io, "wbz")

        # write header
        write(bamout, header(bamin))
        # write records
        for r in bamin
            write(bamout, r)
        end

        # flush and close
        close(bamout)

        # reread
        f = HTSReadWriter(path)
        @test HTSLib.format_description(f) == "BAM version 1 compressed sequence data"
        records = collect(f)
        @test length(records) == 1
        @test refid(records[1]) == 1
        @test rightposition(records[1]) == 102
    end
end

@testset "Regions iterator" begin

    @testset "Constructor" begin

        bam = joinpath(path_of_format("BAM"), "R_12h_D06.uniq.q40.bam")
        bai = joinpath(path_of_format("BAM"), "R_12h_D06.uniq.q40.bam.bai")
        reader = HTSReadWriter(bam; index=bai)

        # Vector{String}
        @test length(collect(HTSRegionsIterator(reader, ["chrM:1-500"]))) == 2
        # String
        @test length(collect(HTSRegionsIterator(reader, "chrM:1-500"))) == 2
        # Tuple
        @test length(collect(HTSRegionsIterator(reader, ("chrM", 1, 500)))) == 2
        @test length(collect(HTSRegionsIterator(reader, (22, 1, 500)))) == 2
        # Vector{Tuple}
        @test length(collect(HTSRegionsIterator(reader, [("chrM", 1, 500)]))) == 2
        @test length(collect(HTSRegionsIterator(reader, [(22, 1, 500)]))) == 2
        # args
        @test length(collect(HTSRegionsIterator(reader, "chrM", 1, 500))) == 2
        @test length(collect(HTSRegionsIterator(reader, 22, 1, 500))) == 2

    end
end


text:
# Load the grammar type.
include("grammar_derivation/grammar.jl")

"""
Provides a definition for a grammar derivation representation.
"""
type GrammarDerivationRepresentationDefinition <: RepresentationDefinition
  """
  The maximum number of times that the genome is allowed to wrap around itself
  if it has run out of codons to use. By default, wrapping is disabled.
  """
  max_wraps::Int

  """
  The grammar that should be used to produce the derivations.
  """
  grammar::Grammar

  GrammarDerivationRepresentationDefinition() = new(0)
end

type GrammarDerivationRepresentation <: Representation
  max_wraps::Int
  grammar::Grammar

  GrammarDerivationRepresentation(max_wraps::Int, grammar::Grammar) =
    new(max_wraps, grammar)
end

chromosome(::GrammarDerivationRepresentation) = AbstractString

"""
TODO: Document grammar derivation representations.
"""
function grammar(def::Function)
  g = GrammarDerivationRepresentationDefinition()
  def(g)
  g
end

"""
Composes a grammar derivation representation from a provided definition.
"""
function compose!(def::GrammarDerivationRepresentationDefinition)
  GrammarDerivationRepresentation(def.max_wraps, def.grammar)
end

"""
Converts from a vector of codons to a grammar derivation.
"""
convert!(
  r_from::IntVectorRepresentation,
  r_to::GrammarDerivationRepresentation,
  c_from::IndividualStage{Vector{Int}},
  c_to::IndividualStage{AbstractString}
) =
  c_to = derive(r_to, c_from.value)

"""
Produces a grammar derivation from a given sequence of codons, according to
a provided grammar derivation representation.
"""
derive(r::GrammarDerivationRepresentation, s::Vector{Int}) =
  derive(r.grammar, s, r.max_wraps)


text:
using Cclib
using Test

@testset "Cclib.jl" begin
    @testset "API surface" begin
        @test Cclib.__version__ == "1.7.1"

        @testset "Checking the exported symbols are available" begin
            for s in [:bridge, :io, :method, :parser]
                @test @isdefined s
            end

            for s in [:ccopen, :ccwrite]
                @test @isdefined s
            end
        end

        @testset "Accessing available parsers" begin
            for p in [:ADF, :DALTON, :FChk, :GAMESS, :GAMESSUK, :Gaussian, :Jaguar,
                :Molcas, :Molpro, :MOPAC, :NWChem, :ORCA, :Psi3, :Psi4,
                :QChem, :Turbomole]
                @test hasproperty(parser, p)
            end
        end

        @testset "Accessing available methods" begin
            for m in [:Bader, :Bickelhaupt, :CDA, :CSPA, :DDEC6, :Density,
                :Electrons, :FragmentAnalysis, :Hirshfeld, :LPA, :MBO,
                :Moments, :MPA, :Nuclear, :OPA, :Orbitals, :Stockholder,
                :Volume]
                @test hasproperty(method, m)
            end
        end

        @testset "Accessing available IO utilities" begin
            for u in [:CJSONReader, :CJSONWriter, :CML, :MOLDEN, :WFXWriter,
                :XYZReader, :XYZWriter]
                @test hasproperty(io, u)
            end

            for u in [:ccframe, :ccopen, :ccread, :ccwrite, :URL_PATTERN]
                @test hasproperty(io, u)
            end
        end
    end

    @testset "Usage" begin
        # An example log file.
        filename = "water.out"

        @testset "Parsing a log file with ccopen and checking some attributes" begin
            parser = ccopen(filename)
            data = parser.parse()

            for a in [:atomcharges, :atomcoords, :atommasses, :atomnos, :charge,
                :coreelectrons, :enthalpy, :entropy, :freeenergy, :geotargets,
                :geovalues, :grads, :homos, :metadata, :moenergies, :moments,
                :mult, :natom, :nbasis, :nmo, :optdone, :pressure, :scfenergies,
                :scftargets, :scfvalues, :temperature, :vibdisps, :vibfreqs,
                :vibirs, :zpve]
                @test hasproperty(data, a)
            end
        end

        @testset "Parsing a log file with ccread and checking some values" begin
            data = io.ccread(filename)

            @test data.natom == 3
            @test data.nmo == 7
            @test data.atomcoords[1, :, :] == [0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
        end
    end
end


text:
module main
using CxxWrap
@wrapmodule(joinpath(@__DIR__, "libs", "libexlib.so"))

function __init__()
    @initcxx
end

function hello()
    println("Hello from julia~")
    println(greet())
end

mutable struct Player
    HP::Int32
    MP::Int32
end

player = Player(10, 22)

function get()::Int32
    println("Hello from julia")
    player.HP += 1
    return player.HP

end

end

max_stars_repo_name: gluaxspeed/embedded-tests

max_issues_repo_name: gluaxspeed/embedded-tests

max_forks_repo_name: gluaxspeed/embedded-tests


text:
using SimpleDifferentialOperators
using Test, LinearAlgebra, DualNumbers

@testset "Operators under reflecting barrier conditions" begin

    # helper function to get differential operators in a more efficient way
    diffusionoperators(x̄, bc) = (L₁₋bc = L₁₋bc(x̄, bc), L₁₊bc = L₁₊bc(x̄, bc), L₂bc = L₂bc(x̄, bc), x = interiornodes(x̄, bc))

    @testset "Accuracy & regression test" begin
        # Uniform grid
        x̄ = 0:4
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Reflecting(), Reflecting()))
        @test @inferred(diffusionoperators(x̄, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = 0.), Mixed(ξ = 0.)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x) # test that the mixed case properly nests the reflecting case
        @test Array(L₁₋bc) == [0. 0. 0.; -1. 1. 0.; 0. -1. 1.]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1. 1.; 0. 0. 0.]
        @test Array(L₂bc) == [-1. 1. 0.; 1. -2. 1.; 0. 1. -1.]
        @test Array(x) == [1; 2; 3]

        # Irregular grid (with Δ_1m = Δ_1p and Δ_Mp = Δ_Mm)
        x̄ = [0.0; 1.0; 2.0; 4.0; 6.0]
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Reflecting(), Reflecting()))
        @test @inferred(diffusionoperators(x̄,(Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = 0.), Mixed(ξ = 0.)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x) # test that the mixed case properly nests the reflecting case
        @test Array(L₁₋bc) == [0. 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1/2 1/2; 0. 0. 0.]
        @test Array(L₂bc) == [-1. 1. 0.; 2/((2.0+1.0)*1.0) -2/(2.0*1.0) 2/((2.0+1.0)*2.0); 0. 1/4 -1/4]
        @test Array(x) == [1.0; 2.0; 4.0]

        # Irregular grid (with Δ_1m != Δ_1p and Δ_Mp != Δ_Mm)
        x̄ = [-1.0; 1.0; 2.0; 4.0; 7.0]

        Δ_1m = x̄[2] - x̄[1]
        Δ_1p = x̄[3] - x̄[2]
        Δ_Mm = x̄[end-1] - x̄[end-2]
        Δ_Mp = x̄[end] - x̄[end-1]
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Reflecting(), Reflecting()))
        @test @inferred(diffusionoperators(x̄,(Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test @inferred(diffusionoperators(x̄,(Mixed(ξ = 0.), Mixed(ξ = 0.)))) == diffusionoperators(x̄,(Reflecting(), Reflecting())) # test that the mixed case properly nests the reflecting case
        @test Array(L₁₋bc) == [0. 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1/2 1/2; 0. 0. 0.]
        Ξ_1 = -2*(1/(Δ_1m*Δ_1p)-1/((Δ_1p+Δ_1m)*(Δ_1m)))
        Ξ_M = -2*(1/(Δ_Mm*Δ_Mp)-1/((Δ_Mp+Δ_Mm)*(Δ_Mp)))

        @test Array(L₂bc) == [Ξ_1 2/(Δ_1p*(Δ_1m+Δ_1p)) 0.;
                            2/((2.0+1.0)*1.0) -2/(2.0*1.0) 2/((2.0+1.0)*2.0);
                            0. 2/(Δ_Mm*(Δ_Mm+Δ_Mp)) Ξ_M]
        @test Array(x) == [1.0; 2.0; 4.0]
    end

    @testset "Consistency" begin
        # Test for consistency
        uniformGrid = range(0.0, 1.0, length = 500)
        irregularGrid = collect(uniformGrid)
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(uniformGrid,(Reflecting(), Reflecting()))
        L₁₋bc_ir, L₁₊_ir, L₂_ir, x_ir = diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))
        @test L₁₋bc ≈ L₁₋bc_ir
        @test L₁₊bc ≈ L₁₊_ir
        @test L₂bc ≈ L₂_ir
    end
end

@testset "Operators under mixed barrier conditions" begin

    # helper function to get differential operators in a more efficient way
    diffusionoperators(x̄, bc) = (L₁₋bc = L₁₋bc(x̄, bc), L₁₊bc = L₁₊bc(x̄, bc), L₂bc = L₂bc(x̄, bc), x = interiornodes(x̄, bc))

    @testset "Accuracy & regression test" begin
        # Uniform grid
        x̄ = 0:4
        ξ_lb, ξ_ub = (-1., -2.)

        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test Array(L₁₋bc) == [1+1/(-1+ξ_lb) 0. 0.; -1. 1. 0.; 0. -1. 1.]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1. 1.; 0. 0. -1+1/(1+ξ_ub)]
        @test Array(L₂bc) == [-2-1/(-1+ξ_lb) 1. 0.; 1. -2. 1.; 0. 1. -2+1/(1+ξ_ub)]
        @test Array(x) == [1; 2; 3]

        # Irregular grid (with Δ_1m = Δ_1p and Δ_Mp = Δ_Mm)
        x̄ = [0.0; 1.0; 2.0; 4.0; 6.0]
        ξ_lb, ξ_ub = (-1., -2.)
        Δ_Mm = x̄[end-1] - x̄[end-2]
        Δ_Mp = x̄[end] - x̄[end-1]

        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test Array(L₁₋bc) == [1+1/(-1+ξ_lb) 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1/2 1/2; 0. 0. (-1+1/(1+ξ_ub*Δ_Mp))/Δ_Mp]
        @test Array(L₂bc) == [-1.5 1. 0.; 2/((2.0+1.0)*1.0) -2/(2.0*1.0) 2/((2.0+1.0)*2.0); 0. 1/4 -2*(1/(Δ_Mm*Δ_Mp)-1/((1+ξ_ub*Δ_Mp)*(Δ_Mp+Δ_Mm)*(Δ_Mp)))]
        @test Array(x) == [1.0; 2.0; 4.0]

        # Irregular grid (with Δ_1m != Δ_1p and Δ_Mp != Δ_Mm)
        x̄ = [-1.0; 1.0; 2.0; 4.0; 7.0]
        ξ_lb, ξ_ub = (-1., -2.)
        Δ_1m = x̄[2] - x̄[1]
        Δ_1p = x̄[3] - x̄[2]
        Δ_Mm = x̄[end-1] - x̄[end-2]
        Δ_Mp = x̄[end] - x̄[end-1]

        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test Array(L₁₋bc) ≈ [1/3 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) ≈ [-1. 1. 0.; 0. -1/2 1/2; 0. 0. -0.4]
        @test Array(L₂bc) ≈ [-8/9 2/3 0.; 2/3 -1. 1/3; 0. 0.2 -0.36]
        @test Array(x) ≈ [1.0; 2.0; 4.0]
    end

    @testset "Accuracy & regression test, for different directions" begin
        # Uniform grid
        x̄ = 0:4
        ξ_lb, ξ_ub = (-1., -2.)

        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub, direction = :forward)))
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub, direction = :forward)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test Array(L₁₋bc) == [1.0 0. 0.; -1. 1. 0.; 0. -1. 1.]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1. 1.; 0. 0. 2.0]
        @test Array(L₂bc) == [-2.0 1. 0.; 1. -2. 1.; 0. 1. 1.0]
        @test Array(x) == [1; 2; 3]

        # Irregular grid (with Δ_1m = Δ_1p and Δ_Mp = Δ_Mm)
        x̄ = [0.0; 1.0; 2.0; 4.0; 6.0]
        ξ_lb, ξ_ub = (-1., -2.)
        Δ_Mm = x̄[end-1] - x̄[end-2]
        Δ_Mp = x̄[end] - x̄[end-1]

        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub, direction = :forward)))
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub, direction = :forward)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test Array(L₁₋bc) == [1.0 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1/2 1/2; 0. 0. 2.0]
        @test Array(L₂bc) == [-2.0 1. 0.; 2/((2.0+1.0)*1.0) -2/(2.0*1.0) 2/((2.0+1.0)*2.0); 0. 1/4 3/4]
        @test Array(x) == [1.0; 2.0; 4.0]

        # Irregular grid (with Δ_1m != Δ_1p and Δ_Mp != Δ_Mm)
        x̄ = [-1.0; 1.0; 2.0; 4.0; 7.0]
        ξ_lb, ξ_ub = (-1., -2.)
        Δ_1m = x̄[2] - x̄[1]
        Δ_1p = x̄[3] - x̄[2]
        Δ_Mm = x̄[end-1] - x̄[end-2]
        Δ_Mp = x̄[end] - x̄[end-1]

        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub, direction = :forward)))
        @test @inferred(diffusionoperators(x̄, (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub, direction = :forward)))) == (L₁₋bc = L₁₋bc, L₁₊bc = L₁₊bc, L₂bc = L₂bc, x = x)
        @test Array(L₁₋bc) ≈ [1.0 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) ≈ [-1. 1. 0.; 0. -1/2 1/2; 0. 0. 2.0]
        @test Array(L₂bc) ≈ [-4/3 2/3 0.; 2/3 -1. 1/3; 0. 0.2 0.6]
        @test Array(x) ≈ [1.0; 2.0; 4.0]
    end

    @testset "Consistency" begin
        ξ_lb, ξ_ub = (1., 2.)
        uniformGrid = range(0.0, 1.0, length = 500)
        irregularGrid = collect(uniformGrid)
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(uniformGrid, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))
        L₁₋bc_ir, L₁₊_ir, L₂_ir, x_ir = diffusionoperators(irregularGrid, (Mixed(ξ = ξ_lb), Mixed(ξ = ξ_ub)))
        @test L₁₋bc ≈ L₁₋bc_ir
        @test L₁₊bc ≈ L₁₊_ir
        @test L₂bc ≈ L₂_ir
        @test x ≈ x_ir
    end
end

@testset "Operators under homogenous absorbing boundary conditions" begin
    # helper function to get differential operators in a more efficient way
    diffusionoperators(x̄, bc) = (L₁₋bc = L₁₋bc(x̄, bc), L₁₊bc = L₁₊bc(x̄, bc), L₂bc = L₂bc(x̄, bc), x = interiornodes(x̄, bc))
    diffusionoperators_without_bc(x̄) = (L₁₋ = L₁₋(x̄), L₁₊ = L₁₊(x̄), L₂ = L₂(x̄), x = interiornodes(x̄))

    @testset "Accuracy & regression test" begin
        # Uniform grid
        x̄ = 0:4
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Absorbing(), Absorbing()))
        L₁₋, L₁₊, L₂, x = diffusionoperators_without_bc(x̄)

        # under absorbing barrier conditions, the corresponding
        # operators on the interiors are identical as interior of extended operators
        @test L₁₋bc == L₁₋[:,2:end-1]
        @test L₁₊bc == L₁₊[:,2:end-1]
        @test L₂bc == L₂[:,2:end-1]

        # check accuracy as regression tests
        @test Array(L₁₋bc) == [1. 0. 0.; -1. 1. 0.; 0. -1. 1.]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1. 1.; 0. 0. -1.]
        @test Array(L₂bc) == [-2. 1. 0.; 1. -2. 1.; 0. 1. -2.]
        @test Array(x) == [1; 2; 3]

        # Reflecting on lb, absorbing on ub
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Reflecting(), Absorbing()))
        L₁₋, L₁₊, L₂, x = diffusionoperators_without_bc(x̄)

        # check accuracy as regression tests
        @test Array(L₁₋bc) == [0. 0. 0.; -1. 1. 0.; 0. -1. 1.]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1. 1.; 0. 0. -1.]
        @test Array(L₂bc) == [-1. 1. 0.; 1. -2. 1.; 0. 1. -2.]
        @test Array(x) == [1; 2; 3]

        # Absorbing on lb, reflecting on ub
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Absorbing(), Reflecting()))
        L₁₋, L₁₊, L₂, x = diffusionoperators_without_bc(x̄)

        # check accuracy as regression tests
        @test Array(L₁₋bc) == [1. 0. 0.; -1. 1. 0.; 0. -1. 1.]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1. 1.; 0. 0. 0.]
        @test Array(L₂bc) == [-2. 1. 0.; 1. -2. 1.; 0. 1. -1.]
        @test Array(x) == [1; 2; 3]

        # Irregular grid (with Δ_1m = Δ_1p and Δ_Mp = Δ_Mm)
        x̄ = [0.0; 1.0; 2.0; 4.0; 6.0]
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Absorbing(), Absorbing()))
        L₁₋, L₁₊, L₂, x = diffusionoperators_without_bc(x̄)

        # under absorbing barrier conditions, the corresponding
        # operators on the interiors are identical as interior of extended operators
        @test L₁₋bc == L₁₋[:,2:end-1]
        @test L₁₊bc == L₁₊[:,2:end-1]
        @test L₂bc == L₂[:,2:end-1]

        # check accuracy as regression tests
        @test Array(L₁₋bc) == [1. 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1/2 1/2; 0. 0. -1/2]
        @test Array(L₂bc) == [-2. 1. 0.; 2/((2.0+1.0)*1.0) -2/(2.0*1.0) 2/((2.0+1.0)*2.0); 0. 1/4 -1/2]
        @test Array(x) == [1.0; 2.0; 4.0]

        # Irregular grid (with Δ_1m != Δ_1p and Δ_Mp != Δ_Mm)
        x̄ = [-1.0; 1.0; 2.0; 4.0; 7.0]

        Δ_1m = x̄[2] - x̄[1]
        Δ_1p = x̄[3] - x̄[2]
        Δ_Mm = x̄[end-1] - x̄[end-2]
        Δ_Mp = x̄[end] - x̄[end-1]
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(x̄, (Absorbing(), Absorbing()))
        L₁₋, L₁₊, L₂, x = diffusionoperators_without_bc(x̄)

        # under absorbing barrier conditions, the corresponding
        # operators on the interiors are identical as interior of extended operators
        @test L₁₋bc == L₁₋[:,2:end-1]
        @test L₁₊bc == L₁₊[:,2:end-1]
        @test L₂bc == L₂[:,2:end-1]

        # check accuracy as regression tests
        @test Array(L₁₋bc) == [1/Δ_1m 0. 0.; -1. 1. 0.; 0. -1/2 1/2]
        @test Array(L₁₊bc) == [-1. 1. 0.; 0. -1/2 1/2; 0. 0. -1/Δ_Mp]
        Ξ_1 = -2*(1/(Δ_1m*Δ_1p)-1/((Δ_1p+Δ_1m)*(Δ_1m)))
        Ξ_M = -2*(1/(Δ_Mm*Δ_Mp)-1/((Δ_Mp+Δ_Mm)*(Δ_Mp)))

        @test Array(L₂bc) == [-2/(Δ_1p*Δ_1m) 2/(Δ_1p*(Δ_1m+Δ_1p)) 0.;
                            2/((2.0+1.0)*1.0) -2/(2.0*1.0) 2/((2.0+1.0)*2.0);
                            0. 2/(Δ_Mm*(Δ_Mm+Δ_Mp)) -2/(Δ_Mm*Δ_Mp) ]
        @test Array(x) == [1.0; 2.0; 4.0]
    end

    @testset "Consistency" begin
        # Test for consistency
        uniformGrid = range(0.0, 1.0, length = 500)
        irregularGrid = collect(uniformGrid)
        L₁₋bc, L₁₊bc, L₂bc, x = diffusionoperators(uniformGrid,(Absorbing(), Absorbing()))
        L₁₋bc_ir, L₁₊_ir, L₂_ir, x_ir = diffusionoperators(irregularGrid, (Absorbing(), Absorbing()))
        @test L₁₋bc ≈ L₁₋bc_ir
        @test L₁₊bc ≈ L₁₊_ir
        @test L₂bc ≈ L₂_ir
    end
end

@testset "Operators under upwind schemes" begin
    f(x) = x^2
    μ(x) = -x # drift depends on state
    σ = 0.1
    ρ = 0.05
    M = 100 # size of grid
    x̄ = range(-1.0, 1.0, length = M+2) # grid
    x = interiornodes(x̄)
    bc = (Reflecting(), Reflecting()) # specify BC (reflecting barrier)
    ## M-vector of drifts stacked according to the states
    μs = μ.(x)

    # Define first order differential operator using upwind scheme
    L_1_upwind = (μs .<= 0) .* L₁₋bc(x̄, bc) + (μs .> 0) .* L₁₊bc(x̄, bc)

    indices_m = findall(μs .<= 0)
    indices_p = findall(μs .> 0)

    @test L_1_upwind[indices_m, indices_m] == Array(L₁₋bc(x̄, bc)[indices_m, indices_m])
    @test L_1_upwind[indices_p, indices_p] == Array(L₁₊bc(x̄, bc)[indices_p, indices_p])
end

@testset "Input Type Variance" begin
    # helper function to get differential operators in a more efficient way
    diffusionoperators(x̄, bc) = (L₁₋bc = L₁₋bc(x̄, bc), L₁₊bc = L₁₊bc(x̄, bc), L₂bc = L₂bc(x̄, bc))

    # BigFloat
    uniformGrid = range(BigFloat(0.0), BigFloat(1.0), length = 100)
    irregularGrid = collect(uniformGrid)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(uniformGrid, (Reflecting(), Reflecting()))
    @test @inferred(diffusionoperators(uniformGrid, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))
    @test @inferred(diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(uniformGrid, (Mixed(ξ = one(BigFloat)), Mixed(ξ = one(BigFloat))))
    @test @inferred(diffusionoperators(uniformGrid, (Mixed(ξ = one(BigFloat)), Mixed(ξ = one(BigFloat))))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(irregularGrid, (Mixed(ξ = one(BigFloat)), Mixed(ξ = one(BigFloat))))
    @test @inferred(diffusionoperators(irregularGrid, (Mixed(ξ = one(BigFloat)), Mixed(ξ = one(BigFloat))))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)

    # Float32
    uniformGrid = range(Float32(0.0), Float32(1.0), length = 100)
    irregularGrid = collect(uniformGrid)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(uniformGrid, (Reflecting(), Reflecting()))
    @test @inferred(diffusionoperators(uniformGrid, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))
    @test @inferred(diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(uniformGrid, (Mixed(ξ = one(Float32)), Mixed(ξ = one(Float32))))
    @test @inferred(diffusionoperators(uniformGrid, (Mixed(ξ = one(Float32)), Mixed(ξ = one(Float32))))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc = L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(irregularGrid, (Mixed(ξ = one(Float32)), Mixed(ξ = one(Float32))))
    @test @inferred(diffusionoperators(irregularGrid, (Mixed(ξ = one(Float32)), Mixed(ξ = one(Float32))))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc = L₁₊_cache, L₂bc = L₂_cache)

    # Duals
    uniformGrid = range(Dual(0.0), Dual(1.0), length = 100)
    irregularGrid = collect(uniformGrid)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(uniformGrid, (Reflecting(), Reflecting()))
    @test @inferred(diffusionoperators(uniformGrid, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))
    @test @inferred(diffusionoperators(irregularGrid, (Reflecting(), Reflecting()))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(uniformGrid, (Mixed(ξ = Dual(1.0)), Mixed(ξ = Dual(1.0))))
    @test @inferred(diffusionoperators(uniformGrid, (Mixed(ξ = Dual(1.0)), Mixed(ξ = Dual(1.0))))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
    L₁₋bc_cache, L₁₊_cache, L₂_cache = diffusionoperators(irregularGrid, (Mixed(ξ = Dual(1.0)), Mixed(ξ = Dual(1.0))))
    @test @inferred(diffusionoperators(irregularGrid, (Mixed(ξ = Dual(1.0)), Mixed(ξ = Dual(1.0))))) == (L₁₋bc = L₁₋bc_cache, L₁₊bc =  L₁₊_cache, L₂bc = L₂_cache)
end


text:
module Fst

include("semiring.jl")
include("wfst.jl")
include("io.jl")
include("compose.jl")
include("topological_sort.jl")
include("trim.jl")

end


text:
# Statistical model for analyzing Radial Velocity observations
# Assumes existance of data in arrays times, obs, sigma_obs
# Assumes existance of function calc_model_rv(theta, time)

export ploglikelihood, plogprior, plogtarget
export set_times, set_obs, set_sigma_obs

# Observational data to compare model to
global times # ::Array{Float64,1}
global obs # ::Array{Float64,1}
global sigma_obs #::Array{Float64,1}

# Functions to set global data within module
function set_times(t::Array{Float64,1}) global times = t   end
function set_obs(o::Array{Float64,1})  global obs = o   end
function set_sigma_obs(so::Array{Float64,1}) global sigma_obs = so  end

function ploglikelihood(p::Vector)
  num_pl = num_planets(p)
  @assert num_pl >= 1
  if !is_valid(p) return -Inf end  # prempt model evaluation
  # Set t, o, and so to point to global arrays with observational data, while enforcing types
  t::Array{Float64,1} = times
  o::Array{Float64,1} = obs
  so::Array{Float64,1} = sigma_obs
  @assert length(times) == length(obs) == length(sigma_obs)
  jitter_sq = RvModelKeplerian.num_jitters >=1 ? extract_jitter(p)^2 : 0.
  chisq = zero(eltype(p))
  log_normalization = -0.5*length(t)*log(2pi)
  for i in 1:length(t)
    model::eltype(p) = calc_model_rv(p,t[i])
    sigma_eff_sq = so[i]^2+jitter_sq
    chisq += abs2(model-o[i])/sigma_eff_sq
    log_normalization -= 0.5*log(sigma_eff_sq)
  end
  return -0.5*(chisq)+log_normalization
end

function plogprior(p::Vector)
  num_pl = num_planets(p)
  @assert num_pl >= 1
  if !is_valid(p) return -Inf end  # prempt model evaluation
  logp = zero(eltype(p))
  logp -= 2*log(2pi)
  const max_period = 10000.0
  const max_amplitude = 10000.0
  for plid in 1:num_pl
    P::eltype(p) = extract_period(p,plid=plid)
    K::eltype(p) = extract_amplitude(p,plid=plid)
    logp += -log((1+P/P0::Float64)*log1p(max_period/P0::Float64)*
                 (1+K/K0::Float64)*log1p(max_amplitude/K0::Float64) )

  end
  if RvModelKeplerian.num_jitters >=1
     const max_jitter = 10000.0
     jitter::eltype(p) = extract_jitter(p)
     logp += -log((1+jitter/Jitter0::Float64)*log1p(max_jitter/Jitter0::Float64))
  end

  return logp::eltype(p)
end

plogtarget(p::Vector) = ploglikelihood(p) + plogprior(p)


text:
module CurrencyData

include("currencies.jl")
include("locale.jl")
include("symbols.jl")

export ISO4217, SHORT_SYMBOL, LONG_SYMBOL, LOCAL_SYMBOL_LOCATION

end  # module CurrencyData


text:
abstract ModelInterface{T}

description(m::ModelInterface) = "$(typeof(m)): $(m.spec), $(m.subspec)"
function Base.show{T<:ModelInterface}(io::IO, m::T)
    @printf io "Model\n"
    @printf io "%s\n" T
    @printf io "description:\n %s\n"          description(m)
    @printf io "data vintage:           %s\n" data_vintage(m)
end

# Interface for I/O settings
spec(m::ModelInterface)         = m.spec
subspec(m::ModelInterface)      = m.subspec
saveroot(m::ModelInterface)     = get_setting(m, :saveroot)
dataroot(m::ModelInterface)     = get_setting(m, :dataroot)

# Interface for data
data_vintage(m::ModelInterface) = get_setting(m, :data_vintage)

# Interface for general computation settings
use_parallel_workers(m::ModelInterface)    = get_setting(m, :use_parallel_workers)

#=
Build paths to where input/output/results data are stored.

Description:
Creates the proper directory structure for input and output files, treating the DSGE/save
    directory as the root of a savepath directory subtree. Specifically, the following
    structure is implemented:

    dataroot/

    savepathroot/
                 output_data/<spec>/<subspec>/log/
                 output_data/<spec>/<subspec>/<out_type>/raw/
                 output_data/<spec>/<subspec>/<out_type>/work/
                 output_data/<spec>/<subspec>/<out_type>/tables/
                 output_data/<spec>/<subspec>/<out_type>/figures/

Note: we refer to the savepathroot/output_data/<spec>/<subspec>/ directory as saveroot.
=#
"""
```
logpath(model)
```
Returns path to log file. Path built as
```
<output root>/output_data/<spec>/<subspec>/log/log_<filestring>.log
```
"""
function logpath(m::ModelInterface)
    return savepath(m, "log", "log.log")
end

strs = [:work, :raw, :tables, :figures]
fns = [symbol(x, "path") for x in strs]
for (str, fn) in zip(strs, fns)
    @eval begin
        # First eval function
        function $fn{T<:AbstractString}(m::ModelInterface,
                                             out_type::T,
                                             file_name::T="",
                                             filestring_addl::Vector{T}=Vector{T}())
            return savepath(m, out_type, $(string(str)), file_name, filestring_addl)
        end

        # Then, add docstring to it
        @doc $(
        """
        ```
        $fn{T<:AbstractString}(m::ModelInterface, out_type::T, file_name::T="")
        ```

        Returns path to specific $str output file, creating containing directory as needed. If
        `file_name` not specified, creates and returns path to containing directory only. Path built
        as
        ```
        <output root>/output_data/<spec>/<subspec>/<out_type>/$str/<file_name>_<filestring>.<ext>
        ```
        """
        ) $fn
    end
end

# Not exposed to user. Actually create path and insert model string to file name.
function savepath{T<:AbstractString}(m::ModelInterface,
                                     out_type::T,
                                     sub_type::T,
                                     file_name::T="",
                                     filestring_addl::Vector{T}=Vector{T}())
    # Containing dir
    path = joinpath(saveroot(m), "output_data", spec(m), subspec(m), out_type, sub_type)
    if !isdir(path)
        mkpath(path)
    end

    # File with model string inserted
    if !isempty(file_name)
        if isempty(filestring_addl)
            myfilestring = filestring(m)
        else
            myfilestring = filestring(m, filestring_addl)
        end
        (base, ext) = splitext(file_name)
        file_name_detail = base * myfilestring * ext
        path = joinpath(path, file_name_detail)
    end

    return path
end

# Input data handled slightly differently, because it is not model-specific.
"""
```
inpath{T<:AbstractString}(m::ModelInterface, in_type::T, file_name::T="")
```

Returns path to specific input data file, creating containing directory as needed. If
`file_name` not specified, creates and returns path to containing directory only. Valid
`in_type` includes:

* `"data"`: recorded data
* `"cond"`: conditional data - nowcasts for the current forecast quarter, or related
* `"user"`: user-supplied data for starting parameter vector, hessian, or related

Path built as
```
<data root>/<in_type>/<file_name>
```
"""
function inpath{T<:AbstractString}(m::ModelInterface, in_type::T, file_name::T="")
    path = dataroot(m)
    # Normal cases.
    if in_type == "data" || in_type == "cond"
        path = joinpath(path, in_type)
    # User-provided inputs. May treat this differently in the future.
    elseif in_type == "user"
        path = joinpath(path, "user")
    else
        error("Invalid in_type: ", in_type)
    end

    # Containing dir
    if !isdir(path)
        mkpath(path)
    end

    # If file_name provided, return full path
    if !isempty(file_name)
        path = joinpath(path, file_name)
    end

    return path
end

filestring(m::ModelInterface) = filestring(m, Vector{AbstractString}())
filestring(m::ModelInterface, d::AbstractString) = filestring(m, [d])
function filestring{T<:AbstractString}(m::ModelInterface,
                                        d::Vector{T})
    if !m.testing
        filestrings = Vector{T}()
        for (skey, sval) in m.settings
            if sval.print
                push!(filestrings, to_filestring(sval))
            end
        end
        append!(filestrings, d)
        sort!(filestrings)
        return "_"*join(filestrings, "_")
    else
        return "_test"
    end
end


text:
# using Distributions

immutable GaussianMixtureModelCommonCovarDiagonal <: Distribution
    mu::Array{Float64,2}
    probs::Vector{Float64}
    covar::Vector{Float64}
    aliastable::Distributions.AliasTable
    function GaussianMixtureModelCommonCovarDiagonal(m::Array{Float64,2}, p::Vector{Float64}, ic::Vector{Float64})
        if size(m,2) != length(p)
            error("means and probs must have the same number of elements")
        end
if( size(m,1) != size(ic,1) )
    error("means and covar matrix not compatible sizes: ",size(m)," vs ",size(ic) )
end
        sump = 0.0
        for i in 1:length(p)
            if p[i] < 0.0
                error("MixtureModel: probabilities must be non-negative")
            end
            sump += p[i]
        end
        table = Distributions.AliasTable(p)
        new(m, p ./ sump, ic, table)
    end
end

function mean(d::GaussianMixtureModelCommonCovarDiagonal)
    np = size(d.mu,2)
    m = zeros(np)
    for i in 1:length(d.probs)
        m += vec(d.mu[:,i]) .* d.probs[i]
    end
    return m
end

function pdf(d::GaussianMixtureModelCommonCovarDiagonal, x::Any)
    p = 0.0
    for i in 1:length(d.probs)
        p += Distributions.pdf(Distributions.MvNormal(sqrt(d.covar)), x .- vec(d.mu[:,i]) ) * d.probs[i]
    end
    return p
end

function rand(d::GaussianMixtureModelCommonCovarDiagonal)
    i = Distributions.rand(d.aliastable)
    return  Distributions.rand(Distributions.MvNormal(vec(d.mu[:,i]),sqrt(d.covar)))
end


text:
export FashionMNIST

"""
Fashion-MNIST

- Authors: Han Xiao, Kashif Rasul, Roland Vollgraf
- Website: https://github.com/zalandoresearch/fashion-mnist

Fashion-MNIST is a dataset of Zalando's article images—consisting
of a training set of 60,000 examples and a test set of 10,000
examples. Each example is a 28x28 grayscale image, associated
with a label from 10 classes. It can serve as a drop-in
replacement for MNIST.

# Interface

- [`FashionMNIST.traintensor`](@ref), [`FashionMNIST.trainlabels`](@ref), [`FashionMNIST.traindata`](@ref)
- [`FashionMNIST.testtensor`](@ref), [`FashionMNIST.testlabels`](@ref), [`FashionMNIST.testdata`](@ref)

# Utilities

- [`FashionMNIST.download`](@ref)

Also, the `FashionMNIST` module is re-exporting [`convert2image`](@ref MNIST.convert2image) from the [`MNIST`](@ref) module.
"""
module FashionMNIST
    using DataDeps
    using FixedPointNumbers
    using ..MLDatasets: bytes_to_type, datafile, download_dep, download_docstring
    import ..MNIST
    using ..MNIST: convert2image
    using ..MNIST.Reader

    export

        classnames,

        traintensor,
        testtensor,

        trainlabels,
        testlabels,

        traindata,
        testdata,

        convert2image,

        download

    @deprecate convert2features reshape

    const DEPNAME = "FashionMNIST"
    const TRAINIMAGES = "train-images-idx3-ubyte.gz"
    const TRAINLABELS = "train-labels-idx1-ubyte.gz"
    const TESTIMAGES  = "t10k-images-idx3-ubyte.gz"
    const TESTLABELS  = "t10k-labels-idx1-ubyte.gz"

    const CLASSES = [
        "T-Shirt",
        "Trouser",
        "Pullover",
        "Dress",
        "Coat",
        "Sandal",
        "Shirt",
        "Sneaker",
        "Bag",
        "Ankle boot"
    ]

    """
        download([dir]; [i_accept_the_terms_of_use])

    Trigger the (interactive) download of the full dataset into
    "`dir`". If no `dir` is provided the dataset will be
    downloaded into "~/.julia/datadeps/$DEPNAME".

    This function will display an interactive dialog unless
    either the keyword parameter `i_accept_the_terms_of_use` or
    the environment variable `DATADEPS_ALWAYS_ACCEPT` is set to
    `true`. Note that using the data responsibly and respecting
    copyright/terms-of-use remains your responsibility.
    """
    download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

    include("interface.jl")

    function __init__()
        register(DataDep(
            DEPNAME,
            """
            Dataset: Fashion-MNIST
            Authors: Han Xiao, Kashif Rasul, Roland Vollgraf
            Website: https://github.com/zalandoresearch/fashion-mnist
            License: MIT

            [Han Xiao et al. 2017]
                Han Xiao, Kashif Rasul, and Roland Vollgraf.
                "Fashion-MNIST: a Novel Image Dataset for Benchmarking Machine Learning Algorithms."
                arXiv:1708.07747

            The files are available for download at the offical
            website linked above. Note that using the data
            responsibly and respecting copyright remains your
            responsibility.
            """,
            "http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/" .* [TRAINIMAGES, TRAINLABELS, TESTIMAGES, TESTLABELS],
            "c916b6e00d3083643332b70f3c5c3543d3941334b802e252976893969ee6af67",
        ))
    end
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"

set -ex
${JULIA} -e 'using Pkg; Pkg.instantiate();
             Pkg.develop(PackageSpec(path=pwd()))'
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

using Documenter, Transducers

makedocs(;
    modules=[Transducers],
    pages=[
        "Home" => "index.md",
        "Manual" => "manual.md",
        "Interface" => "interface.md",
        hide("Internals" => "internals.md"),
    ],
    repo="https://github.com/tkf/Transducers.jl/blob/{commit}{path}#L{line}",
    sitename="Transducers.jl",
    authors="Takafumi Arakaki",
    assets=[],
    strict=true,
)

deploydocs(;
    repo="github.com/tkf/Transducers.jl",
)


text:
# module TestIO

using Test, Dates, Temporal

@testset "Input/Output" begin
    @testset "Sample Data" begin
        filepath = split(pwd(), "/")[end] == "test" ? "$(pwd())/../data/Corn.csv" : "$(pwd())/data/Corn.csv"
        corn = tsread(filepath, header=true, eol='\n')
        @test size(corn, 2) == 8
        @test size(cl(corn), 2) == 1
    end
    # @testset "Web Downloads" begin
    #     crude = quandl("CHRIS/CME_CL1", from="2010-01-01", thru=string(today()))
    #     @test size(crude, 2) == 8
    #     @test op(crude).fields == [:Open]
    # end
end

# end


text:
include("test_highlevel.jl")
include("test_capi.jl")

CUDA_INSTALLED = false
CUDA_FUNCTIONAL = false
try
    import CUDA
    global CUDA_INSTALLED = true
    global CUDA_FUNCTIONAL = CUDA.functional()
catch
end
if CUDA_INSTALLED && CUDA_FUNCTIONAL
    @info """
    Found a working CUDA.jl package, running CUDA tests
    """
    include("test_cuda.jl")
else
    msg = """
    Skipping CUDA tests. Got
    CUDA_INSTALLED = $CUDA_INSTALLED
    CUDA_FUNCTIONAL = $CUDA_FUNCTIONAL
    """
    @warn msg
end


text:
# InflationTotalCPI - Implementación para obtener la medida estándar de ritmo
# inflacionario a través de la variación interanual del IPC

struct InflationTotalCPI <: InflationFunction
end

# Extender el método para obtener el nombre de esta medida
measure_name(::InflationTotalCPI) = "Variación interanual IPC"
measure_tag(::InflationTotalCPI) = "Total"

# Las funciones sobre VarCPIBase deben resumir en variaciones intermensuales

# Método para objetos VarCPIBase cuyo índice base es un escalar
function (inflfn::InflationTotalCPI)(base::VarCPIBase{T, T}) where {T <: AbstractFloat}
    base_ipc = capitalize(base.v, base.baseindex)
    ipc = base_ipc * base.w / base.baseindex
    varinterm!(ipc, ipc, 100)
    ipc
end

# Esta medida sí se comporta diferente de acuerdo a los índices base, por lo que
# se define una versión que toma en cuenta los diferentes índices. Si la medida
# solamente genera resumen de las variaciones intermensuales, no es necesario.
# Método para objetos VarCPIBase cuyos índices base son un vector
function (inflfn::InflationTotalCPI)(base::VarCPIBase{T, B}) where {T <: AbstractFloat, B <: AbstractVector{T}}
    base_ipc = capitalize(base.v, base.baseindex)
    # Obtener índice base y normalizar a 100
    baseindex = base.baseindex' * base.w
    ipc = 100 * (base_ipc * base.w / baseindex)
    varinterm!(ipc, ipc, 100)
    ipc
end


text:
"""

    Interval unions sets of defined by unions of disjoint intervals.
    This file includes constructors, arithmetic (including intervals and scalars)
    and complement functions

    Empty sets and intersecting intervals are appropriately handled in the constructor:

    julia> a = interval(0,2) ∪ interval(3,4)
    [0, 2] ∪ [3, 4]

    julia> b = interval(1,2) ∪ interval(4,5) ∪ ∅
    [1, 2] ∪ [4, 5]

    julia> c = a * b
    [0, 10] ∪ [12, 20]

    julia> complement(c)
    [-∞, 0] ∪ [10, 12] ∪ [20, ∞]

"""

###
#   IntervalUnion constructor. Consists of a vector of intervals
###
struct IntervalU{T<:Real} <: IntervalUnion{T}
    v :: Array{Interval{T}}
end

###
#   Outer constructors
###
function intervalU(x)
    x = IntervalU(x)
    sort!(x.v)
    x = remove_empties(x)
    x = condense(x)
    closeGaps!(x)
    return x
end

intervalU(num :: Real) = IntervalU([interval(num)])

intervalU(lo :: Real, hi :: Real) = IntervalU([interval(lo,hi)])
IntervalU(lo :: Real, hi :: Real) = IntervalU([interval(lo,hi)])

intervalU(x :: Interval) = IntervalU([x])
∪(x :: Interval) = intervalU(x)

∪(x :: Interval, y :: Interval) = intervalU([x; y])
∪(x :: Array{Interval{T}}) where T <:Real = intervalU(x)

intervalU(x :: Interval, y :: IntervalU) = intervalU([x; y.v])
∪(x :: Interval, y :: IntervalU) = intervalU(x,y)

intervalU(x :: IntervalU, y :: Interval) = intervalU([x.v; y])
∪(x :: IntervalU, y :: Interval) = intervalU(x,y)

intervalM(x :: IntervalU, y :: IntervalU) = intervalU([x.v; y.v])
∪(x :: IntervalU, y :: IntervalU) = intervalU(x,y)

# MultiInterval can act like a vector
getindex(x :: IntervalU, ind :: Integer) = getindex(x.v,ind)
getindex(x :: IntervalU, ind :: Array{ <: Integer}) = getindex(x.v,ind)

# Remove ∅ from IntervalUnion
function remove_empties(x :: IntervalU)
    v = x.v
    Vnew = v[v .!= ∅]
    return IntervalU(Vnew)
end

function closeGaps!(x :: IntervalU, maxInts = MAXINTS[1])

    while length(x.v) > maxInts     # Global

        # Complement code
        v = sort(x.v)
        vLo = left.(v)
        vHi = right.(v)
        vLo[1] == -∞ ? popfirst!(vLo) : pushfirst!(vHi, -∞)
        vHi[end] == ∞ ? pop!(vHi) : push!(vLo, ∞)
        xc = interval.(vHi,vLo)

        # Close smallest width of complement
        widths = diam.(xc)
        d, i = findmin(widths)
        merge = hull(x[i-1], x[i])
        popat!(x.v, i)
        popat!(x.v, i-1)
        push!(x.v, merge)
        sort!(x.v)
    end
    return x
end

# Recursively envolpe intervals which intersect.
function condense(x :: IntervalU)

    if iscondensed(x); return x; end

    v = sort(x.v)
    v = unique(v)

    Vnew = Interval{Float64}[]
    for i =1:length(v)

        intersects = intersect.(v[i],v )
        these = findall( intersects .!= ∅)

        push!(Vnew, hull(v[these]))
    end
    return condense( intervalU(Vnew) )
end

function iscondensed(x :: IntervalU)
    v = sort(x.v)
    for i=1:length(v)
        intersects = findall( intersect.(v[i],v[1:end .!= i]) .!= ∅)
        if !isempty(intersects); return false; end
    end
    return true
end

# Recursively envolpe intervals which intersect, except those which touch
function condense_weak(x :: IntervalU)

    if iscondensed(x); return x; end

    v = sort(x.v)
    v = unique(v)

    Vnew = Interval{Float64}[]
    for i =1:length(v)

        intersects = intersect.(v[i],v )

        notempty = intersects .!= ∅
        isItThin = isthin.(intersects)    # Don't hull intervals which touch

        notEmptyOrThin = notempty .* (1 .- isItThin)
        them = findall(notEmptyOrThin .== 1)

        push!(Vnew, hull(v[them]))
    end
    return condense( intervalU(Vnew) )
end

function iscondensed_weak(x :: IntervalU)
    v = sort(x.v)
    for i=1:length(v)
        intersects = intersect.(v[i],v[1:end .!= i])

        notempty = intersects .!= ∅
        isItThin =  isthin.(intersects)

        notEmptyOrThin = notempty .* (1 .- isItThin)

        if sum(notEmptyOrThin) != 0; return false; end
    end
    return true
end


text:
using Test

include("nucleotide-count.jl")

@testset "empty strand" begin
    @test count_nucleotides("") == Dict('A' => 0, 'C' => 0, 'G' => 0, 'T' => 0)
end

@testset "strang with one nucleotide in single-character input" begin
    @test count_nucleotides("G") == Dict('A' => 0, 'C' => 0, 'G' => 1, 'T' => 0)
end

@testset "strand with repeated nucleotide" begin
    @test count_nucleotides("GGGGGGG") == Dict('A' => 0, 'C' => 0, 'G' => 7, 'T' => 0)
end

@testset "strand with multiple nucleotides" begin
    @test count_nucleotides("AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC") == Dict('A' => 20, 'C' => 12, 'G' => 17, 'T' => 21)
end

@testset "strand with invalid nucleotides" begin
    @test_throws DomainError count_nucleotides("AGXXACT")
end


text:
using Plots, StatsPlots
using Plots.PlotMeasures
using DataFrames, Query
using BenchmarkTools

function graphbench_data(suite)
  data = DataFrame(subcat=String[],bench=String[],platform=String[],
                   mt_normalized=Float64[],mediantime=Float64[])
  graphbenches = suite["Graphs"]
  noncatlab_times = Dict{Tuple{String,String},Float64}()
  for (subcat,subsuite) in graphbenches
    for (platform,results) in subsuite
      for (bench,result) in results
        if platform ∈ ["LightGraphs", "MetaGraphs"]
            noncatlab_times[(subcat,bench)] = median(result).time
        end
        new_row = (subcat=subcat,
                    bench=bench,
                    platform=platform,
                    mt_normalized=0.,
                    mediantime=median(result).time)
        push!(data, new_row)
      end
    end
  end
  for i in 1:length(data.subcat)
    key = (data[i,:subcat], data[i,:bench])
    if key ∈ keys(noncatlab_times)
      data[i,:mt_normalized] = data[i,:mediantime] / noncatlab_times[key]
    end
  end
  data
end

function subcat_data(dat,subcat)
  dat |>
    @filter(_.subcat==subcat) |>
    @filter(_.platform == ["Catlab"]) |>
    @orderby((_.bench,_.platform)) |>
    @select(-:subcat) |>
    DataFrame
end

function plot_subcat(dat,subcat,yscale=:linear)
  subcat_data(dat,subcat) |>
    @df groupedbar(:bench,:mt_normalized,group=:platform,
                   xrotation=45,legend=:outerright,bar_width=0.5,
                   yscale=yscale, yguide="Rel. time", bottom_margin=50px)
end

function plot_all_subcats(dat)
  for subcat in unique(dat[!,:subcat])
    yscale = subcat ∈ ["WeightedGraph", "LabeledGraph"] ? :log : :linear
    fig = plot_subcat(dat,subcat,yscale)
    savefig(fig, string("figures/",subcat,".pdf"))
  end
end

max_forks_repo_name: epatters/CompCat


text:
"""
    Root

Object representing a possible root inside a given region. The field `status`
is either `:unknown` or `:unique`. If `status` is `:unique` then we know that
there is a unique root of the function in question inside the given region.

Internally the status may also be `:empty` for region guaranteed to contain no
root, however such `Root`s are discarded by default and thus never returned
by the `roots` function.

# Fields
  - `interval`: a region (either `Interval` of `IntervalBox`) searched for
        roots.
  - `status`: the status of the region, valid values are `:empty`, `unkown` and
        `:unique`.
"""
struct Root{T}
    interval::T
    status::Symbol
end

interval(rt::Root) = rt.interval

"""
    root_status(rt)

Return the status of a `Root`.
"""
root_status(rt::Root) = rt.status  # Use root_status since just `status` is too generic

"""
    root_region(rt)

Return the region associated to a `Root`.
"""
root_region(rt::Root) = rt.interval  # More generic name than `interval`.

"""
    isunique(rt)

Return whether a `Root` is unique.
"""
isunique(rt::Root{T}) where {T} = (rt.status == :unique)

show(io::IO, rt::Root) = print(io, "Root($(rt.interval), :$(rt.status))")

⊆(a::Interval, b::Root) = a ⊆ b.interval   # the Root object has the interval in the first entry
⊆(a::Root, b::Root) = a.interval ⊆ b.interval

big(a::Root) = Root(big(a.interval), a.status)


text:
using Factotum, Statistics, LinearAlgebra, Test
T, n, r = (200, 10, 6)
x = rand(T,n)                  # generate data
fm = Factotum.FactorModel(x, 6; scale = true)      # fit factor model

F = Factotum.factors(fm)
Λ = Factotum.loadings(fm)
σ = Factotum.sdev(fm)
## Check diagonality cov(F)

Σ = cov(F; corrected = false);

@test Σ ≈ diagm(0=>σ.^2)
@test Λ'Λ/n ≈ diagm(0=>ones(r))

## Replicate Ng and Bai (2002)
function simulate_factormodel(r, T, N)
    F = randn(T,r)
    Λ = randn(r,N)
    e = sqrt(r).*randn(T,N)
    F*Λ .+ e
end

X = simulate_factormodel(1, 251, 78)
fm = Factotum.FactorModel(X, 10; scale = true)
Factotum.informationcriteria((IC1, IC2), fm, 10)

X = simulate_factormodel(3, 100, 20)
fm = Factotum.FactorModel(X, 10; scale = true)
ic = Factotum.informationcriteria((IC1, IC2), fm, 10)

findmin(ic)
## Test the proportion of the variance explained

max_stars_repo_name: gragusa/Factotum

max_issues_repo_name: gragusa/Factotum

max_forks_repo_name: gragusa/Factotum


text:
using ChainRules
using ChainRulesCore
using Documenter

@show ENV

makedocs(
    modules=[ChainRules, ChainRulesCore],
    format=Documenter.HTML(prettyurls=false, assets = ["assets/chainrules.css"]),
    sitename="ChainRules",
    authors="Jarrett Revels and other contributors",
    pages=[
        "Introduction" => "index.md",
        "Getting Started" => "getting_started.md",
        "API" => "api.md",
    ],
)

const repo="github.com/JuliaDiff/ChainRules.jl.git"
const PR = get(ENV, "TRAVIS_PULL_REQUEST", "false")
if PR == "false"
    # Normal case, only deply docs if merging to master or release tagged
    deploydocs(repo=repo)
else
    @info "Deploying review docs for PR #$PR"
    # TODO: remove most of this once https://github.com/JuliaDocs/Documenter.jl/issues/1131 is resolved

    # Overwrite Documenter's function for generating the versions.js file
    foreach(Base.delete_method, methods(Documenter.Writers.HTMLWriter.generate_version_file))
    Documenter.Writers.HTMLWriter.generate_version_file(_, _) = nothing
    # Overwrite necessary environment variables to trick Documenter to deploy
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_BRANCH"] = "master"

    deploydocs(
        devurl="preview-PR$(PR)",
        repo=repo,
    )
end


text:
using CUTEst
using CSV
optimal_problems = []
infeasible_problems = []
unbounded_problems = []
max_delta_problems = []
iteration_limit_problems = []
max_time_problems = []
error_problems = []

function cutest_tests()
    @testset "cutest" begin
        @testset "chain" begin
            nlp = CUTEstModel("CHAIN")
            try
				pars = OnePhase.Class_parameters()
		        pars.output_level = 0
	        	iter, status, hist, t, err, timer = OnePhase.one_phase_solve(nlp, pars);
                #@show OnePhase.get_original_x(iter) # gives the primal solution of the solver
                #@show OnePhase.get_y(iter) # gives the dual solution of the solver
                @test true
            catch(e)
                println(e)
                @test false
            finally
                finalize(nlp)
            end
        end
    end
end

function runModelFromProblem(cutest_problem, output_results_to_file)
    nlp = CUTEstModel(cutest_problem)
    try
        println("-----------EXECUTING PROBLEM----------", cutest_problem)
        pars = OnePhase.Class_parameters()
        pars.output_level = 0
        iter, status, hist, t, err, timer = OnePhase.one_phase_solve(nlp, pars);

		if output_results_to_file
			outputResultsToCSVFile(cutest_problem, hist)
		end

        @test true

		if status == :Optimal
			push!(optimal_problems, cutest_problem)
	    elseif status == :primal_infeasible
	        push!(infeasible_problems, cutest_problem)
	    elseif status == :dual_infeasible
	        push!(unbounded_problems, cutest_problem)
	    elseif status == :MAX_IT
	        push!(iteration_limit_problems, cutest_problem)
	    elseif status === :MAX_TIME
	        push!(max_time_problems, cutest_problem)
		else
			push!(max_delta_problems, cutest_problem)
	    end
    catch e
        push!(error_problems, cutest_problem)
        @test false
    finally
        finalize(nlp)
    end
end

function executeCUTEST_Models()
    problems = CUTEst.select(min_con = 1, max_con=5, only_free_var=true, max_var=5)

    for problem in problems
		try
	    	runModelFromProblem(problem, false)
    	catch e
	    	println(e)
    	end
    end

    println("------------RUNNING CUTEST TEST SET SUMMARY------------")
    println("----------------------------------------------Total: ", length(problems))
    println("--------------------------------------------OPTIMAL: ", length(optimal_problems))
    println("--------------------------------------------OPTIMAL: ", optimal_problems)
    println("-----------------------------------------INFEASIBLE: ", length(infeasible_problems))
    println("-----------------------------------------INFEASIBLE: ", infeasible_problems)
    println("------------------------------------------UNBOUNDED: ", length(unbounded_problems))
    println("------------------------------------------UNBOUNDED: ", unbounded_problems)
	println("------------------------------------ITERATION_LIMIT: ", length(iteration_limit_problems))
	println("------------------------------------ITERATION_LIMIT: ", iteration_limit_problems)
	println("-------------------------------------------MAX_TIME: ", length(max_time_problems))
	println("-------------------------------------------MAX_TIME: ", max_time_problems)
	println("------------------------------------------MAX_DELTA: ", length(max_delta_problems))
	println("------------------------------------------MAX_DELTA: ", max_delta_problems)
    println("----------------------------------------------ERROR: ", length(error_problems))
    println("----------------------------------------------ERROR: ", error_problems)
end

function executeCUTEST_Models_benchmark()
    #problems = CUTEst.select(min_con = 1, max_con=100, only_free_var=false, max_var=100)
    #problems = ["HIER13", "CHARDIS1", "HAGER2", "QPCBOEI1", "CHANDHEU", "MPC2", "BLOWEYA", "SINROSNB", "MPC10", "ACOPP57", "ACOPP300", "KISSING2", "LISWET9", "YAO", "VANDERM1", "LUBRIFC", "GMNCASE2", "CHEMRCTA", "CATENARY", "HYDROELM", "CATMIX", "ROTDISC", "LISWET7", "SSNLBEAM", "HAIFAL", "INTEGREQ", "MSQRTB", "DITTERT", "TRAINF", "ROBOTARM", "GPP", "YATP1SS", "QPNBOEI1", "LISWET6", "UBH5", "AGG", "HIER133A", "STEENBRG", "DRCAVTY1", "10FOLDTR", "EIGENA", "ARGLALE", "QR3DBD", "HIER16", "A0NSDSDS", "LEUVEN4", "OPTCDEG2", "YORKNET", "STEENBRA", "EIGMAXC", "MPC1", "TABLE5", "HADAMARD", "STEENBRC", "MPC11", "LEUVEN7", "BLOWEYB", "SOSQP2", "DTOC2", "DALLASM", "TABLE4", "TRO41X9", "CHAIN", "READING8", "QR3D", "QPCBOEI2", "KISSING", "SREADIN3", "LISWET4", "ACOPR118", "OPTCTRL6", "ORTHRGDS", "MPC8", "TRO21X5", "HANGING", "READING7", "CHANDHEQ", "FLOSP2HM", "DTOC4", "POLYGON", "MSS3", "FLOSP2TL", "GROUPING", "HIER133C", "LISWET10", "ORTHREGD", "EIGENACO", "ACOPR57", "FERRISDC", "DTOC1NC", "KSS", "TABLE3", "TABLE1", "VANDERM2", "CBRATU2D", "EIGENB2", "LEUVEN3", "EIGENB", "BROWNALE", "DTOC1ND", "DRCAVTY2", "LEUVEN5", "DTOC3", "SOSQP1", "DIXCHLNV", "STEENBRF", "CORKSCRW", "MSS2", "SPIN", "HAGER1", "MPC15", "EIGMINA", "C-RELOAD", "DRUGDISE", "ZAMB2", "KTMODEL", "HYDROELS", "ARWHDNE", "HAGER3", "LISWET3", "MSQRTA", "STEENBRE", "LEAKNET", "MOSARQP1", "LEUVEN6", "NGONE", "GMNCASE3", "MPC9", "FLOSP2HH", "DALLASL", "LINCONT", "FLOSP2HL", "TOYSARAH", "MPC5", "HIER163A", "CHEMRCTB", "JJTABEL3", "BLOWEYC", "BRATU3D", "ORTHRDS2", "BROYDNBD", "PRIMAL3", "TRAINH", "HVYCRASH", "EIGENBCO", "LEUVEN1", "GMNCASE4", "HIER163E", "DTOC1L", "SPINOP", "CBRATU3D", "EIGMINC", "LISWET12", "YATP1NE", "MOSARQP2", "DTOC1NA", "LISWET8", "MANNE", "ROSEPETAL", "A4X12", "GMNCASE1", "EIGMAXA", "HAGER4", "EIGMAXB", "WOODSNE", "HIE1327D", "HIE1372D", "ROCKET", "LISWET11", "ACOPR300", "HIER133E", "STEERING", "MADSSCHJ", "LISWET5", "ARGLCLE", "VANDERM3", "QPNBOEI2", "EIGENA2", "MPC3", "BROYDN3D", "OPTCDEG3", "LEUVEN2", "EIGENC", "SMMPSF", "EIGENC2", "STEENBRD", "NUFFIELD", "ACOPP118", "ARGTRIG", "LISWET2", "ELEC", "MPC13", "A5NSSNSM", "MPC12", "MINC44", "HIER163D", "LISWET1", "TABLE6", "HIER163C", "HIER133D", "READING3", "EIGENCCO", "DRCAVTY3", "A5NSDSDM", "LUBRIF", "CATENA", "DTOC1NB", "FLOSP2TM", "TABLE7", "MINPERM", "FLOSP2TH", "MPC14", "ZIGZAG", "SVANBERG", "SAWPATH", "ARGLBLE", "SPIN2", "QPNSTAIR", "EIGENAU", "QPCSTAIR", "HELSBY", "HIER133B", "MPC4", "POWELL20", "MPC7", "HIER163B", "OPTCTRL3", "HYDROELL", "ORTHREGF", "STEENBRB", "COSHFUN", "CAMSHAPE", "MPC16", "READING1", "EIGMINB", "SPIN2OP", "MPC6", "ORTHREGC"]
    problems = ["HIER13", "CHARDIS1", "HAGER2", "QPCBOEI1", "CHANDHEU", "MPC2", "BLOWEYA", "SINROSNB", "MPC10", "ACOPP57", "ACOPP300", "KISSING2", "LISWET9", "YAO", "VANDERM1", "LUBRIFC", "GMNCASE2", "CHEMRCTA", "CATENARY", "HYDROELM", "CATMIX", "ROTDISC", "LISWET7", "SSNLBEAM", "HAIFAL", "INTEGREQ", "MSQRTB", "DITTERT", "TRAINF", "ROBOTARM", "GPP", "QPNBOEI1", "LISWET6", "UBH5", "AGG", "HIER133A", "STEENBRG", "DRCAVTY1", "10FOLDTR", "EIGENA", "ARGLALE", "QR3DBD", "HIER16", "A0NSDSDS", "LEUVEN4", "OPTCDEG2", "YORKNET", "STEENBRA", "EIGMAXC", "MPC1", "TABLE5", "HADAMARD", "STEENBRC", "MPC11", "LEUVEN7", "BLOWEYB", "SOSQP2", "DTOC2", "DALLASM", "TABLE4", "TRO41X9", "CHAIN", "READING8", "QR3D", "QPCBOEI2", "KISSING", "SREADIN3", "LISWET4", "ACOPR118", "OPTCTRL6", "ORTHRGDS", "MPC8", "TRO21X5", "HANGING", "READING7", "CHANDHEQ", "FLOSP2HM", "DTOC4", "POLYGON", "MSS3", "FLOSP2TL", "GROUPING", "HIER133C", "LISWET10", "ORTHREGD", "EIGENACO", "ACOPR57", "FERRISDC", "DTOC1NC", "KSS", "TABLE3", "TABLE1", "VANDERM2", "CBRATU2D", "EIGENB2", "LEUVEN3", "EIGENB", "BROWNALE", "DTOC1ND", "DRCAVTY2", "LEUVEN5", "DTOC3", "SOSQP1", "DIXCHLNV", "STEENBRF", "CORKSCRW", "MSS2", "SPIN", "HAGER1", "MPC15", "EIGMINA", "C-RELOAD", "DRUGDISE", "ZAMB2", "KTMODEL", "HYDROELS", "ARWHDNE", "HAGER3", "LISWET3", "MSQRTA", "STEENBRE", "LEAKNET", "MOSARQP1", "LEUVEN6", "NGONE", "GMNCASE3", "MPC9", "FLOSP2HH", "DALLASL", "LINCONT", "FLOSP2HL", "TOYSARAH", "MPC5", "HIER163A", "CHEMRCTB", "JJTABEL3", "BLOWEYC", "BRATU3D", "ORTHRDS2", "BROYDNBD", "PRIMAL3", "TRAINH", "HVYCRASH", "EIGENBCO", "LEUVEN1", "GMNCASE4", "HIER163E", "DTOC1L", "SPINOP", "CBRATU3D", "EIGMINC", "LISWET12", "YATP1NE", "MOSARQP2", "DTOC1NA", "LISWET8", "MANNE", "ROSEPETAL", "A4X12", "GMNCASE1", "EIGMAXA", "HAGER4", "EIGMAXB", "WOODSNE", "HIE1327D", "HIE1372D", "ROCKET", "LISWET11", "ACOPR300", "HIER133E", "STEERING", "MADSSCHJ", "LISWET5", "ARGLCLE", "VANDERM3", "QPNBOEI2", "EIGENA2", "MPC3", "BROYDN3D", "OPTCDEG3", "LEUVEN2", "EIGENC", "SMMPSF", "EIGENC2", "STEENBRD", "NUFFIELD", "ACOPP118", "ARGTRIG", "LISWET2", "ELEC", "MPC13", "A5NSSNSM", "MPC12", "MINC44", "HIER163D", "LISWET1", "TABLE6", "HIER163C", "HIER133D", "READING3", "EIGENCCO", "DRCAVTY3", "A5NSDSDM", "LUBRIF", "CATENA", "DTOC1NB", "FLOSP2TM", "TABLE7", "MINPERM", "FLOSP2TH", "MPC14", "ZIGZAG", "SVANBERG", "SAWPATH", "ARGLBLE", "SPIN2", "QPNSTAIR", "EIGENAU", "QPCSTAIR", "HELSBY", "HIER133B", "MPC4", "POWELL20", "MPC7", "HIER163B", "OPTCTRL3", "HYDROELL", "ORTHREGF", "STEENBRB", "COSHFUN", "CAMSHAPE", "MPC16", "READING1", "EIGMINB", "SPIN2OP", "MPC6", "ORTHREGC"]
    for problem in problems
        runModelFromProblem(problem, true)
    end
	println("------------RUNNING CUTEST TEST SET SUMMARY------------")
    println("----------------------------------------------Total: ", length(problems))
    println("--------------------------------------------OPTIMAL: ", length(optimal_problems))
    println("--------------------------------------------OPTIMAL: ", optimal_problems)
    println("-----------------------------------------INFEASIBLE: ", length(infeasible_problems))
    println("-----------------------------------------INFEASIBLE: ", infeasible_problems)
    println("------------------------------------------UNBOUNDED: ", length(unbounded_problems))
    println("------------------------------------------UNBOUNDED: ", unbounded_problems)
	println("------------------------------------ITERATION_LIMIT: ", length(iteration_limit_problems))
	println("------------------------------------ITERATION_LIMIT: ", iteration_limit_problems)
	println("-------------------------------------------MAX_TIME: ", length(max_time_problems))
	println("-------------------------------------------MAX_TIME: ", max_time_problems)
	println("------------------------------------------MAX_DELTA: ", length(max_delta_problems))
	println("------------------------------------------MAX_DELTA: ", max_delta_problems)
    println("----------------------------------------------ERROR: ", length(error_problems))
    println("----------------------------------------------ERROR: ", error_problems)

    outputTotalResultsToFile(problems, optimal_problems, infeasible_problems, unbounded_problems, error_problems)
end

function outputResultsToCSVFile(cutest_problem, hist)
    current_dir = pwd()
    CSV.write("$current_dir/cutest_outputs/$cutest_problem.csv", hist, header = true)
end

function outputTotalResultsToFile(problems, optimal_problems, infeasible_problems, unbounded_problems, error_problems)
    current_dir = pwd()
    total = length(problems)
    total_optimal = length(optimal_problems)
    total_infeasible = length(infeasible_problems)
    total_unbounded = length(unbounded_problems)
    total_iteration_limit = length(iteration_limit_problems)
	total_max_time = length(max_time_problems)
	total_max_delta = length(max_delta_problems)
	total_erros = length(error_problems)

    open("$current_dir/cutest_outputs/total_results.txt" , "w") do file
	write(file, "----------------------------------------------Total: $total\n")
	write(file, "--------------------------------------------OPTIMAL: $total_optimal\n")
	write(file, "--------------------------------------------OPTIMAL: $optimal_problems\n")
	write(file, "-----------------------------------------INFEASIBLE: $total_infeasible\n")
	write(file, "-----------------------------------------INFEASIBLE: $infeasible_problems\n")
	write(file, "------------------------------------------UNBOUNDED: $total_unbounded\n")
	write(file, "------------------------------------------UNBOUNDED: $unbounded_problems\n")
	write(file, "------------------------------------ITERATION_LIMIT: $total_iteration_limitn\n")
	write(file, "------------------------------------ITERATION_LIMIT: $iteration_limit_problems\n")
	write(file, "-------------------------------------------MAX_TIME: $total_max_time\n")
	write(file, "-------------------------------------------MAX_TIME: $max_time_problems\n")
	write(file, "------------------------------------------MAX_DELTA: $total_max_delta")
	write(file, "------------------------------------------MAX_DELTA: $max_delta_problems\n")
	write(file, "----------------------------------------------ERROR: $total_erros\n")
	write(file, "----------------------------------------------ERROR: $error_problems")
    end
end


text:
############################## Lennard Jones ###################################
struct LennardJones <: EmpiricalPotential
    ϵ::AbstractFloat
    σ::AbstractFloat
    rcutoff::AbstractFloat
end

get_trainable_params(lj::LennardJones) = Parameter{:ϵ,:σ}((lj.ϵ, lj.σ))

get_nontrainable_params(lj::LennardJones) = Parameter{}(())

############################# Energies ##########################################

function potential_energy(R::AbstractFloat, p::LennardJones)
    d = p.σ / R
    4.0 * p.ϵ * (d^12 - d^6)
end

############################### Forces ##########################################

force(r::SVector{3,<:AbstractFloat}, p::LennardJones) = force(r, norm(r), p)

function force(R::AbstractFloat, r::SVector{3,<:AbstractFloat}, p::LennardJones)
    SVector(24.0 * p.ϵ * (2.0 * (p.σ / R)^12 - (p.σ / R)^6) .* r ./ R ./ R)
end
############################## Gradients ########################################
function grad_potential_energy(r::SVector{3,<:AbstractFloat}, p::LennardJones)
    d = p.σ / norm(r)
    (dpdϵ = 4.0 * (d^12 - d^6), dpdσ = 24.0 * p.ϵ / p.σ * (2 * d^12 - d^6))
end

function grad_force(r::SVector{3,<:Real}, p::LennardJones)
    d = norm(r)
    (dfdϵ = 48.0 * ((p.σ / d)^13 - 0.5 * (p.σ / d)^7) .* r ./ d,
        dfdσ = 144.0 * p.ϵ / p.σ * (4.0 * p.σ^12 / d^13 - p.σ^6 / d^7) .* r ./ d)
end

function grad_virial(r::SVector{<:Real}, p::LennardJones)
    df = grad_force(r, p)
    dfde = df[1]
    dfdsig = df[2]

    (dvdϵ = dfde[1] * r[1] + dfde[2] * r[2] + dfde[3] * r[3],
        dvdσ = dfdsig[1] * r[1] + dfdsig[2] * r[2] + dfdsig[3] * r[3])
end


text:
using Compat
using EnglishText

using Test

@testset "Articles" begin
    @test indefinite("pig") == "a"
    @test indefinite("iron") == "an"
    @test indefinite("unicorn") == "a"
    @test indefinite("honest") == "an"
    @test indefinite("honey") == "a"
    @test indefinite("NATO") == "a"
    @test indefinite("NSA") == "an"
end

include("numeric.jl")
include("pluralize.jl")
include("list.jl")
include("text.jl")
include("quantity.jl")


text:
#=
# Tests out using a variety of coefficients in different places
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("coef");

# Try making an optional log
@useLog("coeflog")

# Set up the configuration
@domain(2, SQUARE, UNSTRUCTURED)    # dimension, geometry, decomposition
@solver(CG)                         # DG, CG, etc.
@functionSpace(LEGENDRE, 2)         # function, order (or use testFunction and trialFunction)
@nodes(LOBATTO)                     # elemental node arrangement

# Specify the problem
@mesh(QUADMESH, 20)                 # .msh file or generate our own

@variable(u)                        # same as @variable(u, SCALAR)

@testSymbol(v)                    # sets the symbol for a test function

@boundary(u, 1, DIRICHLET, 0)

# Write the weak form
@coefficient(a, "x")
@coefficient(b, "3")
@coefficient(c, "y")
@coefficient(d, 0.1)
@coefficient(f, "(-5*pi*pi*x*sin(pi*x)*sin(2*pi*y) + pi*cos(pi*x)*sin(2*pi*y)) + y*sin(pi*x)*sin(2*pi*y)")

@weakForm(u, "-a*b*dot(grad(u), grad(v))*d + c*u*d*v*b - d*f*b*v")

solve(u);

# exact solution is sin(pi*x)*sin(2*pi*y)
# check error
maxerr = 0;
exact(x,y) = sin(pi*x)*sin(2*pi*y);

for i=1:size(Femshop.grid_data.allnodes,2)
    x = Femshop.grid_data.allnodes[1,i];
    y = Femshop.grid_data.allnodes[2,i];
    err = abs(u.values[i] - exact(x,y));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

using Plots
pyplot();
display(plot(Femshop.grid_data.allnodes[1,:], Femshop.grid_data.allnodes[1,:], u.values[:], st=:surface))

# check
log_dump_config();
log_dump_prob();

@finalize()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct InstanceActionInfo <: SwaggerModel
    action_type::Union{ Nothing, String } # action_type
    payload::Union{ Nothing, String } # payload

    function InstanceActionInfo(;action_type=nothing, payload=nothing)
        o = new()
        set_field!(o, :action_type, action_type)
        set_field!(o, :payload, payload)
        o
    end
end # type InstanceActionInfo

const _name_map_InstanceActionInfo = Dict{String,Symbol}(["action_type"=>:action_type, "payload"=>:payload])
const _field_map_InstanceActionInfo = Dict{Symbol,String}([:action_type=>"action_type", :payload=>"payload"])
Swagger.name_map(::Type{ InstanceActionInfo }) = _name_map_InstanceActionInfo
Swagger.field_map(::Type{ InstanceActionInfo }) = _field_map_InstanceActionInfo

const _allowed_InstanceActionInfo_action_type = ["BlockDeviceRescan", "FlushMetrics", "InstanceStart", "SendCtrlAltDel"]

function check_required(o::InstanceActionInfo)
    true
end

function validate_field(o::InstanceActionInfo, name::Symbol, val)
    if name === :action_type
        Swagger.validate_param(name, "InstanceActionInfo", :enum, val, _allowed_InstanceActionInfo_action_type)
    end
end


text:
using Revise
using BridgeLandmarks
const BL=BridgeLandmarks
using Random
using Distributions
using DataFrames
using DelimitedFiles
using CSV
using StaticArrays
using LinearAlgebra
using JLD2
using FileIO

Random.seed!(9)

################ set directories ##########################################
workdir = @__DIR__
cd(workdir)
#include("../../outdirpath.jl")
include("outdirpath.jl")

dir1 = joinpath(outdirpath,"exp2/ms")
outdir_ms = mkpath(dir1)

dir2 = joinpath(outdirpath,"exp2/ahs")
outdir_ahs = mkpath(dir2)

################ read data ##########################################
#dat = load("data_exp2.jld2")
dat = load("revised_experiments/exp2/data_exp2.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]
n = dat["n"]

################ settings and mcmc #################################
ups =  [:innov, :rmmala_pos, :parameter]
adaptskip = 100
skip_saveITER = 50#0
printskip = 25
ITER = 2500#0

δpos = (d*n)^(-1/6) * [0.001, 0.0001]
p_ms = Pars_ms(δpos=δpos,  γinit=.1/√n,  ## LET OP
                  aprior=Pareto(1.0, 0.1), η =  n -> 0.0, dt = 0.01,
                adaptskip=adaptskip, skip_saveITER=skip_saveITER, ρlowerbound=0.8, δa=0.02)
p_ahs = Pars_ahs(δpos=δpos,  db=[2.5, 2.5],stdev=.25,γinit=.1, aprior=Pareto(1.0, 0.1), η =  n -> 0.0, dt = 0.001,
                                adaptskip=adaptskip, skip_saveITER=skip_saveITER, ρlowerbound=0.9)
# original experiment
#p_ms = Pars_ms(ρinit = 0.7,covθprop = [0.04 0. 0.; 0. 0.04 0.; 0. 0. 0.04], δpos=0.005, δmom=0.1,cinit=0.2,γinit=2.0  )
#p_ahs = Pars_ahs(ρinit = 0.7,covθprop = [0.04 0. 0.; 0. 0.04 0.; 0. 0. 0.04], δpos=0.01, δmom=0.1,cinit=0.02,γinit=2.0, stdev=0.75, db = [2.5, 2.5]  )

################## deliberately start in a wrong configuration ###################
n = length(xobsT[1])
xinitq = xobsT[1]
θ, ψ =  π/6, 0.25
rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
stretch = SMatrix{2,2}(1.0 + ψ, 0.0, 0.0, 1.0 - ψ)
xinitq_adj = [rot * stretch * xinitq[i] for i in 1:n ]

template_estimation(xobsT; xinitq=xinitq_adj,pars = p_ms,
     ITER=ITER, updatescheme = ups, printskip=printskip, outdir=outdir_ms)

#template_estimation(xobsT; xinitq=xinitq_adj,pars = p_ahs,
#      ITER=ITER, updatescheme = ups, printskip=printskip, outdir=outdir_ahs)
#plottemplate_estimation(outdir)


text:
push!(LOAD_PATH, abspath(@__DIR__,".."))
using LinearAlgebra,Statistics,StatsBase
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs

function onesparsemat(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

##
dt = 1E-3
# two LIF neurons, E and I
myτe = 0.2
myτi = 0.1
vth = 10.
v_r = -5.0
τrefre = 0.3 # refractoriness
τrefri = 0.1
τpcd = 0.2 # post synaptic current decay
e1 = S.NTLIF(myτe,vth,v_r,τrefre,τpcd)
i1 = S.NTLIF(myτi,vth,v_r,τrefri,τpcd)
pse1 = S.PSLIF(e1,1)
psi1 = S.PSLIF(i1,1)

# one static input
h_in_e = onesparsemat(10.1)
h_in_i = onesparsemat(9.0)
in_type = S.InputSimpleOffset()
in_state_e = S.PSSimpleInput(in_type)
in_state_i = S.PSSimpleInput(in_type)

# connect E <-> I , both ways, but no autapses
# i connections should be negative!

conn_ie = S.ConnLIF(onesparsemat(1.0))
conn_ei = S.ConnLIF(onesparsemat(-0.8))
conn_in_e = S.BaseConnection(h_in_e)
conn_in_i = S.BaseConnection(h_in_i)

# connected populations

pop_e = S.Population(pse1,(conn_ei,conn_in_e),(psi1,in_state_e))
pop_i = S.Population(psi1,(conn_ie,conn_in_i),(pse1,in_state_i))

# that's it, let's make the network
myntw = S.RecurrentNetwork(dt,(pop_e,pop_i))

##

Ttot = 10.
times = (0:myntw.dt:Ttot)
nt = length(times)
# initial conditions
pse1.state_now[1] = v_r
psi1.state_now[1] = v_r + 0.95*(vth-v_r)
# things to save
myvse = Vector{Float64}(undef,nt) # voltage
myfiringe = BitVector(undef,nt) # spike raster
myrefre = similar(myfiringe)  # if it is refractory
eicurr = similar(myvse)  # e-i current

myvsi = Vector{Float64}(undef,nt)
myfiringi = BitVector(undef,nt)
myrefri = similar(myfiringe)
iecurr = similar(myvsi)

for (k,t) in enumerate(times)
  S.dynamics_step!(t,myntw)
  myvse[k] = pse1.state_now[1]
  myfiringe[k]=pse1.isfiring[1]
  myrefre[k]=pse1.isrefractory[1]
  myvsi[k] = psi1.state_now[1]
  myfiringi[k]=psi1.isfiring[1]
  myrefri[k]=psi1.isrefractory[1]
  eicurr[k]=conn_ei.post_current[1]
  iecurr[k]=conn_ie.post_current[1]
end

# add spikes for plotting purposes, the eight is set arbitrarily to
# three times the firing threshold
myvse[myfiringe] .= 3 * e1.v_threshold
myvsi[myfiringi] .= 3 * e1.v_threshold;

##

plt=plot(times,myvse;leg=false,linewidth=3,
  ylabel="E (mV)",
  color=colorant"Midnight Blue")

# the green line indicates when the neuron is refractory
plot!(plt,times, 20.0*myrefre; opacity=0.6, color=:green)

plti=plot(times,myvsi;leg=false,linewidth=3,
   ylabel="I (mV)",
   color=colorant"Brick Red")

plot!(plti,times, 20.0*myrefri; opacity=0.6, color=:green)

pltcurr=plot(times, [ iecurr eicurr]; leg=false,
  linewidth=3.0, ylabel="E/I connection curr.",
  color=[colorant"Teal" colorant"Orange"])

plot(plt,plti,pltcurr; layout=(3,1),
  xlabel=["" "" "time (s)"])


text:
function bellmanford(G::DictGraph, s::Int)

end