text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

using SpecialFunctions

import Compat
using Compat.Test

const SF = SpecialFunctions

if isdefined(Base, :MathConstants) && isdefined(Base.MathConstants, :γ)
    import Base.MathConstants: γ
end

# useful test functions for relative error, which differ from isapprox
# in that relerrc separately looks at the real and imaginary parts
relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
≅(a,b) = relerrc(a,b) ≤ 1e-13

@testset "error functions" begin
    @test SF.erf(Float16(1)) ≈ 0.84270079294971486934
    @test SF.erf(1) ≈ 0.84270079294971486934
    @test SF.erfc(1) ≈ 0.15729920705028513066
    @test SF.erfc(Float16(1)) ≈ 0.15729920705028513066
    @test SF.erfcx(1) ≈ 0.42758357615580700442
    @test SF.erfcx(Float32(1)) ≈ 0.42758357615580700442
    @test SF.erfcx(Complex{Float32}(1)) ≈ 0.42758357615580700442
    @test SF.erfi(1) ≈ 1.6504257587975428760
    @test SF.erfinv(0.84270079294971486934) ≈ 1
    @test SF.erfcinv(0.15729920705028513066) ≈ 1
    @test SF.dawson(1) ≈ 0.53807950691276841914

    @test SF.erf(1+2im) ≈ -0.53664356577856503399-5.0491437034470346695im
    @test SF.erfc(1+2im) ≈ 1.5366435657785650340+5.0491437034470346695im
    @test SF.erfcx(1+2im) ≈ 0.14023958136627794370-0.22221344017989910261im
    @test SF.erfi(1+2im) ≈ -0.011259006028815025076+1.0036063427256517509im
    @test SF.dawson(1+2im) ≈ -13.388927316482919244-11.828715103889593303im

    for elty in [Float32,Float64]
        for x in exp10.(Compat.range(-200, stop=-0.01, length=50))
            @test isapprox(SF.erf(SF.erfinv(x)), x, atol=1e-12*x)
            @test isapprox(SF.erf(SF.erfinv(-x)), -x, atol=1e-12*x)
            @test isapprox(SF.erfc(SF.erfcinv(2*x)), 2*x, atol=1e-12*x)
            if x > 1e-20
                xf = Float32(x)
                @test isapprox(SF.erf(SF.erfinv(xf)), xf, atol=1e-5*xf)
                @test isapprox(SF.erf(SF.erfinv(-xf)), -xf, atol=1e-5*xf)
                @test isapprox(SF.erfc(SF.erfcinv(2xf)), 2xf, atol=1e-5*xf)
            end
        end
        @test SF.erfinv(one(elty)) == Inf
        @test SF.erfinv(-one(elty)) == -Inf
        @test_throws DomainError SF.erfinv(convert(elty,2.0))

        @test SF.erfcinv(zero(elty)) == Inf
        @test_throws DomainError SF.erfcinv(-one(elty))
    end

    @test SF.erfinv(one(Int)) == SF.erfinv(1.0)
    @test SF.erfcinv(one(Int)) == SF.erfcinv(1.0)

    @test SF.erfcx(1.8) ≈ SF.erfcx(big(1.8)) rtol=4*eps()
    @test SF.erfcx(1.8e8) ≈ SF.erfcx(big(1.8e8)) rtol=4*eps()
    @test SF.erfcx(1.8e88) ≈ SF.erfcx(big(1.8e88)) rtol=4*eps()
end

@testset "sine and cosine integrals" begin
    # Computed via wolframalpha.com: SinIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]] and CosIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]]
    sinintvals = [0.9460830703671830149, 1.605412976802694849, 1.848652527999468256, 1.75820313894905306, 1.54993124494467414, 1.4246875512805065, 1.4545966142480936, 1.5741868217069421, 1.665040075829602, 1.658347594218874, 1.578306806945727416, 1.504971241526373371, 1.499361722862824564, 1.556211050077665054, 1.618194443708368739, 1.631302268270032886, 1.590136415870701122, 1.536608096861185462, 1.518630031769363932, 1.548241701043439840]
    cosintvals = [0.3374039229009681347, 0.4229808287748649957, 0.119629786008000328, -0.14098169788693041, -0.19002974965664388, -0.06805724389324713, 0.07669527848218452, 0.122433882532010, 0.0553475313331336, -0.045456433004455, -0.08956313549547997948, -0.04978000688411367560, 0.02676412556403455504, 0.06939635592758454727, 0.04627867767436043960, -0.01420019012019002240, -0.05524268226081385053, -0.04347510299950100478, 0.00515037100842612857, 0.04441982084535331654]
    for x in 1:20
        @test SF.sinint(x) ≅ sinintvals[x]
        @test SF.sinint(-x) ≅ -sinintvals[x]
        @test SF.cosint(x) ≅ cosintvals[x]
    end

    @test SF.sinint(1.f0) == Float32(SF.sinint(1.0))
    @test SF.cosint(1.f0) == Float32(SF.cosint(1.0))

    @test SF.sinint(Float16(1.0)) == Float16(SF.sinint(1.0))
    @test SF.cosint(Float16(1.0)) == Float16(SF.cosint(1.0))

    @test SF.sinint(1//2) == SF.sinint(0.5)
    @test SF.cosint(1//2) == SF.cosint(0.5)

    @test SF.sinint(1e300) ≅ π/2
    @test SF.cosint(1e300) ≅ -8.17881912115908554103E-301
    @test SF.sinint(1e-300) ≅ 1.0E-300
    @test SF.cosint(1e-300) ≅ -690.1983122333121

    @test SF.sinint(Inf) == π/2
    @test SF.cosint(Inf) == 0.0
    @test isnan(SF.sinint(NaN))
    @test isnan(SF.cosint(NaN))

    @test_throws ErrorException SF.sinint(big(1.0))
    @test_throws ErrorException SF.cosint(big(1.0))
    @test_throws DomainError SF.cosint(-1.0)
    @test_throws DomainError SF.cosint(Float32(-1.0))
    @test_throws DomainError SF.cosint(Float16(-1.0))
    @test_throws DomainError SF.cosint(-1//2)
end

@testset "airy" begin
    @test_throws SF.AmosException SF.airyai(200im)
    @test_throws SF.AmosException SF.airybi(200)

    for T in [Float32, Float64, Complex{Float32},Complex{Float64}]
        @test SF.airyai(T(1.8)) ≈ 0.0470362168668458052247
        @test SF.airyaiprime(T(1.8)) ≈ -0.0685247801186109345638
        @test SF.airybi(T(1.8)) ≈ 2.595869356743906290060
        @test SF.airybiprime(T(1.8)) ≈ 2.98554005084659907283
    end
    for T in [Complex{Float32}, Complex{Float64}]
        z = convert(T,1.8 + 1.0im)
        @test SF.airyaix(z) ≈ SF.airyai(z) * exp(2/3 * z * sqrt(z))
        @test SF.airyaiprimex(z) ≈ SF.airyaiprime(z) * exp(2/3 * z * sqrt(z))
        @test SF.airybix(z) ≈ SF.airybi(z) * exp(-abs(real(2/3 * z * sqrt(z))))
        @test SF.airybiprimex(z) ≈ SF.airybiprime(z) * exp(-abs(real(2/3 * z * sqrt(z))))
    end
    @test_throws MethodError SF.airyai(complex(big(1.0)))

    for x = -3:3
        @test SF.airyai(x) ≈ SF.airyai(complex(x))
        @test SF.airyaiprime(x) ≈ SF.airyaiprime(complex(x))
        @test SF.airybi(x) ≈ SF.airybi(complex(x))
        @test SF.airybiprime(x) ≈ SF.airybiprime(complex(x))
        if x >= 0
            @test SF.airyaix(x) ≈ SF.airyaix(complex(x))
            @test SF.airyaiprimex(x) ≈ SF.airyaiprimex(complex(x))
        else
            @test_throws DomainError SF.airyaix(x)
            @test_throws DomainError SF.airyaiprimex(x)
        end
        @test SF.airybix(x) ≈ SF.airybix(complex(x))
        @test SF.airybiprimex(x) ≈ SF.airybiprimex(complex(x))
    end
end

@testset "bessel functions" begin
    bessel_funcs = [(SF.bessely0, SF.bessely1, SF.bessely), (SF.besselj0, SF.besselj1, SF.besselj)]
    @testset "$z, $o" for (z, o, f) in bessel_funcs
        @test z(Float32(2.0)) ≈ z(Float64(2.0))
        @test o(Float32(2.0)) ≈ o(Float64(2.0))
        @test z(2) ≈ z(2.0)
        @test o(2) ≈ o(2.0)
        @test z(2.0 + im) ≈ f(0, 2.0 + im)
        @test o(2.0 + im) ≈ f(1, 2.0 + im)
    end
    @testset "besselj error throwing" begin
        @test_throws MethodError SF.besselj(1.2,big(1.0))
        @test_throws MethodError SF.besselj(1,complex(big(1.0)))
        @test_throws MethodError SF.besseljx(1,big(1.0))
        @test_throws MethodError SF.besseljx(1,complex(big(1.0)))
    end
    @testset "besselh" begin
        true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
        @test SF.besselh(3,1,3) ≈ true_h133
        @test SF.besselh(-3,1,3) ≈ -true_h133
        @test SF.besselh(3,2,3) ≈ conj(true_h133)
        @test SF.besselh(-3,2,3) ≈ -conj(true_h133)
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.besselh(1,0)
            @test_throws MethodError SF.besselh(1,big(1.0))
            @test_throws MethodError SF.besselh(1,complex(big(1.0)))
            @test_throws MethodError SF.besselhx(1,big(1.0))
            @test_throws MethodError SF.besselhx(1,complex(big(1.0)))
        end
    end
    @testset "besseli" begin
        true_i33 = 0.95975362949600785698
        @test SF.besseli(3,3) ≈ true_i33
        @test SF.besseli(-3,3) ≈ true_i33
        @test SF.besseli(3,-3) ≈ -true_i33
        @test SF.besseli(-3,-3) ≈ -true_i33
        @test SF.besseli(Float32(-3),Complex{Float32}(-3,0)) ≈ -true_i33
        true_im3p1_3 = 0.84371226532586351965
        @test SF.besseli(-3.1,3) ≈ true_im3p1_3
        for i in [-5 -3 -1 1 3 5]
            @test SF.besseli(i,0) == 0.0
            @test SF.besseli(i,Float32(0)) == 0
            @test SF.besseli(i,Complex{Float32}(0)) == 0
        end
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.besseli(1,1000)
            @test_throws DomainError SF.besseli(0.4,-1.0)
            @test_throws MethodError SF.besseli(1,big(1.0))
            @test_throws MethodError SF.besseli(1,complex(big(1.0)))
            @test_throws MethodError SF.besselix(1,big(1.0))
            @test_throws MethodError SF.besselix(1,complex(big(1.0)))
        end
    end
    @testset "besselj" begin
        @test SF.besselj(0,0) == 1
        for i in [-5 -3 -1 1 3 5]
            @test SF.besselj(i,0) == 0
            @test SF.besselj(i,Float32(0)) == 0
            @test SF.besselj(i,Complex{Float32}(0)) == 0
        end

        j33 = SF.besselj(3,3.)
        @test SF.besselj(3,3) == j33
        @test SF.besselj(-3,-3) == j33
        @test SF.besselj(-3,3) == -j33
        @test SF.besselj(3,-3) == -j33
        @test SF.besselj(3,3f0) ≈ j33
        @test SF.besselj(3,complex(3.)) ≈ j33
        @test SF.besselj(3,complex(3f0)) ≈ j33
        @test SF.besselj(3,complex(3)) ≈ j33

        j43 = SF.besselj(4,3.)
        @test SF.besselj(4,3) == j43
        @test SF.besselj(-4,-3) == j43
        @test SF.besselj(-4,3) == j43
        @test SF.besselj(4,-3) == j43
        @test SF.besselj(4,3f0) ≈ j43
        @test SF.besselj(4,complex(3.)) ≈ j43
        @test SF.besselj(4,complex(3f0)) ≈ j43
        @test SF.besselj(4,complex(3)) ≈ j43

        @test j33 ≈ 0.30906272225525164362
        @test j43 ≈ 0.13203418392461221033
        @test SF.besselj(0.1, complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im
        @test SF.besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im
        @test SF.besselj(1, 3im) ≈ 3.953370217402609396im
        @test SF.besselj(1.0,3im) ≈ SF.besselj(1,3im)

        true_jm3p1_3 = -0.45024252862270713882
        @test SF.besselj(-3.1,3) ≈ true_jm3p1_3

        @testset "Error throwing" begin
            @test_throws DomainError    SF.besselj(0.1, -0.4)
            @test_throws SF.AmosException SF.besselj(20,1000im)
            @test_throws MethodError SF.besselj(big(1.0),3im)
        end
    end

    @testset "besselk" begin
        true_k33 = 0.12217037575718356792
        @test SF.besselk(3,3) ≈ true_k33
        @test SF.besselk(-3,3) ≈ true_k33
        true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
        @test SF.besselk(3,complex(-3)) ≈ true_k3m3
        @test SF.besselk(-3,complex(-3)) ≈ true_k3m3
        # issue #6564
        @test SF.besselk(1.0,0.0) == Inf
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.besselk(200,0.01)
            @test_throws DomainError SF.besselk(3,-3)
            @test_throws MethodError SF.besselk(1,big(1.0))
            @test_throws MethodError SF.besselk(1,complex(big(1.0)))
            @test_throws MethodError SF.besselkx(1,big(1.0))
            @test_throws MethodError SF.besselkx(1,complex(big(1.0)))
        end
    end

    @testset "bessely" begin
        y33 = SF.bessely(3,3.)
        @test SF.bessely(3,3) == y33
        @test SF.bessely(3.,3.) == y33
        @test SF.bessely(3,Float32(3.)) ≈ y33
        @test SF.bessely(-3,3) ≈ -y33
        @test y33 ≈ -0.53854161610503161800
        @test SF.bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.bessely(200.5,0.1)
            @test_throws DomainError SF.bessely(3,-3)
            @test_throws DomainError SF.bessely(0.4,-1.0)
            @test_throws DomainError SF.bessely(0.4,Float32(-1.0))
            @test_throws DomainError SF.bessely(1,Float32(-1.0))
            @test_throws DomainError SF.bessely(0.4,BigFloat(-1.0))
            @test_throws DomainError SF.bessely(1,BigFloat(-1.0))
            @test_throws DomainError SF.bessely(Cint(3),Float32(-3.))
            @test_throws DomainError SF.bessely(Cint(3),Float64(-3.))

            @test_throws MethodError SF.bessely(1.2,big(1.0))
            @test_throws MethodError SF.bessely(1,complex(big(1.0)))
            @test_throws MethodError SF.besselyx(1,big(1.0))
            @test_throws MethodError SF.besselyx(1,complex(big(1.0)))
        end
    end

    @testset "besselhx" begin
        for elty in [Complex{Float32},Complex{Float64}]
            z = convert(elty, 1.0 + 1.9im)
            @test SF.besselhx(1.0, 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
            @test SF.besselhx(Float32(1.0), 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
        end
        @testset "Error throwing" begin
            @test_throws MethodError SF.besselh(1,1,big(1.0))
            @test_throws MethodError SF.besselh(1,1,complex(big(1.0)))
            @test_throws MethodError SF.besselhx(1,1,big(1.0))
            @test_throws MethodError SF.besselhx(1,1,complex(big(1.0)))
        end
    end
    @testset "scaled bessel[ijky] and hankelh[12]" begin
        for x in (1.0, 0.0, -1.0), y in (1.0, 0.0, -1.0), nu in (1.0, 0.0, -1.0)
            z = Complex{Float64}(x + y * im)
            z == zero(z) || @test SF.hankelh1x(nu, z) ≈ SF.hankelh1(nu, z) * exp(-z * im)
            z == zero(z) || @test SF.hankelh2x(nu, z) ≈ SF.hankelh2(nu, z) * exp(z * im)
            (nu < 0 && z == zero(z)) || @test SF.besselix(nu, z) ≈ SF.besseli(nu, z) * exp(-abs(real(z)))
            (nu < 0 && z == zero(z)) || @test SF.besseljx(nu, z) ≈ SF.besselj(nu, z) * exp(-abs(imag(z)))
            z == zero(z) || @test SF.besselkx(nu, z) ≈ SF.besselk(nu, z) * exp(z)
            z == zero(z) || @test SF.besselyx(nu, z) ≈ SF.bessely(nu, z) * exp(-abs(imag(z)))
        end
        @test SF.besselkx(1, 0) == Inf
        for i = [-5 -3 -1 1 3 5]
            @test SF.besseljx(i,0) == 0
            @test SF.besselix(i,0) == 0
            @test SF.besseljx(i,Float32(0)) == 0
            @test SF.besselix(i,Float32(0)) == 0
            @test SF.besseljx(i,Complex{Float32}(0)) == 0
            @test SF.besselix(i,Complex{Float32}(0)) == 0
        end
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.hankelh1x(1, 0)
            @test_throws SF.AmosException SF.hankelh2x(1, 0)
            @test_throws SF.AmosException SF.besselix(-1.01, 0)
            @test_throws SF.AmosException SF.besseljx(-1.01, 0)
            @test_throws SF.AmosException SF.besselyx(1, 0)
            @test_throws DomainError SF.besselix(0.4,-1.0)
            @test_throws DomainError SF.besseljx(0.4, -1.0)
            @test_throws DomainError SF.besselkx(0.4,-1.0)
            @test_throws DomainError SF.besselyx(0.4,-1.0)
        end
    end
    @testset "issue #6653" begin
        @testset "$f" for f in (SF.besselj,SF.bessely,SF.besseli,SF.besselk,SF.hankelh1,SF.hankelh2)
            @test f(0,1) ≈ f(0,Complex{Float64}(1))
            @test f(0,1) ≈ f(0,Complex{Float32}(1))
        end
    end
end

@testset "gamma and friends" begin
    @testset "digamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            @test SF.digamma(convert(elty, 9)) ≈ convert(elty, 2.140641477955609996536345)
            @test SF.digamma(convert(elty, 2.5)) ≈ convert(elty, 0.7031566406452431872257)
            @test SF.digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041107679516822)
            @test SF.digamma(convert(elty, 7e-4)) ≈ convert(elty, -1429.147493371120205005198)
            @test SF.digamma(convert(elty, 7e-5)) ≈ convert(elty, -14286.29138623969227538398)
            @test SF.digamma(convert(elty, 7e-6)) ≈ convert(elty, -142857.7200612932791081972)
            @test SF.digamma(convert(elty, 2e-6)) ≈ convert(elty, -500000.5772123750382073831)
            @test SF.digamma(convert(elty, 1e-6)) ≈ convert(elty, -1000000.577214019968668068)
            @test SF.digamma(convert(elty, 7e-7)) ≈ convert(elty, -1428572.005785942019703646)
            @test SF.digamma(convert(elty, -0.5)) ≈ convert(elty, .03648997397857652055902367)
            @test SF.digamma(convert(elty, -1.1)) ≈ convert(elty,  10.15416395914385769902271)

            @test SF.digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041108)
            @test SF.digamma(convert(elty, 1/2)) ≈ convert(elty, -γ - log(4))
            @test SF.digamma(convert(elty, 1)) ≈ convert(elty, -γ)
            @test SF.digamma(convert(elty, 2)) ≈ convert(elty, 1 - γ)
            @test SF.digamma(convert(elty, 3)) ≈ convert(elty, 3/2 - γ)
            @test SF.digamma(convert(elty, 4)) ≈ convert(elty, 11/6 - γ)
            @test SF.digamma(convert(elty, 5)) ≈ convert(elty, 25/12 - γ)
            @test SF.digamma(convert(elty, 10)) ≈ convert(elty, 7129/2520 - γ)
        end
    end

    @testset "trigamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            @test SF.trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
            @test SF.trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
            @test SF.trigamma(convert(elty, 1/2)) ≈ convert(elty, π^2/2)
            @test SF.trigamma(convert(elty, 1)) ≈ convert(elty, π^2/6)
            @test SF.trigamma(convert(elty, 2)) ≈ convert(elty, π^2/6 - 1)
            @test SF.trigamma(convert(elty, 3)) ≈ convert(elty, π^2/6 - 5/4)
            @test SF.trigamma(convert(elty, 4)) ≈ convert(elty, π^2/6 - 49/36)
            @test SF.trigamma(convert(elty, 5)) ≈ convert(elty, π^2/6 - 205/144)
            @test SF.trigamma(convert(elty, 10)) ≈ convert(elty, π^2/6 - 9778141/6350400)
        end
    end

    @testset "invdigamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            for val in [0.001, 0.01, 0.1, 1.0, 10.0]
                @test abs(SF.invdigamma(SF.digamma(convert(elty, val))) - convert(elty, val)) < 1e-8
            end
        end
        @test abs(SF.invdigamma(2)) == abs(SF.invdigamma(2.))
    end

    @testset "polygamma" begin
        @test SF.polygamma(20, 7.) ≈ -4.644616027240543262561198814998587152547
        @test SF.polygamma(20, Float16(7.)) ≈ -4.644616027240543262561198814998587152547
    end

    @testset "eta" begin
        @test SF.eta(1) ≈ log(2)
        @test SF.eta(2) ≈ pi^2/12
        @test SF.eta(Float32(2)) ≈ SF.eta(2)
        @test SF.eta(Complex{Float32}(2)) ≈ SF.eta(2)
    end
end

@testset "zeta" begin
    @test SF.zeta(0) ≈ -0.5
    @test SF.zeta(2) ≈ pi^2/6
    @test SF.zeta(Complex{Float32}(2)) ≈ SF.zeta(2)
    @test SF.zeta(4) ≈ pi^4/90
    @test SF.zeta(1,Float16(2.)) ≈ SF.zeta(1,2.)
    @test SF.zeta(1.,Float16(2.)) ≈ SF.zeta(1,2.)
    @test SF.zeta(Float16(1.),Float16(2.)) ≈ SF.zeta(1,2.)
    @test isnan(SF.zeta(NaN))
    @test isnan(SF.zeta(1.0e0))
    @test isnan(SF.zeta(1.0f0))
    @test isnan(SF.zeta(complex(0,Inf)))
    @test isnan(SF.zeta(complex(-Inf,0)))
end

#(compared to Wolfram Alpha)
@testset "digamma, trigamma, polygamma & zeta" begin
    for x in -10.2:0.3456:50
        @test 1e-12 > relerr(SF.digamma(x+0im), SF.digamma(x))
    end
    @test SF.digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232
    @test SF.digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im
    @test SF.digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im
    @test SF.trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546
    @test SF.trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im
    @test SF.trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im
    @test SF.polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600
    @test SF.polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im
    @test SF.polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im
    @test 1e-11 > relerr(SF.polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im)
    @test SF.polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im
    @test 1e-11 > relerr(SF.polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291)
    @test SF.zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im
    @test SF.zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im
    @test SF.zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im
    @test 1e-12 > relerrc(SF.zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im)
    @test SF.zeta(3.1,-4.2) ≅ SF.zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295
    @test 1e-15 > relerrc(SF.zeta(3.1+0im,-4.2), SF.zeta(3.1,-4.2+0im))
    @test SF.zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234
    @test SF.zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14
    @test SF.zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9
    @test 1e-12 > relerr(SF.zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14)
    @test SF.polygamma(4, -3.1+Inf*im) == SF.polygamma(4, 3.1+Inf*im) == 0
    @test SF.polygamma(4, -0.0) == Inf == -SF.polygamma(4, +0.0)
    @test SF.zeta(4, +0.0) == SF.zeta(4, -0.0) ≅ pi^4 / 90
    @test SF.zeta(5, +0.0) == SF.zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974
    @test SF.zeta(Inf, 1.) == 1
    @test SF.zeta(Inf, 2.) == 0
    @test isnan(SF.zeta(NaN, 1.))
    @test isa([SF.digamma(x) for x in [1.0]], Vector{Float64})
    @test isa([SF.trigamma(x) for x in [1.0]], Vector{Float64})
    @test isa([SF.polygamma(3,x) for x in [1.0]], Vector{Float64})
    @test SF.zeta(2 + 1im, -1.1) ≅ SF.zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im
    @test SF.polygamma(3,5) ≈ SF.polygamma(3,5.)

    @test SF.zeta(-3.0, 7.0) ≅ -52919/120
    @test SF.zeta(-3.0, -7.0) ≅ 94081/120
    @test SF.zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719
    @test SF.zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576
    @test SF.zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748
    @test SF.zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161
    @test SF.zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im
    @test SF.zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im
    @test SF.zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im
    @test SF.zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im
    @test SF.zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im
    @test SF.zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im
    @test SF.zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im
    @test SF.zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6
    @test SF.zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078
    @test SF.zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977
    @test SF.zeta(-3.2+0.1im, 0.0) ≅ SF.zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im
    @test SF.zeta(-3.2, 0.0) ≅ SF.zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378

    @test 1e-14 > relerr(SF.eta(1+1e-9), 0.693147180719814213126976796937244130533478392539154928250926)
    @test 1e-14 > relerr(SF.eta(1+5e-3), 0.693945708117842473436705502427198307157819636785324430166786)
    @test 1e-13 > relerr(SF.eta(1+7.1e-3), 0.694280602623782381522315484518617968911346216413679911124758)
    @test 1e-13 > relerr(SF.eta(1+8.1e-3), 0.694439974969407464789106040237272613286958025383030083792151)
    @test 1e-13 > relerr(SF.eta(1 - 2.1e-3 + 2e-3 * im), 0.69281144248566007063525513903467244218447562492555491581+0.00032001240133205689782368277733081683574922990400416791019im)
    @test 1e-13 > relerr(SF.eta(1 + 5e-3 + 5e-3 * im), 0.69394652468453741050544512825906295778565788963009705146+0.00079771059614865948716292388790427833787298296229354721960im)
    @test 1e-12 > relerrc(SF.zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im)
    @test 1e-13 > relerrc(SF.zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im)

    # Issue #7169: (TODO: better accuracy should be possible?)
    @test 1e-9  > relerrc(SF.zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
    @test 1e-12 > relerrc(SF.zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im)
    @test 1e-9  > relerrc(SF.zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
    @test 1e-13 > relerrc(SF.zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im)
end

@testset "vectorization of 2-arg functions" begin
    binary_math_functions = [
        SF.besselh, SF.hankelh1, SF.hankelh2, SF.hankelh1x, SF.hankelh2x,
        SF.besseli, SF.besselix, SF.besselj, SF.besseljx, SF.besselk, SF.besselkx, SF.bessely, SF.besselyx,
        SF.polygamma, SF.zeta
    ]
    @testset "$f" for f in binary_math_functions
        x = y = 2
        v = [f(x,y)]
        @test f.([x],y) == v
        @test f.(x,[y]) == v
        @test f.([x],[y]) == v
    end
end

@testset "MPFR" begin
    @testset "bessel functions" begin
        setprecision(53) do
            @test SF.besselj(4, BigFloat(2)) ≈ SF.besselj(4, 2.)
            @test SF.besselj0(BigFloat(2)) ≈ SF.besselj0(2.)
            @test SF.besselj1(BigFloat(2)) ≈ SF.besselj1(2.)
            @test SF.bessely(4, BigFloat(2)) ≈ SF.bessely(4, 2.)
            @test SF.bessely0(BigFloat(2)) ≈ SF.bessely0(2.)
            @test SF.bessely1(BigFloat(2)) ≈ SF.bessely1(2.)
        end
    end

    let err(z, x) = abs(z - x) / abs(x)
        @test 1e-60 > err(SF.eta(parse(BigFloat,"1.005")), parse(BigFloat,"0.693945708117842473436705502427198307157819636785324430166786"))
        @test 1e-60 > err(exp(SF.eta(big(1.0))), 2.0)
    end

    let a = parse(BigInt, "315135")
        @test typeof(SF.erf(a)) == BigFloat
        @test typeof(SF.erfc(a)) == BigFloat
    end
end

@testset "Base Julia issue #17474" begin
    @test SF.f64(complex(1f0,1f0)) === complex(1.0, 1.0)
    @test SF.f64(1f0) === 1.0

    @test typeof(SF.eta(big"2")) == BigFloat
    @test typeof(SF.zeta(big"2")) == BigFloat
    @test typeof(SF.digamma(big"2")) == BigFloat

    @test_throws MethodError SF.trigamma(big"2")
    @test_throws MethodError SF.trigamma(big"2.0")
    @test_throws MethodError SF.invdigamma(big"2")
    @test_throws MethodError SF.invdigamma(big"2.0")

    @test_throws MethodError SF.eta(Complex(big"2"))
    @test_throws MethodError SF.eta(Complex(big"2.0"))
    @test_throws MethodError SF.zeta(Complex(big"2"))
    @test_throws MethodError SF.zeta(Complex(big"2.0"))
    @test_throws MethodError SF.zeta(1.0,big"2")
    @test_throws MethodError SF.zeta(1.0,big"2.0")
    @test_throws MethodError SF.zeta(big"1.0",2.0)
    @test_throws MethodError SF.zeta(big"1",2.0)

    @test typeof(SF.polygamma(3, 0x2)) == Float64
    @test typeof(SF.polygamma(big"3", 2f0)) == Float32
    @test typeof(SF.zeta(1, 2.0)) == Float64
    @test typeof(SF.zeta(1, 2f0)) == Float64 # BitIntegers result in Float64 returns
    @test typeof(SF.zeta(2f0, complex(2f0,0f0))) == Complex{Float32}
    @test typeof(SF.zeta(complex(1,1), 2f0)) == Complex{Float64}
    @test typeof(SF.zeta(complex(1), 2.0)) == Complex{Float64}
end

@test sprint(showerror, SF.AmosException(1)) == "AmosException with id 1: input error."

@testset "gamma and friends" begin
    @testset "gamma, lgamma (complex argument)" begin
        if Base.Math.libm == "libopenlibm"
            @test gamma.(Float64[1:25;]) == gamma.(1:25)
        else
            @test gamma.(Float64[1:25;]) ≈ gamma.(1:25)
        end
        for elty in (Float32, Float64)
            @test gamma(convert(elty,1/2)) ≈ convert(elty,sqrt(π))
            @test gamma(convert(elty,-1/2)) ≈ convert(elty,-2sqrt(π))
            @test lgamma(convert(elty,-1/2)) ≈ convert(elty,log(abs(gamma(-1/2))))
        end
        @test lgamma(1.4+3.7im) ≈ -3.7094025330996841898 + 2.4568090502768651184im
        @test lgamma(1.4+3.7im) ≈ log(gamma(1.4+3.7im))
        @test lgamma(-4.2+0im) ≈ lgamma(-4.2)-5pi*im
        @test factorial(3.0) == gamma(4.0) == factorial(3)
        for x in (3.2, 2+1im, 3//2, 3.2+0.1im)
            @test factorial(x) == gamma(1+x)
        end
        @test lfactorial(0) == lfactorial(1) == 0
        @test lfactorial(2) == lgamma(3)
        # Ensure that the domain of lfactorial matches that of factorial (issue #21318)
        @test_throws DomainError lfactorial(-3)
        @test_throws MethodError lfactorial(1.0)
    end

    # lgamma test cases (from Wolfram Alpha)
    @test lgamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im
    @test lgamma(3.099) ≅ lgamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796
    @test lgamma(1.15) ≅ lgamma(1.15+0im) ≅ -0.06930620867104688224241731415650307100375642207340564554
    @test lgamma(0.89) ≅ lgamma(0.89+0im) ≅ 0.074022173958081423702265889979810658434235008344573396963
    @test lgamma(0.91) ≅ lgamma(0.91+0im) ≅ 0.058922567623832379298241751183907077883592982094770449167
    @test lgamma(0.01) ≅ lgamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419
    @test lgamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im
    @test lgamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im
    @test lgamma(-13.4+0.0im) ≅ conj(lgamma(-13.4-0.0im)) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im
    @test lgamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im
    @test lgamma(1+exp2(-20)) ≅ lgamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7
    @test lgamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im
    @test lgamma(-300+2im) ≅ -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im
    @test lgamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im
    @test lgamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im
    @test lgamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im
    @test lgamma(1+6.5im) ≅ conj(lgamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im
    @test lgamma(1+1im) ≅ conj(lgamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im
    @test lgamma(-pi*1e7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im
    @test lgamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im
    @test lgamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im
    @test lgamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im
    @test lgamma(2.05 + 0.03im) ≅ conj(lgamma(2.05 - 0.03im)) ≅ 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im
    @test lgamma(2+exp2(-20)+exp2(-19)*im) ≅ 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im

    @testset "lgamma for non-finite arguments" begin
        @test lgamma(Inf + 0im) === Inf + 0im
        @test lgamma(Inf - 0.0im) === Inf - 0.0im
        @test lgamma(Inf + 1im) === Inf + Inf*im
        @test lgamma(Inf - 1im) === Inf - Inf*im
        @test lgamma(-Inf + 0.0im) === -Inf - Inf*im
        @test lgamma(-Inf - 0.0im) === -Inf + Inf*im
        @test lgamma(Inf*im) === -Inf + Inf*im
        @test lgamma(-Inf*im) === -Inf - Inf*im
        @test lgamma(Inf + Inf*im) === lgamma(NaN + 0im) === lgamma(NaN*im) === NaN + NaN*im
    end
end

@testset "beta, lbeta" begin
    @test beta(3/2,7/2) ≈ 5π/128
    @test beta(3,5) ≈ 1/105
    @test lbeta(5,4) ≈ log(beta(5,4))
    @test beta(5,4) ≈ beta(4,5)
    @test beta(-1/2, 3) ≈ beta(-1/2 + 0im, 3 + 0im) ≈ -16/3
    @test lbeta(-1/2, 3) ≈ log(16/3)
    @test beta(Float32(5),Float32(4)) == beta(Float32(4),Float32(5))
    @test beta(3,5) ≈ beta(3+0im,5+0im)
    @test(beta(3.2+0.1im,5.3+0.3im) ≈ exp(lbeta(3.2+0.1im,5.3+0.3im)) ≈
          0.00634645247782269506319336871208405439180447035257028310080 -
          0.00169495384841964531409376316336552555952269360134349446910im)

    @test beta(big(1.0),big(1.2)) ≈ beta(1.0,1.2) rtol=4*eps()
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

mutable struct LatticeStruct

    origin::Array{Float64,2};
    orientx::Array{Float64,2};
    orienty::Array{Float64,2};
    orientz::Array{Float64,2};
    spacing::Array{Float64,2};
    a1::Array{Float64,2};
    a2::Array{Float64,2};
    a3::Array{Float64,2};
    basis::Array{Float64,2};
    type::Array{Int64,2};
    style::Int64;
    spaceflag::Int64;
    scale::Float64;
    LatticeStruct() = new();
end

function setlattice(stylein, scale, orientx=[1.0 0.0 0.0], orienty=[0.0 1.0 0.0], orientz=[0.0 0.0 1.0], spacing=[1.0 1.0 1.0], type=nothing, a1=[1.0 0.0 0.0], a2=[0.0 1.0 0.0], a3=[0.0 0.0 1.0], basis=nothing)

lat = LatticeStruct();

lat.scale = scale;
lat.spaceflag = 0;

if (lowercase(stylein) == "none")
    lat.style = 0;
elseif (lowercase(stylein) == "sc")
    lat.style = 1;
elseif (lowercase(stylein) == "bcc")
    lat.style = 2;
elseif (lowercase(stylein) == "fcc")
    lat.style = 3;
elseif (lowercase(stylein) == "hcp")
    lat.style = 4;
elseif (lowercase(stylein) == "diamond")
    lat.style = 5;
elseif (lowercase(stylein) == "sq")
    lat.style = 6;
elseif (lowercase(stylein) == "sq2")
    lat.style = 7;
elseif (lowercase(stylein) == "hex")
    lat.style = 8;
elseif (lowercase(stylein) == "custom")
    lat.style = 9;
else
    error("Invalid lattice style");
end

lat.origin = [0.0 0.0 0.0];
lat.orientx = orientx;
lat.orienty = orienty;
lat.orientz = orientz;
lat.spacing = spacing;

if (lowercase(stylein) == "custom")
    lat.a1 = a1;
    lat.a2 = a2;
    lat.a3 = a3;
    lat.basis = basis;
else
    lat.a1 = [1.0 0.0 0.0];
    lat.a2 = [0.0 1.0 0.0];
    lat.a3 = [0.0 0.0 1.0];

    if (lowercase(stylein) === "hex")
        lat.a2[2] = sqrt(3.0);
    end
    if (lowercase(stylein) === "hcp")
        lat.a2[2] = sqrt(3.0);
        lat.a3[3] = sqrt(8.0/3.0);
    end

    if (lowercase(stylein) === "sc")
        lat.basis = [0 0 0];
    elseif (lowercase(stylein) === "bcc")
        lat.basis = [0 0 0; 0.5 0.5 0.5];
    elseif (lowercase(stylein) === "fcc")
        lat.basis = [0 0 0; 0.5 0.5 0.0; 0.5 0.0 0.5; 0.0 0.5 0.5];
    elseif (lowercase(stylein) === "hcp")
        lat.basis = [0 0 0; 0.5 0.5 0.0; 0.5 5.0/6.0 0.5; 0.0 1.0/3.0 0.5];
    elseif (lowercase(stylein) === "diamond")
        lat.basis = [0 0 0; 0.5 0.5 0.0; 0.5 0.0 0.5; 0.0 0.5 0.5; 0.25 0.25 0.25; 0.25 0.75 0.75; 0.75 0.25 0.75; 0.75 0.75 0.25];
    elseif (lowercase(stylein) === "sq")
        lat.basis = [0 0 0];
    elseif (lowercase(stylein) === "sq2")
        lat.basis = [0 0 0; 0.5 0.5 0.0];
    elseif (lowercase(stylein) === "hex")
        lat.basis = [0 0 0; 0.5 0.5 0.0];
    end
end

nbasis = size(lat.basis)[1];

if type===nothing
    lat.type = Int64.(ones(1,nbasis));
else
    lat.type = type;
end

return lat

end

ext: jl

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module TestSubset

using Test, DataFrames, Statistics, Random

const ≅ = isequal

@testset "subset and subset!" begin
    refdf = DataFrame(x=repeat(Any[true, false], 4),
                      y=repeat([true, false, missing, missing], 2),
                      z=repeat([1, 2, 3, 3], 2),
                      id=1:8)

    for df in (copy(refdf), @view copy(refdf)[1:end-1, :])
        df2 = copy(df)
        @test subset(df, :x) ≅ filter(:x => identity, df)
        @test df ≅ df2
        @test subset(df, :x) isa DataFrame
        @test subset(df, :x, view=true) ≅ filter(:x => identity, df)
        @test subset(df, :x, view=true) isa SubDataFrame
        @test_throws ArgumentError subset(df, :y)
        @test_throws ArgumentError subset(df, :y, :x)
        @test subset(df, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :y, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :y, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :x, :y, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(df, :y, :x, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(df, :x, :y, skipmissing=true, view=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(df, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(df, :x, :y, :id => ByRow(<(4)), skipmissing=true, view=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(df, :x, :id => ByRow(<(4))) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test subset(df, :x, :id => ByRow(<(4)), view=true) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test_throws ArgumentError subset(df)
        @test isempty(subset(df, :x, :x => ByRow(!)))
        @test_throws ArgumentError subset(df, :x => x -> false, :x => x -> missing)
        @test_throws ArgumentError subset(df, :x => x -> true, :x => x -> missing)
        @test_throws ArgumentError subset(df, :x => x -> true, :x => x -> 2)
    end

    for df in (copy(refdf), @view copy(refdf)[1:end-1, :]),
        gdf in (groupby(df, :z), groupby(df, :z)[[3, 2, 1]])
        df2 = copy(df)
        @test subset(gdf, :x) ≅ filter(:x => identity, df)
        @test df ≅ df2
        @test subset(gdf, :x) isa DataFrame
        @test subset(gdf, :x, ungroup=false) ≅
              groupby(filter(:x => identity, df), :z)
        @test subset(gdf, :x, ungroup=false) isa GroupedDataFrame{DataFrame}
        @test subset(gdf, :x, view=true) ≅ filter(:x => identity, df)
        @test subset(gdf, :x, view=true) isa SubDataFrame
        @test subset(gdf, :x, view=true, ungroup=false) ≅
              groupby(filter(:x => identity, df), :z)
        @test subset(gdf, :x, view=true, ungroup=false) isa GroupedDataFrame{<:SubDataFrame}
        @test_throws ArgumentError subset(gdf, :y)
        @test_throws ArgumentError subset(gdf, :y, :x)
        @test subset(gdf, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :y, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :y, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :x, :y, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(gdf, :y, :x, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(gdf, :x, :y, skipmissing=true, view=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true, view=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(gdf, :x, :id => ByRow(<(4))) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test subset(gdf, :x, :id => ByRow(<(4)), view=true) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test_throws ArgumentError subset(gdf)
        @test isempty(subset(gdf, :x, :x => ByRow(!)))
        @test_throws ArgumentError subset(gdf, :x => x -> false, :x => x -> missing)
        @test_throws ArgumentError subset(gdf, :x => x -> true, :x => x -> missing)
        @test_throws ArgumentError subset(gdf, :x => x -> true, :x => x -> 2)
    end

    df = copy(refdf)
    @test subset!(df, :x) === df
    @test subset!(df, :x) ≅ df ≅ filter(:x => identity, refdf)
    df = copy(refdf)
    @test_throws ArgumentError subset!(df, :y)
    @test df ≅ refdf
    df = copy(refdf)
    @test subset!(df, :y, skipmissing=true) === df
    @test subset!(df, :y, skipmissing=true) ≅ df ≅ filter(:y => x -> x === true, refdf)
    df = copy(refdf)
    @test subset!(df, :x, :y, skipmissing=true) === df
    @test subset!(df, :x, :y, skipmissing=true) ≅ df ≅
          filter([:x, :y] => (x, y) -> x && y === true, refdf)
    df = copy(refdf)
    @test subset!(df, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅ df ≅
          filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, refdf)
    df = copy(refdf)
    @test subset!(df, :x, :id => ByRow(<(4))) ≅ df ≅
          filter([:x, :id] => (x, id) -> x && id < 4, refdf)
    df = copy(refdf)
    @test_throws ArgumentError subset!(df)
    df = copy(refdf)
    @test isempty(subset!(df, :x, :x => ByRow(!)))
    @test isempty(df)

    df = copy(refdf)
    @test_throws ArgumentError subset!(df, :x => x -> false, :x => x -> missing)
    @test_throws ArgumentError subset!(df, :x => x -> true, :x => x -> missing)
    @test_throws ArgumentError subset!(df, :x => x -> true, :x => x -> 2)

    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x) === df

    df = copy(refdf)
    gdf = groupby(df, :z)
    gdf2 = subset!(gdf, :x, ungroup=false)
    @test gdf2 isa GroupedDataFrame{DataFrame}
    @test parent(gdf2) === df
    @test gdf2 ≅ groupby(df, :z) ≅ groupby(filter(:x => identity, refdf), :z)

    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x) ≅ df ≅ filter(:x => identity, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test_throws ArgumentError subset!(gdf, :y)
    @test df ≅ refdf
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :y, skipmissing=true) === df
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :y, skipmissing=true) ≅ df ≅ filter(:y => x -> x === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :y, skipmissing=true) === df
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :y, skipmissing=true) ≅ df ≅
          filter([:x, :y] => (x, y) -> x && y === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅ df ≅
          filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :id => ByRow(<(4))) ≅ df ≅
          filter([:x, :id] => (x, id) -> x && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test_throws ArgumentError subset!(gdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test isempty(subset!(gdf, :x, :x => ByRow(!)))
    @test isempty(df)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test_throws ArgumentError subset!(gdf, :x => x -> false, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> 2)

    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x) ≅ df ≅ filter(:x => identity, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test_throws ArgumentError subset!(gdf, :y)
    @test df ≅ refdf
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :y, skipmissing=true) ≅ df ≅ filter(:y => x -> x === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x, :y, skipmissing=true) ≅ df ≅
          filter([:x, :y] => (x, y) -> x && y === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅ df ≅
          filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x, :id => ByRow(<(4))) ≅ df ≅
          filter([:x, :id] => (x, id) -> x && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test_throws ArgumentError subset!(gdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test isempty(subset!(gdf, :x, :x => ByRow(!)))
    @test isempty(df)

    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test_throws ArgumentError subset!(gdf, :x => x -> false, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> 2)

    @test_throws ArgumentError subset!(view(refdf, :, :), :x)
    @test_throws ArgumentError subset!(groupby(view(refdf, :, :), :z), :x)

    df = DataFrame(g=[2, 2, 1, 1, 1, 1, 3, 3, 3], x=1:9)
    @test subset(df, :x => x -> x .< mean(x)) == DataFrame(g=[2, 2, 1, 1], x=1:4)
    @test subset(groupby(df, :g), :x => x -> x .< mean(x)) ==
          DataFrame(g=[2, 1, 1, 3], x=[1, 3, 4, 7])

    @test_throws ArgumentError subset(df, :x => x -> missing)
    @test isempty(subset(df, :x => ByRow(x -> missing), skipmissing=true))
    @test_throws ArgumentError isempty(subset(df, :x => x -> missing, skipmissing=true))
    @test isempty(subset(df, :x => ByRow(x -> false)))
    @test_throws ArgumentError isempty(subset(df, :x => x -> false))
    @test subset(df, :x => ByRow(x -> true)) ≅ df
    @test_throws ArgumentError subset(df, :x => x -> true) ≅ df
    @test_throws ArgumentError subset(df, :x => x -> (a=x,))
    @test_throws ArgumentError subset(df, :x => (x -> (a=x,)) => AsTable)

    @test_throws ArgumentError subset(DataFrame(x=false, y=missing), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=missing, y=false), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=missing, y=false), :x)
    @test_throws ArgumentError subset(DataFrame(x=false, y=missing), :y)
    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :y, :x)
    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :y)

    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :x, :y, skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :x, :y, skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :y, :x, skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :y, skipmissing=true)

    @test_throws ArgumentError DataFrames._and()
    @test_throws ArgumentError DataFrames._and_missing()
end

@testset "subsetting requires passing vector" begin
    @test_throws ArgumentError subset(DataFrame(x=[]), :x => x -> 1)
    @test_throws AssertionError subset(DataFrame(x=[]), :x => ByRow(x -> 1))
    @test_throws ArgumentError subset(DataFrame(x=1:3), [] => () -> true)
    @test subset(DataFrame(x=1:3), [] => ByRow(() -> true)) == DataFrame(x=1:3)
    @test_throws ArgumentError subset(DataFrame(x=[0, 1]), :x => ==(0))
    @test_throws ArgumentError subset(DataFrame(x=[0, missing]), :x => ismissing)
    @test_throws ArgumentError subset(groupby(DataFrame(id=[0, 0, 1, 1],
                                                        x=[-1, 1, 3, 4]), :id),
                                      :x => (x -> sum(x) > 0))
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> fill(missing, length(x)))
    @test isempty(subset(DataFrame(x=1:3), :x => x -> fill(missing, length(x)),
                         skipmissing=true))
    @test subset(DataFrame(x=1:3), :x => x -> fill(true, length(x))) == DataFrame(x=1:3)
    @test subset(DataFrame(x=1:3), :x => x -> trues(length(x))) == DataFrame(x=1:3)
    @test isempty(subset(DataFrame(x=1:3), :x => x -> fill(false, length(x))))
    @test isempty(subset(DataFrame(x=1:3), :x => x -> falses(length(x))))
    @test_throws AssertionError subset(DataFrame(), [] => () -> Union{}[])
    @test_throws AssertionError subset(DataFrame(), [] => () -> Union{}[], skipmissing=true)
    @test subset(DataFrame(x=1:3:15, y=1:5), [:x, :y] => (x, y) -> iseven.(x) .& iseven.(y)) ==
          DataFrame(x=[4, 10], y=[2, 4])
    @test subset(DataFrame(x=1:3:15, y=1:5),
                 AsTable([:x, :y]) => v -> iseven.(v.x) .& iseven.(v.y)) ==
          DataFrame(x=[4, 10], y=[2, 4])

    @test subset(DataFrame(x=1:3), :x => x -> Any[true, false, true]) ==
          DataFrame(x=[1, 3])
    @test subset(DataFrame(x=1:3), :x => x -> view(Any[true, false, true], :)) ==
          DataFrame(x=[1, 3])
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> Any[true, false, missing])
    @test subset(DataFrame(x=1:3), :x => x -> Any[true, false, missing], skipmissing=true) ==
          DataFrame(x=[1])
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> Any[true, false, 1])
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> Any[true, false, 1],
                                      skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> (true for i in 1:3))
end

@testset "multicolumn selectors" begin
    dfref = DataFrame(1:10 .!= permutedims(1:11), :auto)

    for df in (dfref, groupby(dfref, :x2), groupby(dfref, :x11))
        @test subset(df, Not(:x1) .=> identity) ==
            subset(df, Not(:x1) .=> [identity identity]) ==
            subset(df, (names(df, Not(:x1)) .=> identity)...)
        @test subset(df, Not(:x1)) == subset(df, names(df, Not(:x1))...)
        @test subset(df, :) == subset(df, names(df)) == subset(df, names(df)...)
        @test subset(df, Cols(:) .=> ByRow(x -> true)) == subset(df, (names(df, Cols(:)) .=> ByRow(x -> true))...)
        @test subset(df, [], [:x1]) == subset(df, :x1)
        @test subset(df, Matrix{Any}(undef, 0, 0), [:x1]) == subset(df, :x1)
        @test subset(df, [:x1, :x3, :x5]) == subset(df, [1, 3, 5]) == subset(df, :x1, :x3, :x5)
        @test subset(df, [false; fill(true, 10)]) == subset(df, (2:11)...)
        @test_throws ArgumentError subset(df, [:x1 => identity => :x1])
        @test_throws ArgumentError subset(df, Not(:x1) .=> identity .=> Not(:x1))
        @test_throws ArgumentError subset(df, Not(:x1) .=> [identity identity] .=> Not(:x1))
        @test_throws ArgumentError subset(df, [:y])
        @test_throws BoundsError subset(df, [-1])
        @test_throws BoundsError subset(df, [true])

        Random.seed!(12345)
        df = DataFrame(rand(100, 5), :auto)
        columns = names(df, Not(Between(:x3, ncol(df))))
        @test subset(df, vcat.(columns, "x5") .=> ByRow(<)) ==
              subset(df, [:x1, :x5] => ByRow(<), [:x2, :x5] => ByRow(<)) ==
              df[(df.x1 .< df.x5) .& (df.x2 .< df.x5), :]
    end
end

end


text:
lines = readlines("inputs/02.txt")
forward = 0
vert = 0
aim = 0
for x in lines
  num = parse(Int64, split(x,' ')[2])
  if x[1] == 'f'
    forward += num
    vert += aim * num
  elseif x[1] == 'd'
    aim += num
  else
    aim -= num
  end
end
println(forward * vert)


text:
module CachedSubVectorsTests

using Test
using Gridap
using Gridap.CachedSubVectors

v = collect(1:10)
c = CachedSubVector(v)
@test v == c

c = CachedSubVector(v,1,10)
@test v == c

locate!(c,3,7)
@test v[3:7] == c

c[1] = -1
@test c[1] == -1

end # module


text:
"""
    module ExactlyRemoveSingularities - exactly remove singularities from A*v = 0 with A an Integer matrix.

Use this module in the following way:

   using ExactlyRemoveSingularities
   result = removeSingularities(A,ix,iy)    # result = (iya, eqr, ix1, ix2, eqx, A1, A2)
   printRemoveSingularities(A,result,names)

where A is a dense or sparse Int matrix A[:,:] describing a linear system

    A*v = 0

with

  - v[ix]: Variables used in the derivative operator, der(v[ix]), so these are the potential states.

  - v[iy]: Variables that must be solved from ```A*v = 0```, because these variables do not appear in the remaining part of the model equations (e.g. absolute potentials). If variables v[iy] cannot be uniquely solved from ```A*v = 0```, then in argument iya the variables are listed that must be set arbitrarily so that v[iy] can be uniquely solved from the equations. For example, if v[iy] contains the absolute potentials and no ground object is present, then v[iya] are the absolute potentials that must be explicitly set, e.g., by introducing an equation v[iya]=0, in order that the equations have a unique solution. Additionally, a warning message should be printed that v[iya] was set to 0.)

The function return arguments are all of type Int and have the following meaning:

  - iya: variables v[iya] can be arbitrarily set (iya is a sub-set of iy)

  - eqr: equations ```A[eqr,:]*v = 0``` are redundant and should be removed (they do not change the result)

  - ix1: x[ix1] are the dependent states (ix1 is a sub-set of ix)

  - ix2: i[ix2] are the independent states (ix2 is a sub-set of ix)

  - eqx, A1, A2: equations ```A[eqx,:]*v = 0``` should be replaced by the equations ```A1*x1 + A2*x2 = 0```, where A1 is square, upper triangular and has full rank with A1[i,i] != 0. In other words, x1 can be computed explicitely from these equations given x2.

The calling program should do the following:

   - Remove equations ```A[eqr,:]*v = 0```, because these are redundant equations.
   - Replace equations ```A[eqx,:]*v = 0``` by ```A1*x1 + A2*x2 = 0```, in order that index reduction can be performed by a structural algorithm (e.g. Pantelides algorithm).
   - Introduce equations v[iya] = 0 and print a warning message that these variables have been arbitrarily set to zero. This is needed, in order that the equations have a unique solution.
"""
module ExactlyRemoveSingularities

export removeSingularities, upperTrapezoidal
export printobj, printRemoveSingularities

"""
    (Afac, rk, p1, p2, Bfac) = upperTrapezoidal(A,B)
    (Afac, rk, p1, p2)       = upperTrapezoidal(A)

A and B are **Int** matrices with dimensions A[n1,n2], B[n1,:] or B[n1], so B can be a matrix or a vector. Furthermore A and B can be dense or sparse matrices (SparseMatrixCSC). If B is not present, it is interpreted
as zero matrix. Typically A has more columns as rows (n2 > n1). The function transforms the equation

    A*X = B

to

    Afac*X[p2,:] = Bfac

with

- Afac an **upper trapezoidal** Int matrix with Afac[i,i] != 0 for i<=rk; Afac[rk+1:end,:] = 0

- Bfac an Int matrix or vector (not present if B is not present)

- rk the rank of A and of Afac

- p1 the row permutation vector of A

- p2 the column permutation vector of A

The transformation is performed by a fraction-free Gauss elimination with
exact (integer) computation. It is then possible to solve for X by splitting
X[p2,:] = [X1;X2] with X1[rk,rk], providing any values for X2 and computing X1 from

    Afac[1:rk,1:rk]*X1 = Bfac[1:rk,:] - Afac[1:rk,rk+1:end]*X2

Since Afac[1:rk,1:rk] is **upper triangular** with non-zeros on the diagonal,
this is just a forward recursion. Note, the (full) pivoting is performed
in such a way that it is tried to utilize the **upper** part of X in X1.
Therefore, variables that must be solved for should be in the upper part of X.

The permutation vector p1 can be used to determine the equations in A and B that
can be removed without changing the result (provided Bfac[rk+1:end,:] = 0):

    p1red = p1[rg+1:end]
    A[p1red,:]*X = B[p1red,:] are redundant equations that can be removed

This implementation is based on:

  - Bareiss E.H. (1968): Sylvester’s identity and multistep integer-preserving Gaussian elimination. Mathematics of Computation, Vol. 22, No. 103, pp. 565 – 578. Download: http://www.ams.org/journals/mcom/1968-22-103/S0025-5718-1968-0226829-0/

  - Turner P.R. (1995): A simplified fraction-free Integer elimination algorithm. Naval Air Warfare Center. Report NAWCADPAX--96-196-TR. Download: http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA313755

  - Dureisseix D. (2012): Generalized fraction-free LU factorization for singular systems with kernel extraction. Linear Algebra and its Applications, Volume 436, Issue 1, pp. 27-40. Download: http://www.sciencedirect.com/science/article/pii/S0024379511004617#
"""
function upperTrapezoidal(A::Matrix{Int}, B::Matrix{Int})
   assert(size(A,1) == size(B,1))
   (n1,n2) = size(A)
   nb   = size(B,2)
   p1   = collect(1:n1)
   p2   = collect(1:n2)
   Afac = copy(A)
   Bfac = copy(B)
   rk   = n1
   p1k  = 0
   p2k  = 0
   Aik  = 0
   temp = 0
   tempRow    = fill(0,1,nb)
   pivot      = 1
   oldPivot   = 1
   pivotFound = false

   for k = 1:n1
      # Search for a pivot in Afac[k:n1,min(k,n2):n2]
      if k <= n2
         if Afac[k,k] == 0
            pivotFound = false
            for k2 = k:n2
               for k1 = k:n1
                  if Afac[k1,k2] != 0
                     pivotFound = true
                     p1k = k1
                     p2k = k2
                     break
                  end
               end
               pivotFound && break
            end

            if pivotFound
               # Row interchange
               if p1k != k
                  for i = p2k:n2
                     temp = Afac[k,i]
                     Afac[k  ,i] = Afac[p1k,i]
                     Afac[p1k,i] = temp
                  end
                  tempRow = Bfac[k,:]
                  Bfac[k,:]   = Bfac[p1k,:]
                  Bfac[p1k,:] = tempRow

                  temp    = p1[k]
                  p1[k]   = p1[p1k]
                  p1[p1k] = temp
               end

               # Column interchange
               if p2k != k
                  for i = 1:n1
                     temp = Afac[i,k]
                     Afac[i,k  ] = Afac[i,p2k]
                     Afac[i,p2k] = temp
                  end
                  temp    = p2[k]
                  p2[k]   = p2[p2k]
                  p2[p2k] = temp
               end
            else
               # no pivot found: submatrix Afac[k:n1,:] has only zeros
               rk = k-1
               return (Afac, rk, p1, p2, Bfac)
            end
         end
      else
         # submatrix Afac[k:n1,:] has only zeros
         rk = k-1
         return (Afac, rk, p1, p2, Bfac)
      end

      # Afac[k,k] != 0
      pivot = Afac[k,k]
      for i = k+1:n1
         Aik  = Afac[i,k]
         Bfac[i,:] = div.(pivot*Bfac[i,:] - Aik*Bfac[k,:], oldPivot)  # guaranteed no remainder
         for j = k+1:n2
            Afac[i,j] = div(pivot*Afac[i,j] - Aik*Afac[k,j], oldPivot)  # guaranteed no remainder
         end
         Afac[i,k] = 0
      end
      oldPivot = pivot
   end
   return (Afac, rk, p1, p2, Bfac)
end

toVector(result) = (result[1], result[2], result[3], result[4], vec(result[5]))
upperTrapezoidal(A::Matrix{Int}, b::Vector{Int}) = toVector(upperTrapezoidal(A, reshape(b,size(b,1),1)))
upperTrapezoidal(A::Matrix{Int})                 =  upperTrapezoidal(A, fill(0,size(A,1),0))

function upperTrapezoidal(A::SparseMatrixCSC{Int,Int}, B::SparseMatrixCSC{Int,Int})
   assert(size(A,1) == size(B,1))
   (n1,n2) = size(A)
   nb   = size(B,2)
   p1   = collect(1:n1)
   p2   = collect(1:n2)
   Afac = copy(A)
   Bfac = copy(B)

   rk   = n1
   p1k  = 0
   p2k  = 0
   Aik  = 0
   temp = 0
   pivot      = 1
   oldPivot   = 1
   pivotFound = false

   for k = 1:n1
      # Search for a pivot in Afac[k:n1,min(k,n2):n2] columnwise
      if k <= n2
         if Afac[k,k] == 0
            Arows = rowvals(Afac)
            pivotFound = false
            Avals = nonzeros(Afac)
            for k2 = k:n2
                for kk1 in nzrange(Afac, k2)
                   k1  = Arows[kk1]
                   val = Avals[kk1]
                    if k1 >= k && val != 0
                     pivotFound = true
                     p1k = k1
                     p2k = k2
                     break
                  end
               end
               pivotFound && break
            end

            if pivotFound
               # Row interchange
               if p1k != k
                  ii = p2k:n2
                  tempAr       = Afac[k,ii]
                  Afac[k  ,ii] = Afac[p1k,ii]
                  Afac[p1k,ii] = tempAr

                  if nb > 0
                     tempRow = Bfac[k,:]
                     Bfac[k,:]   = Bfac[p1k,:]
                     Bfac[p1k,:] = tempRow
                  end

                  temp    = p1[k]
                  p1[k]   = p1[p1k]
                  p1[p1k] = temp
               end

               # Column interchange
               if p2k != k
                  tempAc      = Afac[:,k]
                  Afac[:,k  ] = Afac[:,p2k]
                  Afac[:,p2k] = tempAc

                  temp    = p2[k]
                  p2[k]   = p2[p2k]
                  p2[p2k] = temp
               end
            else
               # no pivot found: submatrix Afac[k:n1,:] has only zeros
               rk = k-1
               return (Afac, rk, p1, p2, Bfac)
            end
         end
      else
         # no pivot found: submatrix Afac[k:n1,:] has only zeros
         rk = k-1
         return (Afac, rk, p1, p2, Bfac)
      end

      # Subtract row k from row i = k+1:n1
      j = k+1:n2
      for ii in nzrange(Afac, k)    # [i,k], i = k:n1
         Arows = rowvals(Afac)
         Avals = nonzeros(Afac)
         i = Arows[ii]
         if i == k
            # Afac[k,k] != 0
            pivot = Avals[ii]
         elseif i > k
            Aik = Avals[ii]
            Avals[ii] = 0
            if nb > 0
               Bfac[i,:] = div.(pivot*Bfac[i,:] - Aik*Bfac[k,:], oldPivot)  # guaranteed no remainder
            end
            Afac[i,j] = div.(pivot*Afac[i,j] - Aik*Afac[k,j], oldPivot)  # guaranteed no remainder
         end
      end
      oldPivot = pivot
   end
   return (Afac, rk, p1, p2, Bfac)
end

upperTrapezoidal(A::SparseMatrixCSC{Int,Int}) =  upperTrapezoidal(A, sparse(fill(0,size(A,1),0)))

"""-----------------------------------------------------------------------
    (iya, eqr, ix1, ix2, eqx, A1, A2) = removeSingularities(A, ix, iy)

A is a dense or sparse Int matrix A[:,:] describing a linear system

    A*v = 0

where

  - v[ix]: Variables used in the derivative operator, der(v[ix]), so these are the potential states.

  - v[iy]: Variables that must be solved from ```A*v = 0```, because these variables do not appear in the remaining part of the model equations (e.g. absolute potentials). If variables v[iy] cannot be uniquely solved from ```A*v = 0```, then in argument iya the variables are listed that must be set arbitrarily so that v[iy] can be uniquely solved from the equations. For example, if v[iy] contains the absolute potentials and no ground object is present, then v[iya] are the absolute potentials that must be explicitly set, e.g., by introducing an equation v[iya]=0, in order that the equations have a unique solution. Additionally, a warning message should be printed that v[iya] was set to 0.)

The function return arguments are all of type Int and have the following meaning:

  - iya: variables v[iya] can be arbitrarily set (iya is a sub-set of iy)

  - eqr: equations ```A[eqr,:]*v = 0``` are redundant and should be removed (they do not change the result)

  - ix1: x[ix1] are the dependent states (ix1 is a sub-set of ix)

  - ix2: i[ix2] are the independent states (ix2 is a sub-set of ix)

  - eqx, A1, A2: equations ```A[eqx,:]*v = 0``` should be replaced by the equations ```A1*x1 + A2*x2 = 0```, where A1 is square, upper triangular and has full rank with A1[i,i] != 0. In other words, x1 can be computed explicitely from these equations given x2.

The calling program should do the following:

   - Remove equations ```A[eqr,:]*v = 0```, because these are redundant equations.
   - Replace equations ```A[eqx,:]*v = 0``` by ```A1*x1 + A2*x2 = 0```, in order that index reduction can be performed by a structural algorithm (e.g. Pantelides algorithm).
   - Introduce equations v[iya] = 0 and print a warning message that these variables have been arbitrarily set to zero. This is needed, in order that the equations have a unique solution.
"""
function removeSingularities(A::T, ix::Vector{Int}, iy::Vector{Int}) where T <: AbstractMatrix{Int}
  # Split equation system A*v = 0 into V*[vy,vr] = X*(-x), where
  #   x  = v[ix]
  #   vy = v[iy]
  #   vr = the remaining variables
  (ne,nv) = size(A)   # number of equations (ne) and variables (nv)
  nx = length(ix)     # number of potential states
  set_iy = Set(iy)
  ivr = collect( setdiff( setdiff(Set(1:nv), Set(ix)), set_iy ) )  # indices of vr
  ivall = [iy;ivr]
  V = [ A[:,iy] A[:,ivr] ]
  X = A[:,ix];

  # Transform V to upper trapezoidal form
  #    Vfac*v[pv2] = Xfac*(-x)
  (Vfac, rk, pv1, pv2, Xfac) = upperTrapezoidal(V, X)

  # iya are the independent variables that are in set_iy
  ivind = ivall[ pv2[rk+1:end] ]
  iya   = collect( intersect(Set(ivind), set_iy) )

  if nx == 0
     eqr = pv1[rk+1:ne]
     eqx = fill(0,0)
     ix1 = fill(0,0)
     ix2 = fill(0,0)
     A1  = fill(0,0,0)
     A2  = fill(0,0,0)

  else
     # Transform state constraints to upper trapezoidal form
     #    Xc*x = 0     // Xc = Xfac[rk+1:end,:]
     #
     #    [Xcfac1 XcFac2;
     #       0       0  ] * [x1,x2] = 0
     Xc = Xfac[rk+1:end,:]
     (Xcfac, rk2, px1, px2) = upperTrapezoidal( Xc )
     eqx2 = px1[1:rk2]
     eqr2 = px1[rk2+1:end]

     eqx  = pv1[rk+eqx2]
     eqr  = pv1[rk+eqr2]
     ix1  = ix[ px2[1:rk2] ]
     ix2  = ix[ px2[rk2+1:end] ]
     A1   = Xcfac[1:rk2, 1:rk2]
     A2   = Xcfac[1:rk2, rk2+1:end]
  end

  return (iya, eqr, ix1, ix2, eqx, A1, A2)
end

"""---------------------------------------------------------------
    str = displayAsString(x)

Return the result of the Julia built-in command display(x) as string
instead of displaying it on the terminal.
"""
displayAsString(x) = sprint((io,x)->display(TextDisplay(io),x),x)

"""---------------------------------------------------------------
    str = displayAsString2(x)

Return the result of the Julia built-in command display(x) as string
instead of displaying it on the terminal, remove the preceding type
information and shift the whole output to the right
"""
function displayAsString2(obj)
   str  = displayAsString(obj)
   str2 = str[searchindex(str,"\n")+1:end]
   str3 = "   " * replace(str2,"\n","\n   ")
end

"""---------------------------------------------------------------
    printobj(name,obj)

Pretty print an instance obj as "<name> = <display(obj)>".
"""
function printobj(name,obj)
   str = displayAsString2(obj)
   println("\n$name =\n$str")
end

function printEquation(coeff, names::Vector{String})
   assert(length(coeff) == length(names))
   print("   ")
   first = true
   for i=1:length(names)
      if coeff[i] != 0
         if first
            first = false
            if coeff[i] < 0
               print("-")
            end
         else
            if coeff[i] < 0
               print(" - ")
            elseif coeff[i] > 0
               print(" + ")
            end
         end

         c = abs(coeff[i])
         if c != 1
            print(c,"*")
         end
         print(names[i])
      end
   end
   println(" = 0")
end

"""
    printRemoveSingularities(A, r, names)

Prints the result of r=removeSingularities(..), where

  - A: matrix of removeSingularities(..)

  - r: the tuple returned by removeSingularities.

  - names: string vector containing the names of the variables v from A*v = 0.
"""
function printRemoveSingularities(A, r, names)
   (iya, eqr, ix1, ix2, eqx, A1, A2) = r

   if length(iya) > 0
      iya_names = displayAsString2(names[iya])
      println("\n  Variables that can be arbitrarily set:\n", iya_names)
   end

   if length(eqr) > 0
      println("\n  Equations that can be removed, since redundant: ", eqr)
      for i=1:length(eqr)
         printEquation(A[eqr[i],:],names)
      end
   end

   if length(ix1) > 0
      ix1_names = displayAsString2(names[ix1])
      ix2_names = displayAsString2(names[ix2])
      println("\n  Dependent state variables x1:\n"  ,ix1_names)
      println("\n  Independent state variables x2:\n",ix2_names)
      println("\n  Equations: ", eqx)
      for i=1:length(eqx)
         printEquation(A[eqx[i],:],names)
      end
      println("\n  that shall be replaced by A1*x1 + A2*x2 = 0: ")
      printobj("   A1", A1)
      printobj("   A2", A2)
      print("\n")
      for i=1:size(A1,1)
         printEquation([A1 A2][i,:],[names[ix1];names[ix2]])
      end
   end
end

end


text:
## Utilities

function symmetrize!(A::Array)
    for i=1:size(A,1)
        for j=1:i-1
            A[i,j] = A[j,i]
        end
    end
    return A
end


text:
using HypothesisTests, Test
using StatsBase
using HypothesisTests: default_tail

@testset "Power Divergence" begin
#Example 1 in R
#Agresti (2007) p. 39

d = [[762,484] [327,239] [468,477]]

m = PowerDivergenceTest(d)

@test m.theta0 ≈ [0.25523082406125785,0.19670969099133556,0.11593952361049113,0.08935608756107216,0.1935739395970214,0.1491899341788219]
@test m.thetahat ≈ [0.2763873775843308,0.1755531374682626,0.11860718171926006,0.08668842945230323,0.16974972796517954,0.17301414581066377]

c = confint(m)
c0 = [(.25788900979, .29554669435),
      (.15705476968, .19471245423),
      (.10010881393, .13776649848),
      (.06819006166, .10584774622),
      (.15125136017, .18890904473),
      (.15451577802, .19217346257)]

for i = 1:length(c)
    @test c[i][1] ≈ c0[i][1]
    @test c[i][2] ≈ c0[i][2]
end

@test pvalue(m) ≈ 2.9535891832117357e-7
@test default_tail(m) == :right
@test m.stat ≈ 30.070149095754687
@test m.df ≈ 2
@test m.n ≈ 2757
@test m.residuals ≈ reshape([2.198855766015898,-2.504669492560728,0.4113701700566286,-0.46858294710127296,-2.8432397155451494,3.2386734435365825], size(d))
@test m.stdresiduals ≈ reshape([4.502053521086705,-4.502053521086705,0.6994517329844298,-0.6994517329844298,-5.315945542704929,5.315945542704929], size(d))

m = PowerDivergenceTest(d,lambda=0.0)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=-1.0)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=-2.0)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=-0.5)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=2/3)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = ChisqTest(d)
m = MultinomialLRT(d)

confint(m, method = :bootstrap)
confint(m, method = :bootstrap, tail=:left)
confint(m, method = :bootstrap, tail=:right)

confint(m, method = :gold)
confint(m, method = :gold, tail=:left)
confint(m, method = :gold, tail=:right)

confint(m, method = :quesenberry_hurst)
confint(m, method = :quesenberry_hurst, tail=:left)
confint(m, method = :quesenberry_hurst, tail=:right)

confint(m, method = :sison_glaz)
confint(m, method = :sison_glaz, correct=false)
confint(m, method = :sison_glaz, tail=:left)
confint(m, method = :sison_glaz, tail=:right)

@test_throws ArgumentError confint(m, method=:FOO)
@test_throws ArgumentError confint(m, tail=:fox)

#Example 3 in R

d = [ 20, 15, 25 ]
m = PowerDivergenceTest(d)

@test m.theta0 ≈ [0.3333333333333333,0.3333333333333333,0.3333333333333333]
@test m.thetahat ≈ [0.3333333333333333,0.25,0.4166666666666667]

c = confint(m)
c0 = [(.21666666667, .48098082062),
      (.13333333333, .39764748728),
      (.30000000000, .56431415395)]

for i in 1:length(c)
    @test c[i][1] ≈ c0[i][1]
    @test c[i][2] ≈ c0[i][2]
end

@test pvalue(m) ≈ 0.2865047968601901
@test m.stat ≈ 2.5
@test m.df ≈ 2
@test m.n ≈ 60
@test m.residuals ≈ [0.0,-1.118033988749895,1.118033988749895]
@test m.stdresiduals ≈ [0.0,-1.3693063937629153,1.3693063937629153]

m = PowerDivergenceTest(d,lambda=0.0)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=-1.0)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=-2.0)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=-0.5)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = PowerDivergenceTest(d,lambda=2/3)
testname(m)
pvalue(m)
show(IOBuffer(), m)

m = ChisqTest(d)
m = MultinomialLRT(d)

confint(m, method = :bootstrap)
confint(m, method = :bootstrap, tail=:left)
confint(m, method = :bootstrap, tail=:right)

confint(m, method = :gold)
confint(m, method = :gold, tail=:left)
confint(m, method = :gold, tail=:right)

confint(m, method = :quesenberry_hurst)
confint(m, method = :quesenberry_hurst, tail=:left)
confint(m, method = :quesenberry_hurst, tail=:right)

confint(m, method = :sison_glaz)
confint(m, method = :sison_glaz, correct=false)
confint(m, method = :sison_glaz, tail=:left)
confint(m, method = :sison_glaz, tail=:right)

@test_throws ArgumentError confint(m, method=:FOO)
@test_throws ArgumentError confint(m, tail=:fox)

#
x=[1,2,3,1,2,3]
y=[1,1,1,2,2,3]

d = counts(x,y,3)

ChisqTest(d)
MultinomialLRT(d)

PowerDivergenceTest(x,y,3)
PowerDivergenceTest(x,y,(1:3,1:3))

ChisqTest(x,y,3)
ChisqTest(x,y,(1:3,1:3))

MultinomialLRT(x,y,3)
MultinomialLRT(x,y,(1:3,1:3))

# Test that large counts don't cause overflow (issue #43)
d = [113997 1031298
     334453 37471]
PowerDivergenceTest(d)
end


text:
function deserialize_image_summary(summary)
    img = summary.image

    value = _deserialize_png(img.encoded_image_string)

    return value
end

function lookahead_deserialize_image_summary(old_tag, old_val, evs::Summary,
                                                 state_old)
    # prepare the default output (when no action is taken)
    result = old_tag, old_val, state_old
    combined_imgs = [old_val]

    state = state_old + 1
    iter_result = iterate(evs, state)
    while iter_result !== nothing
        ((tag, summary), state) = iter_result
        # iteration body
        typ = summary_type(summary)
        typ !== :image && break
        tag != old_tag && break

        push!(combined_imgs, deserialize_image_summary(summary))

        # end iterate
        iter_result = iterate(evs, state)
    end

    if length(combined_imgs) > 1
        new_val = similar(old_val, size(old_val)..., length(combined_imgs))
        for (i, img)=enumerate(combined_imgs)
            new_val[:, :, i] .= img
        end
        result =  old_tag, new_val, state_old + length(combined_imgs)
    end

    return result
end

# Ugly-ugly hack to work around
# https://github.com/JuliaIO/FileIO.jl/issues/174
# Essentially, ImageMagick.load() cannot handle an in-memory stream
# but only files. Instead, ImageMagick.readblob() should be used.
# On MacOs it's not a problem thanks to QuartzImageIO
if !Sys.isapple()
    did_patch_imagemagick = false

    function _deserialize_png(data)
        if !did_patch_imagemagick
            global did_patch_imagemagick = true

            # First load ImageMagick
            FileIO.checked_import(:ImageMagick)
        end
        # attempt to mimick ImageMagick type reconstruction
        attributes = PNGImage.read_info(data)
        T = PNGImage.png_color_T(attributes)

        # This errors becase we are loading ImageMagick after defining func.
        # data = FileIO.ImageMagick.readblob(data)
        # ERROR: MethodError: no method matching readblob(::Array{UInt8,1})
        # The applicable method may be too new
        _data = Base.invokelatest(FileIO.ImageMagick.readblob, data)
        return reinterpret(T, _data)
    end
else
    _deserialize_png(data) = load(Stream(format"PNG", IOBuffer(data)))
end


text:
using Luxor, Random, Colors

Random.seed!(42)

using Test

function test_arc_sagitta(fname)
    Drawing("A4", fname)
    origin()
    background(Colors.RGB(1, 1, 1))
    setopacity(0.6)

    translate(0, -150)
    p1 = O - (150, 0)
    p2 = O + (150, 0)
    @test_throws ErrorException arc2sagitta(p1, p2, 0)
    @test_throws ErrorException carc2sagitta(p1, p2, 0)

    sethue("red")
    for s in 1:4:150
        @test isarcclockwise(O, p1, p2) == false
        cp, r = arc2sagitta(p1, p2, s, :stroke)
        @test r >= 150
    end
    sethue("green")
    for s in 1:4:150
        cp, r = carc2sagitta(p1, p2, s, :stroke)
        @test r >= 150
    end

    translate(0, 300)

    p1 = O - (0, 150)
    p2 = O + (0, 150)
    @test_throws ErrorException arc2sagitta(p1, p2, 0)
    @test_throws ErrorException carc2sagitta(p1, p2, 0)

    sethue("red")
    for s in 1:4:150
        # the arc is neither clockwise or counterclockwise
        # because points are collinear!
        @test isarcclockwise(O, p2, p1) == false
        @test isarcclockwise(O, p1, p2) == false
        cp, r = arc2sagitta(p1, p2, s, :stroke)
        @test r >= 150
    end
    sethue("green")
    for s in 1:4:150
        cp, r = carc2sagitta(p1, p2, s, :stroke)
        @test r >= 150
    end

    @test finish() == true
end

fname = "arc-sagitta.png"
test_arc_sagitta(fname)
println("...finished arc_sagitta test: output in $(fname)")


text:
module WeberDAQmx
using Weber
using PyCall
import Weber: record, setup

function InitExtension(args...;keys...)
  daq_extension(args...;keys...)
end

type DAQmx <: Weber.Extension
  port::String
  codes::Dict{String,Int}
  reserved::Array{Int}
  wait_len::Float64
  N::Int
  bitarray::BitArray
  pyarray::PyArray
  task::PyObject
end

type EmptyDAQmx <: Weber.Extension
end

const missing_library_error = "Location of niDAQmx library and include file unknown"
function nidaq_missing(e::PyCall.PyError)
  :message in keys(e.val) && startswith(String(e.val[:message]),missing_library_error)
end
nidaq_missing(e) = false

isdaq_error(e::PyCall.PyError) = e.T == pyDAQmx[:DAQmxFunctions][:DAQError]
isdaq_error(e) = false
daq_message(e::PyCall.PyError) = e.val[:mess]

daq_extension(::Void;codes=nothing,eeg_sample_rate=nothing) = EmptyDAQmx()
function daq_extension(port::String;codes=Dict{String,Int}(),
                       eeg_sample_rate=nothing)
  @assert eeg_sample_rate != nothing
  # code clearing waits just a little longer than one eeg sample
  wait_len = 1.25/eeg_sample_rate

  # load pyDAQmx library
  try
    global const pyDAQmx = pyimport_conda("PyDAQmx","PyDAQmx","haberdashPI")
    global const group_by_channel = pyDAQmx[:DAQmx_Val_GroupByChannel]
    global const np = pyimport("numpy")
  catch e
    if nidaq_missing(e)
      error("Could not find NI-DAQ library. Make sure it is installed.")
    else rethrow(e)
    end
  end

  # setup the task
  task = PyObject(nothing)
  try
    task = pyDAQmx[:Task]()
    task[:CreateDOChan](port,"",pyDAQmx[:DAQmx_Val_ChanForAllLines])
    task[:StartTask]()
    atexit(() -> task[:StopTask]())
  catch e
    if isdaq_error(e)
	  error("NI-DAQ Serial port error: "*daq_message(e)*
		    "\n Python Stacktrace:\n"*string(e))
    else rethrow(e) end
  end

  bitarray = BitArray(8)
  pyarray = pycall(np[:zeros],PyArray,8,dtype="uint8")
  reserved = values(codes) |> collect
  push!(reserved,0x00) # 0 is used to clear the line
  DAQmx(port,codes,reserved,wait_len,0,bitarray,pyarray,task)
end

# any other codes are automatically selected by calling `DAQcode(code)`
function daq_code(daq::DAQmx,str::String)
  get!(daq.codes,str) do
    daq.N += 1
    while daq.N in daq.reserved
      daq.N += 1
    end

    if daq.N > 0xff
      error("The code \"$str\" is the 256th code, but only 255 are allowed.")
    end

    return UInt8(daq.N)
  end
end

# DAQwrite sends a code to the NI-DAQ interface

function daq_write(daq::DAQmx,str::String)
  try
    daq.bitarray.chunks[1] = daq_code(daq,str)
    daq.pyarray[:] = daq.bitarray
    daq.task[:WriteDigitalLines](1,1,10.0,group_by_channel,
                                 daq.pyarray,nothing,nothing)
    daq.pyarray[:] = 0
    sleep(daq.wait_len)
    daq.task[:WriteDigitalLines](1,1,10.0,group_by_channel,
                                 daq.pyarray,nothing,nothing)
    sleep(daq.wait_len)
    reinterpret(Int,daq.bitarray.chunks[1])
  catch e
    if isdaq_error(e)
	  error("NI-DAQ Serial port error: "*daq_message(e)*
		    "\n Python Stacktrace:\n"*string(e))
    else rethrow(e) end
  end
end

function setup(fn::Function,e::Weber.ExtendedExperiment{DAQmx})
  setup(next(e)) do
    addcolumn(:daq_code)
    fn()
  end
end

function record(e::ExtendedExperiment{DAQmx},code;kwds...)
  record(next(e),code;daq_code = daq_write(extension(e),code),kwds...)
end

end


text:
"""
comp prob
"""

export fteval!, fteval, ftevalgen!

using DD

"""
eval
"""

function fteval(f::AbstractFaultTreeNode, env::Dict{Symbol,Tx}) where Tx
    top, = bdd(f)
    return fteval(top, env)
end

function fteval(f::AbstractDDNode{Tv,Ti}, env::Dict{Symbol,Tx}) where {Tv,Ti,Tx}
    cache = Dict{AbstractDDNode{Tv,Ti},Tx}()
    return fteval!(f, env, cache)
end

function fteval!(f::DDVariable{Tv,Ti,2}, env::Dict{Symbol,Tx}, cache::Dict{AbstractDDNode{Tv,Ti},Tx}) where {Tv,Ti,Tx}
    get(cache, f) do
        prob = [1-env[f.label], env[f.label]]
        cache[f] = sum([x[1]*fteval!(x[2], env, cache) for x = zip(prob, f.nodes) if !iszero(x[1])])
    end
end

function fteval!(f::DDValue{Tv,Ti}, env::Dict{Symbol,Tx}, cache::Dict{AbstractDDNode{Tv,Ti},Tx}) where {Tv,Ti,Tx}
    (f.val == Tv(0)) ? Tx(0) : Tx(1)
end

"""
deriv1
"""

function fteval(f::AbstractFaultTreeNode, env::Dict{Symbol,Tx}, denv::Dict{Symbol,Tx}) where Tx
    top, = bdd(f)
    return fteval(top, env, denv)
end

function fteval(f::AbstractDDNode{Tv,Ti}, env::Dict{Symbol,Tx}, denv::Dict{Symbol,Tx}) where {Tv,Ti,Tx}
    cache = Dict{AbstractDDNode{Tv,Ti},Tx}()
    dcache = Dict{AbstractDDNode{Tv,Ti},Tx}()
    return fteval!(f, env, denv, cache, dcache)
end

function fteval!(f::DDVariable{Tv,Ti,2}, env::Dict{Symbol,Tx}, denv::Dict{Symbol,Tx},
    cache::Dict{AbstractDDNode{Tv,Ti},Tx}, dcache::Dict{AbstractDDNode{Tv,Ti},Tx}) where {Tv,Ti,Tx}
    get(dcache, f) do
        prob = [1-env[f.label], env[f.label]]
        dprob = [-denv[f.label], denv[f.label]]
        v = [x[1] * fteval!(x[2], env, cache) for x = zip(dprob, f.nodes) if !iszero(x[1])]
        res1 = iszero(v) ? Tx(0) : sum(v)
        res2 = sum([x[1] * fteval!(x[2], env, denv, cache, dcache) for x = zip(prob, f.nodes) if !iszero(x[1])])
        dcache[f] = res1 + res2
    end
end

function fteval!(f::DDValue{Tv,Ti}, env::Dict{Symbol,Tx}, denv::Dict{Symbol,Tx},
    cache::Dict{AbstractDDNode{Tv,Ti},Tx}, dcache::Dict{AbstractDDNode{Tv,Ti},Tx}) where {Tv,Ti,Tx}
    Tx(0)
end

"""
ftevalgen
"""

function ftevalgen!(f::DDVariable{Tv,Ti,2}, env0::Dict{Symbol,Array{Tx,N1}}, env1::Dict{Symbol,Array{Tx,N1}},
    cache::Dict{AbstractDDNode{Tv,Ti},Array{Tx,N2}}) where {Tv,Ti,Tx,N1,N2}
    get(cache, f) do
        v1 = env0[f.label] * ftevalgen!(f.nodes[1], env0, env1, cache)
        v2 = env1[f.label] * ftevalgen!(f.nodes[2], env0, env1, cache)
        cache[f] = v1 + v2
    end
end

function ftevalgen!(f::DDValue{Tv,Ti}, env0::Dict{Symbol,Array{Tx,N1}}, env1::Dict{Symbol,Array{Tx,N1}},
    cache::Dict{AbstractDDNode{Tv,Ti},Array{Tx,N2}}) where {Tv,Ti,Tx,N1,N2}
    cache[f]
end