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

using Test, Distributed, Random

import Logging: Debug, Info, Warn

@testset "@test" begin
    @test true
    @test 1 == 1
    @test 1 != 2
    @test ==(1, 1)
    @test ==((1, 1)...)
    @test 1 ≈ 2 atol=1
    @test strip("\t  hi   \n") == "hi"
    @test strip("\t  this should fail   \n") != "hi"
    @test isequal(1, 1)
    @test isapprox(1, 1, atol=0.1)
    @test isapprox(1, 1; atol=0.1)
    @test isapprox(1, 1; [(:atol, 0)]...)
end
@testset "@test keyword precedence" begin
    # post-semicolon keyword, suffix keyword, pre-semicolon keyword
    @test isapprox(1, 2, atol=0) atol=1
    @test isapprox(1, 3, atol=0; atol=2) atol=1
end
@testset "@test should only evaluate the arguments once" begin
    g = Int[]
    f = (x) -> (push!(g, x); x)
    @test f(1) == 1
    @test g == [1]

    empty!(g)
    @test isequal(f(2), 2)
    @test g == [2]
end

@testset "@test_broken with fail" begin
    @test_broken false
    @test_broken 1 == 2
    @test_broken 1 != 1
    @test_broken strip("\t  hi   \n") != "hi"
    @test_broken strip("\t  this should fail   \n") == "hi"
end
@testset "@test_broken with errors" begin
    @test_broken error()
    @test_broken absolute_nonsense
end
@testset "@test_skip" begin
    @test_skip error()
    @test_skip true
    @test_skip false
    @test_skip gobbeldygook
end
@testset "@test_warn" begin
    @test 1234 === @test_nowarn(1234)
    @test 5678 === @test_warn("WARNING: foo", begin println(STDERR, "WARNING: foo"); 5678; end)
    let a
        @test_throws UndefVarError(:a) a
        @test_nowarn a = 1
        @test a === 1
    end
end
@testset "@test and elements of an array" begin
    a = Array{Float64, 5}(uninitialized, 2, 2, 2, 2, 2)
    a[1, 1, 1, 1, 1] = 10
    @test a[1, 1, 1, 1, 1] == 10
    @test a[1, 1, 1, 1, 1] != 2
end

@test rand() != rand()

@testset "Pass - exception" begin
    @test endswith(sprint(show, @test_throws ErrorException error()),
                   "Thrown: ErrorException")
    @test endswith(sprint(show, @test_throws ErrorException("test") error("test")),
                   "Thrown: ErrorException")
end
# Test printing of Fail results
mutable struct NoThrowTestSet <: Test.AbstractTestSet
    results::Vector
    NoThrowTestSet(desc) = new([])
end
Test.record(ts::NoThrowTestSet, t::Test.Result) = (push!(ts.results, t); t)
Test.finish(ts::NoThrowTestSet) = ts.results
let fails = @testset NoThrowTestSet begin
        # Fail - wrong exception
        @test_throws OverflowError error()
        # Fail - no exception
        @test_throws OverflowError 1 + 1
        # Fail - comparison
        @test 1+1 == 2+2
        # Fail - approximate comparison
        @test 1/1 ≈ 2/1
        # Fail - chained comparison
        @test 1+0 == 2+0 == 3+0
        # Fail - comparison call
        @test ==(1 - 2, 2 - 1)
        # Fail - splatting
        @test ==(1:2...)
        # Fail - isequal
        @test isequal(0 / 0, 1 / 0)
        # Fail - function splatting
        @test isequal(1:2...)
        # Fail - isapprox
        @test isapprox(0 / 1, -1 / 0)
        # Fail - function with keyword
        @test isapprox(1 / 2, 2 / 1, atol=1 / 1)
        @test isapprox(1 - 2, 2 - 1; atol=1 - 1)
        # Fail - function keyword splatting
        k = [(:atol, 0), (:nans, true)]
        @test isapprox(1, 2; k...)
        # Error - unexpected pass
        @test_broken true
        # Error - converting a call into a comparison
        @test ==(1, 1:2...)
    end
    for i in 1:length(fails) - 2
        @test isa(fails[i], Test.Fail)
    end

    let str = sprint(show, fails[1])
        @test contains(str, "Expression: error()")
        @test contains(str, "Thrown: ErrorException")
    end

    let str = sprint(show, fails[2])
        @test contains(str, "Expression: 1 + 1")
        @test contains(str, "No exception thrown")
    end

    let str = sprint(show, fails[3])
        @test contains(str, "Expression: 1 + 1 == 2 + 2")
        @test contains(str, "Evaluated: 2 == 4")
    end

    let str = sprint(show, fails[4])
        @test contains(str, "Expression: 1 / 1 ≈ 2 / 1")
        @test contains(str, "Evaluated: 1.0 ≈ 2.0")
    end

    let str = sprint(show, fails[5])
        @test contains(str, "Expression: 1 + 0 == 2 + 0 == 3 + 0")
        @test contains(str, "Evaluated: 1 == 2 == 3")
    end

    let str = sprint(show, fails[6])
        @test contains(str, "Expression: 1 - 2 == 2 - 1")
        @test contains(str, "Evaluated: -1 == 1")
    end

    let str = sprint(show, fails[7])
        @test contains(str, "Expression: (==)(1:2...)")
        @test !contains(str, "Evaluated")
    end

    let str = sprint(show, fails[8])
        @test contains(str, "Expression: isequal(0 / 0, 1 / 0)")
        @test contains(str, "Evaluated: isequal(NaN, Inf)")
    end

    let str = sprint(show, fails[9])
        @test contains(str, "Expression: isequal(1:2...)")
        @test contains(str, "Evaluated: isequal(1, 2)")
    end

    let str = sprint(show, fails[10])
        @test contains(str, "Expression: isapprox(0 / 1, -1 / 0)")
        @test contains(str, "Evaluated: isapprox(0.0, -Inf)")
    end

    let str = sprint(show, fails[11])
        @test contains(str, "Expression: isapprox(1 / 2, 2 / 1, atol=1 / 1)")
        @test contains(str, "Evaluated: isapprox(0.5, 2.0; atol=1.0)")
    end

    let str = sprint(show, fails[12])
        @test contains(str, "Expression: isapprox(1 - 2, 2 - 1; atol=1 - 1)")
        @test contains(str, "Evaluated: isapprox(-1, 1; atol=0)")
    end

    let str = sprint(show, fails[13])
        @test contains(str, "Expression: isapprox(1, 2; k...)")
        @test contains(str, "Evaluated: isapprox(1, 2; atol=0, nans=true)")
    end

    let str = sprint(show, fails[14])
        @test contains(str, "Unexpected Pass")
        @test contains(str, "Expression: true")
    end

    let str = sprint(show, fails[15])
        @test contains(str, "Expression: ==(1, 1:2...)")
        @test contains(str, "MethodError: no method matching ==(::$Int, ::$Int, ::$Int)")
    end
end

@testset "printing of a TestSetException" begin
    tse_str = sprint(show, Test.TestSetException(1, 2, 3, 4, Vector{Union{Test.Error, Test.Fail}}()))
    @test contains(tse_str, "1 passed")
    @test contains(tse_str, "2 failed")
    @test contains(tse_str, "3 errored")
    @test contains(tse_str, "4 broken")
end

@test Test.finish(Test.FallbackTestSet()) !== nothing

OLD_STDOUT = STDOUT
OLD_STDERR = STDERR
catch_out = IOStream("")
catch_err = IOStream("")
rde, wre = redirect_stderr()
rdo, wro = redirect_stdout()

# test that FallbackTestSet will throw immediately
cmd = `$(Base.julia_cmd()) --startup-file=no --depwarn=error test_exec.jl`
@test !success(pipeline(cmd))

@testset "no errors" begin
    @test true
    @test 1 == 1
end

# Test entirely empty test set
@testset "outer" begin
    @testset "inner" begin
    end
end

@testset "testset types" begin
    ts = @testset "@testset should return the testset" begin
        @test true
    end
    @test typeof(ts) == Test.DefaultTestSet
    @test ts.n_passed == 1
    tss = @testset "@testset/for should return an array of testsets: $i" for i in 1:3
        @test true
    end
    @test length(tss) == 3
    @test typeof(tss[1]) == Test.DefaultTestSet
    @test tss[1].n_passed == 1
end
@testset "accounting" begin
    local ts, fails
    try
        ts = @testset "outer" begin
            @testset "inner1" begin
                @test true
                @test false
                @test 1 == 1
                @test 2 == :foo
                @test 3 == 3
                @testset "d" begin
                    @test 4 == 4
                end
                @testset begin
                    @test :blank != :notblank
                end
            end
            @testset "inner1" begin
                @test 1 == 1
                @test 2 == 2
                @test 3 == :bar
                @test 4 == 4
                @test_throws ErrorException 1+1
                @test_throws ErrorException error()
                @test_throws RemoteException error()
                @testset "errrrr" begin
                    @test "not bool"
                    @test error()
                end

                error("exceptions in testsets should be caught")
                @test 1 == 1 # this test will not be run
            end

            @testset "loop with desc" begin
                @testset "loop1 $T" for T in (Float32, Float64)
                    @test 1 == T(1)
                end
            end
            @testset "loops without desc" begin
                @testset for T in (Float32, Float64)
                    @test 1 == T(1)
                end
                @testset for T in (Float32, Float64), S in (Int32,Int64)
                    @test S(1) == T(1)
                end
            end
            @testset "some loops fail" begin
                @testset for i in 1:5
                    @test i <= 4
                end
                # should add 3 errors and 3 passing tests
                @testset for i in 1:6
                    iseven(i) || error("error outside of test")
                    @test true # only gets run if the above passed
                end
            end
        end
        # These lines shouldn't be called
        error("No exception was thrown!")
    catch ex
        redirect_stdout(OLD_STDOUT)
        redirect_stderr(OLD_STDERR)
        ex
    end
    @testset "ts results" begin
        @test isa(ts, Test.DefaultTestSet)
        passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = Test.get_test_counts(ts)
        total_pass   = passes + c_passes
        total_fail   = fails  + c_fails
        total_error  = errors + c_errors
        total_broken = broken + c_broken
        @test total_pass   == 24
        @test total_fail   == 6
        @test total_error  == 6
        @test total_broken == 0
    end
    ts.anynonpass = false
    deleteat!(Test.get_testset().results,1)
end

@test .1+.1+.1 ≈ .3
@test .1+.1+.1 ≉ .4

ts = @testset "@testset should return the testset" begin
    @test true
end
@test typeof(ts) == Test.DefaultTestSet
@test ts.n_passed == 1

tss = @testset "@testset/for should return an array of testsets: $i" for i in 1:3
    @test true
end
@test length(tss) == 3
@test typeof(tss[1]) == Test.DefaultTestSet
@test tss[1].n_passed == 1

# Issue #17908 (return)
testset_depth17908 = Test.get_testset_depth()
@testset for i in 1:3
    i > 1 && return
    @test i == 1
end
# The return aborts the control flow so the expression above doesn't return a
# value. The only thing we can test is whether the testset is properly popped.
# Do not use `@test` since the issue this is testing will swallow the error.
@assert testset_depth17908 == Test.get_testset_depth()

# Issue #17462 and Issue #17908 (break, continue)
testset_depth17462 = Test.get_testset_depth()
counter_17462_pre = 0
counter_17462_post = 0
tss17462 = @testset for x in [1,2,3,4]
    global counter_17462_pre, counter_17462_post
    counter_17462_pre += 1
    if x == 1
        @test counter_17462_pre == x
        continue
        @test false
    elseif x == 3
        @test counter_17462_pre == x
        break
        @test false
    elseif x == 4
        @test false
    else
        @test counter_17462_pre == x
        @test x == 2
        @test counter_17462_post == 0
    end
    counter_17462_post += 1
end
# Do not use `@test` since the issue this is testing will swallow the error.
# Same for the `@assert` in the for loop below
@assert testset_depth17462 == Test.get_testset_depth()
@assert length(tss17462) == 3
for ts17462 in tss17462
    @assert isa(ts17462, Test.DefaultTestSet)
end
@test counter_17462_pre == 3
@test counter_17462_post == 1

# Issue #21008
ts = try
    @testset "@test_broken and @test_skip should not give an exception" begin
        @test_broken false
        @test_skip true
        @test_skip false
    end
catch
    nothing # Shouldn't get here
end
@test ts isa Test.DefaultTestSet

# now we're done running tests with DefaultTestSet so we can go back to STDOUT
redirect_stdout(OLD_STDOUT)
redirect_stderr(OLD_STDERR)

# import the methods needed for defining our own testset type
import Test: record, finish
using Test: get_testset_depth, get_testset
using Test: AbstractTestSet, Result, Pass, Fail, Error
struct CustomTestSet <: Test.AbstractTestSet
    description::AbstractString
    foo::Int
    results::Vector
    # constructor takes a description string and options keyword arguments
    CustomTestSet(desc; foo=1) = new(desc, foo, [])
end

record(ts::CustomTestSet, child::AbstractTestSet) = push!(ts.results, child)
record(ts::CustomTestSet, res::Result) = push!(ts.results, res)
function finish(ts::CustomTestSet)
    # just record if we're not the top-level parent
    if get_testset_depth() > 0
        record(get_testset(), ts)
    end
    ts
end

ts = @testset CustomTestSet "Testing custom testsets" begin
    # this testset should inherit the parent testset type
    @testset "custom testset inner 1" begin
        @test true
        @test false
        @test error("this error will be reported as an error")
        @test_throws ErrorException nothing
        @test_throws ErrorException error("this error is a success")
    end
    # this testset has its own testset type
    @testset CustomTestSet foo=4 "custom testset inner 2" begin
        # this testset should inherit the type, but not the argument. If a particular
        # testset type wants inheritance behavior they should implement it themselves
        # using get_testset() in the constructor
        @testset "custom testset inner 2 inner 1" begin
            @test true
        end
        # make sure the RHS can use computed values, also tests options without
        # specifying the testset type
        @testset foo=(1+2) "custom testset inner 2 inner 2" begin
            @test true
        end
    end
end

@test typeof(ts) == CustomTestSet
@test ts.foo == 1
@test ts.description == "Testing custom testsets"
@test typeof(ts.results[1]) == CustomTestSet
@test ts.results[1].description == "custom testset inner 1"
@test ts.results[1].foo == 1
@test typeof(ts.results[1].results[1]) == Pass
@test typeof(ts.results[1].results[2]) == Fail
@test typeof(ts.results[1].results[3]) == Error
@test typeof(ts.results[1].results[4]) == Fail
@test typeof(ts.results[1].results[5]) == Pass

@test typeof(ts.results[2]) == CustomTestSet
@test ts.results[2].description == "custom testset inner 2"
@test ts.results[2].foo == 4
@test typeof(ts.results[2].results[1]) == CustomTestSet
@test ts.results[2].results[1].foo == 1
@test typeof(ts.results[2].results[1].results[1]) == Pass
@test typeof(ts.results[2].results[2]) == CustomTestSet
@test ts.results[2].results[2].foo == 3

# test custom testset types on testset/for
tss = @testset CustomTestSet foo=3 "custom testset $i" for i in 1:6
    @testset "inner testset $i-$j" for j in 1:3
        @test iseven(i + j)
    end
    # make sure a testset within a testset/for works
    @testset "inner testset $i" begin
        @test iseven(i)
    end
end

for i in 1:6
    @test typeof(tss[i]) == CustomTestSet
    @test tss[i].foo == 3
    for j in 1:3
        @test typeof(tss[i].results[j]) == CustomTestSet
        @test tss[i].results[j].foo == 1
        @test typeof(tss[i].results[j].results[1]) == (iseven(i+j) ? Pass : Fail)
    end
    @test typeof(tss[i].results[4]) == CustomTestSet
    @test typeof(tss[i].results[4].results[1]) == (iseven(i) ? Pass : Fail)
end

# test @inferred
function uninferrable_function(i)
    q = [1, "1"]
    return q[i]
end

@test_throws ErrorException @inferred(uninferrable_function(1))
@test @inferred(identity(1)) == 1

# Ensure @inferred only evaluates the arguments once
inferred_test_global = 0
function inferred_test_function()
    global inferred_test_global
    inferred_test_global += 1
    true
end
@test @inferred inferred_test_function()
@test inferred_test_global == 1

struct SillyArray <: AbstractArray{Float64,1} end
Base.getindex(a::SillyArray, i) = rand() > 0.5 ? 0 : false
@testset "@inferred works with A[i] expressions" begin
    @test @inferred((1:3)[2]) == 2
    test_result = @test_throws ErrorException @inferred(SillyArray()[2])
    @test contains(test_result.value.msg, "Bool")
end
# Issue #14928
# Make sure abstract error type works.
@test_throws Exception error("")

# Issue #17105
# @inferred with kwargs
function inferrable_kwtest(x; y=1)
    2x
end
function uninferrable_kwtest(x; y=1)
    2x+y
end
@test @inferred(inferrable_kwtest(1)) == 2
@test @inferred(inferrable_kwtest(1; y=1)) == 2
@test @inferred(uninferrable_kwtest(1)) == 3
@test_throws ErrorException @inferred(uninferrable_kwtest(1; y=2)) == 2

@test_throws ErrorException @testset "$(error())" for i in 1:10
end
@test_throws ErrorException @testset "$(error())" begin
end

@testset "backtraces in test errors" begin
    local f = tempname()
    write(f,
    """
    using Test
    @testset begin
        @test 1==2
        @test_throws MethodError 1
    end
    """)
    local msg = read(pipeline(ignorestatus(`$(Base.julia_cmd()) --startup-file=no --color=no $f`), stderr=DevNull), String)
    # NOTE: This test depends on the code generated by @testset getting compiled,
    # to get good backtraces. If it fails, check the implementation of `testset_beginend`.
    @test !contains(msg, "do_test(")
    @test !contains(msg, "include(")
    @test contains(msg, "at " * f * ":3")
    @test contains(msg, "at " * f * ":4")
    rm(f; force=true)
end

let io = IOBuffer()
    exc = Test.TestSetException(1,2,3,4,Vector{Union{Test.Error, Test.Fail}}())
    Base.showerror(io, exc, backtrace())
    @test !contains(String(take!(io)), "backtrace()")
end

@testset "#19750" begin
    io = IOBuffer()
    exc = Test.TestSetException(1,2,3,4,Vector{Union{Test.Error, Test.Fail}}())
    Base.showerror(io, exc, backtrace())
    @test !contains(String(take!(io)), "backtrace()")

    exc = Test.FallbackTestSetException("msg")
    Base.showerror(io, exc, backtrace())
    str = String(take!(io))
    @test contains(str, "msg")
    @test !contains(str, "backtrace()")
end

let msg = read(pipeline(ignorestatus(`$(Base.julia_cmd()) --startup-file=no --color=no -e '
        using Test

        foo(x) = length(x)^2

        @testset "Foo Tests" begin
            @testset "Animals" begin
                @testset "Felines" begin
                    @test foo("cat") == 9
                end
                @testset "Canines" begin
                    @test foo("dog") == 11
                end
            end
            @testset "Arrays" begin
                @test foo(zeros(2)) == 4
                @test foo(fill(1., 4)) == 15
            end
        end'`), stderr=DevNull), String)
    @test contains(msg,
        """
        Test Summary: | Pass  Fail  Total
        Foo Tests     |    2     2      4
          Animals     |    1     1      2
            Felines   |    1            1
            Canines   |          1      1
          Arrays      |    1     1      2
        """)
end

# 20489
let msg = split(read(pipeline(ignorestatus(`$(Base.julia_cmd()) --startup-file=no --color=no -e '
        Test.print_test_results(Test.DefaultTestSet(""))'`), stderr=DevNull), String), "\n")[1]
    @test msg == rstrip(msg)
end

@testset "test guarded srand" begin
    seed = rand(UInt)
    orig = copy(Random.GLOBAL_RNG)
    @test guardsrand(()->rand(), seed) == guardsrand(()->rand(), seed)
    @test guardsrand(()->rand(Int), seed) == guardsrand(()->rand(Int), seed)
    r1, r2 = MersenneTwister(0), MersenneTwister(0)
    a, b = guardsrand(r1) do
        srand(r1, 0)
        rand(r1), rand(r1, Int)
    end::Tuple{Float64,Int}
    c, d = guardsrand(r2) do
        srand(r2, 0)
        rand(r2), rand(r2, Int)
    end::Tuple{Float64,Int}
    @test a == c == rand(r1) == rand(r2)
    @test b == d == rand(r1, Int) == rand(r2, Int)
    @test orig == Random.GLOBAL_RNG
    @test rand(orig) == rand()
end

@testset "file info in test errors" begin
    local f = tempname()

    write(f,
    """
    using Test
    @testset begin
        @test 1==2
        @test_throws UndefVarError 1
        @test_broken 1 == 1
    end
    """)

    local msg = read(pipeline(ignorestatus(`$(Base.julia_cmd()) --startup-file=no --color=no $f`), stderr=DevNull), String)
    @test contains(msg, "at " * f * ":" * "3")
    @test contains(msg, "at " * f * ":" * "4")
    @test contains(msg, "at " * f * ":" * "5")

    rm(f; force=true)
end

# issue #24919
@testset "≈ with atol" begin
    local cmd = `$(Base.julia_cmd()) --startup-file=no --color=no`
    f(src) = read(pipeline(ignorestatus(`$cmd -e $src`), stderr=DevNull), String)

    msg = f("""
    using Test
    x, y = 0.9, 0.1
    @test x ≈ y atol=0.01
    """)
    @test contains(msg, "Evaluated: 0.9 ≈ 0.1 (atol=0.01)")

    msg = f("""
    using Test
    x, y = 0.9, 0.1
    @test x ≈ y nans=true atol=0.01
    """)
    @test contains(msg, "Evaluated: 0.9 ≈ 0.1 (nans=true, atol=0.01)")
end

@testset "@test_logs" begin
    function foo(n)
        @info "Doing foo with n=$n"
        for i=1:n
            @debug "Iteration $i"
        end
    end

    @test_logs (Info,"Doing foo with n=2") foo(2)

    # Log pattern matching
    # Regex
    @test_logs (Info,r"^Doing foo with n=[0-9]+$") foo(10)
    @test_logs (Info,r"^Doing foo with n=[0-9]+$") foo(1)
    # Level symbols
    @test_logs (:debug,) min_level=Debug @debug "foo"
    @test_logs (:info,)  @info  "foo"
    @test_logs (:warn,)  @warn  "foo"
    @test_logs (:error,) @error "foo"

    # Pass through so the value of the expression can also be tested
    @test (@test_logs (Info,"blah") (@info "blah"; 42)) == 42

    # Debug level log collection
    @test_logs (Info,"Doing foo with n=2") (Debug,"Iteration 1") (Debug,"Iteration 2") min_level=Debug foo(2)

    @test_logs (Debug,"Iteration 5") min_level=Debug match_mode=:any foo(10)

    # Test failures
    fails = @testset NoThrowTestSet "check that @test_logs detects bad input" begin
        @test_logs (Warn,) foo(1)
        @test_logs (Warn,) match_mode=:any @info "foo"
        @test_logs (Debug,) @debug "foo"
    end
    @test length(fails) == 3
    @test fails[1] isa Test.LogTestFailure
    @test fails[2] isa Test.LogTestFailure
    @test fails[3] isa Test.LogTestFailure
end

function newfunc()
    42
end
@deprecate oldfunc newfunc

@testset "@test_deprecated" begin
    @test_deprecated oldfunc()

    # Expression passthrough
    if Base.JLOptions().depwarn != 2
        @test (@test_deprecated oldfunc()) == 42

        fails = @testset NoThrowTestSet "check that @test_deprecated detects bad input" begin
            @test_deprecated newfunc()
            @test_deprecated r"Not found in message" oldfunc()
        end
        @test length(fails) == 2
        @test fails[1] isa Test.LogTestFailure
        @test fails[2] isa Test.LogTestFailure
    else
        @warn """Omitting `@test_deprecated` tests which can't yet
                 be tested in --depwarn=error mode"""
    end
end

@testset "@testset preserves GLOBAL_RNG's state, and re-seeds it" begin
    # i.e. it behaves as if it was wrapped in a `guardsrand(GLOBAL_RNG.seed)` block
    seed = rand(UInt128)
    srand(seed)
    a = rand()
    @testset begin
        # global RNG must re-seeded at the beginning of @testset
        @test a == rand()
    end
    @testset for i=1:3
        @test a == rand()
    end
    # the @testset's above must have no consequence for rand() below
    b = rand()
    srand(seed)
    @test a == rand()
    @test b == rand()
end

@testset "non AbstractTestSet as testset" begin
    local f, err = tempname(), tempname()
    write(f,
    """
    using Test
    desc = "description"
    @testset desc begin
        @test 1==1
    end
    """)
    run(pipeline(ignorestatus(`$(Base.julia_cmd()) --startup-file=no --color=no $f`), stderr=err))
    msg = read(err, String)
    @test contains(msg, "Expected `desc` to be an AbstractTestSet, it is a String")
    rm(f; force=true)
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
using dynamicRisk
using Test

@testset "BSDE" begin
    # Write your tests here.
    X0 = [0.0f0];          # initial value of stochastic state
    tspan = (0.0f0,1.0f0);
    g(X) = sum(X);         # needs to be number not [x]                                # terminal condition
    μ_f(X,p,t) = 0.0f0; #Vector d x 1                    # drift of forward
    σ_f(X,p,t) = 1.0f0; #Matrix d x d    # diffusion part of forward --> 0.5*σ^2 is the full factor

    prob2 = SDEProblem(μ_f,σ_f,X0,tspan);

    @test  dynamicRM(prob2) ≈ 1.0 atol=0.1
end

@testset "Tree" begin
    # Write your tests here.
    trr = Tree(402)
    @test AVaR(trr::Tree,alpha::Float64)[1] <= maximum(trr.state)

end

ext: jl


text:
using Test
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
splitln(str) = split(strip(str), "\n")
import ExtensibleEffects: autorun

# Vector/Iterable
# ---------------

program3 = @syntax_eff_noautorun begin
  r0 = [1,4]
  [r0, r0, r0]
end
r3 = runhandlers(Vector, program3)
@test r3 == [1,1,1,4,4,4]
@test autorun(program3) == [1,1,1,4,4,4]

iterables = @syntax_eff begin
  r0 = Iterable([1,4])
  Iterable([r0, r0, r0])
end
@test collect(iterables) == [1,1,1,4,4,4]

program4 = @syntax_eff_noautorun begin
  a = [1,2,3]
  b = iftrue(a % 2 == 0) do
    42
  end
  [b, a+b]
end

r4_1 = runhandlers((Const, Identity, Vector), program4)
r4_2 = runhandlers((Vector, Const, Identity), program4)

@test r4_2 == program4 |>
  eff -> ExtensibleEffects.runhandler(Identity, eff) |>
  eff -> ExtensibleEffects.runhandler(Const, eff) |>
  eff -> ExtensibleEffects.runhandler(Vector, eff) |>
  ExtensibleEffects.runlast

@test r4_1 == Const(nothing)
@test r4_2 == [Const(nothing), Identity(42), Identity(44), Const(nothing)]
@test autorun(program4) == [Const(nothing), Identity(42), Identity(44), Const(nothing)]

# test syntax
wrapper(i::Int) = collect(1:i)
wrapper(any) = any

r5_1 = @runhandlers (Vector, Const, Identity) @syntax_eff_noautorun wrapper begin
  a = 3
  b = iftrue(a % 2 == 0) do
    42
  end
  [b, a+b]
end
@test r5_1 == [Const(nothing), Identity(42), Identity(44), Const(nothing)]

r5_2 = @syntax_eff_noautorun wrapper begin
  a = 3
  b = iftrue(a % 2 == 0) do
    42
  end
  [b, a+b]
end
@test runhandlers((Vector, Const, Identity), r5_2) ==  [Const(nothing), Identity(42), Identity(44), Const(nothing)]

r5_3 = @runhandlers (Vector, Const, Identity) @syntax_eff_noautorun begin
  a = Identity(4)
  b = iftrue(a % 2 == 0) do
    42
  end
  [b, a+b]
end
@test r5_3 == [Identity(42), Identity(46)]

r5_4 = @runhandlers (Const, Identity, Vector) @syntax_eff_noautorun begin
  a = [1,2,3]
  b = iftrue(a % 2 == 0) do
    42
  end
  [b, a+b]
end
@test r5_4 == Const(nothing)

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
using Symbolics

function lotka_volterra!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

@variables t du[1:2] u[1:2] p[1:4]
lotka_volterra!(du, u, p, t)
D = Differential(t)
f = build_function((@. D(u) ~ du), u, p, t, target=Symbolics.JuliaTarget())

Base.remove_linenums!(f[1])


text:
using OmegaCore
using Test

function testvar()
  n = ~ finite(1:10)
  μ = [~ Normal(0, 1) for i = 1:n]
  s = sum(μ)
end

function testmaxmean()
  x = 1 ~ normal(0, 1)
  y = 2 ~ bounded(-4.0, 1.0)
  z(ω) = x(ω) + y(ω)

  # Find y which maximises entropy of x
  argmax(y, mean(z))
end

function test_polynomial()
  function polynomial(ω)
    n = finite(ω, 1:10)
    function (x)
      y = 0.0
      for i = 1:n
        y += ~ normal(ω, 0, 1) * x^i
      end
      y
    end
  end
end

@testset "free-variables" begin
  testvar()
  testmaxmean()
  test_polynomial()
end


text:
using Echidna
using Test

# test for compare_pareto_dominance
nobjs = 2
nvars = 30
zdt1_problem = Problem(
    2, [false for i in 1:nobjs], nvars, [MOGA_Real(0.0, 1.0) for i in 1:30],
    ZDT1
)

objective_vals = [0.,3.], [1., 1.], [2., 2.], [3., 0.]
sols = [ Solution(
    zdt1_problem, [.1 for i in 1:nvars], false, ov, 0.0, 0
) for ov in objective_vals]

expectation = [[0,0,0,0] [0,0,-1,0] [0,1,0,0] [0,0,0,0]]
for i=1:4
    for j=1:4
        @test compare_pareto_dominance(sols[i], sols[j]) == expectation[j, i]
    end
end

a = Archive(compare_pareto_dominance, Vector{Solution}())
insert_solutions!(a, sols)
@test length(a.solutions) == 3

# SWITCH OBJECTIVE GOALS
zdt1_problem = Problem(
    2, [true, false], nvars, [MOGA_Real(0.0, 1.0) for i in 1:30], ZDT1
)
sols = [ Solution(
    zdt1_problem, [.1 for i in 1:nvars], false, ov, 0.0, 0
) for ov in objective_vals]

expectation = [[0,1,1,1] [-1,0,0,1] [-1,0,0,1] [-1,-1,-1,0]]
for i=1:4
    for j=1:4
        @test compare_pareto_dominance(sols[i], sols[j]) == expectation[j, i]
    end
end

a = Archive(compare_pareto_dominance, Vector{Solution}())
insert_solutions!(a, sols)
@test length(a.solutions) == 1

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
module BookTest

using Reexport: @reexport
@reexport begin
using Books:
    build_all,
    gen
using DataFrames:
    DataFrame,
    filter!,
    filter,
    select!,
    select
end # @reexport

using BenchmarkTools
using StaticArrays

include("data.jl")
include("optimizing-serial-code.jl")

export example_dataframe
export inner_rows!
export inner_cols!

"""
    build()

This function is called during CI.
"""
function build()
    println("Building Example Book")
    # To avoid publishing broken websites.
    fail_on_error = true
    gen(; fail_on_error)
    build_all(; fail_on_error)
end

end # module

max_stars_repo_name: jvaverka/BookTest

max_issues_repo_name: jvaverka/BookTest

max_forks_repo_name: jvaverka/BookTest


text:
"Allocate ODBC handles for interacting with the ODBC Driver Manager"
function ODBCAllocHandle(handletype, parenthandle)
    handle = Ref{Ptr{Void}}()
    ODBC.API.SQLAllocHandle(handletype,parenthandle,handle)
    handle = handle[]
    if handletype == ODBC.API.SQL_HANDLE_ENV
        ODBC.API.SQLSetEnvAttr(handle,ODBC.API.SQL_ATTR_ODBC_VERSION,ODBC.API.SQL_OV_ODBC3)
    end
    return handle
end

"Alternative connect function that allows user to create datasources on the fly through opening the ODBC admin"
function ODBCDriverConnect!(dbc::Ptr{Void},conn_string,driver_prompt::UInt16)
    window_handle = C_NULL
    @static if is_windows()
        window_handle = ccall((:GetForegroundWindow, :user32), Ptr{Void}, () )
        driver_prompt = ODBC.API.SQL_DRIVER_PROMPT
    end
    out_conn = Block(ODBC.API.SQLWCHAR,BUFLEN)
    out_buff = Ref{Int16}()
    @CHECK dbc ODBC.API.SQL_HANDLE_DBC ODBC.API.SQLDriverConnect(dbc,window_handle,conn_string,out_conn.ptr,BUFLEN,out_buff,driver_prompt)
    connection_string = String(out_conn,out_buff[])
    # free!(out_conn)
    return connection_string
end

"`ODBC.execute!` is a minimal method for just executing an SQL `query` string. No results are checked for or returned."
function execute!(dsn::DSN, query::AbstractString)
    stmt = dsn.stmt_ptr
    ODBC.ODBCFreeStmt!(stmt)
    ODBC.@CHECK stmt ODBC.API.SQL_HANDLE_STMT ODBC.API.SQLExecDirect(stmt, query)
    return
end

"""
`ODBC.Source` constructs a valid `Data.Source` type that executes an SQL `query` string for the `dsn` ODBC DSN.
Results are checked for and an `ODBC.ResultBlock` is allocated to prepare for fetching the resultset.
"""
function Source(dsn::DSN, query::AbstractString)
    stmt = dsn.stmt_ptr
    ODBC.ODBCFreeStmt!(stmt)
    ODBC.@CHECK stmt ODBC.API.SQL_HANDLE_STMT ODBC.API.SQLExecDirect(stmt, query)
    rows, cols = Ref{Int}(), Ref{Int16}()
    ODBC.API.SQLNumResultCols(stmt,cols)
    ODBC.API.SQLRowCount(stmt,rows)
    rows, cols = rows[], cols[]
    #Allocate arrays to hold each column's metadata
    cnames = Array(String,cols)
    ctypes, csizes = Array(ODBC.API.SQLSMALLINT,cols), Array(ODBC.API.SQLULEN,cols)
    cdigits, cnulls = Array(ODBC.API.SQLSMALLINT,cols), Array(ODBC.API.SQLSMALLINT,cols)
    juliatypes = Array(DataType,cols)
    alloctypes = Array(DataType,cols)
    indcols = Array(Vector{ODBC.API.SQLLEN},cols)
    columns = Array(ODBC.Block,cols)
    longtext = false
    #Allocate space for and fetch the name, type, size, etc. for each column
    len, dt, csize = Ref{ODBC.API.SQLSMALLINT}(), Ref{ODBC.API.SQLSMALLINT}(), Ref{ODBC.API.SQLULEN}()
    digits, null = Ref{ODBC.API.SQLSMALLINT}(), Ref{ODBC.API.SQLSMALLINT}()
    cname = ODBC.Block(ODBC.API.SQLWCHAR, ODBC.BUFLEN)
    for x = 1:cols
        ODBC.API.SQLDescribeCol(stmt, x, cname.ptr, ODBC.BUFLEN, len, dt, csize, digits, null)
        cnames[x] = String(cname, len[])
        t = dt[]
        ctypes[x], csizes[x], cdigits[x], cnulls[x] = t, csize[], digits[], null[]
        alloctypes[x], juliatypes[x], islongtext = ODBC.API.SQL2Julia[t]
        longtext |= islongtext
    end
    # ODBC.free!(cname)
    # Determine fetch strategy
    # rows might be -1 (dbms doesn't return total rows in resultset), 0 (empty resultset), or 1+
    if rows > -1 # known # of rows
        if longtext # longtext column types present in resultset
            rowset = min(1,rows) # in case of rows == 0, empty resultset
            multifetch = rows != 1 # we don't need to multifetch if there's only one row
        else
            if rows < ODBC.API.MAXFETCHSIZE
                rowset = rows # rowset is >= 0 and < MAXFETCHSIZE
                multifetch = false
            else
                rowset = ODBC.API.MAXFETCHSIZE
                multifetch = true
            end
        end
    else
        # unknown # of rows
        multifetch = true
        rowset = longtext ? 1 : ODBC.API.MAXFETCHSIZE
    end
    ODBC.API.SQLSetStmtAttr(stmt, ODBC.API.SQL_ATTR_ROW_ARRAY_SIZE, rowset, ODBC.API.SQL_IS_UINTEGER)
    for x = 1:cols
        block = ODBC.Block(alloctypes[x], rowset, csizes[x]+1)
        ind = Array(ODBC.API.SQLLEN,rowset)
        ODBC.API.SQLBindCols(stmt,x,ODBC.API.SQL2C[ctypes[x]],block.ptr,block.elsize,ind)
        columns[x], indcols[x] = block, ind
    end
    schema = Data.Schema(cnames, juliatypes, rows,
        Dict("types"=>[ODBC.API.SQL_TYPES[c] for c in ctypes], "sizes"=>csizes, "digits"=>cdigits, "nulls"=>cnulls))
    rowsfetched = Ref{ODBC.API.SQLLEN}() # will be populated by call to SQLFetchScroll
    ODBC.API.SQLSetStmtAttr(stmt,ODBC.API.SQL_ATTR_ROWS_FETCHED_PTR,rowsfetched,ODBC.API.SQL_NTS)
    rb = ODBC.ResultBlock(columns,indcols,alloctypes,rowset,rowsfetched)
    return ODBC.Source(schema,cols,dsn,query,rb,ODBC.API.SQLFetchScroll(stmt,ODBC.API.SQL_FETCH_NEXT,0))
end

# Fetch Strategies
  #1. Reinterpet Blocks that were allocated and bound
    # isbits arrays: let Julia take ownership of Block memory; zero out the Block
    # mutable arrays: keep references to Block in `other`
  #2. Copy allocated and bound Blocks to pre-allocated output and reinterpret final output
    # isbits arrays: copy Block memory to dest array, free allocated Blocks
    # mutable arrays: keep references to copied Blocks in `other`; free allocated Blocks
  #3. Grow output by rows fetched, copy allocated and bound Blocks to newly grown output, reinterpret final output
    # isbits arrays: copy Block memory to dest array, free allocated Blocks
    # mutable arrays: keep references to copied Blocks in `other`
# Fetch Strategy Determination
  # if known # of rows and # of rows < MAXFETCHSIZE and no LONGTEXT
    # then set rowset = rows, do one SQLFetchScroll, fetch strategy #1
  # if known # of rows and # of rows > MAXFETCHSIZE and no LONGTEXT
    # then set rowset = MAXFETCHSIZE, do multiple SQLFetchScroll, fetch strategy #2
  # if known # of rows and # of rows < or > MAXFETCHSIZE and LONGTEXT
    # then set rowset = 1, do multiple SQLFetchScroll, fetch strategy #2
  # if unknown # of rows and no LONGTEXT
    # then set rowset = MAXFETCHSIZE, do one or more SQLFetchScroll, fetch strategy #3
  # if unknown # of rows and LONGTEXT
    # then set rowset = 1, do multiple SQLFetchScroll, fetch strategy #3
  # if # of rows = 0 then done

"Checks if an `ODBC.Source` has finished fetching results from an executed query string"
Data.isdone(source::ODBC.Source, x, y) = source.status != ODBC.API.SQL_SUCCESS && source.status != ODBC.API.SQL_SUCCESS_WITH_INFO
# Data.reset!(io::ODBC.Source) = nothing

Data.streamtype{T<:ODBC.Source}(::Type{T}, ::Type{Data.Column}) = true
# Data.streamtype{T<:ODBC.Source}(::Type{T}, ::Type{Data.Field}) = true

# function Data.getfield{T}(source::ODBC.Source, ::Type{T}, row, col)
#     if !isdefined(source.columns, col)
#         source.columns[col] = Data.getcolumn(source, T, col)
#     end
#     return source.columns[col][row]
# end

function booleanize!(ind::Vector{ODBC.API.SQLLEN}, new::Vector{Bool}, len)
    @simd for i = 1:len
        @inbounds new[i] = ind[i] == ODBC.API.SQL_NULL_DATA
    end
    return new
end

# bitstype/immutable types
function Data.getcolumn{T}(source::ODBC.Source, ::Type{T}, i)
    rb = source.rb;
    len = rb.rowsfetched[]
    dest = NullableArray(T, len)
    ccall(:memcpy, Void, (Ptr{T}, Ptr{T}, Csize_t), pointer(dest.values), rb.columns[i].ptr, len * sizeof(T))
    booleanize!(rb.indcols[i], dest.isnull, len)
    if i == source.cols && !Data.isdone(source, 1, 1)
        source.status = ODBC.API.SQLFetchScroll(source.dsn.stmt_ptr, ODBC.API.SQL_FETCH_NEXT, 0)
    end
    return dest
end

# Dec64 type
const DECZERO = Dec64(0)
function cast{T}(::Type{Dec64}, ptr::Ptr{T}, len)
    len <= 0 && return DECZERO
    bytes = unsafe_wrap(Array, ptr, ODBC.bytes2codeunits(T, len))
    str = String(Base.transcode(UInt8, bytes))
    return parse(Dec64, str)
end

"fill a NullableVector by copying the data from a Block that has Dec64 data"
function Data.getcolumn(source::ODBC.Source, ::Type{Dec64}, i)
    rb = source.rb;
    len = rb.rowsfetched[]
    dest = NullableArray(Dec64, len)
    values = dest.values
    isnull = dest.isnull
    cur = rb.columns[i].ptr
    elsize = rb.columns[i].elsize
    for x = 1:len
        @inbounds values[x] = cast(Dec64, cur, rb.indcols[i][x])
        @inbounds isnull[x] = rb.indcols[i][x] == ODBC.API.SQL_NULL_DATA
        cur += elsize
    end
    if i == source.cols && !Data.isdone(source, 1, 1)
        source.status = ODBC.API.SQLFetchScroll(source.dsn.stmt_ptr, ODBC.API.SQL_FETCH_NEXT, 0)
    end
    return dest
end

# String/Vector{UInt8} types
# translate a # of bytes and a code unit type (UInt8, UInt16, UInt32) and return the # of code units; returns 0 if field is null
bytes2codeunits(::Type{UInt8},  bytes) = ifelse(bytes == ODBC.API.SQL_NULL_DATA, 0, Int(bytes))
bytes2codeunits(::Type{UInt16}, bytes) = ifelse(bytes == ODBC.API.SQL_NULL_DATA, 0, Int(bytes >> 1))
bytes2codeunits(::Type{UInt32}, bytes) = ifelse(bytes == ODBC.API.SQL_NULL_DATA, 0, Int(bytes >> 2))

cast{T}(::Type{WeakRefString{T}}, ptr, len, x) = len <= 0 ? WeakRefString{T}(Ptr{T}(0), 0, 0) : WeakRefString{T}(ptr, ODBC.bytes2codeunits(T, len), x)
cast(::Type{Vector{UInt8}}, ptr, len, x) = len <= 0 ? UInt8[] : unsafe_wrap(Array, ptr, len, false)

"fill a NullableVector by copying the data from a Block that has container-type"
function Data.getcolumn{T<:Union{Vector{UInt8},AbstractString}}(source::ODBC.Source, ::Type{T}, i)
    rb = source.rb;
    len = rb.rowsfetched[]
    values = Vector{T}(len)
    isnull = Vector{Bool}(len)
    totalbytes = max(0, len == 1 ? rb.indcols[i][1] : len * rb.columns[i].elsize)
    parent = Vector{UInt8}(totalbytes)
    unsafe_copy!(pointer(parent), convert(Ptr{UInt8}, rb.columns[i].ptr), totalbytes)

    cur = pointer(parent)
    ind = 1
    elsize = rb.columns[i].elsize
    for x = 1:len
        @inbounds values[x] = ODBC.cast(T, cur, rb.indcols[i][x], ind)
        @inbounds isnull[x] = rb.indcols[i][x] == ODBC.API.SQL_NULL_DATA
        cur += elsize
        ind += elsize
    end
    if i == source.cols && !Data.isdone(source, 1, 1)
        source.status = ODBC.API.SQLFetchScroll(source.dsn.stmt_ptr, ODBC.API.SQL_FETCH_NEXT, 0)
    end
    return NullableArray{T,1}(values, isnull, parent)
end

# """
# Stream the results of `source` (if any) to a `DataFrame` type. The `DataFrame` will be allocated according
# to the size of the resultset of the query string.
# """
# function Data.stream!(source::ODBC.Source, ::Type{DataFrame})
#     rb = source.rb;
#     rows, cols = size(source)
#     if rb.fetchsize == rows
#         # fetch strategy #1: reinterpret allocated Blocks
#         data = Array(Any, cols)
#         for col = 1:cols
#             data[col] = NullableArray(rb.jltypes[col],rb.columns[col],rb.indcols[col],rows)
#         end
#         return DataFrame(data, map(Symbol, Data.header(source)))
#     else
#         # fetch strategy #2 or #3
#         df = DataFrame(Data.Schema(Data.header(source), Data.types(source), max(0,rows)))
#         return Data.stream!(source, df)
#     end
# end
#
# "Stream the results of an `ODBC.Source` to a pre-allocated `DataFrame`"
# function Data.stream!(source::ODBC.Source, df::DataFrame)
#     rb = source.rb;
#     rows, cols = size(source)
#     data = df.columns;
#     r = 0
#     while true
#         rowsfetched::ODBC.API.SQLLEN = rb.rowsfetched[]
#         (rowsfetched == 0 || rowsfetched > rb.fetchsize) && break
#         if rows < 0
#             # fetch strategy #3: grow our output and copy Blocks to new output space until done
#             for col = 1:cols
#                 ODBC.append!(rb.jltypes[col],rb.columns[col],rb.indcols[col],data[col],r,rowsfetched)
#             end
#         else
#             # fetch strategy #2: copy allocated Blocks to pre-allocated output
#             for col = 1:cols
#                 #TODO: add a check that we have enough space to copy rowsfetched into data[col] from r
#                 ODBC.copy!(rb.jltypes[col],rb.columns[col],rb.indcols[col],data[col],r,rowsfetched)
#             end
#         end
#         r += rowsfetched
#         Data.isdone(source) && break
#         source.status = ODBC.API.SQLFetchScroll(source.dsn.stmt_ptr,ODBC.API.SQL_FETCH_NEXT,0)
#     end
#     for col = 1:cols
#         free!(rb.columns[col])
#     end
#     source.schema.rows = r
#     return df
# end
#
# function getfield!{T}(jltype::Type{T},block::Block,ind,row,col,cols,sink,null)
#     val = ODBC.getfield(jltype, block, row, ind)
#     CSV.writefield(sink, ind == ODBC.API.SQL_NULL_DATA ? null : val, col, cols)
# end
#
# "Stream the results of an `ODBC.Source` directly to a CSV file, represented by `sink::CSV.Sink`"
# function Data.stream!(source::ODBC.Source, sink::CSV.Sink,header::Bool=true)
#     header && CSV.writeheaders(source,sink)
#     rb = source.rb
#     null = sink.options.null
#     rows, cols = size(source)
#     r = 0
#     while true
#         rowsfetched::ODBC.API.SQLLEN = rb.rowsfetched[]
#         rowsfetched == 0 && break
#         for row = 1:rowsfetched, col = 1:cols
#             ind::ODBC.API.SQLLEN = rb.indcols[col][row]
#             ODBC.getfield!(rb.jltypes[col],rb.columns[col],ind,row,col,cols,sink,null)
#         end
#         r += rowsfetched
#         Data.isdone(source) && break
#         source.status = ODBC.API.SQLFetchScroll(source.dsn.stmt_ptr,ODBC.API.SQL_FETCH_NEXT,0)
#     end
#     for col = 1:cols
#         ODBC.free!(rb.columns[col])
#     end
#     source.schema.rows = r
#     sink.schema = source.schema
#     close(sink)
#     return sink
# end
#
# function getbind!{T}(jltype::Type{T},block::Block,ind,row,col,stmt)
#     val = getfield(jltype, block, row, ind)::T
#     if ind == ODBC.API.SQL_NULL_DATA
#         SQLite.bind!(stmt,col,SQLite.NULL)
#     else
#         SQLite.bind!(stmt,col,val)
#     end
# end
#
# "Stream the results of an `ODBC.Source` directly to an SQLite table, represented by `sink::SQLite.Sink`"
# function Data.stream!(source::ODBC.Source, sink::SQLite.Sink)
#     rb = source.rb
#     rows, cols = size(source)
#     stmt = sink.stmt
#     handle = stmt.handle
#     SQLite.transaction(sink.db) do
#         r = 0
#         while true
#             rowsfetched::ODBC.API.SQLLEN = rb.rowsfetched[]
#             for row = 1:rowsfetched
#                 for col = 1:cols
#                     getbind!(rb.jltypes[col],rb.columns[col],rb.indcols[col][row]::ODBC.API.SQLLEN,row,col,stmt)
#                 end
#                 SQLite.sqlite3_step(handle)
#                 SQLite.sqlite3_reset(handle)
#             end
#             r += rows
#             Data.isdone(source) && break
#             source.status = ODBC.API.SQLFetchScroll(source.dsn.stmt_ptr,ODBC.API.SQL_FETCH_NEXT,0)
#         end
#     end
#     for col = 1:cols
#         free!(rb.columns[col])
#     end
#     SQLite.execute!(sink.db,"analyze $(sink.tablename)")
#     return sink
# end

"""
Convenience method that constructs an `ODBC.Source` and streams the results to `sink`.
`sink` can be any valid `Data.Sink` (`CSV.Sink`,`SQLite.Sink`,etc.) and by default, is a `DataFrame`.
"""
function query(dsn::DSN, querystring::AbstractString, sink=DataFrame)
    source = ODBC.Source(dsn, querystring)
    return Data.stream!(source, sink)
end

"Convenience string macro for executing an SQL statement against a DSN."
macro sql_str(s,dsn)
    query(dsn,s)
end


text:
## Methods for canonical, varfunc and devresid2 for distributions in the exponential family

@doc """
Returns the canonical Link type for an instance of a distribution in the exponential family
""" ->
canonical(::Bernoulli) = LogitLink()
canonical(::Binomial) = LogitLink()
canonical(::Gamma) = InverseLink()
canonical(::Normal) = IdentityLink()
canonical(::Poisson) = LogLink()

@doc """
Evaluate `y*log(y/μ)` with the correct limit as `y` approaches zero from above
"""->
ylogydμ{T<:AbstractFloat}(y::T,μ::T) = y > zero(T) ? y*log(y/μ) : zero(T)

varfunc(::Bernoulli,μ) = μ*(one(μ)-μ)
varfunc(::Binomial,μ) = μ*(one(μ)-μ)
varfunc(::Gamma,μ) = abs2(μ)
varfunc(::Normal,μ) = one(μ)
varfunc(::Poisson,μ) = μ

two(y) = one(y) + one(y)                # equivalent to convert(typeof(y),2)

@doc """
Evaluate the squared deviance residual for a distribution instance and values of `y` and `μ`
"""->
devresid2(::Bernoulli,y,μ) = two(y)*(ylogydμ(y,μ) + ylogydμ(one(y)-y,one(μ)-μ))
devresid2(::Binomial,y,μ) = devresid2(Bernoulli(),y,μ)
devresid2(::Gamma,y,μ) =  two(y)*((y-μ)/μ - (y == zero(y) ? y : log(y/μ)))
devresid2(::Normal,y,μ) = abs2(y-μ)
devresid2(::Poisson,y,μ) = two(y)*(ylogydμ(y,μ)-(y-μ))

@doc """
Initial μ value from the y and wt
""" ->
mustart{T<:AbstractFloat}(::Bernoulli,y::T,wt::T) = (wt*y + convert(T,0.5))/(wt + one(T))
mustart{T<:AbstractFloat}(::Binomial,y::T,wt::T) = (wt*y + convert(T,0.5))/(wt + one(T))
mustart(::Gamma,y,wt) = y
mustart(::Normal,y,wt) = y
mustart{T<:AbstractFloat}(::Poisson,y::T,wt::T) = convert(T,1.1)*y


text:
function DiffEqBase.__solve(prob::QuadraticStateSpaceProblem{uType,uPriorType,tType,P,NP,F,A0Type,
                                                             A1Type,A2Type,BType,C0Type,C1Type,
                                                             C2Type,RType,ObsType,K},
                            alg::DirectIteration, args...;
                            kwargs...) where {uType,uPriorType,tType,P,NP,F,A0Type,A1Type,A2Type,
                                              BType,C0Type,C1Type,C2Type,RType,ObsType,K}
    T = convert(Int64, prob.tspan[2] - prob.tspan[1] + 1)
    noise = get_concrete_noise(prob, prob.noise, prob.B, T - 1)  # concrete noise for simulations as required.
    # checks on bounds
    @assert size(noise, 1) == size(prob.B, 2)
    @assert size(noise, 2) == T - 1
    @assert maybe_check_size(prob.observables, 2, T - 1)

    @unpack A_0, A_1, A_2, B, C_0, C_1, C_2 = prob

    C_2_vec = [C_2[i, :, :] for i in 1:size(C_2, 1)] # should be the native datastructure
    A_2_vec = [A_2[i, :, :] for i in 1:size(A_2, 1)] # should be the native datastructure

    u_f = [zero(prob.u0) for _ in 1:T]  # TODO: move to internal algorithm cache
    u = [zero(prob.u0) for _ in 1:T] # TODO: move to internal algorithm cache
    z = [zero(C_0) for _ in 1:T] # TODO: move to internal algorithm cache

    u[1] .= prob.u0
    u_f[1] .= prob.u0
    z[1] .= C_0
    mul!(z[1], C_1, prob.u0, 1, 1)
    quad_muladd!(z[1], C_2_vec, prob.u0) #z0 .+= quad(C_2, prob.u0)

    loglik = 0.0
    @inbounds @views for t in 2:T
        mul!(u_f[t], A_1, u_f[t - 1])
        mul!(u_f[t], B, view(noise, :, t - 1), 1, 1)

        u[t] .= A_0
        mul!(u[t], A_1, u[t - 1], 1, 1)
        quad_muladd!(u[t], A_2_vec, u_f[t - 1]) # u[t] .+= quad(A_2, u_f[t - 1])
        mul!(u[t], B, view(noise, :, t - 1), 1, 1)

        z[t] .= C_0
        mul!(z[t], C_1, u[t], 1, 1)
        quad_muladd!(z[t], C_2_vec, u_f[t]) # z[t] .+= quad(C_2, u_f[t])
        loglik += maybe_logpdf(prob.observables_noise, prob.observables, t - 1, z[t])
    end

    t_values = prob.tspan[1]:prob.tspan[2]
    return build_solution(prob, alg, t_values, u; W = noise,
                          logpdf = ObsType <: Nothing ? nothing : loglik, retcode = :Success)
end

function ChainRulesCore.rrule(::typeof(DiffEqBase.solve), prob::QuadraticStateSpaceProblem,
                              alg::DirectIteration, args...; kwargs...)
    T = convert(Int64, prob.tspan[2] - prob.tspan[1] + 1)
    noise = get_concrete_noise(prob, prob.noise, prob.B, T - 1)  # concrete noise for simulations as required.
    @assert !isnothing(prob.noise)  # need to have concrete noise for this simple method
    # checks on bounds
    @assert size(noise, 1) == size(prob.B, 2)
    @assert maybe_check_size(prob.observables, 2, T - 1)
    @assert size(noise, 2) == T - 1

    @unpack A_0, A_1, A_2, B, C_0, C_1, C_2 = prob

    C_2_vec = [C_2[i, :, :] for i in 1:size(C_2, 1)] # should be the native datastructure
    A_2_vec = [A_2[i, :, :] for i in 1:size(A_2, 1)] # should be the native datastructure

    u_f = [zero(prob.u0) for _ in 1:T]  # TODO: move to internal algorithm cache
    u = [zero(prob.u0) for _ in 1:T] # TODO: move to internal algorithm cache
    z = [zero(C_0) for _ in 1:T] # TODO: move to internal algorithm cache

    u[1] .= prob.u0
    u_f[1] .= prob.u0
    z[1] .= C_0
    mul!(z[1], C_1, prob.u0, 1, 1)
    quad_muladd!(z[1], C_2_vec, prob.u0) #z0 .+= quad(C_2, prob.u0)

    loglik = 0.0
    @inbounds @views for t in 2:T
        mul!(u_f[t], A_1, u_f[t - 1])
        mul!(u_f[t], B, view(noise, :, t - 1), 1, 1)

        u[t] .= A_0
        mul!(u[t], A_1, u[t - 1], 1, 1)
        quad_muladd!(u[t], A_2_vec, u_f[t - 1]) # u[t] .+= quad(A_2, u_f[t - 1])
        mul!(u[t], B, view(noise, :, t - 1), 1, 1)

        z[t] .= C_0
        mul!(z[t], C_1, u[t], 1, 1)
        quad_muladd!(z[t], C_2_vec, u_f[t]) # z[t] .+= quad(C_2, u_f[t])
        loglik += logpdf(prob.observables_noise, view(prob.observables, :, t - 1) - z[t])
    end
    t_values = prob.tspan[1]:prob.tspan[2]
    sol = build_solution(prob, alg, t_values, u; W = noise, logpdf = loglik, retcode = :Success)

    function solve_pb(Δsol)
        # Currently only changes in the logpdf are supported in the rrule
        @assert Δsol.u == ZeroTangent()
        @assert Δsol.W == ZeroTangent()

        Δlogpdf = Δsol.logpdf
        if iszero(Δlogpdf)
            return (NoTangent(), Tangent{typeof(prob)}(), NoTangent(),
                    map(_ -> NoTangent(), args)...)
        end
        ΔA_0 = zero(A_0)
        ΔA_1 = zero(A_1)
        ΔA_2_vec = [zero(A) for A in A_2_vec] # should be native datastructure
        ΔA_2 = zero(A_2)

        ΔB = zero(B)
        ΔC_0 = zero(C_0)
        ΔC_1 = zero(C_1)
        ΔC_2_vec = [zero(A) for A in C_2_vec] # should be native datastructure
        ΔC_2 = zero(C_2)
        Δu_f_sum = zero(u[1])

        Δnoise = similar(noise)
        Δu = [zero(prob.u0) for _ in 1:T]
        Δu_f = [zero(prob.u0) for _ in 1:T]
        A_2_vec_sum = [(A + A') for A in A_2_vec] # prep the sum since we will use it repeatedly
        C_2_vec_sum = [(A + A') for A in C_2_vec] # prep the sum since we will use it repeatedly

        @views @inbounds for t in T:-1:2
            Δz = Δlogpdf * (view(prob.observables, :, t - 1) - z[t]) ./
                 diag(prob.observables_noise.Σ) # More generally, it should be Σ^-1 * (z_obs - z)

            # inplace adoint of quadratic form with accumulation
            quad_muladd_pb!(ΔC_2_vec, Δu_f[t], Δz, C_2_vec_sum, u_f[t])
            mul!(Δu[t], C_1', Δz, 1, 1)

            quad_muladd_pb!(ΔA_2_vec, Δu_f[t - 1], Δu[t], A_2_vec_sum, u_f[t - 1])
            mul!(Δu[t - 1], A_1', Δu[t])
            mul!(Δu_f[t - 1], A_1', Δu_f[t], 1, 1)

            # Δu_f_sum = Δu[t] .+ Δu_f[t]
            copy!(Δu_f_sum, Δu[t])
            Δu_f_sum .+= Δu_f[t]

            mul!(view(Δnoise, :, t - 1), B', Δu_f_sum)
            # Now, deal with the coefficients
            ΔA_0 += Δu[t]
            mul!(ΔA_1, Δu[t], u[t - 1]', 1, 1)
            mul!(ΔA_1, Δu_f[t], u_f[t - 1]', 1, 1)
            mul!(ΔB, Δu_f_sum, view(noise, :, t - 1)', 1, 1)
            ΔC_0 += Δz
            mul!(ΔC_1, Δz, u[t]', 1, 1)
        end

        # Remove once the vector of matrices or column-major organized 3-tensor is the native datastructure for C_2/A_2
        @views @inbounds for (i, ΔA_2_slice) in enumerate(ΔA_2_vec)
            ΔA_2[i, :, :] .= ΔA_2_slice
        end
        @views @inbounds for (i, ΔC_2_slice) in enumerate(ΔC_2_vec)
            ΔC_2[i, :, :] .= ΔC_2_slice
        end

        return (NoTangent(),
                Tangent{typeof(prob)}(; A_0 = ΔA_0, A_1 = ΔA_1, A_2 = ΔA_2, B = ΔB, C_0 = ΔC_0,
                                      C_1 = ΔC_1, C_2 = ΔC_2, u0 = Δu[1] + Δu_f[1], noise = Δnoise,
                                      observables = NoTangent(), # not implemented
                                      observables_noise = NoTangent()), NoTangent(),
                map(_ -> NoTangent(), args)...)
    end
    return sol, solve_pb
end


text:
using PastaQ
using Random

# Set the random seed so the results are the same
# each time it is run
Random.seed!(1234)

N = 4

# Number to QFT
number = 10
number_bin = reverse(digits(number; base=2, pad=N))
number_bin_st = prod(string.(number_bin))

println("Number of qubits: ", N)
println("Number for QFT: ", number, " -> ", number_bin)
println()

#
# Make state |1010⟩
#

# Start with state |0000⟩
ψ0 = qubits(N)
# Apply gates X_1, X_3
starting_state = [number_bin[n] == 1 ? ("X", n) : ("I", n) for n in 1:N]

ψ0 = runcircuit(ψ0, starting_state)

samples = getsamples(ψ0, 5)
println("Sample from the initial state |$(number_bin_st)⟩:")
display(samples)
println()

#
# Make the QFT circuit
#

circuit = qft(N)

println("QFT circuit gates:")
display(circuit)
println()

println("Running QFT...")
println()
ψ = runcircuit(ψ0, circuit)

println("Sample from QFT|$(number_bin_st)⟩:")
samples = getsamples(ψ, 5)
display(samples)
println()

#
# Run inverse QFT
#

println("Running inverse QFT...")
println()
circuit⁻¹ = qft(N; inverse=true)
ψ = runcircuit(ψ, circuit⁻¹)

println("Sample from QFT⁻¹QFT|$(number_bin_st)⟩:")
samples = getsamples(ψ, 5)
display(samples)


text:
module GWAlpha_ML_module

using Distributed
using SharedArrays
using DelimitedFiles
using Distributions
using LinearAlgebra
using Optim
using ProgressMeter
include("significance_testing_module.jl")

"""
# _____________________________________________________________________________________________________________________________________________
# Negative log-likelihood of the Beta distribution parameters (shape1 and shape2) given the distribution of the allele frequencies across pools

`neg_log_likelihood_cdfbeta(beta::Array{Float64,1}, data_A::Array{Float64,1}, data_B::Array{Float64,1})`

# Inputs
1. beta [Array{Float64,1}]: 4-element vector of the beta distribution shape parameters, 2 for each of the data distributions
2. data_A [Array{Float64,1}; xi ∉ 0]: allele frequency distribution across pools
3. data_B [Array{Float64,1}; xi ∉ 0]: additive inverse of the allele frequency distribution across pools

# Output
1. Negative log-likelihood of the Beta distribution parameters (shape1 and shape2) given the distribution of the allele frequencies across pools

# Examples
```
beta = [1.0, 1.0, 0.5, 1.0]
f = rand(5); bins = repeat([0.2], inner=5); p = 0.5
data_A = cumsum( (f .* bins ./ p) )
data_B = cumsum( ((1 .- f) .* bins ./ (1 .- p)) )
GWAlpha_ML_module.neg_log_likelihood_cdfbeta(beta, data_A, data_B)
```
"""
function neg_log_likelihood_cdfbeta(beta::Array{Float64,1}, data_A::Array{Float64,1}, data_B::Array{Float64,1})
	-sum(
		 log.(10,
		 	 (Distributions.cdf.(Distributions.Beta(beta[1], beta[2]), data_A)) .-
		 	 (Distributions.cdf.(Distributions.Beta(beta[1], beta[2]), append!(zeros(1), data_A[1:(length(data_A)-1)])))
		     )
		 )     -sum(
		 log.(10,
		 	 (Distributions.cdf.(Distributions.Beta(beta[3], beta[4]), data_B)) .-
		 	 (Distributions.cdf.(Distributions.Beta(beta[3], beta[4]), append!(zeros(1), data_B[1:(length(data_B)-1)])))
		     )
		 )
end

"""
# _____________________________________________________________________________________________________________________________________________
# Iterating function per locus for genome-wide additive allelic effects (alpha) association using Pool-seq data

`GWAlpha_ML_iterator(COUNTS::SharedArray{Int64,3}, snp::Int64, BINS::Array{Float64,1}, MAF::Float64, MIN::Float64, MAX::Float64, SD::Float64)`

# Inputs
1. COUNTS [SharedArray{Int64,3}]: 3-dimensional array of allele counts for each locus for each of the 6 SNP alleles (A:T:C:G:N:DEL) for each pool
2. snp [Int64]: index number of the SNP
3. BINS [Array{Float64,1}]: vector of pool sizes
4. MAF [Float64]: minor allele frequency threshold to filter out alleles (MAF <= f > 1-MAF)
5. MIN [Float64]: minimum phenotype value
6. MAX [Float64]: maximum phenotype value
7. SD [Float64]: standard deviation of phenotype value

# Outputs
1. OUT_ALPHA: vector of additive allelic effects of alleles in the locus
2. OUT_allele: vector of alleles passing minor allele frequency (MAF) threshold filtering
3. OUT_snp: vector of repeated locus id for each of the allelic effects or allele
4, OUT_1: vector of ones and zeros indicating the indices of the alleles passing minor allele frequency (MAF) threshold filtering among the 6 SNP alleles per locus (A:T:C:G:N:DEL)
5, OUT_pA: vector of allele frequencies per locus

# Examples
```
NPOOLS=5; NLOCI=20; BINS=repeat([1.0/NPOOLS],inner=NPOOLS); MAF=0.001; MIN=0.00; MAX=10.00; SD=1.00;
COUNTS = SharedArrays.SharedArray{Int64,3}(NPOOLS, 6, NLOCI)
COUNTS[:,:,:] = Int.(round.(rand(NPOOLS, 6, NLOCI)*100))
for snp in 1:NLOCI
	out = GWAlpha_ML_module.GWAlpha_ML_iterator(COUNTS, snp, BINS, MAF, MIN, MAX, SD)
	println(out)
end
```
"""
function GWAlpha_ML_iterator(COUNTS::SharedArray{Int64,3}, snp::Int64, BINS::Array{Float64,1}, MAF::Float64, MIN::Float64, MAX::Float64, SD::Float64)
	OUT_ALPHA = convert(Array{Float64}, zeros(6))
	OUT_allele = zeros(6)
	OUT_snp = zeros(6)
	OUT_1 = zeros(6)
	OUT_pA = convert(Array{Float64}, zeros(6))
	# #parse allele counts from the sync file
	# COUNTS = zeros(Int64, NPOOLS, 6)
	# for i in 1:NPOOLS
	# 	COUNTS[i,:] = [parse(Int64, x) for x in split.(SYNC[snp, 4:(NPOOLS+3)], [':'])[i]]
	# end
	#convert to frequencies per pool
	counts = COUNTS[:,:,snp]
	FREQS = counts ./ ( sum(counts, dims=2) .+ 1e-10 ) #added 1e-10 to the denominator to avoid NAs in pools with no allele counts (zero depth; which should actually have been filtered out after mpileup using awk)
	allele_freqs = sum(FREQS .* BINS, dims=1)
	#iterate across alleles while filtering by MAF
	if (sum(counts) != 0.0)
		if (minimum(allele_freqs[allele_freqs .!= 0.0]) >= MAF) & (maximum(allele_freqs) < (1.0 - MAF)) #locus filtering by mean MAF
			for allele in 1:6
				if (allele_freqs[allele] > 0.0) & (maximum(FREQS[:,allele]) < 0.999999)  #allele filtering remove alleles with no counts and that the number of pools with allele frequency close to one should not occur even once!
				# if (sum(FREQS[:,allele] .== 0.0) < NPOOLS) & (sum(FREQS[:,allele] .> 0.999999) < 1) #filter-out alleles with at least 1 pool fixed for that allele because it causes a failure in the optimization
					freqA = FREQS[:, allele]
					pA = sum(freqA .* BINS)
					pB = 1 - pA
					BINA = (freqA .* BINS) ./ pA
					BINB = ( (1 .- freqA) .* BINS ) ./ (1-pA)
					percA = cumsum(BINA)
					percB = cumsum(BINB)

					### optimize (minimize) -log-likelihood of these major allele frequencies modelled as a beta distribution
					# using Nelder-Mead optimization or Box minimisation (try-catch if one or the other fails with preference to Nelder-Mead)
					lower_limits = [1e-20, 1e-20, 1e-20, 1e-20]
					upper_limits = [1.0, 1.0, 1.0, 1.0]
					initial_values = [0.1, 0.1, 0.1, 0.1]
					BETA = try
						Optim.optimize(beta->neg_log_likelihood_cdfbeta(beta, percA, percB), initial_values, NelderMead())
					catch
						try
							Optim.optimize(beta->neg_log_likelihood_cdfbeta(beta, percA, percB), lower_limits, upper_limits, initial_values)
						catch ### lower limits of 1e-20 to 1e-6 causes beta dist parameter values to shrink to zero somehow - so we're setting lower limits to 1e-5 instead
							lower_limits = [1e-5, 1e-5, 1e-5, 1e-5]
							Optim.optimize(beta->neg_log_likelihood_cdfbeta(beta, percA, percB), lower_limits, upper_limits, initial_values)
						end
					end
					MU_A = MIN + ((MAX-MIN)*BETA.minimizer[1]/(BETA.minimizer[1]+BETA.minimizer[2]))
					MU_B = MIN + ((MAX-MIN)*BETA.minimizer[3]/(BETA.minimizer[3]+BETA.minimizer[4]))

					### compute alpha
					W_PENAL = 2*sqrt(pA*pB)
					ALPHA = W_PENAL*(MU_A - MU_B) / SD
					OUT_ALPHA[allele] =  ALPHA
					OUT_allele[allele] =  allele
					OUT_snp[allele] =  snp
					OUT_1[allele] =  1
					OUT_pA[allele] =  pA
				else
					## for explicit zero effects for null (low to none) frequency alleles
					OUT_ALPHA[allele] =  0.0
					OUT_allele[allele] =  allele
					OUT_snp[allele] =  snp
					OUT_1[allele] =  0
					OUT_pA[allele] =  allele_freqs[allele]
				end
			end
		end
	end
	return([OUT_ALPHA, OUT_allele, OUT_snp, OUT_1, OUT_pA])
end

"""
# _____________________________________________________________________________________________________________________________________________
# Genome-wide additive allelic effects (alpha) association using Pool-seq data

`GWAlpha_ML(filename_sync::String, filename_phen_py::String, MAF::Float64)`

# Inputs
1. *filename_sync* [String]: filename of the synchronized pileup file, i.e. the genotype file from pool sequencing (Pool-seq)
2. *filename_phen_py* [String]; filename of the phenotype specifications (with ".py" extension name for cross-compatibility with [GWAlpha.py](https://github.com/aflevel/GWAlpha))
3. *MAF* [Float64]: minor allele frequency threshold to filter out alleles (MAF <= f > 1-MAF)

# Output
1. Tuple of vectors of the additive allelic affects across loci with the corresponding locus identification and significance test statistics
	- CHROM
	- POS
	- ALLELE
	- FREQ
	- ALPHA
	- PVALUES
	- LOD

# Example
`OUT = GWAlpha_ML_module.GWAlpha_ML(filename_sync="test/test.sync", filename_phen_py="test/test.py", MAF=0.001)`
"""
function GWAlpha_ML(;filename_sync::String, filename_phen_py::String, MAF::Float64)
	### load the sync and phenotype files
	SYNC = DelimitedFiles.readdlm(filename_sync, '\t')
	phen = DelimitedFiles.readdlm(filename_phen_py, '\t')

	### gather phenotype specifications
	NPOOLS = length(split(phen[5], ['=', ',', '[', ']', ';'])) - 3 #less the first leading and trailing elements
	if length(split(phen[1], ['=', '\"'])) < 3
		global NAME = split(phen[1], ['=', '\''])[3]
	else
		global NAME = split(phen[1], ['=', '\"'])[3]
	end
	SD = parse.(Float64, split(phen[2], ['=',';'])[2])
	MIN = parse.(Float64, split(phen[3], ['=',';'])[2])
	MAX = parse.(Float64, split(phen[4], ['=',';'])[2])
	PERC = parse.(Float64, split(phen[5], ['=', ',', '[', ']', ';'])[3:(NPOOLS+1)])
	QUAN = parse.(Float64, split(phen[6], ['=', ',', '[', ']', ';'])[3:(NPOOLS+1)])
	BINS = append!([x for x in PERC], 1) - append!(zeros(1), PERC)

	### gather genotype (allele frequency) specifications
	NSNP = size(SYNC)[1]
	n_pools_sync = size(SYNC)[2] - 3
	if NPOOLS != n_pools_sync
		println("The number of pools with phenotype data does not match the number of pools with allele frequency data!")
		println("Please check you input files :-)")
		println("Remove leading and intervening whitespaces in the phenotype file.")
		exit()
	else
		n_pools_sync = nothing #clear out contents of this redundant n_pools variable
	end

	### creating an allele counts SharedArray from SYNC to minimize memory usage
	### input shared array
	COUNTS = SharedArrays.SharedArray{Int64,3}(NPOOLS, 6, size(SYNC)[1])
	progress_bar = ProgressMeter.Progress(NPOOLS, dt=1, desc="Converting the sync file into a SharedArray of allele counts: ",  barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)
	for i in 1:NPOOLS
		COUNTS[i,:,:] = parse.(Int64, hcat(split.(SYNC[:, 4:(NPOOLS+3)], [':'])[:, i]...))
		ProgressMeter.update!(progress_bar, i)
	end
	### output shared arrays
	alpha_out = SharedArrays.SharedArray{Float64,1}(NSNP*6)
	allele_id = SharedArrays.SharedArray{Int64,1}(NSNP*6)
	locus_id = SharedArrays.SharedArray{Int64,1}(NSNP*6)
	locus_w_eff = SharedArrays.SharedArray{Int64,1}(NSNP*6)
	allele_freq = SharedArrays.SharedArray{Float64,1}(NSNP*6)

	if length(Distributed.procs()) > 1
		println("Performing GWAlpha_ML in parallel...")
	 	@time x = @sync @distributed for snp in 1:NSNP
			# println(snp)
			OUT = GWAlpha_ML_iterator(COUNTS, snp, BINS, MAF, MIN, MAX, SD)
			idx = collect(((6*snp)-5):(6*snp))
			alpha_out[idx] = OUT[1]
			allele_id[idx] = OUT[2]
			locus_id[idx] = OUT[3]
			locus_w_eff[idx] = OUT[4]
			allele_freq[idx] = OUT[5]
		end
	else
		for snp in 1:NSNP
			# println(snp)
			OUT = GWAlpha_ML_iterator(COUNTS, snp, BINS, MAF, MIN, MAX, SD)
			idx = collect(((6*snp)-5):(6*snp))
			alpha_out[idx] = OUT[1]
			allele_id[idx] = OUT[2]
			locus_id[idx] = OUT[3]
			locus_w_eff[idx] = OUT[4]
			allele_freq[idx] = OUT[5]
		end
	end
	ALPHA_OUT = alpha_out[locus_w_eff .== 1]
	ALLELE_ID_INT = allele_id[locus_w_eff .== 1]
	LOCUS_ID = locus_id[locus_w_eff .== 1]
	LOCUS_W_EFF = locus_w_eff[locus_w_eff .== 1]
	ALLELE_FREQ = allele_freq[locus_w_eff .== 1]

	### estimate heuristic p-values
	P_VALUES, LOD = significance_testing_module.estimate_pval_lod(convert(Array{Float64,1}, ALPHA_OUT))
	### output
	ALLELE_ID = repeat(["N"], inner=length(ALLELE_ID_INT))
	for i in 1:length(ALLELE_ID) #convert int allele ID into corresponding A, T, C, G, N, DEL
		if ALLELE_ID_INT[i] == 1; ALLELE_ID[i] = "A"
		elseif ALLELE_ID_INT[i] == 2; ALLELE_ID[i] = "T"
		elseif ALLELE_ID_INT[i] == 3; ALLELE_ID[i] = "C"
		elseif ALLELE_ID_INT[i] == 4; ALLELE_ID[i] = "G"
		elseif ALLELE_ID_INT[i] == 5; ALLELE_ID[i] = "N"
		elseif ALLELE_ID_INT[i] == 6; ALLELE_ID[i] = "DEL"
		end
	end
	OUT = (CHROM=convert(Array{Any,1},SYNC[LOCUS_ID,1]), POS=convert(Array{Int64,1},SYNC[LOCUS_ID,2]), ALLELE=convert(Array{Any,1},ALLELE_ID), FREQ=convert(Array{Float64,1},ALLELE_FREQ), ALPHA=convert(Array{Float64},ALPHA_OUT), PVALUES=convert(Array{Float64},P_VALUES), LOD=convert(Array{Float64},LOD))
	return(OUT)
end

end # end of GWAlpha_ML_module

max_stars_repo_name: jeffersonfparil/GWAlpha


text:
module day12

using StaticArrays, Combinatorics, ProgressMeter

readInput() = [parse.(Int, split(x, r"[^\d|-]", keepempty=false)) for x in readlines(joinpath(@__DIR__, "input.txt"))]

struct Moon
    pos::MVector{3, Int}
    vel::MVector{3, Int}
end
Moon(pos) = Moon(pos, [0,0,0])

compare(a, b) = a == b ? 0 : a < b ? +1 : -1
function applyGravity!(a::Moon, b::Moon)
    Δ = compare.(a.pos, b.pos)
    a.vel .+= Δ
    b.vel .-= Δ
end
applyVelocity!(m::Moon) = m.pos .+= m.vel
function stepUniverse!(moons::Vector{Moon})
    for (a, b) in combinations(moons, 2)
        applyGravity!(a, b)
    end
    applyVelocity!.(moons)
end

potentialEnergy(m::Moon) = sum(abs, m.pos)
kineticEnergy(m::Moon) = sum(abs, m.vel)
energy(m::Moon) = potentialEnergy(m) * kineticEnergy(m)

function part1(input, steps = 1000)
    moons = Moon.(input)
    for _ = 1:steps
        stepUniverse!(moons)
    end
    sum(energy, moons)
end

function part2(input)
    moons = Moon.(input)
    x0, y0, z0 = (Tuple((m.pos[i], m.vel[i]) for m in moons) for i = 1:3)
    cycles = [0, 0, 0]
    iter = 1
    while any(==(0), cycles)
        stepUniverse!(moons)
        x, y, z = (Tuple((m.pos[i], m.vel[i]) for m in moons) for i = 1:3)
        x == x0 && cycles[1] == 0 && (cycles[1] = iter)
        y == y0 && cycles[2] == 0 && (cycles[2] = iter)
        z == z0 && cycles[3] == 0 && (cycles[3] = iter)
        iter += 1
    end
    lcm(cycles)
end

function part2Plot(input, nbIter=100, smooth=10)
    moons = Moon.(input)
    anim = Animation()
    @showprogress for i = 1:nbIter

        x0, y0, z0 = ([m.pos[i] for m in moons] for i = 1:3)
        stepUniverse!(moons)
        x, y, z = ([m.pos[i] for m in moons] for i = 1:3)

        rx = range.(x0, x, length=smooth)
        ry = range.(y0, y, length=smooth)
        rz = range.(z0, z, length=smooth)
        for j = 1:smooth
            xx, yy, zz = (r -> r[j]).(rx), (r -> r[j]).(ry), (r -> r[j]).(rz)
            p = scatter3d(reshape(xx, 1, :), reshape(yy, 1, :), reshape(zz, 1, :), xlims=(-20,20), ylims=(-20,20), zlims=(-20,20),
                label=["Io" "Europa" "Ganymede" "Callisto"], color= [:green :orange :black :purple])
            frame(anim, p)
        end
    end
    anim
end

end


text:
using OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
#using DiffEqBayes
using Stan, Mamba
using Test

# This example uses the local version of stan_inference for testing purposes

include(joinpath(Pkg.dir("Stan"), "Examples", "OtherExamples", "DiffEqBayes",
  "stan_inference.jl"))

  #=
Experimental link to DiffEqBayes examples

Possible issues to investigate:

1. @ode_def_nohes needs to be kept outside local scope block
2. Minimize string replacements in Stan language model to prevent recompilations

=#

g1 = @ode_def_nohes LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ=10.0 ρ=28.0 β=>2.6666

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  #isdir("tmp") && rm("tmp", recursive=true)
  isfile("LorentzExample-summaryplot-1.pdf") &&
    rm("LorentzExample-summaryplot-1.pdf")
  isfile("LorentzExample-summaryplot-2.pdf") &&
    rm("LorentzExample-summaryplot-2.pdf")

  r0 = [0.1;0.0;0.0]
  tspan = (0.0,4.0)
  prob = ODEProblem(g1,r0,tspan)
  sol = solve(prob,Tsit5())
  t = collect(range(0.1, stop=4.0, length=10))
  randomized = VectorOfArray([(sol(t[i]) + .01randn(3)) for i in 1:length(t)])
  data = convert(Array,randomized)
  priors = [Normal(2.6666,1)]

  bayesian_result = stan_inference(prob,t,data,priors;
    num_samples=1500,num_warmup=1500)
  theta1 = bayesian_result.chain_results[:,["theta.1"],:]
  global sim2 = bayesian_result.chain_results[1:end,
    ["sigma.1", "sigma.2", "sigma.3", "theta.1"], :]
  ## Plotting
  p = plot(sim2, [:trace, :mean, :density, :autocor], legend=true);
  draw(p, ncol=4, filename="Lorentz-summaryplot", fmt=:pdf)

  @test mean(theta1.value[:,:,1]) ≈ 2.6666 atol=1e-1

end


text:
# Test linear regression

using Regression
using Base.Test

### auxiliary

function _affmul(A::Matrix, x::VecOrMat, bias::Float64)
	d = size(x,1) - 1
	if ndims(x) == 1
		A * x[1:d] + 2.0 * x[d+1]
	else
		A * x[1:d,:] .+ 2.0 * x[d+1,:]
	end
end

### linear least square

A = randn(8, 3)
y = randn(8)
Y = randn(8, 5)

function verify_llsq(A::Matrix, y::VecOrMat, m::Symbol)
	x = llsq(A, y; method=m)
	u = A * x
	@test_approx_eq_eps 0.0 dot(vec(y - u), vec(u)) 1.0e-10

	At = A'
	x2 = llsq(At, y; method=m, trans=true)
	@test_approx_eq_eps x x2 1.0e-12

	x = llsq(A, y; method=m, bias=2.0)
	u = _affmul(A, x, 2.0)
	@test_approx_eq_eps 0.0 dot(vec(y - u), vec(u)) 1.0e-10

	x2 = llsq(At, y; method=m, bias=2.0, trans=true)
	@test_approx_eq_eps x x2 1.0e-12
end

println("    on llsq")
for m in [:qrlq, :orth, :svd]
	verify_llsq(A, y, m)
	verify_llsq(A, Y, m)
end

### Ridge regression

function verify_ridgereg(A::Matrix, y::VecOrMat, r::Float64)
	d = size(A, 2)
	At = A'

	x = ridgereg(A, y, r)
	u = A * x
	z = A' * (u - y) + r * x
	@test_approx_eq_eps zeros(size(x)) z 1.0e-10

	x2 = ridgereg(At, y, r; trans=true)
	@test_approx_eq_eps x x2 1.0e-12

	x = ridgereg(A, y, r; bias=2.0)
	u = _affmul(A, x, 2.0)
	_x = ndims(x) == 1 ? x[1:d] : x[1:d,:]
	z = A' * (u - y) + r * _x
	@test_approx_eq_eps zeros(size(z)) z 1.0e-10

	x2 = ridgereg(At, y, r; bias=2.0, trans=true)
	@test_approx_eq_eps x x2 1.0e-10
end

println("    on ridgereg")
verify_ridgereg(A, y, 0.6)
verify_ridgereg(A, Y, 0.6)


text:
"Reward function for basic drilling model"
struct DrillReward{R<:AbstractDrillingRevenue,C<:AbstractDrillingCost,E<:AbstractExtensionCost,S<:AbstractScrapValue} <: AbstractStaticPayoff
    revenue::R
    drill::C
    extend::E
    scrap::S
end

DrillReward(r,d,e) = DrillReward(r,d,e,ScrapValue_Zero())

# access components
revenue(x::DrillReward) = x.revenue
drill(  x::DrillReward) = x.drill
extend( x::DrillReward) = x.extend
scrap(  x::DrillReward) = x.scrap
cost(   x::DrillReward) = drill(x)

# -----------------------------------------
# lengths
# -----------------------------------------

@inline _nparms(    x::DrillReward) = (_nparm(cost(x)), _nparm(extend(x)), _nparm(scrap(x)), _nparm(revenue(x)))
@inline _nparm(     x::DrillReward) = sum(_nparms(x))
@inline _nparm_cost_ext_scrap(x::DrillReward) = _nparm(cost(x)) + _nparm(extend(x)) + _nparm(scrap(x))
@inline _nparm_cost_ext(x::DrillReward) = _nparm(cost(x)) + _nparm(extend(x))

idx_cost(   x::DrillReward) = OneTo(_nparm(cost(x)))
idx_extend( x::DrillReward) = OneTo(_nparm(extend(x)))  .+  _nparm(cost(x))
idx_scrap(  x::DrillReward) = OneTo(_nparm(scrap(x)))   .+ _nparm_cost_ext(x)
idx_revenue(x::DrillReward) = OneTo(_nparm(revenue(x))) .+ _nparm_cost_ext_scrap(x)
idx_ρ(      x::DrillReward) = _nparm(x) # idx_ρ(revenue(x), idx_revenue(x)

idx_drill_ρ(x::DrillReward) = idx_ρ(x)

idx_drill_g(x::DrillReward) = _nparm_cost_ext_scrap(x) + idx_g(revenue(x))
idx_drill_ψ(x::DrillReward) = _nparm_cost_ext_scrap(x) + idx_ψ(revenue(x))
idx_drill_t(x::DrillReward) = _nparm_cost_ext_scrap(x) .+ idx_t(revenue(x))
idx_drill_D(x::DrillReward) = _nparm_cost_ext_scrap(x) .+ idx_D(revenue(x))

idx_drill_g(d::DataDrill{<:AbstractDynamicDrillModel}) = idx_drill_g(reward(_model(d)))
idx_drill_ψ(d::DataDrill{<:AbstractDynamicDrillModel}) = idx_drill_ψ(reward(_model(d)))
idx_drill_t(d::DataDrill{<:AbstractDynamicDrillModel}) = idx_drill_t(reward(_model(d)))
idx_drill_D(d::DataDrill{<:AbstractDynamicDrillModel}) = idx_drill_D(reward(_model(d)))

vw_cost(   x::DrillReward, theta) = uview(theta, idx_cost(x))
vw_extend( x::DrillReward, theta) = uview(theta, idx_extend(x))
vw_scrap(  x::DrillReward, theta) = uview(theta, idx_scrap(x))
vw_revenue(x::DrillReward, theta) = uview(theta, idx_revenue(x))

# -----------------------------------------
# flow payoffs
# -----------------------------------------

function flow!(grad, x::DrillReward, d, obs, θ, sim, dograd)
    c = flow!(vw_cost(   x, grad), cost(   x), d, obs, vw_cost(   x, θ), sim, dograd)
    e = flow!(vw_extend( x, grad), extend( x), d, obs, vw_extend( x, θ), sim, dograd)
    s = flow!(vw_scrap(  x, grad), scrap(  x), d, obs, vw_scrap(  x, θ), sim, dograd)
    r = flow!(vw_revenue(x, grad), revenue(x), d, obs, vw_revenue(x, θ), sim, dograd)
    return c+e+s+r
end

function flowdψ(x::DrillReward, d, obs, θ, sim)
    T = eltype(θ)
    if d == 0
        u = flowdψ(extend(x), d, obs, vw_extend(x,θ), sim)
    else
        c = flowdψ(cost(x),    d, obs, vw_cost(x,θ), sim)
        r = flowdψ(revenue(x), d, obs, vw_revenue(x,θ), sim)
        u = r+c
    end
    return u::T
end

function coefnames(x::DrillReward)
    c = coefnames(cost(   x))
    e = coefnames(extend( x))
    s = coefnames(scrap(  x))
    r = coefnames(revenue(x))
    return vcat(c, e, s, r)
end

# -----------------------------------------
# deprecate
# -----------------------------------------

# @deprecate number_of_model_parms(x::DrillReward) _nparm(x)
# @deprecate coef_range_revenue(x)       idx_revenue(x)
# @deprecate coef_range_drillingcost(x)  idx_cost(x)
# @deprecate coef_range_extensioncost(x) idx_extend(x)
# @deprecate flowdθ!(args...) dflow!(args...)
# @deprecate theta_revenue(x::DrillReward, theta) vw_revenue(x, theta)
# @deprecate theta_cost(   x::DrillReward, theta) vw_cost(   x, theta)
# @deprecate theta_extend( x::DrillReward, theta) vw_extend( x, theta)


text:
"""
    VacuumChargeDriftModel <: AbstractChargeDriftModel
"""
struct VacuumChargeDriftModel{T <: SSDFloat} <: AbstractChargeDriftModel{T} end

function get_electron_drift_field(ef::Array{SVector{3,T},3}, ::VacuumChargeDriftModel)::Array{SVector{3,T},3} where {T <: SSDFloat}
    return -ef
end
function get_hole_drift_field(ef::Array{SVector{3,T},3}, ::VacuumChargeDriftModel)::Array{SVector{3,T},3} where {T <: SSDFloat}
    return ef
end


text: include("tsp.jl")


text:
using Documenter, Remarkable

makedocs(;
    modules  = [Remarkable],
    format   = Documenter.HTML(analytics = "UA-68884109-1"),
    pages    = [
        "Home" => "index.md",
    ],
    repo     = "https://github.com/TotalVerb/Remarkable.jl/blob/{commit}{path}#L{line}",
    sitename = "Remarkable.jl",
    authors  = "Fengyang Wang",
)

deploydocs(
    repo = "github.com/TotalVerb/Remarkable.jl",
)


text:
module ChainRulesTestUtils

using ChainRulesCore
using ChainRulesCore: frule, rrule
using Compat: only
using FiniteDifferences
using FiniteDifferences: to_vec
using LinearAlgebra
using Random
using Test

const _fdm = central_fdm(5, 1)

export TestIterator
export check_equal, test_scalar, frule_test, rrule_test, generate_well_conditioned_matrix

include("generate_tangent.jl")
include("data_generation.jl")
include("iterator.jl")
include("check_result.jl")
include("testers.jl")

include("deprecated.jl")
end # module


text:
function make_mask(
    Nside, θpix;
    edge_padding_deg = 2,
    edge_rounding_deg = 1,
    apodization_deg = 1,
    ptsrc_radius_arcmin = 7, # roughly similar to Story et al. 2015
    num_ptsrcs = round(Int, prod(Nside .* (1,1)) * (θpix/60)^2 * 120/100)  # SPT-like density of sources
)

    deg2npix(x) = round(Int, x/θpix*60)
    arcmin2npix(x) = round(Int, x/θpix)

    ptsrc = .!bleed(sim_ptsrcs(Nside, num_ptsrcs), arcmin2npix(ptsrc_radius_arcmin))
    boundary = boundarymask(Nside, deg2npix(edge_padding_deg))
    mask_array = if apodization_deg in (false, 0)
        boundary .& ptsrc
    else
        cos_apod(boundary, deg2npix(apodization_deg), deg2npix(edge_rounding_deg)) .* cos_apod(ptsrc, arcmin2npix(ptsrc_radius_arcmin));
    end

    FlatMap(Float32.(mask_array), θpix=θpix)
end

make_mask(::FlatField{<:Flat{Nside,θpix}}; kwargs...) where {Nside,θpix} = make_mask(Nside, θpix; kwargs...)

# all padding/smoothing/etc... quantities below are in units of numbers of pixels

function boundarymask(Nside, pad)
    m = fill(true, Nside .* (1,1,))
    m[1:pad,:]          .= false
    m[:,1:pad]          .= false
    m[end-pad+1:end,:]  .= false
    m[:,end-pad+1:end]  .= false
    m
end

function bleed(img, w)
    Ny,Nx = size(img)
    nearest = getfield.(@ondemand(ImageMorphology.feature_transform)(img),:I)
    [norm(nearest[i,j] .- [i,j]) < w for i=1:Ny,j=1:Nx]
end

function cos_apod(img, w, smooth_distance=false)
    Ny,Nx = size(img)
    nearest = getfield.(@ondemand(ImageMorphology.feature_transform)(.!img),:I)
    distance = [norm(nearest[i,j] .- [i,j]) for i=1:Ny,j=1:Nx]
    if smooth_distance!=false
        distance = @ondemand(ImageFiltering.imfilter)(distance, @ondemand(ImageFiltering.Kernel.gaussian)(smooth_distance))
    end
    @. (1-cos(min(distance,w)/w*pi))/2
end

function round_edges(img, w)
    .!(@ondemand(ImageFiltering.imfilter)(img, @ondemand(ImageFiltering.Kernel.gaussian)(w)) .< 0.5)
end

function sim_ptsrcs(Nside,nsources)
    Ny, Nx = Nside .* (1,1)
    m = fill(false,Ny,Nx);
    for i=1:nsources
        m[rand(1:Ny),rand(1:Nx)] = true
    end
    m
end


text:
# TODO: deprecate KUsparse -> CudaSparseMatrixCSC
# TODO: reimplement kperceptron in new net language

# (c) Deniz Yuret, June 19, 2015
# This is a standalone (single layer) implementation of the averaged
# kernel perceptron algorithm.
# Based on the k_perceptron_multi_train.m in the DOGMA library
# (http://dogma.sourceforge.net) by Francesco Orabona which references:
#     - Crammer, K., & Singer Y. (2003).
#       Ultraconservative Online Algorithms for Multiclass Problems.
#       Journal of Machine Learning Research 3, (pp. 951-991).
# Averaging optimization trick with w0/w1/w2 based on:
# http://ciml.info/dl/v0_9/ciml-v0_9-ch03.pdf.

type KPerceptron <: Loss
    nclass	# number of output classes
    kernel      # kernel function
    p           # kernel parameter tuple
    s           # support vectors
    x           # input minibatch
    k           # kernel matrix
    y           # model output
    u           # number of training updates
    u2          # the last time w2 was calculated
    w0          # regular weights
    w1          # w2-u*w0 (kept up to date during training)
    w2          # summed weights (computed before prediction)
    dw0         # new weights
    dw1         # new weights
    dj          # indices for new support vectors
    dn          # number of new support vectors
    KPerceptron(nclass,kernel,kparams=nothing)=new(nclass,kernel,kparams)
end

function forw(l::KPerceptron, x; train=true, o...)
    initforw(l, x, train)
    l.k = l.kernel(l.k, l.x, l.s, l.p)          # l.s generally larger, so we will transpose l.x, e.g. k=x'*s
    w = (train ? l.w0 : l.w2)                   # w2 averaged, w0 regular weights
    A_mul_Bt!(l.y, w, l.k)                      # l.y = w * l.k'
    return l.y
end

function back(l::KPerceptron, z; returndx=false, o...)
    returndx && error("KPerceptron does not know how to return dx")
    (y,z) = initback(l,z)
    @inbounds for j=1:size(z,2)
        (cz,cy,ymax,zmax) = (0,0,typemin(eltype(y)),typemin(eltype(z)))
        @inbounds for i=1:l.nclass
            z[i,j] > zmax && ((cz,zmax) = (i,z[i,j])) # find the correct answer
            y[i,j] > ymax && ((cy,ymax) = (i,y[i,j])) # find the model answer
        end
        if cz != cy # if model answer is not correct l.x[:,j] becomes a new support vector
            l.dn += one(l.dn)
            l.dj[l.dn] = j
            u = l.u + j - 1
            l.dw1[cz,j] = -u
            l.dw1[cy,j] = u
            l.dw0[cz,j] = 1
            l.dw0[cy,j] = -1
        end
    end
end

function update!(l::KPerceptron; o...) # 198
    l.u += size(l.x,2)
    dj = sub(l.dj, 1:l.dn)
    l.s  = ccat!(l.s,  l.x,   dj)
    l.w0 = ccat!(l.w0, l.dw0, dj)
    l.w1 = ccat!(l.w1, l.dw1, dj)
    l.k  = resize!(l.k, (size(l.k,1),size(l.s,2)))
end

# hcat!{T}(a::Matrix{T}, b::Matrix{T}, vj=(1:size(b,2)), nj=length(vj))=[a b[:,vj[1:nj]]]
# size!(a::Array, d::Dims)=(size(a)==d ? a : Array(eltype(a),d))

# To preserve the behavior and minimize the space, get rid of everything except:
# nclass, kernel, p, s, u, w0, w1
# Also filter the support vectors to keep only the unique ones.
function strip!(l::KPerceptron)
    l.x=l.k=l.y=l.dw0=l.dw1=l.dj=l.dn=nothing
    (l.s,l.w0,l.w1)=uniq!(l.s,l.w0,l.w1)
    return l
end

# We only support KUsparse/KUdense for x, s; KUdense for all else.

function initforw(l::KPerceptron, x, train)
    @assert isa(x, KUsparse) || isa(x, KUdense)
    if !isdefined(l,:s)                         # first initialization
        similar!(l,:s,x,size(x,1),0)      	# s matches x in location, sparseness, eltype, orientation
        l.w0 = KUdense(atype(x),eltype(x),(l.nclass, 0)) # w matches x in location and eltype but is dense
        l.w1 = copy(l.w0)
        l.w2 = copy(l.w0)
        l.u = l.u2 = 0
    end
    l.x = x                                     # x can be cpu/gpu dense/sparse
    @assert typeof(l.x) == typeof(l.s)
    @assert atype(l.x) == atype(l.w0)
    @assert eltype(l.x) == eltype(l.s) == eltype(l.w0)
    @assert size(l.s) == (size(l.x,1), size(l.w0,2))
    similar!(l,:y,l.w0,(size(l.w0,1),size(l.x,2)))
    similar!(l,:k,l.w0,(size(l.x,2),size(l.s,2)))
    if !train && (l.u2 != l.u)
        copy!(l.w2, l.w1)
        axpy!(l.u, l.w0, l.w2)
        l.u2 = l.u
    end
end

function initback(l::KPerceptron, z::KUdense, y::KUdense=l.y)
    @assert size(z) == size(y)
    z = to_host(z.arr)
    y = to_host(y.arr)
    similar!(l,:dw0,y)
    similar!(l,:dw1,y)
    similar!(l,:dj,y,Int32,size(y,2))
    l.dn = 0
    fill!(l.dw0,0)
    fill!(l.dw1,0)
    return (y,z)
end

# Some common kernels
# http://crsouza.com/2010/03/kernel-functions-for-machine-learning-applications

klinear0(k, x, s, p)=(x.' * s)
kpoly0(k, x, s, p)=((x.' * s + p[1]) .^ p[2])
kgauss0(k, x, s, p)=exp(-p[1] * broadcast(+, sum(x.^2,1).', broadcast(+, sum(s.^2,1), -2*(x.' * s))))

# More efficient implementations:

klinear(k, x, s, p)=At_mul_B!(k, x, s)          # k=x'*s

function kpoly(k, x, s, p)
    k = klinear(k, x, s, p)                                               # 1670
    kpolymap(k, p[1], p[2])
    return k
end

function kpolymap(k::Array, c, d)
    @inbounds @simd for i=1:length(k)
        k[i] = (k[i] + c).^d
    end
    return k
end

kpolymap(k::KUdense, c, d)=(kpolymap(k.arr, c, d); k)

# This is a more efficient implementation that covers kpoly + kpolymap
# Do we need gpu kpolymap if we use this?
kpoly{A<:CudaArray}(k::KUdense{A}, x::KUsparse{A}, s::KUsparse{A}, p)=(kpoly(k.arr, convert(KUsparse,x), convert(KUsparse,s), p); k)

kgauss(k::KUdense, x::KUdense, s::KUdense, p)=(kgauss(k.arr, x.arr, s.arr, p); k)
kgauss{A<:BaseArray}(k::KUdense{A}, x::KUsparse{A}, s::KUsparse{A}, p)=(kgauss(k.arr, convert(KUsparse, x), convert(KUsparse, s), p); k)
kgauss(k::Array, x::KUsparse{Array}, s::KUsparse{Array}, p)=(kgauss(k, convert(SparseMatrixCSC,x), convert(SparseMatrixCSC,s), p); k)

function kgauss(k::AbstractArray, x::AbstractArray, s::AbstractArray, p)         # 2582
    k = klinear(k, x, s, p) # 1741
    xx = sum(x.^2,1) # 10
    ss = sum(s.^2,1) # 419 Can be cached
    # return exp(-p[1] * broadcast!(+, k, xx', broadcast!(+, k, ss, -2*k)))
    k1,k2 = size(k); p1 = p[1]
    @inbounds @simd for i=1:k1
        @inbounds @simd for j=1:k2
            k[i,j] = exp(-p1 * (xx[i] + ss[j] - 2*k[i,j]))
        end
    end
    return k
end

if GPU

function kgauss(k::CudaArray{Float32}, x::CudaArray{Float32}, s::CudaArray{Float32}, p)
    @assert size(x,1)==size(s,1)
    @assert size(k)==(size(x,2),size(s,2))
    k = klinear(k, x, s, p)
    x2 = CudaArray(Float32, size(x,2))
    s2 = CudaArray(Float32, size(s,2))
    ccall((:kgauss32sum,libkunet),Void,(Cint,Cint,Ptr{Cfloat},Ptr{Cfloat}),size(x,1),size(x,2),x,x2)
    ccall((:kgauss32sum,libkunet),Void,(Cint,Cint,Ptr{Cfloat},Ptr{Cfloat}),size(s,1),size(s,2),s,s2)
    ccall((:kgauss32map,libkunet),Void,(Cint,Cint,Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat},Cfloat),
          size(x,2),size(s,2),x2,s2,k,p[1])
    gpusync()
    free(x2); free(s2)
    return k
end

function kgauss(k::CudaArray{Float64}, x::CudaArray{Float64}, s::CudaArray{Float64}, p)
    @assert size(x,1)==size(s,1)
    @assert size(k)==(size(x,2),size(s,2))
    k = klinear(k, x, s, p)
    x2 = CudaArray(Float64, size(x,2))
    s2 = CudaArray(Float64, size(s,2))
    ccall((:kgauss64sum,libkunet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),size(x,1),size(x,2),x,x2)
    ccall((:kgauss64sum,libkunet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),size(s,1),size(s,2),s,s2)
    ccall((:kgauss64map,libkunet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Cdouble),
          size(x,2),size(s,2),x2,s2,k,p[1])
    free(x2); free(s2)
    gpusync()
    return k
end

function kgauss(k::CudaArray{Float32}, x::KUsparse{CudaArray,Float32}, s::KUsparse{CudaArray,Float32}, p)
    @assert size(k)==(size(x,2),size(s,2))
    ccall((:kgauss32,libkunet),Void,
          (Cint,Cint,Ptr{Cfloat},Ptr{Cint},Ptr{Cint},Ptr{Cfloat},Ptr{Cint},Ptr{Cint},Ptr{Cfloat},Cfloat),
          size(x,2),size(s,2),x.nzval,x.rowval,x.colptr,s.nzval,s.rowval,s.colptr,k,p[1])
    gpusync()
    return k
end

function kgauss(k::CudaArray{Float64}, x::KUsparse{CudaArray,Float64}, s::KUsparse{CudaArray,Float64}, p)
    @assert size(k)==(size(x,2),size(s,2))
    ccall((:kgauss64,libkunet),Void,
          (Cint,Cint,Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Cdouble),
          size(x,2),size(s,2),x.nzval,x.rowval,x.colptr,s.nzval,s.rowval,s.colptr,k,p[1])
    gpusync()
    return k
end

function kpolymap(k::CudaArray{Float32}, c, d)
    ccall((:kpolymap32,libkunet),Void,
          (Cint,Ptr{Cfloat},Cfloat,Cfloat),
          length(k),k,c,d)
    gpusync()
    return k
end

function kpolymap(k::CudaArray{Float64}, c, d)
    ccall((:kpolymap64,libkunet),Void,
          (Cint,Ptr{Cdouble},Cdouble,Cdouble),
          length(k),k,c,d)
    gpusync()
    return k
end

function kpoly(k::CudaArray{Float32}, x::KUsparse{CudaArray,Float32}, s::KUsparse{CudaArray,Float32}, p)
    @assert size(k)==(size(x,2),size(s,2))
    ccall((:kpoly32,libkunet),Void,
          (Cint,Cint,Ptr{Cfloat},Ptr{Cint},Ptr{Cint},Ptr{Cfloat},Ptr{Cint},Ptr{Cint},Ptr{Cfloat},Cfloat,Cfloat),
          size(x,2),size(s,2),x.nzval,x.rowval,x.colptr,s.nzval,s.rowval,s.colptr,k,p[1],p[2])
    gpusync()
    return k
end

function kpoly(k::CudaArray{Float64}, x::KUsparse{CudaArray,Float64}, s::KUsparse{CudaArray,Float64}, p)
    @assert size(k)==(size(x,2),size(s,2))
    ccall((:kpoly64,libkunet),Void,
          (Cint,Cint,Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Cdouble,Cdouble),
          size(x,2),size(s,2),x.nzval,x.rowval,x.colptr,s.nzval,s.rowval,s.colptr,k,p[1],p[2])
    gpusync()
    return k
end

end # if GPU

# Failed optimization experiments:

# using NumericExtensions # but this slows kgauss down!
# using InplaceOps # maybe later for more readable code...

# # Why is kpoly slower than kgauss?  This is not faster either:
# function kpoly1(k, x, s, p)
#     k = klinear(k, x, s, p)
#     return (k + p[1]).^p[2]
# end

# function kgauss2(x::SparseMatrixCSC, s::SparseMatrixCSC, p, k)         # 2582
#     k = klinear(k, x, s, p) # 1741
#     xx = sum(x.^2,1) # 10
#     ss = sum(s.^2,1) # 419 Can be cached
#     k = broadcast!(+, k, xx', broadcast!(+, k, ss, -2*k))
#     g = -p[1]
#     @in1! k .* g
#     return exp!(k)
# end

# # This is much slower than kgauss and kgauss0

# function kgauss1(x::SparseMatrixCSC, s::SparseMatrixCSC, p, k)
#     k = klinear(k, x, s, p) # 1741
#     xx = sum(x.^2,1) # 10
#     ss = sum(s.^2,1) # 419 Can be cached
#     # return exp(-p[1] * broadcast(+, xx', broadcast(+, ss, -2*k))) # 412
#     g = p[1]
#     @inbounds @simd for i=1:size(k,1)
#         @inbounds @simd for j=1:size(k,2)
#             k[i,j] = exp(-g * (xx[i] + ss[j] - 2*k[i,j]))
#         end
#     end
#     return k
# end

# function kgauss2(k, x, s, p)                    # buggy: does not take into account cells where one matrix is 0
#     k2 = kgauss(x, s, p, copy(k))
#     x = x'
#     @assert size(k)==(size(x,1), size(s,2))
#     fill!(k, zero(eltype(k)))
#     @inbounds @simd for scol=1:size(s,2)
#         @inbounds @simd for sp=s.colptr[scol]:(s.colptr[scol+1]-1)
#             srow = s.rowval[sp]
#             sval = s.nzval[sp]  # 133
#             @inbounds @simd for xp=x.colptr[srow]:(x.colptr[srow+1]-1)
#                 xrow = x.rowval[xp] # 63
#                 xval = x.nzval[xp]  # 217
#                 kinc = (xval - sval)^2
#                 k[xrow,scol] += kinc
#             end
#         end
#     end
#     g = p[1]
#     @inbounds @simd for i=1:length(k); k[i] = exp(-g*k[i]); end
#     isempty(k) || (@show maximum(abs(k-k2)))
#     return k
# end

# buggy for the same reason
# function kgauss1(x::Sparse{CudaArray,Float32}, s::Sparse{CudaArray,Float32}, p::Vector{Float32}, k::CudaArray{Float32})
#     t = x' # do this somewhere else?
#     @assert size(k)==(size(t,1), size(s,2))
#     isempty(k) && return k
#     fill!(k, zero(eltype(k)))
#     ccall((:kgauss32,libkunet),Void,
#           (Cint,Cint,Ptr{Float32},Ptr{Cint},Ptr{Cint},Ptr{Float32},Ptr{Cint},Ptr{Cint},Float32,Ptr{Float32}),
#           size(t,1),size(s,2),t.nzval,t.rowval,t.colptr,s.nzval,s.rowval,s.colptr,p[1],k)
#     k1 = cpucopy(k)
#     k2 = kgauss1(cpucopy(x),cpucopy(s),p,cpucopy(k))
#     isempty(k1) || (@show maximum(abs(k1-k2)))
#     return k
#     #copy!(k, kgauss(cpucopy(x),cpucopy(s),p,cpucopy(k)))
# end

# function kgauss3(x,s,p,k)  # Too slow
#     @assert size(k)==(size(x,2),size(s,2))
#     for i=1:size(x,2)
#         for j=1:size(s,2)
#             k[i,j] = zero(eltype(k)) # calculate k[i,j]
#             x1 = x.colptr[i]; x2 = x.colptr[i+1]
#             s1 = s.colptr[j]; s2 = s.colptr[j+1]
#             xr = x.rowval[x1]; sr = s.rowval[s1]
#             xv = x.nzval[x1]; sv = s.nzval[s1]
#             while (x1 < x2 || s1 < s2)
#                 xr = (x1 < x2 ? x.rowval[x1] : size(x,1)+1)
#                 sr = (s1 < s2 ? s.rowval[s1] : size(s,1)+1)
#                 if xr == sr
#                     xv = x.nzval[x1]; sv = s.nzval[s1]
#                     k[i,j] += (xv-sv)^2
#                     x1 += 1; s1 += 1
#                 elseif xr < sr
#                     xv = x.nzval[x1]
#                     k[i,j] += xv^2
#                     x1 += 1
#                 elseif sr < xr
#                     sv = s.nzval[s1]
#                     k[i,j] += sv^2
#                     s1 += 1
#                 end
#             end
#         end
#     end
#     g = p[1]
#     for i=1:length(k); k[i] = exp(-g*k[i]); end
#     return k
# end

# function klinear3(x,s,p,k)  # Too slow on cpu, ok on gpu
#     @assert size(k)==(size(x,2),size(s,2))
#     @inbounds @simd for i=1:size(x,2)
#         @inbounds @simd for j=1:size(s,2)
#             kij = zero(Float64)
#             x1 = x.colptr[i]; x2 = x.colptr[i+1]
#             s1 = s.colptr[j]; s2 = s.colptr[j+1]
#             while ((x1 < x2) && (s1 < s2))
#                 xr = x.rowval[x1]; sr = s.rowval[s1]
#                 (xr < sr ? x1+=1 :
#                  sr < xr ? s1+=1 :
#                  (kij += x.nzval[x1] * s.nzval[s1]; x1+=1; s1+=1))
#             end
#             k[i,j] = kij
#         end
#     end
#     return k
# end

# klinear4(x,s,p,k)=A_mul_B!(k,x.',s)


text:
include("example_one.jl")
include("example_two.jl")
include("example_three.jl")
include("example_four.jl")
include("example_five.jl")


text:
@testset "GGGenome" begin
    @testset "gggsearch" begin
        query = "TTCATTGACAACATT"
        res = gggsearch(query, db="felCat5", format = "html")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/html")
        @test occursin(Regex("query="*query), String(res.body))

        res = gggsearch(query, db="felCat5", format = "txt", output = "toString")
        @test typeof(res) == String

        res = gggsearch(query, db="felCat5", format = "txt", output = "extractTopHit")
        @test typeof(res) == String
    end

    @testset "gggdbs" begin
       @test  length(gggdbs()) > 0
    end

end


text:
for (func, typ) in [(:clblasSsyr2k, cl.CL_float),
                    (:clblasDsyr2k, cl.CL_double),
                    (:clblasCsyr2k, CL_float2),
                    (:clblasZsyr2k, CL_double2)]

    @eval @blasfun $func(order::clblasOrder, uplo::clblasUplo,
                              transAB::clblasTranspose,
                              N::Csize_t, K::Csize_t,
                              alpha::$typ,
                              A::cl.CL_mem, offA::Csize_t, lda::Csize_t,
                              B::cl.CL_mem, offB::Csize_t, ldb::Csize_t,
                              beta::$typ,
                              C::cl.CL_mem, offC::Csize_t, ldc::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder, uplo::clblasUplo,
                               transAB::clblasTranspose,
                               N::Csize_t, K::Csize_t,
                               alpha::$typ,
                               A::cl.CL_mem, offA::Csize_t, lda::Csize_t,
                               B::cl.CL_mem, offB::Csize_t, ldb::Csize_t,
                               beta::$typ,
                               C::cl.CL_mem, offC::Csize_t, ldc::Csize_t)

end