text:
@testset "5.5.2 Inverse secant functions" begin
    (a, b, c, d, n, x, ) = @variables a b c d n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Inverse*Secants=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*asec(a*x^n)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*asec(a*x^n)=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [asec(a*x^5)/x, x, 7, (1/10)*im*asec(a*x^5)^2 - (1/5)*asec(a*x^5)*log(1 + ℯ^(2*im*asec(a*x^5))) + (1/10)*im*PolyLog(2, -ℯ^(2*im*asec(a*x^5)))]

    @test_int [x^3*asec(sqrt(x)), x, 4, (-(1/4))*sqrt(-1 + x) - (1/4)*(-1 + x)^(3/2) - (3/20)*(-1 + x)^(5/2) - (1/28)*(-1 + x)^(7/2) + (1/4)*x^4*asec(sqrt(x))]
    @test_int [x^2*asec(sqrt(x)), x, 4, (-(1/3))*sqrt(-1 + x) - (2/9)*(-1 + x)^(3/2) - (1/15)*(-1 + x)^(5/2) + (1/3)*x^3*asec(sqrt(x))]
    @test_int [x^1*asec(sqrt(x)), x, 4, (-(1/2))*sqrt(-1 + x) - (1/6)*(-1 + x)^(3/2) + (1/2)*x^2*asec(sqrt(x))]
    @test_int [x^0*asec(sqrt(x)), x, 3, -sqrt(-1 + x) + x*asec(sqrt(x))]
    @test_int [asec(sqrt(x))/x^1, x, 7, im*asec(sqrt(x))^2 - 2*asec(sqrt(x))*log(1 + ℯ^(2*im*asec(sqrt(x)))) + im*PolyLog(2, -ℯ^(2*im*asec(sqrt(x))))]
    @test_int [asec(sqrt(x))/x^2, x, 5, sqrt(-1 + x)/(2*x) - asec(sqrt(x))/x + (1/2)*atan(sqrt(-1 + x))]
    @test_int [asec(sqrt(x))/x^3, x, 6, sqrt(-1 + x)/(8*x^2) + (3*sqrt(-1 + x))/(16*x) - asec(sqrt(x))/(2*x^2) + (3/16)*atan(sqrt(-1 + x))]
    @test_int [asec(sqrt(x))/x^4, x, 7, sqrt(-1 + x)/(18*x^3) + (5*sqrt(-1 + x))/(72*x^2) + (5*sqrt(-1 + x))/(48*x) - asec(sqrt(x))/(3*x^3) + (5/48)*atan(sqrt(-1 + x))]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [x^2*asec(a/x), x, 5, (-(1/3))*a^3*sqrt(1 - x^2/a^2) + (1/9)*a^3*(1 - x^2/a^2)^(3/2) + (1/3)*x^3*acos(x/a)]
    @test_int [x^1*asec(a/x), x, 4, (-(1/4))*a*x*sqrt(1 - x^2/a^2) + (1/2)*x^2*acos(x/a) + (1/4)*a^2*asin(x/a)]
    @test_int [x^0*asec(a/x), x, 3, (-a)*sqrt(1 - x^2/a^2) + x*acos(x/a)]
    @test_int [asec(a/x)/x^1, x, 6, (-(1/2))*im*acos(x/a)^2 + acos(x/a)*log(1 + ℯ^(2*im*acos(x/a))) - (1/2)*im*PolyLog(2, -ℯ^(2*im*acos(x/a)))]
    @test_int [asec(a/x)/x^2, x, 5, -(acos(x/a)/x) + atanh(sqrt(1 - x^2/a^2))/a]
    @test_int [asec(a/x)/x^3, x, 3, sqrt(1 - x^2/a^2)/(2*a*x) - acos(x/a)/(2*x^2)]
    @test_int [asec(a/x)/x^4, x, 6, sqrt(1 - x^2/a^2)/(6*a*x^2) - acos(x/a)/(3*x^3) + atanh(sqrt(1 - x^2/a^2))/(6*a^3)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*asec(a*x^n)*when*n*symbolic=#

    @test_int [asec(a*x^n)/x, x, 7, (im*asec(a*x^n)^2)/(2*n) - (asec(a*x^n)*log(1 + ℯ^(2*im*asec(a*x^n))))/n + (im*PolyLog(2, -ℯ^(2*im*asec(a*x^n))))/(2*n)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*asec(a+b*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^4*asec(a + b*x), x, 9, (a*(20 + 53*a^2)*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(30*b^5) + (11*a*x^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(60*b^3) - (x^3*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(20*b^2) - ((9 + 58*a^2)*(a + b*x)^2*sqrt(1 - 1/(a + b*x)^2))/(120*b^5) + (a^5*asec(a + b*x))/(5*b^5) + (1/5)*x^5*asec(a + b*x) - ((3 + 40*a^2 + 40*a^4)*atanh(sqrt(1 - 1/(a + b*x)^2)))/(40*b^5)]
    @test_int [x^3*asec(a + b*x), x, 8, -(((2 + 17*a^2)*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(12*b^4)) - (x^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(12*b^2) + (a*(a + b*x)^2*sqrt(1 - 1/(a + b*x)^2))/(3*b^4) - (a^4*asec(a + b*x))/(4*b^4) + (1/4)*x^4*asec(a + b*x) + (a*(1 + 2*a^2)*atanh(sqrt(1 - 1/(a + b*x)^2)))/(2*b^4)]
    @test_int [x^2*asec(a + b*x), x, 7, (5*a*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*b^3) - (x*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*b^2) + (a^3*asec(a + b*x))/(3*b^3) + (1/3)*x^3*asec(a + b*x) - ((1 + 6*a^2)*atanh(sqrt(1 - 1/(a + b*x)^2)))/(6*b^3)]
    @test_int [x^1*asec(a + b*x), x, 6, -(((a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(2*b^2)) - (a^2*asec(a + b*x))/(2*b^2) + (1/2)*x^2*asec(a + b*x) + (a*atanh(sqrt(1 - 1/(a + b*x)^2)))/b^2]
    @test_int [x^0*asec(a + b*x), x, 5, ((a + b*x)*asec(a + b*x))/b - atanh(sqrt(1 - 1/(a + b*x)^2))/b]
    @test_int [asec(a + b*x)/x^1, x, 14, asec(a + b*x)*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) + asec(a + b*x)*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) - asec(a + b*x)*log(1 + ℯ^(2*im*asec(a + b*x))) - im*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) - im*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) + (1/2)*im*PolyLog(2, -ℯ^(2*im*asec(a + b*x)))]
    @test_int [asec(a + b*x)/x^2, x, 5, -((b*asec(a + b*x))/a) - asec(a + b*x)/x + (2*b*atan((sqrt(1 + a)*tan((1/2)*asec(a + b*x)))/sqrt(1 - a)))/(a*sqrt(1 - a^2))]
    @test_int [asec(a + b*x)/x^3, x, 7, (b*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(2*a*(1 - a^2)*x) + (b^2*asec(a + b*x))/(2*a^2) - asec(a + b*x)/(2*x^2) - ((1 - 2*a^2)*b^2*atan((sqrt(1 + a)*tan((1/2)*asec(a + b*x)))/sqrt(1 - a)))/(a^2*(1 - a^2)^(3/2))]
    @test_int [asec(a + b*x)/x^4, x, 8, (b*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*a*(1 - a^2)*x^2) - ((2 - 5*a^2)*b^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*a^2*(1 - a^2)^2*x) - (b^3*asec(a + b*x))/(3*a^3) - asec(a + b*x)/(3*x^3) + ((2 - 5*a^2 + 6*a^4)*b^3*atan((sqrt(1 + a)*tan((1/2)*asec(a + b*x)))/sqrt(1 - a)))/(3*a^3*(1 - a^2)^(5/2))]

    @test_int [x^3*asec(a + b*x)^2, x, 20, -((a*x)/b^3) + (a + b*x)^2/(12*b^4) - ((a + b*x)*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/(3*b^4) - (3*a^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/b^4 + (a*(a + b*x)^2*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/b^4 - ((a + b*x)^3*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/(6*b^4) - (a^4*asec(a + b*x)^2)/(4*b^4) + (1/4)*x^4*asec(a + b*x)^2 - (2*im*a*asec(a + b*x)*atan(ℯ^(im*asec(a + b*x))))/b^4 - (4*im*a^3*asec(a + b*x)*atan(ℯ^(im*asec(a + b*x))))/b^4 + log(a + b*x)/(3*b^4) + (3*a^2*log(a + b*x))/b^4 + (im*a*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^4 + (2*im*a^3*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^4 - (im*a*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^4 - (2*im*a^3*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^4]
    @test_int [x^2*asec(a + b*x)^2, x, 17, x/(3*b^2) + (2*a*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/b^3 - ((a + b*x)^2*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/(3*b^3) + (a^3*asec(a + b*x)^2)/(3*b^3) + (1/3)*x^3*asec(a + b*x)^2 + (2*im*asec(a + b*x)*atan(ℯ^(im*asec(a + b*x))))/(3*b^3) + (4*im*a^2*asec(a + b*x)*atan(ℯ^(im*asec(a + b*x))))/b^3 - (2*a*log(a + b*x))/b^3 - (im*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/(3*b^3) - (2*im*a^2*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^3 + (im*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/(3*b^3) + (2*im*a^2*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^3]
    @test_int [x^1*asec(a + b*x)^2, x, 11, -(((a + b*x)*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x))/b^2) - (a^2*asec(a + b*x)^2)/(2*b^2) + (1/2)*x^2*asec(a + b*x)^2 - (4*im*a*asec(a + b*x)*atan(ℯ^(im*asec(a + b*x))))/b^2 + log(a + b*x)/b^2 + (2*im*a*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^2 - (2*im*a*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^2]
    @test_int [x^0*asec(a + b*x)^2, x, 8, ((a + b*x)*asec(a + b*x)^2)/b + (4*im*asec(a + b*x)*atan(ℯ^(im*asec(a + b*x))))/b - (2*im*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b + (2*im*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b]
    @test_int [asec(a + b*x)^2/x^1, x, 17, (-asec(a + b*x)^2)*log(1 + ℯ^(2*im*asec(a + b*x))) + im*asec(a + b*x)*PolyLog(2, -ℯ^(2*im*asec(a + b*x))) - (1/2)*PolyLog(3, -ℯ^(2*im*asec(a + b*x))) + asec(a + b*x)^2*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) + asec(a + b*x)^2*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) - 2*im*asec(a + b*x)*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) - 2*im*asec(a + b*x)*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) + 2*PolyLog(3, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) + 2*PolyLog(3, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2)))]
    @test_int [asec(a + b*x)^2/x^2, x, 12, -((b*asec(a + b*x)^2)/a) - asec(a + b*x)^2/x - (2*im*b*asec(a + b*x)*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (2*im*b*asec(a + b*x)*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) - (2*b*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (2*b*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2))]

    @test_int [x^2*asec(a + b*x)^3, x, 25, ((a + b*x)*asec(a + b*x))/b^3 - (3*im*a*asec(a + b*x)^2)/b^3 + (3*a*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x)^2)/b^3 - ((a + b*x)^2*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x)^2)/(2*b^3) + (a^3*asec(a + b*x)^3)/(3*b^3) + (1/3)*x^3*asec(a + b*x)^3 + (im*asec(a + b*x)^2*atan(ℯ^(im*asec(a + b*x))))/b^3 + (6*im*a^2*asec(a + b*x)^2*atan(ℯ^(im*asec(a + b*x))))/b^3 - atanh(sqrt(1 - 1/(a + b*x)^2))/b^3 + (6*a*asec(a + b*x)*log(1 + ℯ^(2*im*asec(a + b*x))))/b^3 - (im*asec(a + b*x)*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^3 - (6*im*a^2*asec(a + b*x)*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^3 + (im*asec(a + b*x)*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^3 + (6*im*a^2*asec(a + b*x)*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^3 - (3*im*a*PolyLog(2, -ℯ^(2*im*asec(a + b*x))))/b^3 + PolyLog(3, (-im)*ℯ^(im*asec(a + b*x)))/b^3 + (6*a^2*PolyLog(3, (-im)*ℯ^(im*asec(a + b*x))))/b^3 - PolyLog(3, im*ℯ^(im*asec(a + b*x)))/b^3 - (6*a^2*PolyLog(3, im*ℯ^(im*asec(a + b*x))))/b^3]
    @test_int [x^1*asec(a + b*x)^3, x, 16, (3*im*asec(a + b*x)^2)/(2*b^2) - (3*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*asec(a + b*x)^2)/(2*b^2) - (a^2*asec(a + b*x)^3)/(2*b^2) + (1/2)*x^2*asec(a + b*x)^3 - (6*im*a*asec(a + b*x)^2*atan(ℯ^(im*asec(a + b*x))))/b^2 - (3*asec(a + b*x)*log(1 + ℯ^(2*im*asec(a + b*x))))/b^2 + (6*im*a*asec(a + b*x)*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b^2 - (6*im*a*asec(a + b*x)*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b^2 + (3*im*PolyLog(2, -ℯ^(2*im*asec(a + b*x))))/(2*b^2) - (6*a*PolyLog(3, (-im)*ℯ^(im*asec(a + b*x))))/b^2 + (6*a*PolyLog(3, im*ℯ^(im*asec(a + b*x))))/b^2]
    @test_int [x^0*asec(a + b*x)^3, x, 10, ((a + b*x)*asec(a + b*x)^3)/b + (6*im*asec(a + b*x)^2*atan(ℯ^(im*asec(a + b*x))))/b - (6*im*asec(a + b*x)*PolyLog(2, (-im)*ℯ^(im*asec(a + b*x))))/b + (6*im*asec(a + b*x)*PolyLog(2, im*ℯ^(im*asec(a + b*x))))/b + (6*PolyLog(3, (-im)*ℯ^(im*asec(a + b*x))))/b - (6*PolyLog(3, im*ℯ^(im*asec(a + b*x))))/b]
    @test_int [asec(a + b*x)^3/x^1, x, 20, asec(a + b*x)^3*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) + asec(a + b*x)^3*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) - asec(a + b*x)^3*log(1 + ℯ^(2*im*asec(a + b*x))) - 3*im*asec(a + b*x)^2*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) - 3*im*asec(a + b*x)^2*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) + (3/2)*im*asec(a + b*x)^2*PolyLog(2, -ℯ^(2*im*asec(a + b*x))) + 6*asec(a + b*x)*PolyLog(3, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) + 6*asec(a + b*x)*PolyLog(3, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) - (3/2)*asec(a + b*x)*PolyLog(3, -ℯ^(2*im*asec(a + b*x))) + 6*im*PolyLog(4, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))) + 6*im*PolyLog(4, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))) - (3/4)*im*PolyLog(4, -ℯ^(2*im*asec(a + b*x)))]
    @test_int [asec(a + b*x)^3/x^2, x, 14, -((b*asec(a + b*x)^3)/a) - asec(a + b*x)^3/x - (3*im*b*asec(a + b*x)^2*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (3*im*b*asec(a + b*x)^2*log(1 - (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) - (6*b*asec(a + b*x)*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (6*b*asec(a + b*x)*PolyLog(2, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) - (6*im*b*PolyLog(3, (a*ℯ^(im*asec(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (6*im*b*PolyLog(3, (a*ℯ^(im*asec(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2))]

    #= ::Subsubsection:: =#
    #=n<0=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*asec(c+d*x^n))=#

    @test_int [x^1*(a + b*asec(c + d*x^2)), x, 7, (a*x^2)/2 + (b*(c + d*x^2)*asec(c + d*x^2))/(2*d) - (b*atanh(sqrt(1 - 1/(c + d*x^2)^2)))/(2*d)]

    @test_int [x^2*(a + b*asec(c + d*x^3)), x, 7, (a*x^3)/3 + (b*(c + d*x^3)*asec(c + d*x^3))/(3*d) - (b*atanh(sqrt(1 - 1/(c + d*x^3)^2)))/(3*d)]

    @test_int [x^3*(a + b*asec(c + d*x^4)), x, 7, (a*x^4)/4 + (b*(c + d*x^4)*asec(c + d*x^4))/(4*d) - (b*atanh(sqrt(1 - 1/(c + d*x^4)^2)))/(4*d)]

    @test_int [x^(n-1)*asec(a + b*x^n), x, 6, ((a + b*x^n)*asec(a + b*x^n))/(b*n) - atanh(sqrt(1 - 1/(a + b*x^n)^2))/(b*n)]

    #= ::Section::Closed:: =#
    #=Integrands*involving*inverse*secant*functions*of*exponentials=#

    @test_int [asec(c*ℯ^(a + b*x)), x, 7, (im*asec(c*ℯ^(a + b*x))^2)/(2*b) - (asec(c*ℯ^(a + b*x))*log(1 + ℯ^(2*im*asec(c*ℯ^(a + b*x)))))/b + (im*PolyLog(2, -ℯ^(2*im*asec(c*ℯ^(a + b*x)))))/(2*b)]

    #= ::Section::Closed:: =#
    #=Integrands*involving*exponentials*of*inverse*secant*functions=#

    @test_int [ℯ^asec(a*x)*x^2, x, 6, -(((12/5 + (4*im)/5)*ℯ^((1 + 3*im)*asec(a*x))*HypergeometricFunctions._₂F₁(3/2 - im/2, 3, 5/2 - im/2, -ℯ^(2*im*asec(a*x))))/a^3) + ((24/5 + (8*im)/5)*ℯ^((1 + 3*im)*asec(a*x))*HypergeometricFunctions._₂F₁(3/2 - im/2, 4, 5/2 - im/2, -ℯ^(2*im*asec(a*x))))/a^3]
    @test_int [ℯ^asec(a*x)*x^1, x, 6, -(((8/5 + (4*im)/5)*ℯ^((1 + 2*im)*asec(a*x))*HypergeometricFunctions._₂F₁(1 - im/2, 2, 2 - im/2, -ℯ^(2*im*asec(a*x))))/a^2) + ((16/5 + (8*im)/5)*ℯ^((1 + 2*im)*asec(a*x))*HypergeometricFunctions._₂F₁(1 - im/2, 3, 2 - im/2, -ℯ^(2*im*asec(a*x))))/a^2]
    @test_int [ℯ^asec(a*x)*x^0, x, 5, -(((1 + im)*ℯ^((1 + im)*asec(a*x))*HypergeometricFunctions._₂F₁(1/2 - im/2, 1, 3/2 - im/2, -ℯ^(2*im*asec(a*x))))/a) + ((2 + 2*im)*ℯ^((1 + im)*asec(a*x))*HypergeometricFunctions._₂F₁(1/2 - im/2, 2, 3/2 - im/2, -ℯ^(2*im*asec(a*x))))/a]
    @test_int [ℯ^asec(a*x)/x^1, x, 6, (-im)*ℯ^asec(a*x) + 2*im*ℯ^asec(a*x)*HypergeometricFunctions._₂F₁(-(im/2), 1, 1 - im/2, -ℯ^(2*im*asec(a*x)))]
    @test_int [ℯ^asec(a*x)/x^2, x, 3, (1/2)*a*ℯ^asec(a*x)*sqrt(1 - 1/(a^2*x^2)) - ℯ^asec(a*x)/(2*x)]
    @test_int [ℯ^asec(a*x)/x^3, x, 5, (-(1/5))*a^2*ℯ^asec(a*x)*cos(2*asec(a*x)) + (1/10)*a^2*ℯ^asec(a*x)*sin(2*asec(a*x))]
    @test_int [ℯ^asec(a*x)/x^4, x, 6, (1/8)*a^3*ℯ^asec(a*x)*sqrt(1 - 1/(a^2*x^2)) - (a^2*ℯ^asec(a*x))/(8*x) - (3/40)*a^3*ℯ^asec(a*x)*cos(3*asec(a*x)) + (1/40)*a^3*ℯ^asec(a*x)*sin(3*asec(a*x))]

    #= ::Section::Closed:: =#
    #=Miscellaneous*integrands*involving*inverse*secants=#

    @test_int [asec(a + b*x)/((a*d)/b + d*x), x, 8, (im*asec(a + b*x)^2)/(2*d) - (asec(a + b*x)*log(1 + ℯ^(2*im*asec(a + b*x))))/d + (im*PolyLog(2, -ℯ^(2*im*asec(a + b*x))))/(2*d)]
    end

max_stars_repo_path: test/5 Inverse trig functions/5.5 Inverse secant/5.5.2 Inverse secant functions.jl


text:
abstract type Results end

get_results(result::Results) = nothing
struct OperationsProblemResults <: Results
    variables::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
end

get_duals(result::OperationsProblemResults) = nothing
"""This function creates the correct results struct for the context"""
function _make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::DataFrames.DataFrame)
    return OperationsProblemResults(variables, total_cost, optimizer_log, time_stamp)
end
"""This function creates the correct results struct for the context"""
function _make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::Array)
    time_stamp = DataFrames.DataFrame(Range = time_stamp)
    return OperationsProblemResults(variables, total_cost, optimizer_log, time_stamp)
end
"""This function creates the correct results struct for the context"""
function _make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::Array,
                      duals::Dict)
    time_stamp = DataFrames.DataFrame(Range = time_stamp)
    return AggregatedResults(variables, total_cost, optimizer_log, time_stamp, duals)
end
function get_variable(res_model::OperationsProblemResults, key::Symbol)
        try
            !isnothing(results.variables)
            return get(results.variables, key, nothing)
        catch
            throw(ArgumentError("No variable with key $(key) has been found."))
        end
end

function get_optimizer_log(results::OperationsProblemResults)
    return results.optimizer_log
end

function get_time_stamps(results::OperationsProblemResults, key::Symbol)
    return results.time_stamp
end

"""
    results = load_operation_results(path)

This function can be used to load results from a folder
of results from a single-step problem, or for a single foulder
within a simulation.

# Arguments
- `path::AbstractString = folder path`
- `directory::AbstractString = "2019-10-03T09-18-00"`: the foulder name that contains
feather files of the results.

# Example
```julia
results = load_operation_results("/Users/test/2019-10-03T09-18-00")
```
"""
function load_operation_results(folder_path::AbstractString)

    if isfile(folder_path)
        throw(ArgumentError("Not a folder path."))
    end
    files_in_folder = collect(readdir(folder_path))
    variable_list = setdiff(files_in_folder, ["time_stamp.feather", "optimizer_log.json", "check_sum.json"])
    variables = Dict{Symbol, DataFrames.DataFrame}()
    for name in variable_list
        variable_name = splitext(name)[1]
        file_path = joinpath(folder_path, name)
        variables[Symbol(variable_name)] = Feather.read(file_path)
    end
    optimizer = JSON.parse(open(joinpath(folder_path, "optimizer_log.json")))
    time_stamp = Feather.read(joinpath(folder_path, "time_stamp.feather"))
    if size(time_stamp, 1) > find_var_length(variables, variable_list)
        time_stamp = shorten_time_stamp(time_stamp)
    end
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer["obj_value"])
    if any(x -> x == "check_sum.json", files_in_folder)
        check_sum = JSON.parse(open(joinpath(folder_path, "check_sum.json")))
        results = _make_results(variables, obj_value, optimizer, time_stamp, check_sum)
    else
        results = _make_results(variables, obj_value, optimizer, time_stamp)
    end
    return results
end

# this ensures that the time_stamp is not double shortened
function find_var_length(variables::Dict, variable_list::Array)
    return size(variables[Symbol(splitext(variable_list[1])[1])], 1)
end

function shorten_time_stamp(time::DataFrames.DataFrame)
    time = time[1:(size(time, 1)-1),:]
    return time
end

ext: jl


text:
@testset "allapproxequal()" begin
    @test nsi.allapproxequal([])                     == true
    @test nsi.allapproxequal([1.0])                  == true
    @test nsi.allapproxequal([1.00000001, 1.0, 1.0]) == true
    # tolerance is 1e-8          1234567
    @test nsi.allapproxequal([ 1.0000001, 1.0, 1.0]) == false
    @test nsi.allapproxequal([ 1, 1, 1])             == true
    @test nsi.allapproxequal([ 2, 1, 1])             == false
    @test_throws MethodError nsi.allapproxequal([ "A", "A"])
end

@testset "verify_same_length()" begin
    @test nsi.verify_same_length([], [])                 === nothing
    @test nsi.verify_same_length(["A"], ["B"])           === nothing
    @test nsi.verify_same_length(["A", "B"], ["C", "D"]) === nothing
    @test nsi.verify_same_length([1], [2])               === nothing
    @test nsi.verify_same_length([1, 2], [3, 4])         === nothing
    @test nsi.verify_same_length([1 2; 3 4], [5 6; 7 8]) === nothing
    @test_throws ArgumentError nsi.verify_same_length([1], [])
    @test_throws ArgumentError nsi.verify_same_length(["A"], [])
    @test_throws ArgumentError nsi.verify_same_length(["A", "B"], ["C"])
    @test_throws ArgumentError nsi.verify_same_length([1], [])
    @test_throws ArgumentError nsi.verify_same_length([1, 2], [3])
    @test_throws ArgumentError nsi.verify_same_length([1 2; 3 4], [5 6;])
end

@testset "left_right_from_peak()" begin
    x = collect(Float64, 1:20)
    y = zeros(Float64, 20)
    y[10] = 1
    @test nsi.left_right_from_peak(x, y, 10, 1) == [9.5, 10.5]
    @test nsi.left_right_from_peak(x, y, 10, 0.7) == [9.65, 10.35]
    y[10] = 0
    y[1]  = 1
    @test nsi.left_right_from_peak(x, y, 1, 1) == [0.5, 1.5]
    y[1] = 0
    y[20]  = 1
    @test nsi.left_right_from_peak(x, y, 20, 1) == [19.5, 20.5]
end

@testset "trapz()" begin
    # Test integration function trapz using analytic integrals;
    # integral functions (F) of functions f from
    # https://en.wikipedia.org/wiki/Lists_of_integrals

    bounds = [(-1, 4), (1, 4), (-1.4, 7.9), (11.6, 17.4)]

    test_data = [
        (f = x->  exp(3x), F = x->        1/3*exp(3x)),
        (f = x->  1.2^(x), F = x->   1.2^(x)/log(1.2)),
        (f = x->   sin(x), F = x->            -cos(x)),
        (f = x-> 1/(2x+3), F = x-> 1/2*log(abs(2x+3)))
    ]
    x = collect(-20:0.0001:20)
    for t in test_data
        for (l, r) in bounds
            f, F = t
            ref_1 = F(r) - F(l)
            ref_2 = ref_1 - 1/2*(f(r) + f(l))*(r - l) # 1/2*(f(r) + f(l))*(r - l) term: subtracting baseline
            @test nsi.trapz(x, f.(x), l, r) ≈ ref_1 rtol=1e-7
        end
    end

    begin
        x = collect(Float64, 1:10)
        y = collect(Float64, 1:10)
        @test nsi.trapz(x, y, 1, 3) == 4.0
        @test nsi.trapz(x, y, 1.1, 2.1) ≈ 1.6
        @test nsi.trapz(x, y, 1.1, 1.2) ≈ 1.1*0.1 + 0.1*0.1*0.5
        @test nsi.trapz(x, y, 1.0, 1.1) ≈ 0.1 + 0.1*0.1*0.5
        @test nsi.trapz(x, y, 3.9, 4.0) ≈ 3.9*0.1 + 0.1*0.1*0.5
    end

    begin
        x = [1, 2, 3, 4, 6, 7.0]
        y = [1, 2, 4, 2.5, 3, 2]
        @test nsi.trapz(x, y, 1, 7) == 15.75
        @test nsi.trapz(x, y, 0, 10) == 15.75
        @test nsi.trapz(x, y, 1, 2) == 1.5
        @test nsi.trapz(x, y, 3, 4) == 3.25
        @test nsi.trapz(x, y, 1, 6) == 13.25
        @test nsi.trapz(x, y, 6, 1) == 13.25
        @test nsi.trapz(x, y, 1.5, 7) == 15.125
    end
end

@testset "@samples" begin
    seed!(1)
    @testset "Normal" begin
        ref = Particles(10_000, Normal(1, 1))
        out = @samples 10_000 1 ± 1
        @test mean(ref) ≈ mean(out)
        @test std(ref) ≈ std(out)
        # pass expression to macro
        out = @samples 10_000 (2 / 2) ± (3 - 2)
        @test mean(ref) ≈ mean(out)
        @test std(ref) ≈ std(out)
        # pass variables to macro
        out = let
            a = 1
            b = 1
            @samples 10_000 a ± b
        end
        @test mean(ref) ≈ mean(out)
        @test std(ref) ≈ std(out)
    end
    @testset "Uniform" begin
        ref = Particles(10_000, Uniform(1, 2))
        out = @samples 10_000 1 .. 2
        @test mean(ref) ≈ mean(out)
        @test std(ref) ≈ std(out)
        # pass expression to macro
        out = @samples 10_000 (2 / 2) .. (3 - 1)
        @test mean(ref) ≈ mean(out)
        @test std(ref) ≈ std(out)
        # pass variables to macro
        out = let
            a = 1
            b = 2
            @samples 10_000 a .. b
        end
        @test mean(ref) ≈ mean(out)
        @test std(ref) ≈ std(out)
    end
    @test_throws ArgumentError (@samples 1 1 + 1)
    @test_throws ArgumentError (@samples 1 1 + 1 + 1)
end

@testset "lininterp" begin
    xs = 0:0.5:10 |> collect
    ys = 0:0.5:10 |> collect
    @test NoisySignalIntegration.lininterp(2.7, xs, ys) == 2.7
    ys = ys.^2
    @test NoisySignalIntegration.lininterp(2.7, xs, ys) ≈ 7.35
    @test NoisySignalIntegration.lininterp(2.5, xs, ys) == 6.25
    @test NoisySignalIntegration.lininterp(10.0, xs, ys) == 100.0
    @test_throws ErrorException NoisySignalIntegration.lininterp(0.0, xs, ys)
    @test_throws ErrorException NoisySignalIntegration.lininterp(10.1, xs, ys)
end


text:
using BenchmarkTools
using Fed
using Random
using SHA

const suite = BenchmarkGroup()

suite["payload_serde"] = BenchmarkGroup(["serialize", "deserialize"])

# config
const SERVERURL = "http://127.0.0.1:8080"
const NUM_COMM_ROUNDS = 100
const FRACTION_CLIENTS = 0.1f0
const NUM_TOTAL_CLIENTS = 100

const base_config = Fed.Config.BaseConfig(
    SERVERURL,
    NUM_COMM_ROUNDS,
    FRACTION_CLIENTS,
    NUM_TOTAL_CLIENTS
)

weights = rand(Float32, Int(1024*1024/4)) # 1MB
chunksize = 256

# Vanilla
vanilla_config = Fed.Config.VanillaConfig(base_config)
initialize_stats(vanilla_config, length(weights))
vanilla_payload = Fed.Serde.serialize_payload(vanilla_config.payload_serde, weights)

suite["payload_serde"]["vanilla_serialize"] = @benchmarkable Fed.Serde.serialize_payload(vanilla_config.payload_serde, weights)
suite["payload_serde"]["vanilla_deserialize"] = @benchmarkable Fed.Serde.deserialize_payload(vanilla_config.payload_serde, vanilla_payload, "")

# quantized
quantized_config = Fed.Config.QuantizedConfig{UInt8}(base_config)
initialize_stats(quantized_config, length(weights))
quantized_payload = Fed.Serde.serialize_payload(quantized_config.payload_serde, weights)

suite["payload_serde"]["quantized_serialize"] = @benchmarkable Fed.Serde.serialize_payload(quantized_config.payload_serde, weights)
suite["payload_serde"]["quantized_deserialize"] = @benchmarkable Fed.Serde.deserialize_payload(quantized_config.payload_serde, quantized_payload, "")

# qdiff
qdiff_config = Fed.Config.QDiffConfig{UInt8}(base_config, 256, false)
initialize_stats(qdiff_config, length(weights))
old_weights = deepcopy(weights)
old_weights[1:131072]  = rand(Float32, 131072)
old_weights = shuffle(old_weights)
qdiff_payload = Fed.Serde.serialize_payload(qdiff_config.payload_serde, old_weights)

qdiff_config.payload_serde.is_client = true
suite["payload_serde"]["qdiff_serialize"] = @benchmarkable Fed.Serde.serialize_payload(qdiff_config.payload_serde, weights)

qdiff_config.payload_serde.is_client = false
suite["payload_serde"]["qdiff_deserialize"] = @benchmarkable Fed.Serde.deserialize_payload(qdiff_config.payload_serde, qdiff_payload, "")

# gd
gd_config = Fed.Config.GDConfig{UInt8}(
    base_config,
    chunksize,
    sha1,
    # hash_crc32,
    round(UInt8, 0.6 * sizeof(UInt8) * 8),
    "127.0.0.1",
    9090,
    false
)

initialize_stats(gd_config, length(weights))
old_weights = deepcopy(weights)
old_weights[1:131072]  = rand(Float32, 131072)
gd_payload = Fed.Serde.serialize_payload(gd_config.payload_serde, old_weights)

gd_config.payload_serde.is_client = true
suite["payload_serde"]["gd_serialize"] = @benchmarkable Fed.Serde.serialize_payload(gd_config.payload_serde, weights)

gd_config.payload_serde.is_client = false
suite["payload_serde"]["gd_deserialize"] = @benchmarkable Fed.Serde.deserialize_payload(gd_config.payload_serde, gd_payload, "http://127.0.0.1:9090")

# Saving params

paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end


text:
abstract type KernelType end;

abstract type ARDKernel <: KernelType end;

abstract type IsoKernel <: KernelType end;

abstract type KernelCombination <: KernelType end;

"Macro giving common fields for all kernels"
mutable struct KernelFields{T,KT}
    name::String #Name of the kernel function
    variance::HyperParameter{T} #variance of the kernel
    lengthscales::HyperParameters{T} #lengthscale variables
    Ndim::Int64 #Number of lengthscales
    metric::PreMetric
end

# params::HyperParameters{T} #Parameters of the kernel
# Nparam::Int64 #Number of parameters
"""Abstract type for all kernels"""
abstract type Kernel{T<:Real, KT<:KernelType} end;

"""Return the metric of a kernel"""
function getmetric(k::Kernel{T,KT}) where {T<:Real,KT<:KernelType}
    return k.fields.metric
end

"""Return the variance of the kernel"""
function getvariance(k::Kernel{T,KT}) where {T<:Real,KT<:KernelType}
    return getvalue(k.fields.variance)
end

"""Return the lengthscale of the IsoKernel"""
function getlengthscales(k::Kernel{T,IsoKernel}) where {T<:Real}
    return getvalue(k.fields.lengthscales[1])
end

"""Return the lengthscales of the ARD Kernel"""
function getlengthscales(k::Kernel{T,ARDKernel}) where {T<:Real}
    return getvalue(k.fields.lengthscales)
end

"""Set optimizers of all parameters of the kernel"""
function setoptimizer!(k::Kernel{T},opt::Optimizer) where {T<:Real}
    setparamoptimizer!(k.fields.variance,opt)
    setparamoptimizer!(k.fields.lengthscales,opt)
end

"""Get optimizer for variance of kernel"""
function getvarianceoptimizer(k::Kernel{T}) where {T<:Real}
    k.fields.variance.opt
end

isARD(::Kernel{T,KT}) where {T<:Real,KT<:ARDKernel} = true
isARD(::Kernel{T,KT}) where {T<:Real,KT<:IsoKernel} = false
isIso(::Kernel{T,KT}) where {T<:Real,KT<:ARDKernel} = false
isIso(::Kernel{T,KT}) where {T<:Real,KT<:IsoKernel} = true

function Base.show(io::IO,k::Kernel{T,KT}) where {T,KT}
    print(io,"$(k.fields.name)"*(isARD(k) ? " ARD" : "")*" kernel, with variance $(getvariance(k)) and lengthscales $(getlengthscales(k))")
end

"""Standard conversion when giving scalar and not vectors"""
function compute(k::Kernel{T,KT},X1::T,X2::T) where {T<:Real,KT<:KernelType}
    compute(k,[X1],[X2])
end

"""Function to determine most adapted type between a selection"""
function floattype(T_i::DataType...)
    T_max = promote_type(T_i...)
    T_max <: Real ? T_max : Float64
end

function convert(::Type{K},x::Kernel{T2,KT}) where {K<:Kernel{T} where T<:Real,T2<:Real,KT<:KernelType}

end

function apply_gradients_lengthscale!(kernel::Kernel{T,IsoKernel},gradient::T) where {T}
    update!(kernel.fields.lengthscales,[gradient])
end

function apply_gradients_lengthscale!(kernel::Kernel{T,ARDKernel},gradients::Vector{T}) where {T}
    update!(kernel.fields.lengthscales,gradients)
    updateweights!(kernel,getlengthscales(kernel))
end

function apply_gradients_variance!(kernel::Kernel{T,KT},gradient::T) where {T,KT}
    update!(kernel.fields.variance,gradient)
end
#
# function apply_gradients!(kernel::KernelSum,gradients,variance::Bool=true)
#     for i in 1:kernel.Nkernels
#         apply_gradients!(kernel.kernel_array[i],gradients[i],true)
#     end
# end
#
# function apply_gradients!(kernel::KernelProduct,gradients,variance::Bool=true)
#     for i in 1:kernel.Nkernels
#         apply_gradients!(kernel.kernel_array[i],gradients[i],false);
#     end
#     if variance
#         update!(kernel.variance,gradients[end][1]);
#     end
# end


text:
########################
# Completion Functions #
########################
function complete_activate(options, partial, i1, i2)
    possible = String[]
    shared = get(options, :shared, false)
    if shared
        for depot in Base.DEPOT_PATH
            envdir = joinpath(depot, "environments")
            isdir(envdir) || continue
            append!(possible, readdir(envdir))
        end
        return possible
    else
        return complete_local_dir(partial, i1, i2)
    end
end

function complete_local_dir(s, i1, i2)
    if !isempty(s) && s[1] == '~'
        return String[expanduser(s)], i1:i2, true
    end

    cmp = REPL.REPLCompletions.complete_path(s, i2)
    completions = [REPL.REPLCompletions.completion_text(p) for p in cmp[1]]
    completions = filter!(x -> isdir(s[1:prevind(s, first(cmp[2])-i1+1)]*x), completions)
    return completions, cmp[2], !isempty(completions)
end

function complete_remote_package(partial)
    cmp = String[]
    julia_version = VERSION
    for reg in Types.collect_registries()
        data = Types.read_registry(joinpath(reg.path, "Registry.toml"))
        for (uuid, pkginfo) in data["packages"]
            name = pkginfo["name"]
            if startswith(name, partial)
                compat_data = Operations.load_package_data_raw(
                    VersionSpec, joinpath(reg.path, pkginfo["path"], "Compat.toml"))
                supported_julia_versions = VersionSpec()
                found_julia_compat = false
                for (ver_range, compats) in compat_data
                    for (compat, v) in compats
                        if compat == "julia"
                            found_julia_compat = true
                            union!(supported_julia_versions, VersionSpec(v))
                        end
                    end
                end
                if VERSION in supported_julia_versions || !found_julia_compat
                    push!(cmp, name)
                end
            end
        end
    end
    return cmp
end

function complete_help(options, partial)
    names = String[]
    for cmds in values(super_specs)
         append!(names, [spec.canonical_name for spec in values(cmds)])
    end
    return sort!(unique!(append!(names, collect(keys(super_specs)))))
end

function complete_installed_packages(options, partial)
    env = try EnvCache()
    catch err
        err isa PkgError || rethrow()
        return String[]
    end
    mode = get(options, :mode, PKGMODE_PROJECT)
    return mode == PKGMODE_PROJECT ?
        collect(keys(env.project.deps)) :
        unique!([entry.name for (uuid, entry) in env.manifest])
end

function complete_add_dev(options, partial, i1, i2)
    comps, idx, _ = complete_local_dir(partial, i1, i2)
    if occursin(Base.Filesystem.path_separator_re, partial)
        return comps, idx, !isempty(comps)
    end
    comps = vcat(comps, complete_remote_package(partial))
    comps = vcat(comps, filter(x->startswith(x,partial) && !(x in comps),
                               collect(values(Types.stdlib()))))
    return comps, idx, !isempty(comps)
end

########################
# COMPLETION INTERFACE #
########################
function default_commands()
    names = collect(keys(super_specs))
    append!(names, map(x -> getproperty(x, :canonical_name), values(super_specs["package"])))
    return sort(unique(names))
end

function complete_command(statement::Statement, final::Bool, on_sub::Bool)
    if statement.super !== nothing
        if (!on_sub && final) || (on_sub && !final)
            # last thing determined was the super -> complete canonical names of subcommands
            specs = super_specs[statement.super]
            names = map(x -> getproperty(x, :canonical_name), values(specs))
            return sort(unique(names))
        end
    end
    # complete default names
    return default_commands()
end

complete_opt(opt_specs) =
    unique(sort(map(wrap_option,
                    map(x -> getproperty(x, :name),
                        collect(values(opt_specs))))))

function complete_argument(spec::CommandSpec, options::Vector{String},
                           partial::AbstractString, offset::Int,
                           index::Int)
    spec.completions === nothing && return String[]
    # finish parsing opts
    opts = APIOptions(map(parse_option, options), spec.option_specs)
    return applicable(spec.completions, opts, partial, offset, index) ?
        spec.completions(opts, partial, offset, index) :
        spec.completions(opts, partial)
end

function _completions(input, final, offset, index)
    words = tokenize(input)[end]
    word_count = length(words)
    statement, partial = core_parse(words)
    final && (partial = "") # last token is finalized -> no partial
    # number of tokens which specify the command
    command_size = count([statement.preview, statement.super !== nothing, true])
    command_is_focused() = !((word_count == command_size && final) || word_count > command_size)

    if statement.spec === nothing # spec not determined -> complete command
        !command_is_focused() && return String[], 0:-1, false
        x = complete_command(statement, final, word_count == (statement.preview ? 3 : 2))
    else
        command_is_focused() && return String[], 0:-1, false

        if final # complete arg by default
            x = complete_argument(statement.spec, statement.options, partial, offset, index)
        else # complete arg or opt depending on last token
            x = is_opt(partial) ?
                complete_opt(statement.spec.option_specs) :
                complete_argument(statement.spec, statement.options, partial, offset, index)
        end
    end

    # In the case where the completion function wants to deal with indices, it will return a fully
    # computed completion tuple, just return it
    # Else, the completions function will just deal with strings and will return a Vector{String}
    if isa(x, Tuple)
        return x
    else
        possible = filter(possible -> startswith(possible, partial), x)
        return possible, offset:index, !isempty(possible)
    end
end

function completions(full, index)::Tuple{Vector{String},UnitRange{Int},Bool}
    pre = full[1:index]
    isempty(pre) && return default_commands(), 0:-1, false # empty input -> complete commands
    last   = split(pre, ' ', keepempty=true)[end]
    offset = isempty(last) ? index+1 : last.offset+1
    final  = isempty(last) # is the cursor still attached to the final token?
    return _completions(pre, final, offset, index)
end


text:
using BinaryBuilder

name = "Arrow"
version = v"0.13.0"
sources = [
    "http://mirror.cogentco.com/pub/apache/arrow/arrow-0.13.0/apache-arrow-0.13.0.tar.gz" =>
        "ac2a77dd9168e9892e432c474611e86ded0be6dfe15f689c948751d37f81391a"
]

script = raw"""
    cd ${WORKSPACE}/srcdir/apache-arrow-*

    echo 'Building C++ library (required for C library)...'
    cd cpp
    mkdir release
    cd release
    cmake .. \
        -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DCMAKE_TOOLCHAIN_FILE=/opt/${target}/${target}.toolchain \
        -DARROW_DEPENDENCY_SOURCE=BUNDLED \
        -DARROW_COMPUTE=OFF \
        -DARROW_IPC=OFF \
        -DARROW_PARQUET=ON \
        -DBOOST_ROOT=${WORKSPACE}/destdir
    make -j${nproc}
    make -j${nproc} parquet
    make install

    echo 'Building C library...'
    cd ../../c_glib
    ./configure
        --prefix=${prefix} \
        --build=${target} \
        --with-sysroot=/opt/${target}/${target}/sys-root \
        --disable-static \
        --disable-glibtest \
        --disable-gtk-doc
    make -j${nproc}
    make install
    """

platforms = [Linux(:x86_64)]

products(prefix) = [
    LibraryProduct(prefix, "libarrow-glib", :libarrow),
    LibraryProduct(prefix, "libparquet-glib", :libparquet),
]

dependencies = [
    "https://github.com/staticfloat/GlibBuilder/releases/download/v2.54.2-2/build.jl",
    "https://github.com/twadleigh/BoostBuilder/releases/download/v1.68.0-4/build_Boost.v1.68.0.jl",
    "https://github.com/BenJWard/ZLibBuilder/releases/download/v1.2.11/build_ZLibBuilder.v1.2.11.jl",
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: ararslan/ArrowBuilder

max_issues_repo_name: ararslan/ArrowBuilder

max_forks_repo_name: ararslan/ArrowBuilder


text:
module TestSparseRollCall
    using Base.Test
    import RollCallDataIO

    encoding = RollCallDataIO.RollCallEncoding()
    votes = [1 0 0; 0 6 6;]
    roll_call = RollCallDataIO.RollCall(votes, encoding)

    new_roll_call = convert(RollCallDataIO.SparseRollCall, roll_call)
    @test isa(new_roll_call, RollCallDataIO.SparseRollCall)

    @test new_roll_call.n_legislators === 2
    @test new_roll_call.n_bills === 3
    @test new_roll_call.legislators == Int[1, 2, 2]
    @test new_roll_call.bills == Int[1, 2, 3]
    @test new_roll_call.votes == Int[1, 0, 0]
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
# This is deprecated and will be removed soon.
# WebIO+Blink integration is now implemented in Blink.jl.
# https://github.com/JunoLab/Blink.jl/pull/201

using AssetRegistry
using Base64: stringmime

using Sockets

struct BlinkConnection <: WebIO.AbstractConnection
    page::Blink.Page
end

Blink.body!(p::Blink.Page, x::AbstractWidget) = Blink.body!(p, Widgets.render(x))
function Blink.body!(p::Blink.Page, x::Union{Node, Scope})
    wait(p)

    bs = AssetRegistry.register(BLINK_BUNDLE_PATH)
    Blink.loadjs!(p, bs)

    conn = BlinkConnection(p)
    Blink.handle(p, "webio") do msg
        WebIO.dispatch(conn, msg)
    end

    Blink.body!(p, stringmime(MIME"text/html"(), x))
end

Blink.body!(p::Blink.Window, x::AbstractWidget) = Blink.body!(p, Widgets.render(x))
function Blink.body!(p::Blink.Window, x::Union{Node, Scope})
    Blink.body!(p.content, x)
end

function Sockets.send(b::BlinkConnection, data)
    Blink.msg(b.page, Dict(:type=>"webio", :data=>data))
end

Base.isopen(b::BlinkConnection) = Blink.active(b.page)

function WebIO.register_renderable(T::Type, ::Val{:blink})
    eval(:(Blink.body!(p::Union{Blink.Window, Blink.Page}, x::$T) =
           Blink.body!(p, WebIO.render(x))))
end

WebIO.setup_provider(::Val{:blink}) = nothing  # blink setup has no side-effects
WebIO.setup(:blink)


text:
using LinearAlgebra, SnpArrays, SparseArrays, Test

const EUR = SnpArray(SnpArrays.datadir("EUR_subset.bed")) # no missing genotypes
const mouse = SnpArray(SnpArrays.datadir("mouse.bed")) # has missing genotypes

@testset "size" begin
@test size(EUR) == (379, 54051)
@test size(EUR, 1) == 379
@test size(EUR, 2) == 54051
@test size(mouse) == (1940, 10150)
@test size(mouse, 3) == 1
@test_throws ErrorException size(mouse, 0)
@test length(mouse) == 1940 * 10150
@test eltype(mouse) == UInt8
end

@testset "counts" begin
cc = counts(EUR, dims=1)
@test counts(EUR) == sum(cc, dims=2)
@test size(cc) == (4, size(EUR, 2))
@test view(cc, :, 1) == [2, 0, 70, 307]
@test all(sum(cc, dims = 1) .== size(EUR, 1))
@test all(iszero, view(cc, 2, :))
rc = counts(EUR, dims=2)
@test size(rc) == (4, size(EUR, 1))
@test view(rc, :, 1) == [2997, 0, 13143, 37911]
@test all(sum(rc, dims = 1) .== size(EUR, 2))
@test all(iszero, view(rc, 2, :))
@test_throws ArgumentError counts(EUR, dims=3)
end

@testset "means" begin
# ADDITIVE_MODEL
cmns = mean(EUR, dims = 1)
@test size(cmns) == (1, size(EUR, 2))
@test cmns[1] ≈ 1.804749340369393
@test mean(cmns) ≈ mean(EUR)
@test minimum(cmns) ≈ 1.0
@test maximum(cmns) ≈ 1.9788918205804749
rmns = mean(EUR, dims = 2)
@test size(rmns) == (size(EUR, 1), 1)
@test rmns[1] ≈ 1.6459454959205195
@test mean(rmns) ≈ mean(EUR)
rmnextrma = extrema(rmns)
@test rmnextrma[1] ≈ 1.637749532848606
@test rmnextrma[2] ≈ 1.6627259440158368
@test_throws ArgumentError mean(EUR, dims=3)
# DOMINANT_MODEL
cmns = mean(EUR, dims = 1, model = DOMINANT_MODEL)
@test size(cmns) == (1, size(EUR, 2))
@test cmns[1] ≈ 0.9947229551451188
@test mean(cmns) ≈ mean(EUR, model = DOMINANT_MODEL)
@test minimum(cmns) ≈ 0.5989445910290238
@test maximum(cmns) ≈ 1.0
rmns = mean(EUR, dims = 2, model = DOMINANT_MODEL)
@test size(rmns) == (size(EUR, 1), 1)
@test rmns[1] ≈ 0.9445523672087472
@test mean(rmns) ≈ mean(EUR, model = DOMINANT_MODEL)
rmnextrma = extrema(rmns)
@test rmnextrma[1] ≈ 0.9325266877578583
@test rmnextrma[2] ≈ 0.9530073449149877
@test_throws ArgumentError mean(EUR, dims=3, model=DOMINANT_MODEL)
# RECESSIVE_MODEL
cmns = mean(EUR, dims = 1, model = RECESSIVE_MODEL)
@test size(cmns) == (1, size(EUR, 2))
@test cmns[1] ≈ 0.8100263852242744
@test mean(cmns) ≈ mean(EUR, model = RECESSIVE_MODEL)
@test minimum(cmns) ≈ 0.0
@test maximum(cmns) ≈ 0.9868073878627969
rmns = mean(EUR, dims = 2, model = RECESSIVE_MODEL)
@test size(rmns) == (size(EUR, 1), 1)
@test rmns[1] ≈ 0.7013931287117722
@test mean(rmns) ≈ mean(EUR, model = RECESSIVE_MODEL)
rmnextrma = extrema(rmns)
@test rmnextrma[1] ≈ 0.6932341677304767
@test rmnextrma[2] ≈ 0.7191171301178516
@test_throws ArgumentError mean(EUR, dims=3, model=RECESSIVE_MODEL)
end

@testset "var/std" begin
cvars = var(EUR, dims = 1)
@test size(cvars) == (1, size(EUR, 2))
@test cvars[1] ≈ 0.16812553224162724
@test iszero(minimum(cvars))
@test maximum(cvars) ≈ 0.8705727966941688
@test std(EUR, dims = 1) ≈ sqrt.(cvars)
rvars = var(EUR, dims = 2)
@test size(rvars) == (size(EUR, 1), 1)
@test rvars[1] ≈ 0.33960146078503206
@test minimum(rvars) ≈ 0.31972707244268267
@test maximum(rvars) ≈ 0.365901927927013
@test std(EUR, dims = 2) ≈ sqrt.(rvars)
end

@testset "missingpos" begin
mp = missingpos(mouse)
cc = counts(mouse, dims=1)
@test isa(mp, SparseMatrixCSC)
@test sum(mp, dims = 1) == view(cc, 2:2, :)
@test sum(mp, dims = 2) == view(counts(mouse, dims=2), 2:2, :)'
end

@testset "create bed" begin
tmpbf = SnpArray("tmp.bed", 5, 3)
@test isfile("tmp.bed")
@test all(tmpbf .== 0x00)
fill!(tmpbf, 0x02)
tmpbf2 = SnpArray("tmp.bed", 5)
@test all(tmpbf2 .== 0x02)
Sys.iswindows() || rm("tmp.bed", force=true)
tmpbf2 = SnpArray("tmp2.bed", SnpArray(undef, 5, 3))
fill!(tmpbf2, 0x01)
@test all(tmpbf2 .== 0x01)
Sys.iswindows() || rm("tmp2.bed", force=true)
end

@testset "convert" begin
# check convert coding
tmpbf = SnpArray(undef, 4, 1)
tmpbf[1] = 0x00
tmpbf[2] = 0x01
tmpbf[3] = 0x02
tmpbf[4] = 0x03
# additive model
v = convert(Matrix{Float64}, tmpbf, model=ADDITIVE_MODEL)
@test all([v[1], v[3], v[4]] .== [0.0, 1.0, 2.0])
@test isnan(v[2])
# dominant model
v = convert(Matrix{Float64}, tmpbf, model=DOMINANT_MODEL)
@test all([v[1], v[3], v[4]] .== [0.0, 1.0, 1.0])
@test isnan(v[2])
# recessive model
v = convert(Matrix{Float64}, tmpbf, model=RECESSIVE_MODEL)
@test all([v[1], v[3], v[4]] .== [0.0, 0.0, 1.0])
@test isnan(v[2])
# mouse test data
cmns = mean(mouse, dims = 1)
cvrs = var(mouse, dims = 1)
cstd = std(mouse, dims = 1)
v = convert(Vector{Float64}, @view(mouse[:, 1]))
@test isnan(mean(v))
@test mean(filter(!isnan, v)) ≈ cmns[1]
@test isnan(var(v))
@test var(filter(!isnan, v)) ≈ cvrs[1]
@test isnan(std(v))
@test std(filter(!isnan, v)) ≈ cstd[1]
end

@testset "convert (ADMIXTURE)" begin
# convert A2 allele frequency to dosage
@test convert(Float64, 0.5, ADDITIVE_MODEL)  == 1.0
@test convert(Float64, 0.5, DOMINANT_MODEL)  == 0.75
@test convert(Float64, 0.5, RECESSIVE_MODEL) == 0.25
end

@testset "maf" begin
cc = counts(mouse, dims=1)
@test all(0 .≤ maf(mouse) .≤ 0.5)
@test all(minorallele(mouse) .== (cc[1, :] .> cc[4, :]))
cc = counts(EUR, dims=1)
@test all(0 .≤ maf(EUR) .≤ 0.5)
@test all(minorallele(EUR) .== (cc[1, :] .> cc[4, :]))
end

@testset "grm" begin
Φgrm = grm(EUR, method=:GRM)
@test size(Φgrm) == (size(EUR, 1), size(EUR, 1))
@test issymmetric(Φgrm)
@test eigmin(Φgrm) > -1e-8
Φmom = grm(EUR, method=:MoM)
@test size(Φmom) == (size(EUR, 1), size(EUR, 1))
@test issymmetric(Φmom)
Φrbs = grm(EUR, method=:Robust)
@test size(Φrbs) == (size(EUR, 1), size(EUR, 1))
@test issymmetric(Φrbs)
@test eigmin(Φrbs) > -1e-8
m, n, K = size(EUR, 1), size(EUR, 2), 3
P = rand(n, K)
Q = rand(m, K)
Q ./= sum(Q, dims=2)
Φreap = grm_admixture(EUR, transpose(P), transpose(Q))
@test size(Φreap) == (size(EUR, 1), size(EUR, 1))
@test issymmetric(Φreap)
@test eigmin(Φrbs) > -1e-8
end

@testset "filter" begin
rowmask, colmask =  SnpArrays.filter(mouse, min_success_rate_per_row=0.99,
    min_success_rate_per_col=0.99)
SnpArrays.filter(SnpArrays.datadir("mouse"), rowmask, colmask; des="tmp")
tmpbf = SnpArray("tmp.bed")
@test size(tmpbf) == (1912, 9997)
@test all(missingrate(tmpbf, 1) .≤ 0.01)
@test all(missingrate(tmpbf, 2) .≤ 0.01)
Sys.iswindows() || rm("tmp.bed", force=true)
rm("tmp.bim", force=true)
rm("tmp.fam", force=true)
rowmask, colmask =  SnpArrays.filter(mouse, min_success_rate_per_row=0.99,
    min_success_rate_per_col=0.99, min_maf=0.01, min_hwe_pval=1e-8)
@test (count(rowmask), count(colmask)) == (1911, 9510)
compress_plink(SnpArrays.datadir("mouse"), "gz")
SnpArrays.filter(SnpArrays.datadir("mouse.bed.gz"), SnpArrays.datadir("mouse.bim.gz"),
SnpArrays.datadir("mouse.fam.gz"), 1:5, 1:3; des="tmp")
tmpbf = SnpArray("tmp.bed")
@test size(tmpbf) == (5, 3)
@test isfile("tmp.bed")
@test isfile("tmp.fam")
@test isfile("tmp.bim")
rm(SnpArrays.datadir("mouse.bed.gz"), force=true)
rm(SnpArrays.datadir("mouse.fam.gz"), force=true)
rm(SnpArrays.datadir("mouse.bim.gz"), force=true)
Sys.iswindows() || rm("tmp.bed", force=true)
rm("tmp.bim", force=true)
rm("tmp.fam", force=true)
end

@testset "SnpBitMatrix-vector multiplication" begin
reltol = 5e-4
for t in [Float32, Float64]
    v1 = randn(t, size(EUR, 1))
    v2 = randn(t, size(EUR, 2))
    for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
        @test norm(SnpBitMatrix{t}(EUR, model=model) * v2 -
            convert(Matrix{t}, EUR, model=model) * v2) /
            norm(convert(Matrix{t}, EUR, model=model) * v2) < reltol
        @test norm(SnpBitMatrix{t}(EUR, center=true, model=model) * v2 -
            convert(Matrix{t}, EUR, center=true, model=model) * v2) /
            norm(convert(Matrix{t}, EUR, center=true, model=model) * v2) < reltol
        @test norm(SnpBitMatrix{t}(EUR, scale=true, model=model) * v2 -
            convert(Matrix{t}, EUR, scale=true, model=model) * v2) /
            norm(convert(Matrix{t}, EUR, scale=true, model=model) * v2) < reltol
        @test norm(SnpBitMatrix{t}(EUR, center=true, scale=true, model=model) * v2 -
            convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) /
            norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) < reltol
        @test norm(transpose(SnpBitMatrix{t}(EUR, model=model)) * v1 -
            transpose(convert(Matrix{t}, EUR, model=model)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, model=model)) * v1) < reltol
        @test norm(transpose(SnpBitMatrix{t}(EUR, center=true, model=model)) * v1 -
            transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) < reltol
        @test norm(transpose(SnpBitMatrix{t}(EUR, scale=true, model=model)) * v1 -
            transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) < reltol
        @test norm(transpose(SnpBitMatrix{t}(EUR, center=true, scale=true, model=model)) * v1 -
            transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) < reltol
    end
end
end

@testset "copyto bitmatrix" begin
for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
    mousebm = SnpBitMatrix{t}(mouse, model=ADDITIVE_MODEL, center=false, scale=false)
    v = copyto!(zeros(1), @view(mousebm[702]))
    @test v[1] == 0
    mousebm = SnpBitMatrix{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true)
    v = copyto!(zeros(1), @view(mousebm[702]))
    @test v[1] ≈ -1 * mousebm.μ[1] * mousebm.σinv[1]
    EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
    EURtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
    EURtest = copyto!(zeros(t, size(EUR)), EURbm)
    @test all(EURtest .≈ EURtrue)
    EURtest = copyto!(zeros(10, 5), @view(EURbm[11:20, 1:2:10]))
    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
    EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
    EURtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
    EURtest = copyto!(zeros(t, size(EUR)), EURbm)
    @test all(EURtest .≈ EURtrue)
    EURtest = copyto!(zeros(10, 5), @view(EURbm[11:20, 1:2:10]))
    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
end
end

@testset "convert (SnpBitMatrix)" begin
for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
    EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
    xtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
    xtest = convert(Matrix{t}, EURbm)
    @test all(xtrue .≈ xtest)
    EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
    xtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
    xtest = convert(Matrix{t}, EURbm)
    @test all(xtrue .≈ xtest)
    EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=true, scale=true)
    xtest = convert(Matrix{t}, @view(EURbm[:, 2:2:10]))
    @test all(xtrue .≈ xtest)
    EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=false, scale=false)
    xtest = convert(Matrix{t}, @view(EURbm[:, 2:2:10]))
    @test all(xtrue .≈ xtest)
end
end

@testset "copyto SnpLinAlg" begin
for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
    # imputing missing data
    mousela = SnpLinAlg{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true, impute=true)
    v = copyto!(zeros(1), @view(mousela[702])) # missing entry
    @test isapprox(v[1], 1.113003134727478, atol=1e-6)
    # not imputing missing data
    mousela = SnpLinAlg{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true, impute=false)
    v = copyto!(zeros(1), @view(mousela[702])) # missing entry
    @test isnan(v[1])
    # no missing data
    EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
    EURtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
    EURtest = copyto!(zeros(t, size(EUR)), EURla)
    @test all(EURtest .≈ EURtrue)
    EURtest = copyto!(zeros(10, 5), @view(EURla[11:20, 1:2:10]))
    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
    EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
    EURtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
    EURtest = copyto!(zeros(t, size(EUR)), EURla)
    @test all(EURtest .≈ EURtrue)
    EURtest = copyto!(zeros(10, 5), @view(EURla[11:20, 1:2:10]))
    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
end
end

@testset "convert (SnpLinAlg)" begin
for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
    EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
    xtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
    xtest = convert(Matrix{t}, EURla)
    @test all(xtrue .≈ xtest)
    EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
    xtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
    xtest = convert(Matrix{t}, EURla)
    @test all(xtrue .≈ xtest)
    EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=true, scale=true)
    xtest = convert(Matrix{t}, @view(EURla[:, 2:2:10]))
    @test all(xtrue .≈ xtest)
    EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=false, scale=false)
    xtest = convert(Matrix{t}, @view(EURla[:, 2:2:10]))
    @test all(xtrue .≈ xtest)
end
end

@testset "SnpLinAlg-vector multiplication (zeroimpute)" begin
reltol = 5e-4
for t in [Float32, Float64]
    v1 = randn(t, size(EUR, 1))
    v2 = randn(t, size(EUR, 2))
    for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
        @test norm(SnpLinAlg{t}(EUR, model=model, impute=false) * v2 -
            convert(Matrix{t}, EUR, model=model, impute=false) * v2) /
            norm(convert(Matrix{t}, EUR, model=model, impute=false) * v2) < reltol
        @test norm(SnpLinAlg{t}(EUR, center=true, model=model, impute=false) * v2 -
            convert(Matrix{t}, EUR, center=true, model=model, impute=false) * v2) /
            norm(convert(Matrix{t}, EUR, center=true, model=model, impute=false) * v2) < reltol
        @test norm(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false) * v2 -
            convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * v2) /
            norm(convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * v2) < reltol
        @test norm(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * v2 -
            convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * v2) /
            norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * v2) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, model=model, impute=false)) * v1 -
            transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * v1) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=false)) * v1 -
            transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * v1) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false)) * v1 -
            transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * v1) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false)) * v1 -
            transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * v1) < reltol
    end
end
end

@testset "SnpLinAlg-matrix multiplication zeroimpute" begin
    for t in [Float32, Float64] # Y = A*X
        X1 = randn(t, size(EUR, 1), 3)
        X2 = randn(t, size(EUR, 2), 3)
        for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
            @test all(isapprox(SnpLinAlg{t}(EUR, model=model, impute=false) * X2,
                convert(Matrix{t}, EUR, model=model, impute=false) * X2))
            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, model=model, impute=false) * X2,
                convert(Matrix{t}, EUR, center=true, model=model, impute=false) * X2))
            @test all(isapprox(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false) * X2,
                convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * X2))
            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * X2,
                convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * X2))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, model=model, impute=false)) * X1,
                transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * X1))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=false)) * X1,
                transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * X1))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false)) * X1,
                transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * X1))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false)) * X1,
                transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * X1))
        end
    end
end

@testset "SnpLinAlg-vector multiplication (Miter > 0)" begin
    EUR11 = [EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR]
    EUR11la = SnpLinAlg{Float64}(EUR11, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
    v = rand(size(EUR11la, 2))
    vtest = EUR11la * v
    vtrue = convert(Matrix{Float64}, EUR11, model=ADDITIVE_MODEL, impute=true, center=true, scale=true) * v
    @test norm(vtest - vtrue) < 5e-4
end

if get(ENV,"JULIA_SNPARRAYS_TEST_CUDA","") == "true"
    using CUDA, Adapt
    @testset "lin. alg. cuda" begin
    reltol = 5e-4
    for t in [Float32, Float64]
        v1 = randn(t, size(EUR, 1))
        v2 = randn(t, size(EUR, 2))
        v1_d = adapt(CuVector{t}, v1)
        v2_d = adapt(CuVector{t}, v2)
        for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
            @test norm(collect(CuSnpArray{t}(EUR, model=model) * v2_d) -
                convert(Matrix{t}, EUR, model=model) * v2) /
                norm(convert(Matrix{t}, EUR, model=model) * v2) < reltol
            @test norm(collect(CuSnpArray{t}(EUR, center=true, model=model) * v2_d) -
                convert(Matrix{t}, EUR, center=true, model=model) * v2) /
                norm(convert(Matrix{t}, EUR, center=true, model=model) * v2) < reltol
            @test norm(collect(CuSnpArray{t}(EUR, scale=true, model=model) * v2_d) -
                convert(Matrix{t}, EUR, scale=true, model=model) * v2) /
                norm(convert(Matrix{t}, EUR, scale=true, model=model) * v2) < reltol
            @test norm(collect(CuSnpArray{t}(EUR, center=true, scale=true, model=model) * v2_d) -
                convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) /
                norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) < reltol
            @test norm(collect(transpose(CuSnpArray{t}(EUR, model=model)) * v1_d) -
                transpose(convert(Matrix{t}, EUR, model=model)) * v1) /
                norm(transpose(convert(Matrix{t}, EUR, model=model)) * v1) < reltol
            @test norm(collect(transpose(CuSnpArray{t}(EUR, center=true, model=model)) * v1_d) -
                transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) /
                norm(transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) < reltol
            @test norm(collect(transpose(CuSnpArray{t}(EUR, scale=true, model=model)) * v1_d) -
                transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) /
                norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) < reltol
            @test norm(collect(transpose(CuSnpArray{t}(EUR, center=true, scale=true, model=model)) * v1_d) -
                transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) /
                norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) < reltol
        end
    end
    end
end

@testset "SnpLinAlg-vector multiplication meanimpute" begin
reltol = 5e-4
for t in [Float32, Float64]
    v1 = randn(t, size(EUR, 1))
    v2 = randn(t, size(EUR, 2))
    for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
        @test norm(SnpLinAlg{t}(EUR, model=model, impute=true) * v2 -
            convert(Matrix{t}, EUR, model=model, impute=true) * v2) /
            norm(convert(Matrix{t}, EUR, model=model, impute=true) * v2) < reltol
        @test norm(SnpLinAlg{t}(EUR, center=true, model=model, impute=true) * v2 -
            convert(Matrix{t}, EUR, center=true, model=model, impute=true) * v2) /
            norm(convert(Matrix{t}, EUR, center=true, model=model, impute=true) * v2) < reltol
        @test norm(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true) * v2 -
            convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * v2) /
            norm(convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * v2) < reltol
        @test norm(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true) * v2 -
            convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * v2) /
            norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * v2) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, model=model, impute=true)) * v1 -
            transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * v1) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=true)) * v1 -
            transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * v1) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true)) * v1 -
            transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * v1) < reltol
        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true)) * v1 -
            transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * v1) /
            norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * v1) < reltol
    end
end
end

@testset "SnpLinAlg-matrix multiplication meanimpute" begin
    for t in [Float32, Float64] # Y = A*X
        X1 = randn(t, size(EUR, 1), 3)
        X2 = randn(t, size(EUR, 2), 3)
        for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
            @test all(isapprox(SnpLinAlg{t}(EUR, model=model, impute=true) * X2,
                convert(Matrix{t}, EUR, model=model, impute=true) * X2))
            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, model=model, impute=true) * X2,
                convert(Matrix{t}, EUR, center=true, model=model, impute=true) * X2))
            @test all(isapprox(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true) * X2,
                convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * X2))
            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * X2,
                convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * X2))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, model=model, impute=true)) * X1,
                transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * X1))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=true)) * X1,
                transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * X1))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true)) * X1,
                transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * X1))
            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true)) * X1,
                transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * X1))
        end
    end
end

@testset "subarrays" begin
@test all(@view(EUR[1:2:10, 1:2:10]) .==
[[0x03 0x03 0x02 0x02 0x03];
[0x03 0x03 0x03 0x03 0x03];
[0x03 0x03 0x03 0x03 0x03];
[0x02 0x03 0x02 0x00 0x02];
[0x03 0x03 0x02 0x02 0x03]
])
@test all(convert(Matrix{Float64}, @view(EUR[1:2:10, 1:2:10])) .≈
[[2.0 2.0 1.0 1.0 2.0];
[2.0 2.0 2.0 2.0 2.0];
[2.0 2.0 2.0 2.0 2.0];
[1.0 2.0 1.0 0.0 1.0];
[2.0 2.0 1.0 1.0 2.0]
])

EURsub = @view EUR[1:2:100, 1:2:100]
EURsubbm = SnpBitMatrix{Float64}(EURsub, model=ADDITIVE_MODEL, center=true, scale=true) # BitMatrix from the SubArray
EURsubfm = convert(Matrix{Float64}, EURsub, model=ADDITIVE_MODEL, center=true, scale=true) # Float64 Matrix from the SubArray

v1 = randn(size(EURsub, 1))
v2 = randn(size(EURsub, 2))

@test isapprox(EURsubbm * v2, EURsubfm * v2)
@test isapprox(EURsubbm' * v1, EURsubfm' * v1)
end

@testset "split-merge-readwrite" begin
EUR_data = SnpData(SnpArrays.datadir("EUR_subset"))
# a small subset of EUR_data
SnpArrays.filter(EUR_data, collect(1:10), collect(1:20000); des="EUR_subset.small")
EUR_small = SnpData("EUR_subset.small")

# filter
chr17 = SnpArrays.filter(SnpArrays.datadir("EUR_subset"); des="tmp.filter.chr17", f_snp = x -> String(x[:chromosome]) == "17")
@test chr17.snps == 11041
@test chr17.people == 379
male = SnpArrays.filter(EUR_data; des="tmp.filter.male", f_person = x -> String(x[:sex]) == "1")
@test male.snps == 54051
@test male.people == 178
chr17_male = SnpArrays.filter(EUR_data; des="tmp.filter.chr17.male", f_snp = x -> String(x[:chromosome]) == "17", f_person = x -> String(x[:sex]) == "1")
@test chr17_male.snps == 11041
@test chr17_male.people == 178
@test size(chr17_male.snparray) == (178, 11041)
# cleanup
for ft in ["bim", "fam", "bed"]
    ft == "bed" && Sys.iswindows() && continue
    rm("tmp.filter.chr17." * ft, force=true)
    rm("tmp.filter.male." * ft, force=true)
    rm("tmp.filter.chr17.male." * ft, force=true)
end

# split
splitted = SnpArrays.split_plink("EUR_subset.small"; prefix="tmp.split.chr.")
piece = splitted["17"]
@test piece.people == 10
@test piece.snps == 11041
@test size(piece.person_info) == (10, 6)
@test size(piece.snp_info) == (11041, 6)
@test size(piece.snparray) == (10, 11041)
@test size(piece.snparray.columncounts) == (4, 11041)
@test piece.snparray.m == 10

splitted_bysex = SnpArrays.split_plink(EUR_data, :sex; prefix="tmp.split.sex.")
@test splitted_bysex["1"].people == 178
@test splitted_bysex["2"].people == 201

# merge
#@time merged = SnpArrays.merge_plink("tmp.merged", splitted) # write_plink is included here
#@test EUR_data.people == merged.people
#@test EUR_data.snps == merged.snps
#@test EUR_data.person_info == merged.person_info
#@test EUR_data.snp_info == merged.snp_info # note: the ordering of merged data might be different on other dataset b/c sorted order of chromosomes

#output = SnpData("tmp.merged")
#@test EUR_data.people == output.people
#@test EUR_data.snps == output.snps
#@test EUR_data.person_info == output.person_info
#@test EUR_data.snp_info == output.snp_info

# cleanup
#rm("tmp.merged.bim", force=true)
#rm("tmp.merged.fam", force=true)
#Sys.iswindows() || rm("tmp.merged.bed", force=true)

# merge from splitted files
@time merged_from_splitted_files = merge_plink("tmp.split.chr"; des = "tmp2.merged")
@test EUR_small.people == merged_from_splitted_files.people
@test EUR_small.snps == merged_from_splitted_files.snps
@test EUR_small.person_info == merged_from_splitted_files.person_info
@test EUR_small.snp_info == merged_from_splitted_files.snp_info

# cleanup
rm("tmp2.merged.bim", force=true)
rm("tmp2.merged.fam", force=true)
Sys.iswindows() || rm("tmp2.merged.bed", force=true)

rm("EUR_subset.small.bim", force=true)
rm("EUR_subset.small.fam", force=true)
Sys.iswindows() || rm("EUR_subset.small.bed", force=true)

for k in keys(splitted)
    for ft in ["bim", "fam", "bed"]
        ft == "bed" && Sys.iswindows() && continue
        rm("tmp.split.chr.$(k)." * ft, force=true)
    end
end
for k in keys(splitted_bysex)
    for ft in ["bim", "fam", "bed"]
        ft == "bed" && Sys.iswindows() && continue
        rm("tmp.split.sex.$(k)." * ft, force=true)
    end
end
end

@testset "(de)compress" begin
for format in SnpArrays.ALLOWED_FORMAT
    # compress mouse Plink files
    compress_plink(SnpArrays.datadir("mouse"), format)
    @test isfile(SnpArrays.datadir("mouse.bed." * format))
    @test isfile(SnpArrays.datadir("mouse.fam." * format))
    @test isfile(SnpArrays.datadir("mouse.bim." * format))
    # read in compressed Plink files
    mouse_zip = SnpArray(SnpArrays.datadir("mouse.bed." * format))
    fill!(mouse.rowcounts, 0)
    fill!(mouse.columncounts, 0)
    @time begin
        SnpArrays._counts(mouse, 1)
        SnpArrays._counts(mouse, 2)
    end
    fill!(mouse_zip.rowcounts, 0)
    fill!(mouse_zip.columncounts, 0)
    @time begin
        SnpArrays._counts(mouse_zip, 1)
        SnpArrays._counts(mouse_zip, 2)
    end
    @test mouse.m == mouse_zip.m
    @test norm(mouse.rowcounts - mouse_zip.rowcounts) < 1e-8
    @test norm(mouse.columncounts - mouse_zip.columncounts) < 1e-8
    @test norm(mouse.data - mouse_zip.data) < 1e-8
    # Decompress Plink files
    decompress_plink(SnpArrays.datadir("mouse"), format, SnpArrays.datadir("mouse2"))
    @test stat(SnpArrays.datadir("mouse") * ".bed").size ==
    stat(SnpArrays.datadir("mouse2") * ".bed").size
    @test stat(SnpArrays.datadir("mouse") * ".fam").size ==
    stat(SnpArrays.datadir("mouse2") * ".fam").size
    @test stat(SnpArrays.datadir("mouse") * ".bim").size ==
    stat(SnpArrays.datadir("mouse2") * ".bim").size
    # clean up
    rm(SnpArrays.datadir("mouse.bed." * format), force=true)
    rm(SnpArrays.datadir("mouse.fam." * format), force=true)
    rm(SnpArrays.datadir("mouse.bim." * format), force=true)
    rm(SnpArrays.datadir("mouse2.bed"), force=true)
    rm(SnpArrays.datadir("mouse2.fam"), force=true)
    rm(SnpArrays.datadir("mouse2.bim"), force=true)
end
@test_throws ArgumentError SnpArray(SnpArrays.datadir("mouse.bed.zip"))
end

@testset "indexin" begin
    b = ['1', '2', '3', '4']
    for a in [['4', '3', '1', '2'], ['3', '2', '6'], ['b', 'c', 'd', 'c']]
        aind, bmask = SnpArrays.indexin_general(a, b)
        @test all(a[aind] .== b[bmask])
    end
end

@testset "concat" begin
s = SnpArrays.filter(SnpArrays.datadir("mouse"), 1:2, 1:3; des="mouse.tmp.filtered")
sd = SnpData("mouse.tmp.filtered")
sd_hcat = hcat(sd, sd, sd; des="mouse.hcat")
sd_vcat = vcat(sd, sd, sd; des="mouse.vcat")
sd_hvcat = hvcat((2, 2), sd, sd, sd, sd; des="mouse.hvcat")

ref = [[0x02 0x02 0x02];
[0x02 0x02 0x03]]

@test all(sd_hcat.snparray .== [ref ref ref])
@test all(sd_vcat.snparray .== [ref; ref; ref])
@test all(sd_hvcat.snparray .== [ref ref; ref ref])

@test sd_hcat.people == 2
@test sd_hcat.snps == 9

@test sd_vcat.people == 6
@test sd_vcat.snps == 3

@test sd_hvcat.people == 4
@test sd_hvcat.snps == 6

for ft in ["bim", "fam", "bed"]
    ft == "bed" && Sys.iswindows() && continue
    rm("mouse.tmp.filtered." * ft, force=true)
    rm("mouse.hcat." * ft, force=true)
    rm("mouse.vcat." * ft, force=true)
    rm("mouse.hvcat." * ft, force=true)
end
end

@testset "reorder" begin
mouse_prefix = SnpArrays.datadir("mouse")
run(`cp $(mouse_prefix * ".bed") mouse_testreorder.bed`)
run(`cp $(mouse_prefix * ".bim") mouse_testreorder.bim`)
run(`cp $(mouse_prefix * ".fam") mouse_testreorder.fam`)

mouse_data = SnpData(mouse_prefix)
mouse_toreorder = SnpData("mouse_testreorder", "r+")
m, n = size(mouse_toreorder.snparray)
using Random
ind = randperm(m);
SnpArrays.reorder!(mouse_toreorder, ind)

# reread the file
mouse_toreorder_read = SnpData("mouse_testreorder"; famnm="mouse_testreorder.reordered.fam")
@test all(mouse_data.snparray[ind, :] .== mouse_toreorder_read.snparray[:, :])
@test all(mouse_data.person_info[ind, 1] .== mouse_toreorder_read.person_info[:, 1])

# cleanup
rm("mouse_testreorder.bed", force=true)
rm("mouse_testreorder.bim", force=true)
rm("mouse_testreorder.fam", force=true)
rm("mouse_testreorder.reordered.fam", force=true)
end

@testset "vcf2plink" begin
# Download an example VCF file
isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
    joinpath(pwd(), "test.08Jun17.d8b.vcf.gz"));
s = vcf2plink("test.08Jun17.d8b.vcf.gz", "test.08Jun17.d8b")
@test size(s, 1) == 191
@test size(s, 2) == 1354
@test s[1, 1] == 0x00
@test s[1, 5] == 0x02
@test s[38, 5] == 0x03
rm("test.08Jun17.d8b.vcf.gz", force=true)
rm("test.08Jun17.d8b.bed", force=true)
rm("test.08Jun17.d8b.bim", force=true)
rm("test.08Jun17.d8b.fam", force=true)
end

@testset "kinship_pruning" begin
g = grm(mouse)
@test count(kinship_pruning(g; method=:gcta)) == 68
@test count(kinship_pruning(g; method=:top_down)) == 125
@test count(kinship_pruning(g; method=:bottom_up)) == 132
@test count(kinship_pruning(g; method=:plink)) == 126
end


text:
# The following was inspired by AtariAlgos.jl
# Modifications by Dennis Wilson @d9w

using Colors
using ImageCore
using ImageTransformations
export
    Game,
    close!,
    draw,
    get_inputs,
    get_rgb

struct Game
    ale::ALEPtr
    width::Int
    height::Int
    actions::Array{Int32}
end

function Game(romfile::String, seed::Int)
    ale = ALE_new()
    setInt(ale, "random_seed", Cint(seed))
    loadROM(ale, romfile)
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    actions = getMinimalActionSet(ale)
    Game(ale, w, h, actions)
end

function close!(game::Game)
    ALE_del(game.ale)
end

function draw(game::Game)
    rawscreen = getScreenRGB(game.ale)
    colorview(RGB, Float64.(reshape(rawscreen/256.,
                                    (3, game.width, game.height))))';
end

function get_inputs(game::Game)
    screen = getScreen(game.ale)/(0xff*1.0)
    screen = reshape(screen, (game.width, game.height))'
    # imresize(screen, (42, 32))/256.
    screen
end

function get_rgb(game::Game)
    # TODO: speed this up
    # rawscreen = Array{Cuchar}(undef, game.width * game.height * 3)
    rawscreen = getScreenRGB(game.ale)
    rgb = Float64.(reshape(rawscreen/256., (3, game.width, game.height)));
    [Array{Float64}(rgb[i,:,:]) for i in 1:3]
end

function get_ram(game::Game)
    getRAM(game.ale) ./ typemax(UInt8)
end


text:
# map: apply function to elts of iterator
#
# PROBLEM: Type inference in NamedTuples does not work for Union{Missing, T} types

export map

struct Map{I,F,M} <: AbstractGorIter{I}
    rows::I
    func::F
end

"""
    map(rows, func)
    rows |> map(func)

Apply function `func` to elements of genome ordered stream `rows`.
The function `func` should return a `NamedTuple`.

NOTE: Julia cannot infer the type of NamedTuples with `Union{Missing,T}`.
This means that if the input stream `rows` has columns of type `Union{Missing,T}`,
the pipeline probably fails.
"""
map(rows, func) = Map{typeof(rows), typeof(func),
                         Base.return_types(func, (eltype(rows),))[1]}(rows, func)
map(func) = rows -> map(rows, func)

#transform = map

#mutate(rows, func) = gormap(rows, row -> merge(row, func(row)))
#mutate(func) = rows -> mutate(rows, func)

Base.IteratorEltype(::Type{Map{I,F,M}}) where {I,F,M} = Base.HasEltype()
Base.eltype(::Type{Map{I,F,M}}) where {I,F,M} = M
Base.IteratorSize(::Type{Map{I,F,M}}) where {I,F,M} = Base.IteratorSize(I)
Base.length(m::Map{I,F,M}) where {I,F,M} = Base.length(m.rows)
#Base.size(m::Map{I,F,M}) where {I,F,M} = Base.size(m.rows)

Tables.istable(m::Map{I,F,M}) where {I,F,M} = true
Tables.rowaccess(m::Map{I,F,M}) where {I,F,M} = true
Tables.schema(m::Map{I,F,M}) where {I,F,M} = Tables.Schema(M)

function Base.iterate(m::Map{I,F,M}) where {I,F,M}
    elt_s = iterate(m.rows)
    state = elt_s
    iterate(m, state)
end

function Base.iterate(m::Map{I,F,M}, state) where {I,F,M}
    elt_s = state
    elt_s == nothing && return nothing
    M(m.func(elt_s[1])), iterate(m.rows, elt_s[2])
end


text:
function send(url::AbstractString, to::Vector, from::AbstractString, body::IO,
              options::SendOptions = SendOptions())
  if options.blocking
    rd::ReadData = ReadData()

    rd.typ = :io
    rd.src = body
    seekend(body)
    rd.sz = position(body)
    seekstart(body)

    _do_send(url, to, from, options, rd)
  else
    remotecall(send, myid(), url, to, from, body, set_opt_blocking(options))
  end
end

function _do_send(url::AbstractString, to::Vector, from::AbstractString,
                  options::SendOptions, rd::ReadData)
  ctxt = false
  slist::Ptr{Cvoid} = C_NULL
  try
    ctxt = setup_easy_handle(url, options)
    ctxt.rd = rd

    @ce_curl curl_easy_setopt ctxt.curl CURLOPT_READFUNCTION c_curl_read_cb
    @ce_curl curl_easy_setopt ctxt.curl CURLOPT_READDATA ctxt

    for tos in to
      slist = curl_slist_append(slist, tos)
    end
    @ce_curl curl_easy_setopt ctxt.curl CURLOPT_MAIL_RCPT slist
    @ce_curl curl_easy_setopt ctxt.curl CURLOPT_MAIL_FROM from

    @ce_curl curl_easy_perform ctxt.curl
    process_response(ctxt)
    return ctxt.resp
  finally
    if (slist != C_NULL)
      curl_slist_free_all(slist)
    end
    cleanup_easy_context(ctxt)
  end
end

function exec_as_multi(ctxt)
  curl = ctxt.curl
  curlm = curl_multi_init()

  if curlm == C_NULL
    error("Unable to initialize curl_multi_init()")
  end

  try
    if isa(ctxt.options.callback, Function)
      ctxt.options.callback(curl)
    end

    @ce_curlm curl_multi_add_handle curl

    n_active = Array{Int}(1)
    n_active[1] = 1

    no_to = 30 * 24 * 3600.0
    request_timeout = 0.001 + (ctxt.options.request_timeout == 0.0 ? no_to : ctxt.options.request_timeout)

    started_at = time()
    time_left = request_timeout

    # poll_fd is unreliable when multiple parallel fds are active, hence using curl_multi_perform
    cmc = curl_multi_perform(curlm, n_active);
    while (n_active[1] > 0) &&  (time_left > 0)
      nb1 = ctxt.bytes_recd
      cmc = curl_multi_perform(curlm, n_active);
      if(cmc != CURLM_OK) error("curl_multi_perform() failed: " * string(curl_multi_strerror(cmc))) end

      nb2 = ctxt.bytes_recd

      if (nb2 > nb1)
        yield() # Just yield to other tasks
      else
        sleep(0.005) # Just to prevent unnecessary CPU spinning
      end

      time_left = request_timeout - (time() - started_at)
    end
    # END OF curl_multi_perform

    if (n_active[1] == 0)
      msgs_in_queue = Array{Cint}(1)
      p_msg::Ptr{CURLMsg2} = curl_multi_info_read(curlm, msgs_in_queue)

      while (p_msg != C_NULL)
        # println("Messages left in Q : " * string(msgs_in_queue[1]))
        msg = unsafe_load(p_msg)

        if (msg.msg == CURLMSG_DONE)
          ec = convert(Int, msg.data)
          if (ec != CURLE_OK)
            # println("Result of transfer: " * string(msg.data))
            throw("Error executing request : " * string(curl_easy_strerror(ec)))
          else
            process_response(ctxt)
          end
        end

        p_msg = curl_multi_info_read(curlm, msgs_in_queue)
      end
    else
      error("request timed out")
    end

  finally
    curl_multi_remove_handle(curlm, curl)
    curl_multi_cleanup(curlm)
  end

  ctxt.resp
end


text:
using StingerWrapper
using BenchmarkTools
using GraphBenchmarks
import GraphBenchmarks: construct, preparebench, picksources

import StingerWrapper: bfs

export StingerGraph, SerialBFS, LevelSynchronousBFS, StingerBenchmarks

immutable StingerGraph <: GraphType
    s::Stinger
    nv::Int64
end

immutable StingerBenchmarks <: GraphBenchmarkSpec end

immutable SerialBFS <: AbstractBFS end
immutable LevelSynchronousBFS <: AbstractBFS end

function construct(t::Type{StingerGraph}, edges::Array{Int64, 2})
    s = Stinger()
    for i=1:size(edges,2)
        insert_edge!(s, 0, edges[1, i]-1, edges[2, i]-1, 0, 0)
    end
    StingerGraph(s, maximum(edges)) #max number of active vertices
end

function picksources(::StingerBenchmarks, nv::Int64)
    sources = collect(1:1000)
end

function preparebench(benchmark::GraphBenchmarkSpec, alg::SerialBFS, g::StingerGraph)
    sources = picksources(benchmark, g.nv) - 1
    bfsbench = @benchmarkable begin
        for src in $sources
            StingerWrapper.bfs($(g.s), src, $(g.nv))
        end
    end seconds=6000 samples=10
    bfsbench
end

function preparebench(benchmark::GraphBenchmarkSpec, alg::LevelSynchronousBFS, g::StingerGraph)
    sources = picksources(benchmark, g.nv) - 1
    bfsbench = @benchmarkable begin
        for src in $sources
            StingerWrapper.bfs(StingerWrapper.LevelSynchronous(), $(g.s), src, $(g.nv))
        end
    end seconds=6000 samples=10
    bfsbench
end

function bfsbenchsuite(T::Type{StingerGraph}, bench::GraphBenchmarkSpec, scale::Int, edgefactor::Int)
    srand(0)
    gen = Kronecker(scale, edgefactor)
    suite = BenchmarkGroup(["bfs", "stinger"])
    suite["serial"] = BenchmarkGroup(["serialbfs"])
    suite["serial"]["stingerserialbfs"] = preparebenchmark(bench, T, SerialBFS(), gen)
    suite["parallel"] = BenchmarkGroup(["parallelbfs"])
    suite["parallel"]["levelsyncbfs"] = preparebenchmark(bench, T, LevelSynchronousBFS(), gen)
    suite
end

function comparebfs(T::Type{StingerGraph}, bench::GraphBenchmarkSpec, scale::Int, edgefactor::Int)
    suite = bfsbenchsuite(T, bench, scale, edgefactor)
    run(suite, verbose=true)
end


text:
using JLD2

export obs_env
export AbstractLattice, SquareLattice
abstract type AbstractLattice end
struct SquareLattice <: AbstractLattice end

export VUMPSRuntime, SquareVUMPSRuntime
# NOTE: should be renamed to more explicit names
"""
    VUMPSRuntime{LT}

a struct to hold the tensors during the `vumps` algorithm, each is a `Ni` x `Nj` Matrix, containing
- `d × d × d × d'` `M[i,j]` tensor
- `D × d' × D` `AL[i,j]` tensor
- `D × D`     `C[i,j]` tensor
- `D × d' × D` `AR[i,j]` tensor
- `D × d' × D` `FL[i,j]` tensor
- `D × d' × D` `FR[i,j]` tensor
and `LT` is a AbstractLattice to define the lattice type.
"""
struct VUMPSRuntime{LT,T,N,AT <: AbstractArray{<:AbstractArray,2},CT,ET1,ET2}
    M::AT
    AL::ET1
    C::CT
    AR::ET1
    FL::ET2
    FR::ET2
    function VUMPSRuntime{LT}(M::AT, AL::ET1, C::CT, AR::ET1, FL::ET2, FR::ET2) where {LT <: AbstractLattice,AT <: AbstractArray{<:AbstractArray,2}, CT <: AbstractArray{<:AbstractArray,2}, ET1 <: AbstractArray{<:AbstractArray,2}, ET2 <: AbstractArray{<:AbstractArray,2}}
        T, N = eltype(M[1,1]), ndims(M[1,1])
        new{LT,T,N,AT,CT,ET1,ET2}(M, AL, C, AR, FL, FR)
    end
end

const SquareVUMPSRuntime{T,AT} = VUMPSRuntime{SquareLattice,T,4,AT}
function SquareVUMPSRuntime(M::AT, AL, C, AR, FL, FR) where {AT <: AbstractArray{<:AbstractArray,2}}
    ndims(M[1,1]) == 4 || throw(DimensionMismatch("M dimensions error, should be `4`, got $(ndims(M[1,1]))."))
    VUMPSRuntime{SquareLattice}(M, AL, C, AR, FL, FR)
end

@doc raw"
    SquareVUMPSRuntime(M::AbstractArray{T,4}, env::Val, χ::Int)

create a `SquareVUMPSRuntime` with M-tensor `M`. The `NixNj` `AL,C,AR,FL,FR`
tensors are initialized according to `env`. If `env = Val(:random)`,
the `A[i,j]` is initialized as a random `D×d×D` tensor,and `AL[i,j],C[i,j],AR[i,j]` are the corresponding
canonical form. `FL,FR` is the left and right environment.

# example

```jldoctest; setup = :(using VUMPS)
julia> Ni, Nj = 2, 2;

julia> M = Array{Array{ComplexF64,3},2}(undef, Ni, Nj);

julia> for j = 1:Nj, i = 1:Ni
           M[i,j] = rand(2,2,2,2)
       end

julia> rt = SquareVUMPSRuntime(M, Val(:random), 4);

julia> size(rt.AL) == (2,2)
true

julia> size(rt.AL[1,1]) == (4,2,4)
true
```
"
function SquareVUMPSRuntime(M::AbstractArray{<:AbstractArray,2}, env, χ::Int; verbose=false)
    return SquareVUMPSRuntime(M, _initializect_square(M, env, χ; verbose=verbose)...)
end

function _initializect_square(M::AbstractArray{<:AbstractArray,2}, env::Val{:random}, χ::Int; verbose=false)
    A = initialA(M, χ)
    AL, L = leftorth(A)
    R, AR = rightorth(AL)
    _, FL = leftenv(AL, AL, M)
    _, FR = rightenv(AR, AR, M)
    C = LRtoC(L,R)
    Ni, Nj = size(M)
    verbose && print("random initial $(Ni)×$(Nj) vumps_χ$(χ) environment-> ")
    AL, C, AR, FL, FR
end

function _initializect_square(M::AbstractArray{<:AbstractArray,2}, chkp_file::String, χ::Int; verbose=false)
    env = load(chkp_file)["env"]
    Ni, Nj = size(M)
    atype = _arraytype(M[1,1])
    verbose && print("vumps $(Ni)×$(Nj) environment load from $(chkp_file) -> ")
    AL, C, AR, FL, FR = env.AL, env.C, env.AR, env.FL, env.FR
    Zygote.@ignore begin
        AL, C, AR, FL, FR = Array{atype{ComplexF64,3},2}(env.AL), Array{atype{ComplexF64,2},2}(env.C), Array{atype{ComplexF64,3},2}(env.AR), Array{atype{ComplexF64,3},2}(env.FL), Array{atype{ComplexF64,3},2}(env.FR)
    end
    AL, C, AR, FL, FR
end

function vumps(rt::VUMPSRuntime; tol::Real=1e-10, maxiter::Int=10, miniter::Int=1, verbose=false, show_every = Inf)
    # initialize
    olderror = Inf
    vumps_counting = show_every_count(show_every)

    stopfun = StopFunction(olderror, -1, tol, maxiter, miniter)
    rt, err = fixedpoint(res -> vumpstep(res...;show_counting=vumps_counting), (rt, olderror), stopfun)
    verbose && println("vumps done@step: $(stopfun.counter), error=$(err)")
    return rt
end

function show_every_count(n::Number)
    i = 0
    counting() = (i += 1; mod(i,n)==0 ? i : 0)
    return counting
end

function vumpstep(rt::VUMPSRuntime, err; show_counting = show_every_count(Inf))
    temp = show_counting()
    temp != 0 && println("vumps@step: $(temp), error=$(err)")
    M, AL, C, AR, FL, FR = rt.M, rt.AL, rt.C, rt.AR, rt.FL, rt.FR
    AC = ALCtoAC(AL,C)
    _, ACp = ACenv(AC, FL, M, FR)
    _, Cp = Cenv(C, FL, FR)
    ALp, ARp, _, _ = ACCtoALAR(ACp, Cp)
    _, FL = leftenv(AL, ALp, M, FL)
    _, FR = rightenv(AR, ARp, M, FR)
    _, ACp = ACenv(ACp, FL, M, FR)
    _, Cp = Cenv(Cp, FL, FR)
    ALp, ARp, errL, errR = ACCtoALAR(ACp, Cp)
    erroverlap = error(ALp, Cp, ARp, FL, M, FR)
    err = erroverlap + errL + errR
    # @show error(ALp, Cp, ARp, FL, M, FR)
    # err = errL + errR
    err > 1e-8 && temp >= 10 && println("errL=$errL, errR=$errR, erroverlap=$erroverlap")
    return SquareVUMPSRuntime(M, ALp, Cp, ARp, FL, FR), err
end

"""
    uptodown(i,Ni,Nj)

````
i -> (i,j) -> (Nj +1 - i,j) -> ir
````
"""
function uptodown(i,Ni,Nj)
    Liner = LinearIndices((1:Ni,1:Nj))
    Cart = CartesianIndices((1:Ni,1:Nj))
    Index = Cart[i]
    i,j = Index[1],Index[2]
    ir = Ni + 1 - i
    Liner[ir,j]
end

"""
    env = vumps_env(M::AbstractArray; χ::Int, tol::Real=1e-10, maxiter::Int=10, miniter::Int=1, verbose = false, savefile = false, infolder::String="./data/", outfolder::String="./data/", direction::String= "up", downfromup = false, show_every = Inf)

sometimes the finally observable is symetric, so we can use the same up and down environment.
"""
function vumps_env(M::AbstractArray; χ::Int, tol::Real=1e-10, maxiter::Int=10, miniter::Int=1, verbose = false, savefile = false, infolder::String="./data/", outfolder::String="./data/", direction::String= "up", downfromup = false, show_every = Inf)
    verbose && (direction == "up" ? print("↑ ") : print("↓ "))
    downfromup && direction == "down" && (direction = "up")

    D = size(M[1,1],1)
    savefile && mkpath(outfolder)
    in_chkp_file = infolder*"/$(direction)_D$(D)_χ$(χ).jld2"

    if isfile(in_chkp_file)
        rtup = SquareVUMPSRuntime(M, in_chkp_file, χ; verbose = verbose)
    else
        rtup = SquareVUMPSRuntime(M, Val(:random), χ; verbose = verbose)
    end
    env = vumps(rtup; tol=tol, maxiter=maxiter, miniter=miniter, verbose = verbose, show_every = show_every)

    Zygote.@ignore savefile && begin
        out_chkp_file = outfolder*"/$(direction)_D$(D)_χ$(χ).jld2"
        ALs, Cs, ARs, FLs, FRs = Array{Array{ComplexF64,3},2}(env.AL), Array{Array{ComplexF64,2},2}(env.C), Array{Array{ComplexF64,3},2}(env.AR), Array{Array{ComplexF64,3},2}(env.FL), Array{Array{ComplexF64,3},2}(env.FR)
        envsave = SquareVUMPSRuntime(M, ALs, Cs, ARs, FLs, FRs)
        save(out_chkp_file, "env", envsave)
    end
    env
end

"""
    M, ALu, Cu, ARu, ALd, Cd, ARd, FL, FR, envup.FL, envup.FR = obs_env(M::AbstractArray; χ::Int, tol::Real=1e-10, maxiter::Int=10, miniter::Int=1, verbose=false, savefile= false, infolder::String="./data/", outfolder::String="./data/", updown = true, downfromup = false, show_every = Inf)

If `Ni,Nj>1` and `Mij` are different bulk tensor, the up and down environment are different. So to calculate observable, we must get ACup and ACdown, which is easy to get by overturning the `Mij`. Then be cautious to get the new `FL` and `FR` environment.
"""
function obs_env(M::AbstractArray; χ::Int, tol::Real=1e-10, maxiter::Int=10, miniter::Int=1, verbose=false, savefile= false, infolder::String="./data/", outfolder::String="./data/", updown = true, downfromup = false, show_every = Inf)
    M /= norm(M)
    envup = vumps_env(M; χ=χ, tol=tol, maxiter=maxiter, miniter=miniter, verbose=verbose, savefile=savefile, infolder=infolder,outfolder=outfolder, direction="up", downfromup=downfromup, show_every = show_every)
    ALu,ARu,Cu = envup.AL,envup.AR,envup.C

    D = size(M[1,1],1)
    atype = _arraytype(M[1,1])
    in_chkp_file_obs = infolder*"/obs_D$(D)_χ$(χ).jld2"
    if isfile(in_chkp_file_obs)
        verbose && println("←→ observable environment load from $(in_chkp_file_obs)")
        FL, FR = load(in_chkp_file_obs)["env"]
        Zygote.@ignore begin
            FL, FR = Array{atype{ComplexF64,3},2}(FL), Array{atype{ComplexF64,3},2}(FR)
        end
    else
        FL, FR = envup.FL, envup.FR
    end

    if updown
        Ni, Nj = size(ALu)
        Md = [permutedims(M[uptodown(i,Ni,Nj)], (1,4,3,2)) for i = 1:Ni*Nj]
        Md = reshape(Md, Ni, Nj)

        envdown = vumps_env(Md; χ=χ, tol=tol, maxiter=maxiter, miniter=miniter, verbose=verbose, savefile=savefile, infolder=infolder, outfolder=outfolder, direction="down", downfromup=downfromup, show_every = show_every)
        ALd, ARd, Cd = envdown.AL, envdown.AR, envdown.C
    else
        ALd, ARd, Cd = ALu, ARu, Cu
    end

    _, FL = obs_FL(ALu, ALd, M, FL)
    _, FR = obs_FR(ARu, ARd, M, FR)
    Zygote.@ignore savefile && begin
        out_chkp_file_obs = outfolder*"/obs_D$(D)_χ$(χ).jld2"
        envsave = (Array{Array{ComplexF64,3},2}(FL), Array{Array{ComplexF64,3},2}(FR))
        save(out_chkp_file_obs, "env", envsave)
    end
    return M, ALu, Cu, ARu, ALd, Cd, ARd, FL, FR, envup.FL, envup.FR
end


text:
# # Baroclinic double gyre
#
# This script simulates the spin-up of a three-dimensional,
# density-stratified baroclinic double gyre in a closed basin
# using `ClimateMachine.Ocean.HydrostaticBoussinesqSuperModel`.

using Printf
using Plots
using StaticArrays
using ClimateMachine

ClimateMachine.init()

using ClimateMachine.Ocean
using ClimateMachine.Ocean.Domains
using ClimateMachine.Ocean.Fields

using ClimateMachine.GenericCallbacks: EveryXSimulationTime
using ClimateMachine.Ocean: current_step, Δt, current_time
using CLIMAParameters: AbstractEarthParameterSet, Planet

# We begin by specifying a few parameters:

wind_stress_magnitude = 2e-4 # m² s⁻²
Lx = 4e6 # m
Ly = 6e6 # m
H = 3000 # m
N² = 1e-5 # s⁻²
hour = 3600.0 # comes in handy

# Next we specify the domain

domain = RectangularDomain(
    Ne = (16, 24, 3),
    Np = 4,
    x = (0, Lx),
    y = (0, Ly),
    z = (-H, 0),
    periodicity = (false, false, false),
)

# # Boundary conditions
#
# The ``x``-component of the surface stress varies in ``y``:

wind_stress(y) = @SVector [wind_stress_magnitude * cos(2π * y / Ly), 0]

top_bc = OceanBC(Penetrable(KinematicStress(wind_stress)), Insulating())

# while we use no-slip at side walls and free-slip at the bottom,

side_wall_bc = OceanBC(Impenetrable(NoSlip()), Insulating())
bottom_bc = OceanBC(Impenetrable(FreeSlip()), Insulating())

boundary_conditions = (side_wall_bc, bottom_bc, top_bc)

# The indices of `boundary_conditions` correspond to the integers
# specified above in the kwarg `boundary` passed to `RectangularDomain`.
#
# We use parameters such that `θ` represents buoyancy, and initialize
# `θ` with a linear gradient corresponding to a linear density stratification:

struct EarthParameters <: AbstractEarthParameterSet end
Planet.grav(::EarthParameters) = 0.1

linear_gradient(x, y, z) = N² * z
initial_conditions = InitialConditions(θ = linear_gradient)

# We're now ready to build the model. We use large viscosity and diffusion
# coefficients, and Coriolis parameters appropriate for mid-latitudes on Earth.

model = Ocean.HydrostaticBoussinesqSuperModel(
    domain = domain,
    time_step = hour / 4,
    initial_conditions = initial_conditions,
    parameters = EarthParameters(),
    turbulence_closure = (νʰ = 5e3, νᶻ = 5e-3, κʰ = 1e3, κᶻ = 1e-4),
    rusanov_wave_speeds = (cʰ = 1, cᶻ = 1e-3),
    buoyancy = (αᵀ = 1 / Planet.grav(EarthParameters()),),
    coriolis = (f₀ = 1e-4, β = 1e-11),
    boundary_tags = ((1, 1), (1, 1), (2, 3)),
    boundary_conditions = boundary_conditions,
)

# We prepare a callback that periodically fetches the horizontal velocity and
# tracer concentration for later animation,

u, v, η, θ = model.fields
fetched_states = []

data_fetcher = EveryXSimulationTime(24hour) do
    @info "Step: $(current_step(model)), hours: $(current_time(model) / hour), max|u|: $(maximum(abs, u))"

    push!(
        fetched_states,
        (
            u = assemble(u),
            θ = assemble(θ),
            η = assemble(η),
            time = current_time(model),
        ),
    )
end

# and then run the simulation.

day = 24hour

model.solver_configuration.timeend = 365day

result = ClimateMachine.invoke!(
    model.solver_configuration;
    user_callbacks = [data_fetcher],
)

# Finally, we make an animation of the double gyre spinning up.

animation = @animate for (i, state) in enumerate(fetched_states)
    @info "Plotting frame $i of $(length(fetched_states))..."

    j = round(Int, domain.Ne.y / 2)

    xy_kwargs =
        (xlim = domain.x, ylim = domain.y, linewidth = 0, aspectratio = 1)
    xz_kwargs = (
        xlim = domain.x,
        ylim = domain.z,
        linewidth = 0,
        aspectratio = 1000,
    )

    x, y, z = state.u.x, state.u.y, state.u.z

    u_xy_plot = contourf(
        x,
        y,
        state.u.data[:, :, 1]';
        color = :balance,
        xy_kwargs...,
    )
    η_xy_plot = contourf(
        x,
        y,
        state.η.data[:, :, 1]';
        color = :balance,
        xy_kwargs...,
    )
    θ_xy_plot = contourf(
        x,
        y,
        state.θ.data[:, :, 1]';
        color = :thermal,
        xy_kwargs...,
    )

    u_xz_plot = contourf(
        x,
        z,
        state.u.data[:, j, :]';
        color = :balance,
        xz_kwargs...,
    )
    η_xz_plot = contourf(
        x,
        z,
        state.η.data[:, j, :]';
        color = :balance,
        xz_kwargs...,
    )
    θ_xz_plot = contourf(
        x,
        z,
        state.θ.data[:, j, :]';
        color = :thermal,
        xz_kwargs...,
    )

    u_title = @sprintf("u (m s⁻¹) at t = %.2f days", state.time / day)
    η_title = @sprintf("η (m) at t = %.2f days", state.time / day)
    θ_title = @sprintf("b (m s⁻²) at t = %.2f days", state.time / day)

    plot(
        u_xy_plot,
        η_xy_plot,
        θ_xy_plot,
        u_xz_plot,
        η_xz_plot,
        θ_xz_plot,
        title = [u_title η_title θ_title],
        size = (1200, 500),
        layout = (2, 3),
    )
end

gif(animation, "baroclinic_double_gyre.gif", fps = 8)


text:
export ActuatorSet

mutable struct ActuatorSet{ActuatorList} <: AbstractActuator
    devname::String
    actuators::ActuatorList
end

"""
`ActuatorSet(dname, acts::ActuatorList)`

Create a meta actuator from individual actuators.
"""
function ActuatorSet(dname, acts::ActuatorList) where {ActuatorList}
    return ActuatorSet(dname, actuators)
end

function ActuatorSet(dname, acts...)

    return ActuatorSet(dname, acts)
end

import Base.getindex
getindex(act::ActuatorSet, i) = act.actuators[i]

function getindex(act::ActuatorSet, dname::AbstractString)
    for dev in act.actuators
        if devname(dev) == dname
            return dev
        end
    end

    # If we got here, this the actuator doesn't exist
    # throw an exception:
    throw(KeyError(dname))
end

import Base.*

*(a1::AbstractActuator, a2::AbstractActuator) = ActuatorSet("actuators", (a1, a2))

"Number of axes (degrees of freedom) of the actuator set"
numaxes(actuators::ActuatorSet) = sum(numaxes(act) for act in actuators.actuators)

"Return the axes names of the `ActuatorSet`"
function axesnames(actuators::ActuatorSet)

    nn = String[]

    for act in actuators.actuators
        axes = axesnames(act)
        for ax in axes
            push!(nn, ax)
        end
    end
    return nn
end

"Get the actuators to `move` to a point given by the components of x"
function moveto(actuators::ActuatorSet, x)

    naxes = 0

    for act in actuators.actuators
        nx = numaxes(act)
        moveto(act, x[naxes+1:naxes+nx])
        naxes += nx
    end
    return

end