text:
using LinearAlgebra, Statistics, SpecialFunctions
using ColorVectorSpace, Colors, FixedPointNumbers

using Test

n8sum(x,y) = Float64(N0f8(x)) + Float64(N0f8(y))

macro test_colortype_approx_eq(a, b)
    :(test_colortype_approx_eq($(esc(a)), $(esc(b)), $(string(a)), $(string(b))))
end

function test_colortype_approx_eq(a::Colorant, b::Colorant, astr, bstr)
    @test typeof(a) == typeof(b)
    n = length(fieldnames(typeof(a)))
    for i = 1:n
        @test getfield(a, i) ≈ getfield(b,i)
    end
end

struct RatRGB <: AbstractRGB{Rational{Int}}
    r::Rational{Int}
    g::Rational{Int}
    b::Rational{Int}
end
ColorTypes.red(c::RatRGB)   = c.r
ColorTypes.green(c::RatRGB) = c.g
ColorTypes.blue(c::RatRGB)  = c.b

struct RGBA32 <: AbstractRGBA{RGB24, N0f8}
    color::UInt32
    RGBA32(c::UInt32, ::Type{Val{true}}) = new(c)
end
function RGBA32(r, g, b, alpha=1N0f8)
    u32 = reinterpret(UInt32, ARGB32(r, g, b, alpha))
    RGBA32((u32 << 0x8) | (u32 >> 0x18), Val{true})
end
ColorTypes.red(  c::RGBA32) = reinterpret(N0f8, (c.color >> 0x18) % UInt8)
ColorTypes.green(c::RGBA32) = reinterpret(N0f8, (c.color >> 0x10) % UInt8)
ColorTypes.blue( c::RGBA32) = reinterpret(N0f8, (c.color >> 0x08) % UInt8)
ColorTypes.alpha(c::RGBA32) = reinterpret(N0f8, c.color % UInt8)
ColorTypes.comp4(c::RGBA32) = alpha(c)

struct GrayA32 <: AbstractGrayA{Gray24, N0f8}
    color::UInt32
    GrayA32(c::UInt32, ::Type{Val{true}}) = new(c)
end
GrayA32(g, alpha=1N0f8) = GrayA32(bswap(reinterpret(UInt32, AGray32(g, alpha))), Val{true})
ColorTypes.gray( c::GrayA32) = reinterpret(N0f8, (c.color >> 0x18) % UInt8)
ColorTypes.alpha(c::GrayA32) = reinterpret(N0f8, c.color % UInt8)
ColorTypes.comp2(c::RGBA32) = alpha(c)

@testset "Colortypes" begin
    @testset "ambiguities" begin
        @test isempty(detect_ambiguities(ColorVectorSpace))
    end

    @testset "convert" begin
        for x in (0.5, 0.5f0, NaN, NaN32, N0f8(0.5))
            @test @inferred(convert(Gray{typeof(x)}, x))  === @inferred(convert(Gray, x))  === Gray(x)
            @test @inferred(convert(RGB{typeof(x)}, x))   === @inferred(convert(RGB, x))   === RGB(x, x, x)
            # These should be fixed by a future release of ColorTypes
            @test_broken @inferred(convert(AGray{typeof(x)}, x)) === @inferred(convert(AGray, x)) === AGray(x, 1)
            @test_broken @inferred(convert(ARGB{typeof(x)}, x))  === @inferred(convert(ARGB, x))  === ARGB(x, x, x, 1)
            @test_broken @inferred(convert(GrayA{typeof(x)}, x)) === @inferred(convert(GrayA, x)) === GrayA(x, 1)
            @test_broken @inferred(convert(RGBA{typeof(x)}, x))  === @inferred(convert(RGBA, x))  === RGBA(x, x, x, 1)
        end
    end

    @testset "nan" begin
        function make_checked_nan(::Type{T}) where T
            x = nan(T)
            isa(x, T) && mapreducec(isnan, &, true, x)
        end
        for S in (Float32, Float64)
            @test make_checked_nan(S)
            @test make_checked_nan(Gray{S})
            @test make_checked_nan(AGray{S})
            @test make_checked_nan(GrayA{S})
            @test make_checked_nan(RGB{S})
            @test make_checked_nan(ARGB{S})
            @test make_checked_nan(RGBA{S})
        end
    end

    @testset "traits" begin
        @test floattype(Gray{N0f8}) === Gray{float(N0f8)}
    end

    @testset "Arithmetic with Gray" begin
        cf = Gray{Float32}(0.1)
        @test @inferred(+cf) === cf
        @test @inferred(-cf) === Gray(-0.1f0)
        @test @inferred(one(cf)*cf) === cf
        @test oneunit(cf) === Gray(1.0f0)
        ccmp = Gray{Float32}(0.2)
        @test @inferred(2*cf) === cf*2 === 2.0f0*cf === cf*2.0f0 === ccmp
        @test @inferred(ccmp/2) === cf
        @test @inferred(cf*cf) === Gray{Float32}(0.1f0*0.1f0)
        @test @inferred(Gray{N0f32}(0.5)*Gray(0.5f0)) === Gray(Float64(N0f32(0.5)) * 0.5)
        @test @inferred(cf^2 ) === Gray{Float32}(0.1f0*0.1f0)
        @test @inferred(cf^3.0f0) === Gray{Float32}(0.1f0^3.0f0)
        @test @inferred(2.0*cf) === cf*2.0 === Gray(2.0*0.1f0)
        cf64 = Gray(0.2)
        @test cf / cf64 === Gray(0.1f0/0.2)
        @test_throws MethodError cf ÷ cf
        @test cf + 0.1     === 0.1 + cf        === Gray(Float64(0.1f0) + 0.1)
        @test cf64 - 0.1f0 === -(0.1f0 - cf64) === Gray( 0.2 - Float64(0.1f0))
        @test_throws MethodError abs2(ccmp)
        @test norm(cf) == norm(cf, 2) == norm(gray(cf))
        @test norm(cf, 1)   == norm(gray(cf), 1)
        @test norm(cf, Inf) == norm(gray(cf), Inf)
        @test @inferred(abs(cf)) === Gray(0.1f0)
        cu = Gray{N0f8}(0.1)
        @test @inferred(2*cu) === cu*2 === Gray(2*gray(cu))
        @test @inferred(2.0f0*cu) === cu*2.0f0 === Gray(2.0f0*gray(cu))
        f = N0f8(0.5)
        @test @inferred(gray(f*cu)) === gray(cu*f) ===f*gray(cu)
        @test @inferred(cf/2.0f0) === Gray{Float32}(0.05)
        @test @inferred(cu/2) === Gray(cu.val/2)
        @test @inferred(cu/0.5f0) === Gray(cu.val/0.5f0)
        @test @inferred(cf+cf) === ccmp
        @test isfinite(cf)
        @test isfinite(Gray(true))
        @test !isinf(cf)
        @test !isinf(Gray(f))
        @test !isnan(cf)
        @test !isfinite(Gray(NaN))
        @test !isinf(Gray(NaN))
        @test isnan(Gray(NaN))
        @test !isfinite(Gray(Inf))
        @test Gray(Inf) |> isinf
        @test !isnan(Gray(Inf))
        @test abs(Gray(0.1)) === Gray(0.1)
        @test eps(Gray{N0f8}) === Gray(eps(N0f8))  # #282
        @test atan(Gray(0.1), Gray(0.2)) == atan(0.1, 0.2)
        @test hypot(Gray(0.2), Gray(0.3)) === hypot(0.2, 0.3)
        # Multiplication
        @test cf ⋅ cf   === gray(cf)^2
        @test cf ⋅ cf64 === gray(cf) * gray(cf64)
        @test cf ⊙ cf   === Gray(gray(cf)^2)
        @test cf ⊙ cf64 === Gray(gray(cf) * gray(cf64))
        @test cf ⊗ cf   === Gray(gray(cf)^2)
        @test cf ⊗ cf64 === Gray(gray(cf) * gray(cf64))

        acu = Gray{N0f8}[cu]
        acf = Gray{Float32}[cf]
        @test @inferred(acu./trues(1)) == acu
        @test typeof(acu./trues(1)) == Vector{typeof(cu/true)}
        @test @inferred(ones(Int, 1)./acu) == [1/cu]
        @test typeof(ones(Int, 1)./acu) == Vector{typeof(1/cu)}
        @test @inferred(acu./acu) == [1]
        @test typeof(acu./acu) == Vector{typeof(cu/cu)}
        @test typeof(acu+acf) == Vector{Gray{Float32}}
        @test typeof(acu-acf) == Vector{Gray{Float32}}
        @test typeof(acu.+acf) == Vector{Gray{Float32}}
        @test typeof(acu.-acf) == Vector{Gray{Float32}}
        @test typeof(acu.+cf) == Vector{Gray{Float32}}
        @test typeof(acu.-cf) == Vector{Gray{Float32}}
        @test typeof(2*acf) == Vector{Gray{Float32}}
        @test typeof(2 .* acf) == Vector{Gray{Float32}}
        @test typeof(0x02*acu) == Vector{Gray{Float32}}
        @test typeof(acu/2) == Vector{Gray{typeof(N0f8(0.5)/2)}}
        @test typeof(acf.^2) == Vector{Gray{Float32}}
        @test (acu./Gray{N0f8}(0.5))[1] == gray(acu[1])/N0f8(0.5)
        @test (acf./Gray{Float32}(2))[1] ≈ 0.05f0
        @test (acu/2)[1] == Gray(gray(acu[1])/2)
        @test (acf/2)[1] ≈ Gray{Float32}(0.05f0)

        @test gray(0.8) === 0.8

        a = Gray{N0f8}[0.8,0.7]
        @test a == a
        @test a === a
        @test isapprox(a, a)
        @test sum(a) == Gray(n8sum(0.8,0.7))
        @test sum(a[1:1]) == a[1]
        @test abs( varmult(*, a) - (a[1]-a[2])^2 / 2 ) <= 0.001

        @test real(Gray{Float32}) <: Real
        @test zero(ColorTypes.Gray) == 0
        @test oneunit(ColorTypes.Gray) == 1

        @test typeof(float(Gray{N0f16}(0.5))) <: AbstractFloat
        @test quantile( Gray{N0f16}[0.0,0.5,1.0], 0.1) ≈ 0.1 atol=eps(N0f16)
        @test middle(Gray(0.2)) === Gray(0.2)
        @test middle(Gray(0.2), Gray(0.4)) === Gray((0.2+0.4)/2)

        # issue #56
        @test Gray24(0.8)*N0f8(0.5) === Gray{N0f8}(0.4)
        @test Gray24(0.8)*0.5 === Gray(0.4)
        @test Gray24(0.8)/2   === Gray(0.5f0*N0f8(0.8))
        @test Gray24(0.8)/2.0 === Gray(0.4)

        # issue #133
        @test Gray24(0.2) + Gray24(0.4) === Gray24(0.6)
        @test Gray24(1)   - Gray24(0.2) === Gray24(0.8)
        @test Gray24(1)   * Gray24(0.2) === Gray24(0.2)
    end

    @testset "Comparisons with Gray" begin
        g1 = Gray{N0f8}(0.2)
        g2 = Gray{N0f8}(0.3)
        @test isless(g1, g2)
        @test !(isless(g2, g1))
        @test g1 < g2
        @test !(g2 < g1)
        @test isless(g1, 0.5)
        @test !(isless(0.5, g1))
        @test g1 < 0.5
        @test !(0.5 < g1)
        @test @inferred(max(g1, g2)) === g2
        @test @inferred(max(g1, Gray(0.3))) === Gray(0.3)
        @test max(g1, 0.1) === max(0.1, g1) === Float64(gray(g1))
        @test (@inferred(min(g1, g2)) ) == g1
        @test min(g1, 0.1) === min(0.1, g1) === 0.1
        a = Gray{Float64}(0.9999999999999999)
        b = Gray{Float64}(1.0)

        @test (Gray(0.3) < Gray(NaN)) == (0.3 < NaN)
        @test (Gray(NaN) < Gray(0.3)) == (NaN < 0.3)
        @test isless(Gray(0.3), Gray(NaN)) == isless(0.3, NaN)
        @test isless(Gray(NaN), Gray(0.3)) == isless(NaN, 0.3)
        @test isless(Gray(0.3), NaN) == isless(0.3, NaN)
        @test isless(Gray(NaN), 0.3) == isless(NaN, 0.3)
        @test isless(0.3, Gray(NaN)) == isless(0.3, NaN)
        @test isless(NaN, Gray(0.3)) == isless(NaN, 0.3)

        @test isapprox(a, b)
        a = Gray{Float64}(0.99)
        @test !(isapprox(a, b, rtol = 0.01))
        @test isapprox(a, b, rtol = 0.1)
    end

    @testset "Unary operations with Gray" begin
        ntested = 0
        for g in (Gray(0.4), Gray{N0f8}(0.4))
            @test @inferred(zero(g)) === typeof(g)(0)
            @test @inferred(oneunit(g)) === typeof(g)(1)
            for opgroup in (ColorVectorSpace.unaryOps, (:trunc, :floor, :round, :ceil, :eps, :bswap))
                for op in opgroup
                    op ∈ (:frexp, :exponent, :modf, :lfact) && continue
                    op === :~ && eltype(g) === Float64 && continue
                    op === :significand && eltype(g) === N0f8 && continue
                    try
                        v = @eval $op(gray($g))  # if this fails, don't bother with the next test
                        @test @eval($op($g)) === Gray(v)
                        ntested += 1
                    catch ex
                        @test ex isa Union{DomainError,ArgumentError}
                    end
                end
            end
        end
        @test ntested > 130
        @test logabsgamma(Gray(0.2)) == (Gray(logabsgamma(0.2)[1]), 1)
        for g in (Gray{N0f8}(0.4), Gray{N0f8}(0.6))
            for op in (:trunc, :floor, :round, :ceil)
                v = @eval $op(Bool, gray($g))
                @test @eval($op(Bool, $g)) === Gray(v)
            end
        end
        for (g1, g2) in ((Gray(0.4), Gray(0.3)), (Gray(N0f8(0.4)), Gray(N0f8(0.3))))
            for op in (:mod, :rem, :mod1)
                v = @eval $op(gray($g1), gray($g2))
                @test @eval($op($g1, $g2)) === Gray(v)
            end
        end
        u = N0f8(0.4)
        @test ~Gray(u) == Gray(~u)
        @test -Gray(u) == Gray(-u)
    end

    @testset "Arithmetic with TransparentGray" begin
        p1 = GrayA{Float32}(Gray(0.8), 0.2)
        @test @inferred(zero(p1)) === GrayA{Float32}(0,0)
        @test @inferred(oneunit(p1)) === GrayA{Float32}(1,1)
        @test +p1 == p1
        @test -p1 == GrayA(-0.8f0, -0.2f0)
        p2 = GrayA{Float32}(Gray(0.6), 0.3)
        @test_colortype_approx_eq p1+p2 GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq (p1+p2)/2 GrayA{Float32}(Gray(0.7),0.25)
        @test_colortype_approx_eq 0.4f0*p1+0.6f0*p2 GrayA{Float32}(Gray(0.68),0.26)
        @test_colortype_approx_eq ([p1]+[p2])[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1].+[p2])[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1].+p2)[1] GrayA{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1]-[p2])[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1].-[p2])[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1].-p2)[1] GrayA{Float32}(Gray(0.2),-0.1)
        @test_colortype_approx_eq ([p1]/2)[1] GrayA{Float32}(Gray(0.4),0.1)
        @test_colortype_approx_eq (0.4f0*[p1]+0.6f0*[p2])[1] GrayA{Float32}(Gray(0.68),0.26)

        a = GrayA{N0f8}[GrayA(0.8,0.7), GrayA(0.5,0.2)]
        @test sum(a) == GrayA(n8sum(0.8,0.5), n8sum(0.7,0.2))
        @test isapprox(a, a)
        a = AGray{Float64}(1.0, 0.9999999999999999)
        b = AGray{Float64}(1.0, 1.0)

        @test a ≈ b
        a = AGray{Float64}(1.0, 0.99)
        @test !isapprox(a, b, rtol = 0.01)
        @test isapprox(a, b, rtol = 0.1)

        # issue #56
        @test AGray32(0.8,0.2)*N0f8(0.5) === AGray{N0f8}(0.4,0.1)
        @test AGray32(0.8,0.2)*0.5 === AGray(0.4,0.1)
        @test AGray32(0.8,0.2)/2   === AGray(0.5f0*N0f8(0.8),0.5f0*N0f8(0.2))
        @test AGray32(0.8,0.2)/2.0 === AGray(0.4,0.1)

        # issue #133
        @test AGray32(1, 0.4) - AGray32(0.2, 0.2) === AGray32(0.8, 0.2)

        # issue #146
        @test @inferred(GrayA32(0.8,0.2)*N0f8(0.5)) === GrayA{N0f8}(0.4,0.1)
        @test @inferred(GrayA32(0.8,0.2)*0.5) === GrayA(0.4,0.1)
        @test @inferred(GrayA32(0.8,0.2)/2)   === GrayA(0.5f0*N0f8(0.8),0.5f0*N0f8(0.2))
        @test @inferred(GrayA32(0.8,0.2)/2.0) === GrayA(0.4,0.1)
        @test @inferred(GrayA32(1, 0.4) - GrayA32(0.2, 0.2)) === GrayA32(0.8, 0.2)

        # Multiplication
        cf = AGray{Float32}(0.8, 0.2)
        @test_throws MethodError cf * cf
        @test_throws MethodError cf ⋅ cf
        @test_throws MethodError cf ⊗ cf
        cf64 = mapc(Float64, cf)
        @test @inferred(cf ⊙ cf)   === AGray{Float32}(0.8f0^2, 0.2f0^2)
        @test @inferred(cf ⊙ cf64) === AGray{Float64}(0.8f0*(0.8f0*1.0), 0.2f0*(0.2f0*1.0))
    end

    @testset "Arithemtic with RGB" begin
        cf = RGB{Float32}(0.1,0.2,0.3)
        @test @inferred(zero(cf)) === RGB{Float32}(0,0,0)
        @test @inferred(oneunit(cf)) === RGB{Float32}(1,1,1)
        @test +cf == cf
        @test -cf == RGB(-0.1f0, -0.2f0, -0.3f0)
        ccmp = RGB{Float32}(0.2,0.4,0.6)
        @test 2*cf == ccmp
        @test cf*2 == ccmp
        @test ccmp/2 == cf
        @test 2.0f0*cf == ccmp
        @test eltype(2.0*cf) == Float64
        cu = RGB{N0f8}(0.1,0.2,0.3)
        @test_colortype_approx_eq 2*cu RGB(2*cu.r, 2*cu.g, 2*cu.b)
        @test_colortype_approx_eq 2.0f0*cu RGB(2.0f0*cu.r, 2.0f0*cu.g, 2.0f0*cu.b)
        f = N0f8(0.5)
        @test (f*cu).r ≈ f*cu.r
        @test cf/2.0f0 == RGB{Float32}(0.05,0.1,0.15)
        @test cu/2 ≈ RGB(cu.r/2,cu.g/2,cu.b/2)
        @test cu/0.5f0 ≈ RGB(cu.r/0.5f0, cu.g/0.5f0, cu.b/0.5f0)
        @test cf+cf == ccmp
        @test cu * 1//2 == mapc(x->Float64(Rational(x)/2), cu)
        @test_colortype_approx_eq (cf.*[0.8f0])[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq ([0.8f0].*cf)[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test isfinite(cf)
        @test !isinf(cf)
        @test !isnan(cf)
        @test !isfinite(RGB(NaN, 1, 0.5))
        @test !isinf(RGB(NaN, 1, 0.5))
        @test isnan(RGB(NaN, 1, 0.5))
        @test !isfinite(RGB(1, Inf, 0.5))
        @test isinf(RGB(1, Inf, 0.5))
        @test !isnan(RGB(1, Inf, 0.5))
        @test abs(RGB(0.1,0.2,0.3)) == RGB(0.1,0.2,0.3)
        @test_throws MethodError abs2(RGB(0.1,0.2,0.3))
        @test_throws MethodError sum(abs2, RGB(0.1,0.2,0.3))
        @test norm(RGB(0.1,0.2,0.3)) ≈ sqrt(0.14)/sqrt(3)

        acu = RGB{N0f8}[cu]
        acf = RGB{Float32}[cf]
        @test typeof(acu+acf) == Vector{RGB{Float32}}
        @test typeof(acu-acf) == Vector{RGB{Float32}}
        @test typeof(acu.+acf) == Vector{RGB{Float32}}
        @test typeof(acu.-acf) == Vector{RGB{Float32}}
        @test typeof(acu.+cf) == Vector{RGB{Float32}}
        @test typeof(acu.-cf) == Vector{RGB{Float32}}
        @test typeof(2*acf) == Vector{RGB{Float32}}
        @test typeof(convert(UInt8, 2)*acu) == Vector{RGB{Float32}}
        @test typeof(acu/2) == Vector{RGB{typeof(N0f8(0.5)/2)}}
        rcu = rand(RGB{N0f8}, 3, 5)
        @test @inferred(rcu./trues(3, 5)) == rcu
        @test typeof(rcu./trues(3, 5)) == Matrix{typeof(cu/true)}

        a = RGB{N0f8}[RGB(1,0,0), RGB(1,0.8,0)]
        @test sum(a) == RGB(2.0,0.8,0)
        @test sum(typeof(a)()) == RGB(0.0,0.0,0)
        @test isapprox(a, a)
        a = RGB{Float64}(1.0, 1.0, 0.9999999999999999)
        b = RGB{Float64}(1.0, 1.0, 1.0)

        @test isapprox(a, b)
        a = RGB{Float64}(1.0, 1.0, 0.99)
        @test !(isapprox(a, b, rtol = 0.01))
        @test isapprox(a, b, rtol = 0.1)
        # issue #56
        @test RGB24(1,0,0)*N0f8(0.5) === RGB{N0f8}(0.5,0,0)
        @test RGB24(1,0,0)*0.5 === RGB(0.5,0,0)
        @test RGB24(1,0,0)/2   === RGB(0.5f0,0,0)
        @test RGB24(1,0,0)/2.0 === RGB(0.5,0,0)
        # issue #133
        @test RGB24(1, 0, 0) + RGB24(0, 0, 1) === RGB24(1, 0, 1)

        # Multiplication
        @test_throws MethodError cf*cf
        cf64 = mapc(Float64, cf)
        @test cf ⋅ cf   === (red(cf)^2 + green(cf)^2 + blue(cf)^2)/3
        @test cf ⋅ cf64 === (red(cf)*red(cf64) + green(cf)*green(cf64) + blue(cf)*blue(cf64))/3
        @test cf ⊙ cf   === RGB(red(cf)^2, green(cf)^2, blue(cf)^2)
        @test cf ⊙ cf64 === RGB(red(cf)*red(cf64), green(cf)*green(cf64), blue(cf)*blue(cf64))
        c2 = rand(RGB{Float64})
        rr = cf ⊗ c2
        matrix_rr = Matrix(rr)
        @test matrix_rr  == [red(cf)*red(c2)   red(cf)*green(c2)   red(cf)*blue(c2);
                             green(cf)*red(c2) green(cf)*green(c2) green(cf)*blue(c2);
                             blue(cf)*red(c2)  blue(cf)*green(c2)  blue(cf)*blue(c2)]
        @test RGBRGB(matrix_rr) == rr
        @test RGBRGB{N0f8}(zeros(3, 3)) == zero(RGBRGB{N0f8})
        @test_throws MethodError RGBRGB(rand(9))
        @test_throws DimensionMismatch RGBRGB(rand(9, 1))
        @test +rr === rr
        @test -rr === RGBRGB(-rr.rr, -rr.gr, -rr.br, -rr.rg, -rr.gg, -rr.bg, -rr.rb, -rr.gb, -rr.bb)
        @test rr + rr == 2*rr == rr*2
        @test rr - rr == zero(rr)
        io = IOBuffer()
        print(io, N0f8)
        Tstr = String(take!(io))
        cfn = RGB{N0f8}(0.1, 0.2, 0.3)
        show(io, cfn ⊗ cfn)
        nsuf = string(0.0N0f8)[4:end] # FixedPointNumbers issue #241
        @test String(take!(io)) == "RGBRGB{$Tstr}([0.012$nsuf 0.02$nsuf 0.031$nsuf; " *
                                                  "0.02$nsuf 0.039$nsuf 0.059$nsuf; " *
                                                  "0.031$nsuf 0.059$nsuf 0.09$nsuf])"
        show(io, "text/plain", cfn ⊗ cfn)
        spstr = Base.VERSION >= v"1.5" ? "" : " "
        @test String(take!(io)) == "RGBRGB{$Tstr}:\n 0.012  0.02   0.031\n 0.02   0.039  0.059\n 0.031  0.059  0.09$spstr"
        show(IOContext(io, :compact => false), "text/plain", cf64 ⊗ cf64)
        regex = Regex(raw"^RGBRGB{Float64}:\n" *
                      raw" 0.010\d+\s+0.020\d+\s+0.030\d+\s*\n" *
                      raw" 0.020\d+\s+0.040\d+\s+0.060\d+\s*\n" *
                      raw" 0.030\d+\s+0.060\d+\s+0.090\d+\s*$")
        @test occursin(regex, String(take!(io)))
    end

    @testset "Arithemtic with TransparentRGB" begin
        cf = RGBA{Float32}(0.1,0.2,0.3,0.4)
        @test @inferred(zero(cf)) === RGBA{Float32}(0,0,0,0)
        @test @inferred(oneunit(cf)) === RGBA{Float32}(1,1,1,1)
        @test +cf == cf
        @test -cf == RGBA(-0.1f0, -0.2f0, -0.3f0, -0.4f0)
        ccmp = RGBA{Float32}(0.2,0.4,0.6,0.8)
        @test 2*cf == ccmp
        @test cf*2 == ccmp
        @test ccmp/2 == cf
        @test 2.0f0*cf == ccmp
        @test eltype(2.0*cf) == Float64
        cu = RGBA{N0f8}(0.1,0.2,0.3,0.4)
        @test_colortype_approx_eq 2*cu RGBA(2*cu.r, 2*cu.g, 2*cu.b, 2*cu.alpha)
        @test_colortype_approx_eq 2.0f0*cu RGBA(2.0f0*cu.r, 2.0f0*cu.g, 2.0f0*cu.b, 2.0f0*cu.alpha)
        f = N0f8(0.5)
        @test (f*cu).r ≈ f*cu.r
        @test cf/2.0f0 == RGBA{Float32}(0.05,0.1,0.15,0.2)
        @test cu/2 == RGBA(cu.r/2,cu.g/2,cu.b/2,cu.alpha/2)
        @test cu/0.5f0 == RGBA(cu.r/0.5f0, cu.g/0.5f0, cu.b/0.5f0, cu.alpha/0.5f0)
        @test cf+cf == ccmp
        @test_colortype_approx_eq (cf.*[0.8f0])[1] RGBA{Float32}(0.8*0.1,0.8*0.2,0.8*0.3,0.8*0.4)
        @test_colortype_approx_eq ([0.8f0].*cf)[1] RGBA{Float32}(0.8*0.1,0.8*0.2,0.8*0.3,0.8*0.4)
        @test isfinite(cf)
        @test !isinf(cf)
        @test !isnan(cf)
        @test isnan(RGBA(NaN, 1, 0.5, 0.8))
        @test !isinf(RGBA(NaN, 1, 0.5))
        @test isnan(RGBA(NaN, 1, 0.5))
        @test !isfinite(RGBA(1, Inf, 0.5))
        @test RGBA(1, Inf, 0.5) |> isinf
        @test !isnan(RGBA(1, Inf, 0.5))
        @test !isfinite(RGBA(0.2, 1, 0.5, NaN))
        @test !isinf(RGBA(0.2, 1, 0.5, NaN))
        @test isnan(RGBA(0.2, 1, 0.5, NaN))
        @test !isfinite(RGBA(0.2, 1, 0.5, Inf))
        @test RGBA(0.2, 1, 0.5, Inf) |> isinf
        @test !isnan(RGBA(0.2, 1, 0.5, Inf))
        @test abs(RGBA(0.1,0.2,0.3,0.2)) === RGBA(0.1,0.2,0.3,0.2)

        acu = RGBA{N0f8}[cu]
        acf = RGBA{Float32}[cf]
        @test typeof(acu+acf) == Vector{RGBA{Float32}}
        @test typeof(acu-acf) == Vector{RGBA{Float32}}
        @test typeof(acu.+acf) == Vector{RGBA{Float32}}
        @test typeof(acu.-acf) == Vector{RGBA{Float32}}
        @test typeof(acu.+cf) == Vector{RGBA{Float32}}
        @test typeof(acu.-cf) == Vector{RGBA{Float32}}
        @test typeof(2*acf) == Vector{RGBA{Float32}}
        @test typeof(convert(UInt8, 2)*acu) == Vector{RGBA{Float32}}
        @test typeof(acu/2) == Vector{RGBA{typeof(N0f8(0.5)/2)}}

        a = RGBA{N0f8}[RGBA(1,0,0,0.8), RGBA(0.7,0.8,0,0.9)]
        @test sum(a) == RGBA(n8sum(1,0.7),0.8,0,n8sum(0.8,0.9))
        @test isapprox(a, a)
        a = ARGB{Float64}(1.0, 1.0, 1.0, 0.9999999999999999)
        b = ARGB{Float64}(1.0, 1.0, 1.0, 1.0)

        @test isapprox(a, b)
        a = ARGB{Float64}(1.0, 1.0, 1.0, 0.99)
        @test !(isapprox(a, b, rtol = 0.01))
        @test isapprox(a, b, rtol = 0.1)
        # issue #56
        @test ARGB32(1,0,0,0.8)*N0f8(0.5) === ARGB{N0f8}(0.5,0,0,0.4)
        @test ARGB32(1,0,0,0.8)*0.5 === ARGB(0.5,0,0,0.4)
        @test ARGB32(1,0,0,0.8)/2   === ARGB(0.5f0,0,0,0.5f0*N0f8(0.8))
        @test ARGB32(1,0,0,0.8)/2.0 === ARGB(0.5,0,0,0.4)
        # issue #133
        @test ARGB32(1, 0, 0, 0.2) + ARGB32(0, 0, 1, 0.2) === ARGB32(1, 0, 1, 0.4)

        # issue #146
        @test @inferred(RGBA32(1,0,0,0.8)*N0f8(0.5)) === RGBA{N0f8}(0.5,0,0,0.4)
        @test @inferred(RGBA32(1,0,0,0.8)*0.5) === RGBA(0.5,0,0,0.4)
        @test @inferred(RGBA32(1,0,0,0.8)/2)   === RGBA(0.5f0,0,0,0.5f0*N0f8(0.8))
        @test @inferred(RGBA32(1,0,0,0.8)/2.0) === RGBA(0.5,0,0,0.4)
        @test @inferred(RGBA32(1, 0, 0, 0.2) + RGBA32(0, 0, 1, 0.2)) === RGBA32(1, 0, 1, 0.4)

        # Multiplication
        @test_throws MethodError cf * cf
        @test_throws MethodError cf ⋅ cf
        @test_throws MethodError cf ⊗ cf
        cf64 = mapc(Float64, cf)
        @test @inferred(cf ⊙ cf)   === RGBA{Float32}(0.1f0^2, 0.2f0^2, 0.3f0^2, 0.4f0^2)
        @test @inferred(cf ⊙ cf64) === RGBA{Float64}(0.1f0*(0.1f0*1.0), 0.2f0*(0.2f0*1.0),
                                                     0.3f0*(0.3f0*1.0), 0.4f0*(0.4f0*1.0))
    end

    @testset "Mixed-type arithmetic" begin
        # issue 155
        @test @inferred(Gray(0.2f0) + Gray24(0.2)) === Gray{Float32}(0.2 + 0.2N0f8)
        @test @inferred(RGBX(0, 0, 1) + XRGB(1, 0, 0)) === XRGB{N0f8}(1, 0, 1)
        @test @inferred(BGR(0, 0, 1) + RGB24(1, 0, 0)) === RGB{N0f8}(1, 0, 1)
        @test_throws Exception HSV(100, 0.2, 0.4) + Gray(0.2)

        @test AGray32(0.2, 0.4) + Gray24(0.2) === AGray32(0.4, 0.4N0f8+1N0f8)
        @test AGray32(0.2, 0.4) + Gray(0.2f0) === AGray{Float32}(0.2+0.2N0f8, 0.4N0f8+1)
        @test RGB(1, 0, 0)      + Gray(0.2f0) === RGB{Float32}(1.2, 0.2, 0.2)
        @test RGB(1, 0, 0)      - Gray(0.2f0) === RGB{Float32}(0.8, -0.2, -0.2)
        @test RGB24(1, 0, 0)    + Gray(0.2f0) === RGB{Float32}(1.2, 0.2, 0.2)
        @test RGB24(1, 0, 0)    - Gray(0.2f0) === RGB{Float32}(0.8, -0.2, -0.2)
        @test RGB(1.0f0, 0, 0)  + Gray24(0.2) === RGB{Float32}(1.2, 0.2, 0.2)
        @test RGB(1.0f0, 0, 0)  - Gray24(0.2) === RGB{Float32}(0.8, -0.2, -0.2)
        @test RGB24(1, 0, 0)    + Gray24(0.2) === RGB24(1N0f8+0.2N0f8, 0.2, 0.2)
        @test RGB24(0.4, 0, 0.2)   + AGray32(0.4, 1)   === ARGB32(0.8, 0.4, 0.6, 1N0f8+1N0f8)
        @test RGB24(0.4, 0.6, 0.5) - AGray32(0.4, 0.2) === ARGB32(0, 0.2, 0.1, 0.8)
        @test ARGB32(0.4, 0, 0.2, 0.5) + Gray24(0.4)   === ARGB32(0.8, 0.4, 0.6, 0.5N0f8+1N0f8)
        @test ARGB32(0.4, 0, 0.2, 0.5) + AGray32(0.4, 0.2) === ARGB32(0.8, 0.4, 0.6, 0.5N0f8+0.2N0f8)
        @test ARGB32(0.4, 0, 0.2, 0.5) + RGB(0.4f0, 0, 0) === ARGB{Float32}(0.4N0f8+0.4, 0, 0.2N0f8, 0.5N0f8+1)

        g, rgb = Gray{Float32}(0.2), RGB{Float64}(0.1, 0.2, 0.3)
        ag, argb = AGray{Float64}(0.2, 0.8), ARGB{Float32}(0.1, 0.2, 0.3, 0.4)
        @test g ⋅ rgb === rgb ⋅ g === 0.2f0*(0.1 + 0.2 + 0.3)/3
        @test_throws MethodError g ⋅ ag
        @test_throws MethodError g ⋅ argb
        @test_throws MethodError ag ⋅ rgb
        @test_throws MethodError ag ⋅ argb
        @test_throws MethodError rgb ⋅ argb
        @test g ⊙ rgb === rgb ⊙ g === RGB{Float64}(0.2f0*0.1, 0.2f0*0.2, 0.2f0*0.3)
        @test g ⊙ ag === ag ⊙ g === AGray{Float64}(0.2f0*0.2, 1.0f0*0.8)
        @test g ⊙ argb === argb ⊙ g === ARGB{Float32}(0.2f0*0.1f0, 0.2f0*0.2f0, 0.2f0*0.3f0, 1.0f0*0.4f0)
        @test ag ⊙ rgb === rgb ⊙ ag === ARGB{Float64}(0.2*0.1, 0.2*0.2, 0.2*0.3, 0.8*1.0)
        @test ag ⊙ argb === argb ⊙ ag === ARGB{Float64}(0.2*0.1f0, 0.2*0.2f0, 0.2*0.3f0, 0.8*0.4f0)
        @test rgb ⊙ argb === argb ⊙ rgb === ARGB{Float64}(0.1*0.1f0, 0.2*0.2f0, 0.3*0.3f0, 1.0*0.4f0)
        @test g ⊗ rgb === RGB(g) ⊗ rgb
        @test rgb ⊗ g === rgb ⊗ RGB(g)
        @test Matrix(g ⊗ rgb) == Matrix(rgb ⊗ g)'
        @test_throws MethodError g ⊗ ag
        @test_throws MethodError g ⊗ argb
        @test_throws MethodError ag ⊗ rgb
        @test_throws MethodError ag ⊗ argb
        @test_throws MethodError rgb ⊗ argb
    end

    @testset "Custom RGB arithmetic" begin # see also the `RGBA32` cases above
        cf = RatRGB(1//10, 2//10, 3//10)
        @test cf ⋅ cf   === (Float64(red(cf))^2 + Float64(green(cf))^2 + Float64(blue(cf))^2)/3
    end

    @testset "arithmetic with Bool" begin
        cb = Gray{Bool}(1)
        @test @inferred(+cb) === cb
        @test @inferred(-cb) === cb # wrapped around
        @test_broken @inferred(one(cb) * cb) === cb
        @test oneunit(cb) === Gray(true)
        # TODO: add more tests (cf. issue #148)
    end

    @testset "Complement" begin
        @test complement(Gray(0.2)) === Gray(0.8)
        @test complement(AGray(0.2f0, 0.7f0)) === AGray(0.8f0, 0.7f0)
        @test complement(GrayA{N0f8}(0.2, 0.7)) === GrayA{N0f8}(0.8, 0.7)
        @test complement(Gray24(0.2)) === Gray24(0.8)
        @test complement(AGray32(0.2, 0.7)) === AGray32(0.8, 0.7)

        @test complement(RGB(0, 0.3, 1)) === RGB(1, 0.7, 0)
        @test complement(ARGB(0, 0.3f0, 1, 0.7f0)) === ARGB(1, 0.7f0, 0, 0.7f0)
        @test complement(RGBA{N0f8}(0, 0.6, 1, 0.7)) === RGBA{N0f8}(1, 0.4, 0.0, 0.7)
        @test complement(RGB24(0, 0.6, 1)) === RGB24(1, 0.4, 0.0)
        @test complement(ARGB32(0, 0.6, 1, 0.7)) === ARGB32(1, 0.4, 0.0, 0.7)
    end

    @testset "dotc" begin
        @test dotc(0.2, 0.2) == 0.2^2
        @test dotc(Int8(3), Int16(6)) === 18
        @test dotc(0.2, 0.3f0) == 0.2*0.3f0
        @test dotc(N0f8(0.2), N0f8(0.3)) == Float32(N0f8(0.2))*Float32(N0f8(0.3))
        @test dotc(Gray{N0f8}(0.2), Gray24(0.3)) == Float32(N0f8(0.2))*Float32(N0f8(0.3))
        xc, yc = RGB(0.2,0.2,0.2), RGB{N0f8}(0.3,0.3,0.3)
        @test isapprox(dotc(xc, yc) , dotc(convert(Gray, xc), convert(Gray, yc)), atol=1e-6)
        @test dotc(RGB(1,0,0), RGB(0,1,1)) == 0
    end

    @testset "typemin/max" begin
        for T in (Normed{UInt8,8}, Normed{UInt8,6}, Normed{UInt16,16}, Normed{UInt16,14}, Float32, Float64)
            @test typemin(Gray{T}) === Gray{T}(typemin(T))
            @test typemax(Gray{T}) === Gray{T}(typemax(T))
            @test typemin(Gray{T}(0.5)) === Gray{T}(typemin(T))
            @test typemax(Gray{T}(0.5)) === Gray{T}(typemax(T))
            A = maximum(Gray{T}.([1 0 0; 0 1 0]); dims=1)  # see PR#44 discussion
            @test isa(A, Matrix{Gray{T}})
            @test size(A) == (1,3)
        end
    end

    @testset "Colors issue #326" begin
        A = rand(RGB{N0f8}, 2, 2)
        @test @inferred(mean(A)) == mean(map(c->mapc(FixedPointNumbers.Treduce, c), A))
    end

    @testset "Equivalence" begin
        x = 0.4
        g = Gray(x)
        c = RGB(g)
        for p in (0, 1, 2, Inf)
            @test norm(x, p) == norm(g, p) ≈ norm(c, p)
        end
        @test dot(x, x) == dot(g, g) ≈ dot(c, c)
        @test_throws MethodError mapreduce(x->x^2, +, c)   # this risks breaking equivalence & noniterability
    end

    @testset "varmult" begin
        cs = [RGB(0.2, 0.3, 0.4), RGB(0.5, 0.3, 0.2)]
        @test varmult(⋅, cs) ≈ 2*(0.15^2 + 0.1^2)/3    # the /3 is for the 3 color channels, i.e., equivalence
        @test varmult(⋅, cs; corrected=false) ≈ (0.15^2 + 0.1^2)/3
        @test varmult(⋅, cs; mean=RGB(0, 0, 0)) ≈ (0.2^2+0.3^2+0.4^2 + 0.5^2+0.3^2+0.2^2)/3
        @test varmult(⊙, cs) ≈ 2*RGB(0.15^2, 0, 0.1^2)
        @test Matrix(varmult(⊗, cs)) ≈ 2*[0.15^2 0 -0.1*0.15; 0 0 0; -0.1*0.15 0 0.1^2]

        cs = [RGB(0.1, 0.2,  0.3)  RGB(0.3, 0.5, 0.3);
              RGB(0.2, 0.21, 0.33) RGB(0.4, 0.51, 0.33);
              RGB(0.3, 0.22, 0.36) RGB(0.5, 0.52, 0.36)]
        v1 = RGB(0.1^2, 0.15^2, 0)
        @test varmult(⊙, cs, dims=2) ≈ 2*[v1, v1, v1]
        v2 = RGB(0.1^2, 0.01^2, 0.03^2)
        @test varmult(⊙, cs, dims=1) ≈ [v2 v2]
    end

    @testset "copy" begin
        g = Gray{N0f8}(0.2)
        @test copy(g) === g
        c = RGB(0.1, 0.2, 0.3)
        @test copy(c) === c
    end

end


text:
const NTypes{N} = NTuple{N, Val}

valueof(::Val{x}) where {x} = x

# Not exactly `Base.aligned_sizeof`
Base.@pure function sizeof_aligned(T::Type)
    if isbitstype(T)
        al = Base.datatype_alignment(T)
        return (Core.sizeof(T) + al - 1) & -al
    else
        return nothing
    end
end

@inline foldlargs(op, x) = x
@inline foldlargs(op, x1, x2, xs...) = foldlargs(op, @return_if_reduced(op(x1, x2)), xs...)

terminating_foldlargs(op, fallback) = fallback()
@inline function terminating_foldlargs(op, fallback::F, x1, x2, xs...) where {F}
    acc = op(x1, x2)
    acc isa Reduced && return unreduced(acc)
    xs isa Tuple{} && return fallback()
    return terminating_foldlargs(op, fallback, acc, xs...)
end

@inline foldrunion(op, ::Type{T}, init) where {T} =
    if T isa Union
        acc = @return_if_reduced foldrunion(op, T.b, init)
        foldrunion(op, T.a, acc)
    else
        op(T, init)
    end

@generated asntypes(::Type{T}) where {T} =
    QuoteNode(foldrunion((S, types) -> (Val(S), types...), T, ()))

struct UnionValue{T <: NTypes,NBytes}
    types::T
    data::NTuple{NBytes,UInt32}
    typeid::UInt8
end

@noinline unreachable() = error("unreachable")

@inline function unionvalue(::Type{T}, v::T) where {T}
    if T isa Union
        types = asntypes(T)
        nbytes = foldrunion(T, 0) do S, n
            Base.@_inline_meta
            max(sizeof_aligned(S), n)
        end
        dest = Ref(ntuple(_ -> UInt32(0), Val(nbytes)))
        GC.@preserve dest begin
            unsafe_store!(Ptr{typeof(v)}(pointer_from_objref(dest)), v)
        end
        @inline function searchid((v, id), t)
            if v isa valueof(t)
                Reduced(id)
            else
                (v, id + 1)
            end
        end
        typeid = terminating_foldlargs(searchid, unreachable, (v, 1), types...)
        return UnionValue(types, dest[], UInt8(typeid))
    else
        return v
    end
end

@noinline invalid_typeid() = error("invalid typeid")

interpret(x) = x
@inline function interpret(uv::UnionValue)
    data = uv.data
    typeid = uv.typeid
    @inline function _get(id, t)
        if id == typeid
            T = valueof(t)
            ref = Ref(data)
            GC.@preserve ref begin
                v = unsafe_load(Ptr{T}(pointer_from_objref(ref)))
            end
            return Reduced(v)
        else
            id + 1
        end
    end
    return terminating_foldlargs(_get, invalid_typeid, 1, uv.types...)
end

@inline function CUDA.shfl_recurse(op, uv::UnionValue)
    data = map(op, uv.data)
    typeid = op(uv.typeid)
    return UnionValue(uv.types, data, typeid)
end

ext: jl


text:
using CGT
using Test

@testset "CGT.jl" begin
    # Write your tests here.
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    GrahamScan

Compute the convex hull of a set of points or geometries using the
Graham's scan method. See [https://en.wikipedia.org/wiki/Graham_scan]
(https://www.pnas.org/content/117/52/33711).

The method has complexity `O(n*log(n))` where `n` is the number of points.

## References

* Cormen et al. 2009. [Introduction to Algorithms]
  (https://mitpress.mit.edu/books/introduction-algorithms-third-edition)
"""
struct GrahamScan <: HullMethod end

function hull(pset::PointSet{2,T}, ::GrahamScan) where {T}
  # remove duplicates
  Q = coordinates.(pset) |> unique

  # sort by y then by x
  sort!(Q, by=reverse)

  # corner cases
  n = length(Q)
  n == 1 && return Point(Q[1])
  n == 2 && return Segment(Point(Q[1]), Point(Q[2]))
  if n == 3
    p₀, p₁, p₂ = Q
    θ = ∠(Point(p₁), Point(p₀), Point(p₂))
    if isapprox(θ, zero(T), atol=atol(T))
      return Segment(Point(p₀), Point(p₂))
    else
      c = Chain(Point(p₀), Point(p₁), Point(p₂), Point(p₀))
      c = orientation(c) == :CCW ? c : reverse(c)
      return PolyArea(c)
    end
  end

  # sort by polar angle
  p₀ = Point(Q[1])
  p  = Point.(Q[2:n])
  x  = p₀ + Vec{2,T}(1, 0)
  θ  = [∠(x, p₀, pᵢ) for pᵢ in p]
  p  = p[sortperm(θ)]

  # rotational sweep
  c = [p₀, p[1], p[2]]
  for pᵢ in p[3:end]
    while ∠(c[end-1], c[end], pᵢ) > atol(T)
      pop!(c)
    end
    push!(c, pᵢ)
  end

  # close chain
  push!(c, c[1])

  PolyArea(c)
end


text:
# Add parallel workers and make package available to workers
using Distributed, CSV, GLM, Distributions, Random, SparseArrays
import HurdleDMR
addprocs(Sys.CPU_THREADS-2)
@everywhere using HurdleDMR

pminmax(x, xmin, xmax) = max(xmin, min(xmax, x))

function bdata(n, p, d; seed=13, qmin=1e-4, qmax=1.0-1e-4)
  Random.seed!(seed)

  λ0 = 500
  m = 1 .+ rand(Poisson(λ0*log(d)),n)
  covars = rand(n,p)
  ηfn(vi) = exp.([0 + λ0/d*log(1+i)*sum(vi) for i=1:d])

  # raw category probabilities
  q = [ηfn(covars[i,:]) for i=1:n]

  # rescale once and chop extremes
  for i=1:n
    q[i] = pminmax.(q[i]/sum(q[i]), qmin, qmax)
  end

  # rescale again to make sure sums to 1
  for i=1:n
    q[i] ./= sum(q[i])
  end

  counts = convert(SparseMatrixCSC{Float64,Int},hcat(broadcast((qi,mi)->rand(Multinomial(mi, qi)),q,m)...)')
  # projdir = p

  covars, counts
end


text:
module SfRocketry

import PhysicalConstants.CODATA2018: c_0, g_n
using Unitful

import ..SfUnits: Angle, to_angle, Acceleration, Speed

export rocket_propulsive_efficiency, brachistochrone_transit_time, brachistochrone_acceleration, brachistochrone_delta_v,
	boost_coast_transit_time, boost_coast_thrust_time, beam_core_mass_ratio, mass_ratio, delta_v, rocket_thrust, rocket_power,
	mass_flow, exhaust_velocity, shkadov_thrust

rocket_propulsive_efficiency(ship_velocity::Speed, exhaust_velocity::Speed) =
	2(ship_velocity / exhaust_velocity) / (1+(ship_velocity / exhaust_velocity)^2)

"""
    brachistochrone_transit_time(d::Unitful.Length, a::Acceleration)

Time to travel a distance of `d` given a constant acceleration of `a` and a final velocity of 0 relative to the start.
"""
brachistochrone_transit_time(d::Unitful.Length, a::Acceleration) = upreferred(2sqrt(d/a))

"""
	brachistochrone_acceleration(d::Unitful.Length, t::Unitful.Time)

Required constant acceleration to travel a distance of `d` in time `t` with a final velocity of 0 relative to the start.
"""
brachistochrone_acceleration(d::Unitful.Length, t::Unitful.Time) = upreferred(4d/t^2)

"""
	brachistochrone_delta_v(d::Unitful.Length, a::Acceleration)

Required ΔV to travel a distance `d` with constant acceleration `a` and attain a final velocity of 0 relative to the start.
"""
brachistochrone_delta_v(d::Unitful.Length, a::Acceleration) = 2sqrt(d*a) |> u"km/s"

"""
    boost_coast_transit_time(d::Unitful.Length, a::Acceleration, boost::Unitful.Time)

Total time to travel a distance of `d`, given speedup and slowdown phases of length `boost` undergoing acceleration `a`.
"""
boost_coast_transit_time(d::Unitful.Length, a::Acceleration, boost::Unitful.Time) = ((d - (a * boost^2)) / (a * boost)) + 2boost

"""
    boost_coast_thrust_time(dist::Unitful.Length, acc::Acceleration, transit::Unitful.Time)

Combined duration of boost and brake phases required to cross a distance of `dist` in time `transit` given an acceleration under thrust of `acc`.
"""
function boost_coast_thrust_time(dist::Unitful.Length, acc::Acceleration, transit::Unitful.Time)
	# Tt - t^2 - d/a
	# ax^2 + bx + c

	a = -1
	b = transit
	c = -dist / acc

	discr = b^2 - 4*a*c

	t1 = (-b + sqrt(discr))/(2a)
	t2 = (-b - sqrt(discr))/(2a)

	if t1 < 0u"s"
		return t2
	elseif t2 < 0u"s"
		return t1
	else
		return min(t1, t2)
	end
end

"""
	beam_core_mass_ratio(a, dv::Speed, ve::Speed)

Mass ratio of a rocket where mass is lost without contributing to thrust, eg. a beam core rocket.

# Arguments
 - `a`: Fraction of propellant mass remaining after reaction (.22 for Frisbee rocket).
 - `dv::Speed`: mission delta-V (.25c for Frisbee rocket).
 - `ve::Speed`: exhaust velocity (.33c for Frisbee rocket).

 Derivation in Robert Frisbee's [How to build an antimatter rocket for interstellar missions](https://web.archive.org/web/20060601234257/http://www.aiaa.org/Participate/Uploads/2003-4676.pdf).
"""
function beam_core_mass_ratio(a, dv::Speed, ve::Speed)
	# un-tidied form of the equation looks more like this:
	# bcmr(a, ∆V, Isp) = ( ( (-2Isp*∆V/c_0^2+(1-a)-((1-a)^2+4a*Isp^2/c_0^2)^0.5)*(1-a+((1-a)^2+4a*Isp^2/c_0^2)^0.5) ) /( (-2Isp*∆V/c_0^2+(1-a)+((1-a)^2+4a*Isp^2/c_0^2)^0.5)*(1-a-((1-a)^2+4a*Isp^2/c_0^2)^0.5) ) )^(1 / ((1-a)^2+4a*Isp^2/c_0^2)^0.5)

	k1 = sqrt((1-a)^2+4a*ve^2/c_0^2)
	k2 = (-2ve * dv / c_0^2) + 1 - a

	p1 = k2 - k1
	p2 = 1 - a + k1
	p3 = k2 + k1
	p4 = 1 - a - k1
	p5 = 1 / k1

	((p1 * p2) / (p3 * p4))^p5
end

"""
    mass_ratio(dv::Speed, ve::Speed)

Compute the mass ratio for a rocket with total delta-V `dv` and exhaust velocity `ve`.
"""
mass_ratio(dv::Speed, ve::Speed) = exp(dv/ve)

delta_v(ve::Speed, mr::Real) = ve * log(mr)
delta_v(isp::Unitful.Time, mr::Real) = isp * g_n * log(mr)

rocket_thrust(ve::Speed, mdot::Unitful.MassFlow) = Unitful.uconvert(Unitful.N, mdot * ve)
rocket_thrust(ve::Speed, fp::Unitful.Power) = Unitful.uconvert(Unitful.N, (2 * fp) / ve)

rocket_power(ve::Speed, mdot::Unitful.MassFlow) = Unitful.uconvert(Unitful.W, (mdot * ve^2) / 2)
rocket_power(ve::Speed, thrust::Unitful.Force) = Unitful.uconvert(Unitful.W, (thrust * ve) / 2)

mass_flow(ve::Speed, thrust::Unitful.Force) = Unitful.uconvert(u"kg/s", thrust / ve)

"""
    exhaust_velocity(isp::Unitful.Time)

Convert specific impulse `isp` in seconds to exhaust velocity.
"""
exhaust_velocity(isp::Unitful.Time) = isp * g_n

"""
    exhaust_velocity(thrust::Unitful.Force, mdot::Unitful.MassFlow)

Compute exhaust velocity required for the given mass flow `mdot` to develop `thrust`.
"""
exhaust_velocity(thrust::Unitful.Force, mdot::Unitful.MassFlow) = Unitful.uconvert(u"m/s", thrust / mdot)

"""
    exhaust_velocity(fp::Unitful.Power, mdot::Unitful.MassFlow)

Compute exhaust velocity required for the given mass flow `mdot` and thrust power `fp`.
"""
exhaust_velocity(fp::Unitful.Power, mdot::Unitful.MassFlow) = Unitful.uconvert(u"m/s", sqrt((2 * fp) / mdot))

"""
    shkadov_thrust(stellar_luminosity::Unitful.Power, rim_angle::Angle)

Compute thrust of a Shkadov thruster with given rim angle (90° for max thrust_ about a star with the given `stellar_luminosity` in watts.

Note that this is a simple photon thruster.
"""
shkadov_thrust(stellar_luminosity::Unitful.Power, rim_angle::Angle) = (stellar_luminosity / 2c_0) * (1 - cos(rim_angle)) |> u"N"

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
#!/usr/local/julia/bin/julia

# @license BSD-3 https://opensource.org/licenses/BSD-3-Clause
# Copyright (c) 2020, Tim Redick
# All rights reserved.

using Pkg
Pkg.add("IJulia")

max_stars_repo_name: Tuebel/BayesionMethodsForHackers

max_issues_repo_name: Tuebel/BayesionMethodsForHackers

max_forks_repo_name: Tuebel/BayesianMethodsForHackers


text:
import CMBMissionSim
const CMBSim = CMBMissionSim
using Test

@testset "Quaternions" begin
    include("test_quaternions.jl")
end

@testset "Pointing generation" begin
    include("test_genpointings.jl")
end

@testset "Dipole analysis" begin
    include("test_dipole.jl")
end

@testset "Beam functions" begin
    include("test_beams.jl")
end


text:
include("../envs/RLEnvironments.jl")
using .RLEnvironments
using Plots

function act(Q, s, ϵ)
    rand() <= ϵ && return rand(1:size(Q, 1))
    a = @view Q[:, s]
    return rand(findall(a .== maximum(a)))
end

function run!(Q, env; epochs=20, ϵ=0.1, α=0.5, γ=0.99)
    rewards = zeros(epochs)
    steps = zeros(Int64, epochs)

    for i in 1:epochs
        done = false
        s1 = reset!(env)
        while !done
            a = act(Q, s1, ϵ)
            s2, r, done = step!(env, a)
            r2 = done ? 0.0 : maximum(view(Q, :, s2))
            Q[a, s1] += α * (r + γ * r2 - Q[a, s1])
            s1 = s2

            rewards[i] += r
            steps[i] += 1
        end
        ϵ = max(ϵ * 0.995, 0.01)
    end
    return rewards, steps
end

env = SimpleRooms{Int64, Int64, Float64}()
Q = zeros(length(env.actionspace), length(env.observationspace))

r, s = run!(Q, env)
plt = plot([r s], labels=["Reward" "Steps"], color=[:red :blue], layout=(2, 1))
display(plt)
#savefig(plt, "qlearning")


text:
"""
    DoNotCare{T}

A type that is equal to all things that are of type <:T.
For internal use.
Will interact weirdly with hash based dicts
"""
struct DoNotCare{T}
end

function Base.isequal(::DoNotCare{T}, ::A) where {T,A}
    A <: T
end

function Base.isequal(::DoNotCare{T1}, ::DoNotCare{T2}) where {T1,T2}
    T1 == T2
end

Base.isequal(a, d::DoNotCare)=isequal(d,a)
Base.:(==)(a, d::DoNotCare)=isequal(d,a)
Base.:(==)(d::DoNotCare, a)=isequal(d,a)
Base.:(==)(d::DoNotCare, a::DoNotCare)=isequal(d,a)

Base.isless(a, d::DoNotCare)=false
Base.isless(d::DoNotCare, a)=false
Base.isless(d::DoNotCare, a::DoNotCare)=false


text:
let md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
    dummy = CuFunction(md, "dummy")

    active_blocks(dummy, 1)
    active_blocks(dummy, 1; shmem=64)

    occupancy(dummy, 1)
    occupancy(dummy, 1; shmem=64)

    launch_configuration(dummy)
    launch_configuration(dummy; shmem=64)
    launch_configuration(dummy; shmem=64, max_threads=64)

    let cb_calls = 0
        launch_configuration(dummy; shmem=threads->(cb_calls += 1; 0))
        @test cb_calls > 0
    end
end


text:
"""
`DPGopts(m;σβ=1.,αΘ=0.0001,αw=0.001,αv=0.001,αu=0.001,γ=0.99,τ=0.001,iters=20_000, critic_update=:gradient,λrls=0.999,stepreduce_interval=1000,stepreduce_factor=0.995,hold_actor=1000)`

Structure with options to the DMP

# Fields

`σβ` Exploration noise covariance\n
`αΘ::Float64` Actor step size\n
`αw::Float64` Q-function step size 1\n
`αv::Float64` Q-function step size 2\n
`αu::Float64` Currently not used\n
`γ::Float64` Discount factor\n
`τ::Float64` Tracking factor between target and training networks\n
`λ::Float64` Regularization parameter for leastsquares Q-function fitting\n
`iters::Int64` Number of iterations to run\n
`m::Int64` Action dimension\n
`critic_update::Symbol` How to update the critic, can be chosen as `:gradient`, `:rls`, `:kalman`, `:leastsquares`\n
`λrls::Float64` If rls is used for critic update, use this forgetting factor\n
`stepreduce_interval::Int` The stepsize is reduced with this interval\n
`stepreduce_factor::Float64` The stepsize is reduced with this factor\n
`hold_actor::Int` Keep the actor from being updated for a few iterations in the beginning to allow the critic to obtain reasonable values\n
`experience_replay::Int` How many transitions to store in the replay buffer\n
`experience_ratio::Int` How many stored transitions to train with for every new experience\n
`momentum::Float64` Momentum term in actor gradient ascent (0.9)\n
`rmsprop::Bool` Use rmsprop for gradient training\n
`rmspropfactor::Float64` How much rmsprop, (0.9)\n
`critic_update_interval::Int` If critic is trained using leastsquares, only relearn the critic between this many mini-batches\n
`eval_interval::Int` With this interval the policy is evaluated without noise, parameters saved and stepsizes are reduced if divergence is detected.\n
`divergence_threshold::Float64` If the cost is higher than `divergence_threshold*bestcost`, the parameters are reset to their previously best values and stepsizes are reduced.\n
See example file or the paper Deterministic Policy Gradient Algorithms, Silver et al. 2014
"""
immutable DPGopts
    σβ
    αΘ::Float64
    αw::Float64
    αv::Float64
    αu::Float64
    γ::Float64
    τ::Float64
    λ::Float64
    iters::Int64
    m::Int64
    critic_update::Symbol
    λrls::Float64
    stepreduce_interval::Int
    stepreduce_factor::Float64
    hold_actor::Int
    experience_replay::Int
    experience_ratio::Int
    momentum::Float64
    rmsprop::Bool
    rmspropfactor::Float64
    critic_update_interval::Int
    eval_interval::Int
    divergence_threshold::Float64
end

DPGopts(m;
σβ=1.,
αΘ=0.001,
αw=0.01,
αv=0.01,
αu=0.001,
γ=0.99,
τ=0.001,
λ=1e-4,
iters=20_000,
critic_update=:gradient,
λrls=0.999,
stepreduce_interval=1000,
stepreduce_factor=0.995,
hold_actor=1000,
experience_replay=0,
experience_ratio=10,
momentum=0.9,
rmsprop=true,
rmspropfactor=0.9,
critic_update_interval=1,
eval_interval=10,
divergence_threshold=1.5) =
DPGopts(σβ,αΘ,αw,αv,αu,γ,τ,λ,iters,m,critic_update,λrls,stepreduce_interval,stepreduce_factor,hold_actor,experience_replay,experience_ratio,momentum,rmsprop,rmspropfactor,critic_update_interval,eval_interval,divergence_threshold)

"""
Structure with functions to pass to the DMP

# Fields

`μ,Q,gradients,simulate,exploration,reward`

See example file or the paper by Silver et al. 2014
"""
immutable DPGfuns
    μ::Function
    Q::Function
    gradients::Function
    simulate::Function
    exploration::Function
    reward::Function
end

"""
Structure which contains the parameters of the DPG optimization problem\n
`Θ` parameters in the actor\n
`w` parameters in the Q-function\n
`v` parameters in the Q-function\n
All parameters should be a subtype of AbstractVector
A typical Q-function looks like `Q = (∇μ(s)*(a-μ(s)))'w + V(s,v)`

"""
type DPGstate{T1<:AbstractVector,T2<:AbstractVector,T3<:AbstractVector}
    Θ::T1
    w::T2
    v::T3
end

function J(x,a,r)
    cost = @parallel (+) for t = 1:size(x,1)
    r(x[t,:][:],a[t,:][:],t)
end
-cost
end

# """
# `cost, Θ, w, v = dpg(opts, funs, state0, x0)`
#
# Main function.
#
# # Arguments
# `opts::DPGopts` structure with options and parameters\n
# `funs::DPGfuns` structure with functions\n
# `state0::DPGstate` initial parameters
# `x0` initial system state
# """
function dpg(opts, funs, state0, x0,C, progressfun = (Θ,w,v,C,i,x,uout,cost)->0)
    println("=== Deterministic Policy Gradient ===")
    # Expand input structs
    σβ          = opts.σβ
    αΘ          = opts.αΘ
    αw          = opts.αw
    αv          = opts.αv
    αu          = opts.αu
    γ           = opts.γ
    τ           = opts.τ
    iters       = opts.iters
    m           = opts.m
    n           = length(x0)
    critic_update= opts.critic_update
    λrls        = opts.λrls
    mem         = SequentialReplayMemory(opts.experience_replay)
    μ           = funs.μ
    Q           = funs.Q
    gradients   = funs.gradients
    simulate    = funs.simulate
    exploration = funs.exploration
    r           = funs.reward
    println("Training using $critic_update")

    # Initialize parameters
    Θ           = deepcopy(state0.Θ) # Weights
    w           = deepcopy(state0.w)
    v           = deepcopy(state0.v)
    Θt          = deepcopy(Θ) # Tracking weights
    wt          = deepcopy(w)
    vt          = deepcopy(v)
    Pw          = size(Θ,1)
    Pv          = size(v,1)
    Θb          = deepcopy(Θ) # Best weights
    wb          = deepcopy(w)
    vb          = deepcopy(v)
    dΘs         = 1000ones(Pw) # Weight gradient states
    dΘs2        = 0*ones(Pw)
    dws         = 100ones(Pw)
    dvs         = 100ones(Pv)
    cost        = zeros(iters)
    bestcost    = Inf
    x,uout      = simulate(Θ, x0)
    T           = size(x,1)

    # TODO: Make the parameters below part of the options
    if critic_update == :rls
        Pvw = 10eye(Pw+Pv)
    elseif critic_update == :kalman
        Pk = 1000eye(Pw+Pv)
        R2 = 1
        R12 = zeros(Pw+Pv)
    end

    s = zeros(n)

    function train!(tvec, update_actor,update_critic=true)
        dΘ,dw,dv = zeros(Θ),zeros(w),zeros(v)
        if critic_update ==:leastsquares && update_critic
            basis = sample_uniform!(mem,Pv)
            for (it,t) in enumerate(basis)
                C[:,it] = t.s
            end
            push!(mem,basis)
            A = Matrix{Float64}(length(mem),Pv+Pw)
            targets = Vector{Float64}(length(mem))
            for (it,t) in enumerate(mem)
                a1          = μ(t.s1,Θ,t.t)
                _, ∇wQ, ∇vQ, _ = gradients(t.s1,t.s,a1,t.a,Θ,w,v,t.t,C)
                y           = t.r + γ * Q(t.s1,a1,vt,wt,Θt,t.t,C)
                targets[it] = y
                # t.δ         = (y - Q(t.s,t.a,v,w,Θ,t.t))[1] # Not needed for batch learning
                A[it,1:Pv]  = ∇vQ
                A[it,Pv+1:end]  = ∇wQ
            end
            vw = [A; opts.λ*eye(Pv+Pw)]\[targets+100;zeros(Pv+Pw)] #TODO magic numer 100?
            v[:],w[:] = vw[1:Pv],vw[Pv+1:end]
        end
        for t in tvec
            a1          = μ(t.s1,Θt,t.t)
            ∇aQ, ∇wQ,∇vQ, ∇μ = gradients(t.s1,t.s,a1,t.a,Θ,w,v,t.t,C)
            y           = t.r + γ * Q(t.s1,a1,vt,wt,Θt,t.t,C)
            t.δ         = (y - Q(t.s,t.a,v,w,Θ,t.t,C))[1]
            if critic_update == :rls
                vw,Pvw      = RLS([v;w], y, [∇vQ;∇wQ], Pvw, λrls)
                v[:],w[:]   = vw[1:Pv],vw[Pv+1:end]
            elseif critic_update == :kalman

                Φ = [∇vQ;∇wQ]
                # R1 = ΦΦ', to only update covariance in the direction of incoming data
                vw,Pk = kalman(Φ*Φ',R2,R12,[v;w], y, Φ, Pk)
                v[:],w[:] = vw[1:Pv],vw[Pv+1:end]
            elseif critic_update ==:gradient
                dw += t.δ * ∇wQ  #- γ * ϕ(s1,a1) * ϕu
                dv += t.δ * ∇vQ   #- γ * ϕ(s1) * ϕu
            end
            dΘ += ∇μ*∇aQ
        end
        dΘ./= length(tvec)
        # RMS prop update parameters (gradient divided by running average of RMS gradient, see. http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf
        dΘs[:] = opts.rmspropfactor*dΘs + (1-opts.rmspropfactor)*dΘ.^2
        if update_actor
            ΔΘ = αΘ * dΘ
            if opts.rmsprop
                ΔΘ ./= sqrt(dΘs+1e-10) # RMSprop
            end
            # ΔΘ = dΘ./sqrt(dΘs+1e-10).*sqrt(dΘs2) # ADAdelta
            # dΘs2[:] = 0.9dΘs2 + 0.1ΔΘ.^2 # ADAdelta
            dΘs2 = opts.momentum*dΘs2 + ΔΘ # Momentum + RMSProp
            Θ[:]  += dΘs2
        end

        if critic_update == :gradient
            Δw = αw/T * dw
            Δv = αv/T * dv
            dws[:] = opts.rmspropfactor*dws + (1-opts.rmspropfactor)*dw.^2
            dvs[:] = opts.rmspropfactor*dvs + (1-opts.rmspropfactor)*dv.^2
            if opts.rmsprop
                Δw ./= sqrt(dws)+1e-10
                Δv ./= sqrt(dvs)+1e-10
            end
            w += Δw
            v += Δv
        end

        # Update tracking networks
        Θt[:], wt[:], vt[:] = τ*Θ + (1-τ)*Θt, τ*w + (1-τ)*wt, τ*v + (1-τ)*vt
        nothing
    end

    # Main loop ================================================================
    for i = 1:iters
        x0i         = x0 #+ 2randn(n) # TODO: this should not be hard coded
        noise       = exploration(σβ)
        x,uout      = simulate(Θ, x0i, noise)
        batch = Vector{Transition}(T-1)
        for ti = 1:T-1
            s1          = x[ti+1,:][:]
            s           = x[ti,:][:]
            a           = uout[ti,:][:]
            ri          = r(s1,a,ti)
            cost[i]    -= ri
            trans       = Transition(s,s1,a,ri,0.,ti,i)
            batch[ti]   = trans
            opts.experience_replay > 0 && push!(mem, trans)
        end
        train!(batch, i > opts.hold_actor)

        if i % opts.stepreduce_interval == 0
            αΘ  *= opts.stepreduce_factor
            αw  *= opts.stepreduce_factor
            αv  *= opts.stepreduce_factor
        end

        if opts.experience_replay > 0 && i > 10
            for ei = 1:opts.experience_ratio
                trans = sample_uniform!(mem,T-1)
                train!(trans, i > opts.hold_actor, ei-1 % opts.critic_update_interval == 0)
                push!(mem, trans)
            end
            ((i % 50) == 0) && sort!(mem)
        end

        if (i-1) % opts.eval_interval == 0 # Simulate without noise and evaluate cost
            x,uout = simulate(Θ, x0) # Evaluate using latest parameters and possibly revert back to tracking parameters
            cost[i] = J(x,uout,r)
            progressfun(Θ,w,v,C,i,x,uout, cost)
            if critic_update == :gradient
                println(i, ", cost: ", cost[i] |> r5, " norm ∇Θ: ", Σ½(dΘs) |> r5, " norm ∇w: ", Σ½(dws) |> r5, " norm ∇v: ", Σ½(dvs) |> r5)#, " trace(P): ", trace(Pvw) |> r5)
            else
                println(i, ", cost: ", cost[i] |> r5, " norm ∇Θ: ", Σ½(dΘs) |> r5)
            end
            if cost[i] < bestcost
                bestcost = cost[i]
                # TODO: changed to saveing tracking networks, to be more likely to escape local minima
                Θb = deepcopy(Θ)
                wb = deepcopy(w)
                vb = deepcopy(v)

            elseif cost[i] > opts.divergence_threshold*bestcost
                print_with_color(:orange,"Reducing stepsizes due to divergence\n")
                αΘ  /= 5
                αw  /= 5
                αv  /= 5
                αu  /= 5
                σβ  /= 2
                Θ, w, v = deepcopy(Θb), deepcopy(wb), deepcopy(vb) # reset parameters
            end
        end

    end
    println("Done. Minimum cost: $(minimum(cost[1:opts.eval_interval:end])), ($(minimum(cost)))")
    return cost, Θb, wb, vb, mem # Select the parameters with lowest cost
end


text:
function gasPlasma(Ip::Float64, npts::Int64)
    ns = sqrt(13.59844*e/Ip)
    wp = Ip/hbar
    cn2 = 2^(2*ns)/(ns*gamma(ns+1)*gamma(ns))
    return GasPlasma(ns=ns,wp=wp,cn2=cn2, npts=npts)
end

function ionRateAdk(Et::Array{Float64,1}, WE::Array{Float64,1}, Ip::Float64, ns::Float64, wp::Float64, cn2::Float64)
    @. WE = e*abs(Et)/sqrt(2*m_e*Ip)
    # for i in eachindex(WE) the same as above
    #     WE[i] = e*abs(Et[i])/sqrt(2*m_e*Ip)
    # end
    @. WE =  wp*cn2*(4wp/WE)^(2ns-1)*exp(-4.0/3.0*wp/WE)
    # for i in eachindex(WE) The same as above
    #     WE[i] = wp*cn2*(4wp/WE[i])^(2ns-1)*exp(-4.0/3.0*wp/WE[i])
    # end
    # @time  WE[notfinite.(WE)] .= 0.0
    for i in eachindex(WE)
        if ~isfinite(WE[i])
            WE[i] = 0.0
        end
    end
end

function adk(Et::Array{Float64,1}, τ::Array{Float64,1}, gas::Gas)
    adk(Et, τ, gas.gp.Na, gas.gp.Ne, gas.gp.WE, gas)
end

function adk(Et::Array{Float64,1}, τ::Array{Float64,1}, Na::Array{Float64,1},
            Ne::Array{Float64,1}, WE::Array{Float64,1}, gas::Gas)
    ionRateAdk(Et, gas.gp.WE, gas.Ip, gas.gp.ns, gas.gp.wp, gas.gp.cn2)
    # @time ionRateAdk(Et, gas)

    # @time Na .= exp.(-cumul_integrate(τ, WE))
    cumtrapz(τ, WE, gas.gp.temp1)

    @. Na = exp(-gas.gp.temp1)
    # @time for i in eachindex(Na) the above is better
    #     Na[i] = exp(-gas.gp.temp1[i])
    # end

    # Ne[notfinite.(Ne) .| (Ne .< 0.0)] .= 0.0
    # Ne .= gas.p*gas.pp*N0.*(1.0.-Na)
    for i in eachindex(Ne)
        Ne[i] = gas.p*gas.pp*N0*(1.0-Na[i])

        if ~isfinite(Ne[i]) | (Ne[i] < 0.0)
            Ne[i] = 0.0
        end
    end

    # Na .= gas.p*gas.pp*N0.*Na
    # @time Na[notfinite.(Na) .| (Na .< 0.0)] .= 1.0
    for i in eachindex(Na)
        Na[i] = gas.p*gas.pp*N0*Na[i]

        if ~isfinite(Na[i]) | (Na[i] < 0.0)
            Na[i] = 0.0
        end
    end
end

function geisslerPlasmaUPPE(Et::Array{Float64,1}, τ::Array{Float64,1}, dτ::Float64, gas::Gas)
    adk(Et, τ, gas)
    @. gas.gp.temp2 = gas.gp.Na*gas.gp.WE/Et
    #the above is better. something wrong with passing gas.gp.element to for loop
    # (seems to be a type problem. I need to specify all types to git rid of the problem)
    # @time for i in eachindex(gas.gp.temp1)
    #     gas.gp.temp2[i] = gas.gp.Na[i]*gas.gp.WE[i]/Et[i]
    # end
    cumtrapz(τ, gas.gp.temp2, gas.gp.Pe_l)
    @. gas.gp.Pe_l .= gas.Ip .* gas.gp.Pe_l

    # gas.gp.Pe_f .= cumul_integrate(τ, gas.gp.Ne.*Et)
    # gas.gp.Pe_f .= e^2/m_e.*cumul_integrate(τ, gas.gp.Pe_f)
    @. gas.gp.temp3 = gas.gp.Ne.*Et
    cumtrapz(τ, gas.gp.temp3, gas.gp.temp4)
    cumtrapz(τ, gas.gp.temp4, gas.gp.Pe_f)

    @. gas.gp.Pe = e^2/m_e*gas.gp.Pe_f + gas.gp.Pe_l
    @. gas.gp.Pe[notfinite.(gas.gp.Pe)] = 0.0
end

struct IonizationRatePPT
    ip::Float64           # ionization potential [J]
    Ui_a::Float64         # ionization potential [au]
    E0_a::Float64         # atomic field strength [au]
    ns::Float64           # effective principle quantum number
    ls::Float64           # effective azimuthal quantum number
    Cnl2::Float64         #
    w_a::Float64          # frequency [au]
    l::Int                # azimuthal quantum number
end

function IonizationRatePPT(Ip, n, w; Z=1)
    l = n-1
    ip = Ip*e
    Ui_a = ip/eh
    E0_a = (2Ui_a)^(3/2)
    ns = Z/sqrt(2Ui_a)
    ls = ns - 1
    Cnl2 = 2^(2ns)/(ns*gamma(ns + ls + 1)*gamma(ns - ls))
    w_a = w*a_t
    IonizationRatePPT(ip, Ui_a, E0_a, ns, ls, Cnl2, w_a, l)
end

function ppt(Et, ir::IonizationRatePPT; tol::Float64 = 1e-3, rcycle::Bool = true)
    E_a = abs(E)/a_E
    Eratio = E_a/ir.E0_a
    #if Eratio > 0.5
    #    error("electric field intensity beyond PPT validity: $Eratio")
    #end
    g = ir.w_a*sqrt(2ir.Ui_a)/E_a
    g2 = g*g
    beta = 2g/sqrt(1 + g2)
    alpha = 2(asinh(g) - g/sqrt(1 + g2))
    Up_a = E_a*E_a/(4ir.w_a*ir.w_a)
    Uit_a = ir.Ui_a + Up_a
    v = Uit_a/ir.w_a
    G = 3/(2g)*((1 + 1/(2g2))*asinh(g) - sqrt(1 + g2)/(2g))
    ret = 0.0
    div = 0.0
    for m=-ir.l:ir.l
        div += 1
        am = abs(m)
        flm = ((2ir.l + 1)*factorial(ir.l + am)
               / (2^am*factorial(am)*factorial(ir.l - am)))
        Am = 4/(sqrt(3pi)*factorial(am))*g2/(1 + g2)
        lret = sqrt(3/(2pi))*ir.Cnl2*flm*ir.Ui_a
        lret *= (2ir.E0_a/(E_a*sqrt(1 + g2)))^(2ir.ns - am - 3/2)
        lret *= Am*exp(-2ir.E0_a*G/(3E_a))
        # remove cycle average # TODO: check this
        if rcycle
            lret *= sqrt(pi*ir.E0_a/(3E_a))
        end
        k = ceil(v)
        sum = 0.0
        dsum = 0.0
        while true
            diff = k - v
            dsum = exp(-alpha*diff)*phi(m, sqrt(beta*diff))
            sum += dsum
            k += 1.0
            dsum/sum < tol && break
        end
        lret *= sum
        ret += lret
    end
    ret/(a_t*div)
end

#=
for gnlse
    function ionRateAdk(Et, Ip)
        Ip *= e
        ns = sqrt(13.59844*e/Ip)
        wp = Ip/hbar
        cn2 = 2^(2*ns)/(ns*gamma(ns+1)*gamma(ns))
        wt = @. e*abs(Et)/sqrt(2*m_e*Ip)
        WE = @. wp*cn2*(4wp/wt)^(2ns-1)*exp(-4.0/3.0*wp/wt)
        WE[notfinite.(WE)] .= 0.0
        WE
    end

    #for gnlse
    function adk(Et, τ, Ip, p)
        WE = ionRateAdk(Et, Ip)
        Na = exp.(-cumul_integrate(τ, WE))
        Ne = p*N0.*(1.0.-Na)
        Na = p*N0.*Na
        Na[notfinite.(Na) .| (Na .< 0.0)] .= 1.0
        # for i in eachindex(Na)
        Ne[notfinite.(Ne) .| (Ne .< 0.0)] .= 0.0
        [Na, Ne, WE]
    end

    #for gnlse
    function plasmaTerms(Et, τ, U, Ip, p, λ0, diameter)
        ω0 = 2pi*c/λ0
        dτ = τ[2] - τ[1]
        Na, Ne, WE = adk(Et, τ, Ip, p)
        frequancy_term = @. e^2.0*Ne/(2ϵ0*m_e*c*ω0)
        Aeff = effectiveArea(diameter)
        loss_term = @. Na*WE*Aeff*Ip*e/(2.0*abs2(U))
        loss_term[notfinite.(loss_term)] .= 0.0
        Pl = -1im*shifted_ifft(U .* frequancy_term)  .- shifted_ifft(U .* loss_term)
        return Pl
    end

   # old plasma gas mixture for air
    # if fiber.gas == :Air
    #
    #     O2_Ip = 12.0697*e
    #     N2_Ip = 15.581*e
    #
    #     nsO2 = sqrt(13.59844*e/O2_Ip)
    #     wpO2 = O2_Ip/hbar
    #     cn2O2 = 2^(2*nsO2)/(nsO2*gamma(nsO2+1)*gamma(nsO2))
    #
    #     nsN2 = sqrt(13.59844*e/N2_Ip)
    #     wpN2 = N2_Ip/hbar
    #     cn2N2 = 2^(2*nsN2)/(nsN2*gamma(nsN2+1)*gamma(nsN2))
    #
    #     plconstantsO2 = [O2_Ip, nsO2, wpO2, cn2O2, fiber.p]
    #     plconstantsN2 = [N2_Ip, nsN2, wpN2, cn2N2, fiber.p[1]]
    #
    #     GePlO2 = geisslerPlasmaUPPE{Float64}(npts = size(τ)[1], p = fiber.p[1], plconstants = plconstantsO2)
    #     GePlN2 = geisslerPlasmaUPPE{Float64}(npts = size(τ)[1], p = fiber.p[1], plconstants = plconstantsN2)
    #     GePl = [GePlO2, GePlN2]
    # else
    #     Ip = fiber.Ip*e
    #     ns = sqrt(13.59844*e/Ip)
    #     wp = Ip/hbar
    #     cn2 = 2^(2*ns)/(ns*gamma(ns+1)*gamma(ns))
    #     plconstants = [Ip, ns, wp, cn2, fiber.p[1]]
    #     ir = IonizationRatePPT(fiber.Ip, 3, ω0; Z=1)
    #     GePl = [geisslerPlasmaUPPE{Float64}(npts = size(τ)[1],plconstants = plconstants, pltype=pltype, p=fiber.p[1])]
    # end
    # if plasma == true
    #     Ip = fiber.Ip*e
    #     ns = sqrt(13.59844*e/Ip)
    #     wp = Ip/hbar
    #     cn2 = 2^(2*ns)/(ns*gamma(ns+1)*gamma(ns))
    #     plconstants = [Ip, ns, wp, cn2, fiber.p[1]]
    #     ir = IonizationRatePPT(fiber.Ip, 3, g.ω0; Z=1)
    #     GePl = [geisslerPlasmaUPPE{Float64}(npts = size(τ)[1],plconstants = plconstants, pltype=pltype, p=fiber.p[1])]
    # else
    #     GePl =   [geisslerPlasmaUPPE{Float64}( npts=g.npts, plconstants=[],  p=gases[1].p)]
    # end

# function geisslerPlasmaUPPE(Et, τ, gas, pltype)
#     gePl = gas.gasplasma

#     Na, Ne, WE = adk(Et, τ, gas.Ip, gas.p*gas.pp)
#     gePl.WE .= WE
#     gePl.Na .= Na
#     gePl.Ne .= Ne
#     gePl.Pe_l .= gas.Ip .* cumul_integrate(τ, gePl.Na.*gePl.WE./Et)
#     gePl.Pe_f .= cumul_integrate(τ, gePl.Ne.*Et)
#     gePl.Pe_f .= e^2/m_e.*cumul_integrate(τ, gePl.Pe_f)
#     gePl.Pe .= gePl.Pe_f .+ gePl.Pe_l
#     gePl.Pe[notfinite.(gePl.Pe)] .= 0.0

#     # @. gePl.temp3 = gePl.Na*gePl.WE/Et

#     # # @time for i  in eachindex(temp3)
#     # #     temp3[i] = Na[i]*WE[i]/Et[i]
#     # # end
#     # cumtrapz(τ, gePl.temp3, gePl.Pe_l)
#     # println(gePl.Pe_l)

#     # @. gePl.Pe_l = gas.Ip*gePl.Pe_l

#     # # @time for i in eachindex(Pe_l)
#     # #     Pe_l[i] = Ip*Pe_l[i]
#     # # end
#     # @. gePl.temp4 = gePl.Ne*Et
#     # # @time for i in eachindex(temp4)
#     # #     temp4[i] = Ne[i]*Et[i]
#     # # end
#     # cumtrapz(τ, gePl.temp4, gePl.temp1)
#     # cumtrapz(τ, gePl.temp1, gePl.Pe_f)
#     # for i in eachindex(gePl.Pe)
#     #     gePl.Pe[i] = e^2/m_e*gePl.Pe_f[i] + gePl.Pe_l[i]
#     #     if ~isfinite(gePl.Pe[i])
#     #         gePl.Pe[i] = 0.0
#     #     end
#     # end
# end

#=
all below is for uppe
=#
# function ionRateAdk(Et, Ip, gePl)
#     @. gePl.WE = e*abs(Et)/sqrt(2*m_e*Ip)
#     @. gePl.WE =  gePl.wp*gePl.cn2*(4*gePl.wp/gePl.WE)^(2*gePl.ns-1)*exp(-4.0/3.0*gePl.wp/gePl.WE)
#     # @time  WE[notfinite.(WE)] .= 0.0
#     for i in eachindex(gePl.WE)
#         if ~isfinite(gePl.WE[i])
#             gePl.WE[i] = 0.0
#         end
#     end
# end

# function ionRateAdk(Et, WE, plconstants)
#     Ip, ns, wp, cn2, p = plconstants
#     @. WE = e*abs(Et)/sqrt(2*m_e*Ip)
#     @. WE =  wp*cn2*(4wp/WE)^(2ns-1)*exp(-4.0/3.0*wp/WE)
#     # @time  WE[notfinite.(WE)] .= 0.0
#     for i in eachindex(WE)
#         if ~isfinite(WE[i])
#             WE[i] = 0.0
#         end
#     end
# end

# initial gas mixture try (did not work)
# function ionRate(Et, τ, gePl, gas,pltype)
#     if pltype == :adk
#         ionRateAdk(Et, gas.Ip, gePl)
#     elseif pltype == :ppt
#         ppt.(Et, gePl.ir::IonizationRatePPT)
#     end
#     # p = gePl.plconstants[5]
#     cumtrapz(τ, gePl.WE, gePl.temp2)
#     # Na .= exp.(-temp2)
#     for i in eachindex(gePl.Na)
#         gePl.Na[i] = exp(-gePl.temp2[i])
#     end
#     for i in eachindex(gePl.Ne)
#         # gePl.Ne[i] = p*N0*(1.0-gePl.Na[i]
#         gePl.Ne[i] = gas.p*gas.pp*N0*(1.0-gePl.Na[i])

#         if ~isfinite(gePl.Ne[i]) | (gePl.Ne[i] < 0.0)
#             gePl.Ne[i] = 0.0
#         end
#     end
#     for i in eachindex(gePl.Na)
#         gePl.Na[i] = gas.p*gas.pp*N0*gePl.Na[i]
#         if ~isfinite(gePl.Na[i]) | (gePl.Na[i] < 0.0)
#             gePl.Na[i] = 0.0
#         end
#     end
# end

# function geisslerPlasmaUPPE(Et, τ, dτ, Ip, p, Pe0, Na, Ne, WE, Pe_l, Pe_f)
#     Na, Ne, WE = adk(Et, τ, Ip, p)
#     Pe_l .= Ip .* cumul_integrate(τ, Na.*WE./Et)
#     Pe_f .= cumul_integrate(τ, Ne.*Et)
#     Pe_f .= e^2/m_e.*cumul_integrate(τ, Pe_f)
#     Pe0 .= Pe_f .+ Pe_l
#     Pe0[notfinite.(Pe0)] .= 0.0
#     return Pe0
# end

# function ionRateAdk(Et, Ip)

#     ns = sqrt(13.59844*e/Ip)
#     wp = Ip/hbar
#     cn2 = 2^(2*ns)/(ns*gamma(ns+1)*gamma(ns))
#     wt = @. e*abs(Et)/sqrt(2*m_e*Ip)
#     WE = @. wp*cn2*(4wp/wt)^(2ns-1)*exp(-4.0/3.0*wp/wt)
#     WE[notfinite.(WE)] .= 0.0
#     WE
# end

# function adk(Et, τ, Ip, p)
#     WE = ionRateAdk(Et, Ip)
#     Na = exp.(-cumul_integrate(τ, WE))
#     Ne = p*N0.*(1.0.-Na)
#     Na = p*N0.*Na
#     Na[notfinite.(Na) .| (Na .< 0.0)] .= 1.0
#     Ne[notfinite.(Ne) .| (Ne .< 0.0)] .= 0.0
#     [Na, Ne, WE]
# end

# function ionRateAdk(Et, gas)
#     @. gas.gp.wt =  e*abs(Et)/sqrt(2*m_e*gas.Ip)
#     @. gas.gp.WE =  gas.gp.wp*gas.gp.cn2*(4*gas.gp.wp/gas.gp.wt)^(2*gas.gp.ns-1)*exp(-4.0/3.0*gas.gp.wp/gas.gp.wt)
#     gas.gp.WE[notfinite.(gas.gp.WE)] .= 0.0
# end
=#

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev