text:
using ApproxFun, SingularIntegralEquations, DualNumbers, RiemannHilbert, LinearAlgebra, FastTransforms, SpecialFunctions, Test
import ApproxFunBase: ArraySpace, pieces, dotu, interlace
import RiemannHilbert: RiemannDual, LogNumber, fpstieltjesmatrix!, fpstieltjesmatrix, orientedleftendpoint, orientedrightendpoint, finitepart, fpcauchymatrix, collocationvalues, collocationpoints
import SingularIntegralEquations: stieltjesmoment, stieltjesmoment!, undirected, Directed, ⁺, ⁻, istieltjes
import SingularIntegralEquations.HypergeometricFunctions: speciallog

@testset "RiemannDual" begin
    for h in (0.1,0.01), a in (2exp(0.1im),1.1)
        @test log(RiemannDual(0,a))(h) ≈ log(h*a)
        @test log(RiemannDual(Inf,a))(h) ≈ log(a/h)
    end

    for h in (0.1,0.01), a in (2exp(0.1im),1.1)
        @test log1p(RiemannDual(-1,a))(h) ≈ log(h*a)
        @test log1p(RiemannDual(Inf,a))(h) ≈ log(a/h)
    end

    h = 0.0000001
    for z in (RiemannDual(-1,-1), RiemannDual(1,1), RiemannDual(-1,2exp(0.1im)), RiemannDual(1,2exp(0.1im))),
            k = 0:1
        l = stieltjesjacobimoment(0,0,k,z)
        @test l(h) ≈ stieltjesjacobimoment(0,0,k,realpart(z)+epsilon(z)h) atol=1E-5
    end

    h=0.0001
    for z in (RiemannDual(1,3exp(0.2im)), RiemannDual(1,0.5exp(-1.3im)),
                RiemannDual(-1,3exp(0.2im)), RiemannDual(-1,0.5exp(-1.3im)),
                RiemannDual(-1,1), RiemannDual(1,-1))
        @test atanh(z)(h) ≈  atanh(realpart(z)+epsilon(z)h) atol = 1E-4
    end

    z = RiemannDual(1,-0.25)
    h = 0.0000001
    @test speciallog(z)(h) ≈ speciallog(realpart(z)+epsilon(z)h) atol=1E-4

    h = 0.00001
    for z in (RiemannDual(-1,-1), RiemannDual(-1,exp(0.1im)), RiemannDual(-1,exp(-0.1im)))
        @test stieltjesjacobimoment(0.5,0,0,z)(h) ≈ stieltjesjacobimoment(0.5,0,0,realpart(z)+epsilon(z)h) atol=1E-4
    end
end

@testset "Legendre Cauchy" begin
    f = Fun(exp,Legendre())

    h = 0.00001
    for z in  (RiemannDual(-1,-1), RiemannDual(-1,1+im), RiemannDual(-1,1-im))
        @test cauchy(f, z)(h) ≈ cauchy(f, realpart(z) + epsilon(z)h) atol=1E-4
    end
end

@testset "Directed and RiemannDual" begin
    @test undirected(Directed{false}(RiemannDual(0,-1))) == 0

    @test real(LogNumber(2im,im+1)) == LogNumber(0,1)
    @test imag(LogNumber(2im,im+1)) == LogNumber(2,1)
    @test conj(LogNumber(2im,im+1)) == LogNumber(-2im,1-im)

    @test log(Directed{false}(RiemannDual(0,-1))) == LogNumber(1,π*im)
    @test log(Directed{true}(RiemannDual(0,-1))) == LogNumber(1,-π*im)

    @test log(Directed{false}(RiemannDual(0,-1-eps()*im))) == LogNumber(1,π*im)
    @test log(Directed{false}(RiemannDual(0,-1+eps()*im))) == LogNumber(1,π*im)

    @test log(Directed{true}(RiemannDual(0,-1-eps()*im))) == LogNumber(1,-π*im)
    @test log(Directed{true}(RiemannDual(0,-1+eps()*im))) == LogNumber(1,-π*im)

    z = Directed{false}(RiemannDual(1,-2))

    for k=0:1, s=(false,true)
        z = Directed{s}(RiemannDual(-1,2))
        l = stieltjesmoment(Legendre(),k,z)
        h = 0.00000001
        @test l(h) ≈ stieltjesmoment(Legendre(),k,-1 + epsilon(z.x)h + (s ? 1 : -1)*eps()*im) atol=1E-5
    end

    @test RiemannHilbert.orientedleftendpoint(ChebyshevInterval()) == RiemannDual(-1.0,1)
    @test RiemannHilbert.orientedrightendpoint(ChebyshevInterval()) == RiemannDual(1.0,-1)
end

@testset "Interval FPStieltjes" begin
    Γ = ChebyshevInterval()
    f = Fun(x->exp(-40(x-0.1)^2), Legendre())
    C = Array{ComplexF64}(undef, ncoefficients(f), ncoefficients(f))
    d = Segment(im,2im)

    fpstieltjesmatrix!(C, space(f), d)
    c = Fun(d, chebyshevtransform(C*coefficients(f); kind=2))
    @test c(1.5im) ≈ stieltjes(f,1.5im)

    d = Segment(-1,-1+im)
    fpstieltjesmatrix!(C, space(f), d)
    @test norm(C) ≤ 100
    c = Fun(d, chebyshevtransform(C*coefficients(f); kind=2))
    @test c(-1+0.5im) ≈ stieltjes(f,-1+0.5im)

    d = Segment(1,1+im)
    fpstieltjesmatrix!(C, space(f), d)
    @test norm(C) ≤ 200
    c = Fun(d, chebyshevtransform(C*coefficients(f); kind=2))
    @test c(1+0.5im) ≈ stieltjes(f,1+0.5im)

    d = ChebyshevInterval()
    fpstieltjesmatrix!(C, space(f), d)
    @test norm(C) ≤ 200
    c = Fun(d, chebyshevtransform(C*coefficients(f); kind=2))
    @test c(0.5) ≈ stieltjes(f,0.5⁻)

    d = Segment(0,1)
    f = Fun(x->exp(-200(x-0.6)^2), Legendre(d))
    C = fpstieltjesmatrix(space(f), ncoefficients(f), ncoefficients(f))
    @test norm(C) ≤ 200
    c = Fun(d, chebyshevtransform(C*coefficients(f); kind=2))
    @test c(0.5) ≈ stieltjes(f,0.5⁻)
end

@testset "finitepart stieltjes" begin
    f = Fun(exp,Legendre())
    f1 = Fun(exp,Legendre(-1..0))
    f2 = Fun(exp,Legendre(0..1))
    fp = f1+f2

    @test stieltjes(f,0.0⁻) ≈ finitepart(stieltjes(f1,RiemannDual(0.0,-im)) + stieltjes(f2,RiemannDual(0.0,-im)))
    @test stieltjes(f,0.0⁻) ≈ finitepart(stieltjes(f1,RiemannDual(0.0,exp(-0.1im))) + stieltjes(f2,RiemannDual(0.0,exp(-0.1im))))
    @test stieltjes(f,0.0⁻) ≈ finitepart(stieltjes(f1,Directed{false}(RiemannDual(0.0,-1.0))) + stieltjes(f2,RiemannDual(0.0,-1.0)))

    @test stieltjes(fp,RiemannDual(0.0,-im)) ≈ stieltjes(f1,RiemannDual(0.0,-im)) + stieltjes(f2,RiemannDual(0.0,-im))
    @test stieltjes(f,0.0⁻) ≈ finitepart(stieltjes(fp,RiemannDual(0.0,-im)))
end

@testset "Two interval" begin
    @testset "-1..0 and 0..1" begin
        sp = Legendre(-1 .. 0) ⊕ Legendre(0 .. 1)
        f = Fun(x->exp(-40(x-0.1)^2), sp)
        v = components(f)
        ns = ncoefficients.(v)
        C11 = fpstieltjesmatrix(space(v[1]), ncoefficients(v[1]), ncoefficients(v[1]))
        c_vals11 = C11*v[1].coefficients
        h = 0.00000001; @test stieltjes(v[1],h*im) ≈ stieltjes(v[1],RiemannDual(0.0,im))(h) atol=1E-6
        @test finitepart(stieltjes(v[1],Directed{false}(RiemannDual(0.0,-1)))) ≈ c_vals11[1]

        C22 = fpstieltjesmatrix(space(v[2]), ncoefficients(v[2]), ncoefficients(v[2]))
        c_vals22 = C22*v[2].coefficients
        h = 0.00000001; @test stieltjes(v[2],h*im) ≈ stieltjes(v[2],RiemannDual(0.0,im))(h) atol=1E-6
        @test finitepart(stieltjes(v[2],Directed{false}(RiemannDual(0.0,1)))) ≈ c_vals22[end]

        C12 = fpstieltjesmatrix(space(v[2]), domain(v[1]), ncoefficients(v[1]), ncoefficients(v[2]))
        c_vals12 = C12*v[2].coefficients
        h = 0.00000001; @test stieltjes(v[2],-h) ≈ stieltjes(v[2],RiemannDual(0.0,-1))(h) atol=1E-6
        @test finitepart(stieltjes(v[2],RiemannDual(0.0,-1))) ≈ c_vals12[1]

        C = fpstieltjesmatrix(space(f), ns, ns)
        @test norm(C) ≤ 200

        @test C[1:44,1:2:end] == C11

        c_vals = C*coefficients(f)
        pts = RiemannHilbert.collocationpoints(space(f), ns)

        f_ex = Fun(x->exp(-40(x-0.1)^2), Legendre())
        @test stieltjes(f_ex, Directed{false}(0.0)) ≈ finitepart(stieltjes(f,Directed{false}(RiemannDual(0.0,-im))))
        @test finitepart(stieltjes(v[1],RiemannDual(0.0,-im))+stieltjes(v[2],RiemannDual(0.0,-im))) ≈ stieltjes(f_ex,Directed{false}(0.0))

        @test c_vals[1] ≈ finitepart(stieltjes(f,Directed{false}(RiemannDual(0.0,-im))))
        @test c_vals[1] ≈ finitepart(stieltjes(f,Directed{false}(RiemannDual(0.0,exp(-0.1im)))))
        @test c_vals[1] ≈ finitepart(stieltjes(f,RiemannDual(0.0,-im)))
        @test c_vals[2:ns[1]-1] ≈  stieltjes.(f,pts[2:ns[1]-1]⁻)
        @test c_vals[ns[1]] ≈ finitepart(stieltjes(f,RiemannDual(-1.0,-1.0)))

        @test c_vals[ns[1]+1] ≈ finitepart(stieltjes(f,RiemannDual(1.0,1.0)))
        @test c_vals[ns[1]+2:end-1] ≈  stieltjes.(f,pts[ns[1]+2:end-1]⁻)
        @test c_vals[end] ≈ finitepart(stieltjes(f,RiemannDual(0.0,-im)))

        h =0.00001
        @test stieltjes(v[1], Directed{false}(RiemannDual(0.0,-1.0))) ≈ stieltjes(v[1], RiemannDual(0.0,-1.0-eps()*im))

        @test finitepart(stieltjes(v[1], Directed{false}(RiemannDual(0.0,-1.0)))+ stieltjes(v[2], RiemannDual(0.0,-1.0))) ≈
            stieltjes(v[1], -0.00000000001im)+stieltjes(v[2], -0.00000000001im)

        @test finitepart(stieltjes(v[1], Directed{false}(RiemannDual(0.0,-1.0)))+ stieltjes(v[2], RiemannDual(0.0,-1.0))) ≈
            stieltjes(f, -0.00000000001im)

        @test c_vals[1] ≈ stieltjes(f, -0.0000000001im)
        @test c_vals[end] ≈ stieltjes(f, -0.0000000001im)

        @test finitepart(stieltjes(f, RiemannDual(0.0,-im))) ≈
            finitepart(stieltjes(v[1], RiemannDual(0.0,-im))+stieltjes(v[2], RiemannDual(0.0,-im)))

        C11 = fpstieltjesmatrix(space(v[1]), ncoefficients(v[1]), ncoefficients(v[1]))
        C12 = fpstieltjesmatrix(space(v[2]), domain(v[1]), ncoefficients(v[1]), ncoefficients(v[2]))

        @test finitepart(stieltjes(v[1], Directed{false}(RiemannDual(0.0,-1.0)))) ≈ (C11*coefficients(v[1]))[1]
        @test finitepart(stieltjes(v[2], RiemannDual(0.0,-1.0))) ≈ dotu(stieltjesmoment!(Array{ComplexF64}(undef,ncoefficients(v[2])), space(v[2]), orientedrightendpoint(domain(v[1])), finitepart),
                    coefficients(v[2]))

        @test C12[1,:] ≈ stieltjesmoment!(Array{ComplexF64}(undef,ncoefficients(v[2])), space(v[2]), orientedrightendpoint(domain(v[1])), finitepart)

        @test finitepart(stieltjes(v[2], RiemannDual(0.0,-1.0))) ≈ (C12*coefficients(v[2]))[1]
        @test C[1,:] ≈ interlace(C11[1,:], C12[1,:])
    end

    @testset "0..1 and 0..1" begin
        sp = Legendre(Segment(0 , -1)) ⊕ Legendre(0 .. 1)

        f = Fun(x->sign(x)*exp(-40(x-0.1)^2), sp)
        v = components(f)
        ns = ncoefficients.(v)
        C = fpstieltjesmatrix(space(f), ns, ns)
        @test norm(C) ≤ 200
        c_vals = C*coefficients(f)
        pts = RiemannHilbert.collocationpoints(space(f), ns)

        @test stieltjes(f,-1-eps()) ≈ stieltjes(f,RiemannDual(-1.0,-1.0)).c
        @test c_vals[1] ≈ finitepart(stieltjes(f,RiemannDual(-1.0,-1.0)))
        @test c_vals[2:ns[1]-1] ≈ stieltjes.(f,pts[2:ns[1]-1]⁻)
        @test c_vals[ns[1]] ≈ finitepart(stieltjes(f,RiemannDual(0.0,+im)))
        @test c_vals[ns[1]+1] ≈ finitepart(stieltjes(f,RiemannDual(1.0,1.0)))
        @test c_vals[ns[1]+2:end-1] ≈  stieltjes.(f,pts[ns[1]+2:end-1]⁻)
        @test c_vals[end] ≈ finitepart(stieltjes(f,RiemannDual(0.0,-im)))

        h =0.00001
        @test stieltjes(v[1], Directed{false}(RiemannDual(0.0,-1.0))) ≈ stieltjes(v[1], RiemannDual(0.0,-1.0+eps()*im))
        @test stieltjes(v[1], Directed{true}(RiemannDual(0.0,-1.0))) ≈ stieltjes(v[1], RiemannDual(0.0,-1.0-eps()*im))

        @test stieltjes(v[2], RiemannDual(0.0,-1.0))(h) ≈ stieltjes(v[2], -h) atol=1E-3

        @test stieltjes(v[1], Directed{true}(RiemannDual(0.0,-1.0)))(h) ≈ stieltjes(v[1], -h-eps()*im) atol=1E-3

        @test finitepart(stieltjes(v[1], Directed{true}(RiemannDual(0.0,-1.0)))+ stieltjes(v[2], RiemannDual(0.0,-1.0))) ≈
            stieltjes(v[1], -0.00000000001im)+stieltjes(v[2], -0.00000000001im)

        @test finitepart(stieltjes(v[1], Directed{true}(RiemannDual(0.0,-1.0)))+ stieltjes(v[2], RiemannDual(0.0,-1.0))) ≈
            stieltjes(f, -0.00000000001im)

        @test c_vals[ns[1]] ≈ stieltjes(f, +0.0000000001im)
        @test c_vals[end] ≈ stieltjes(f, -0.0000000001im)
    end

    @testset "-1..0 and 0..1" begin
        sp = Legendre(-1 .. 0) ⊕ Legendre(0 .. 1)
        f = Fun(x->exp(-40(x-0.1)^2), sp)
        v = components(f)
        ns = ncoefficients.(v)
        C = fpstieltjesmatrix(space(f), ns, ns)
        @test norm(C) ≤ 200

        c_vals = C*coefficients(f)
        pts = RiemannHilbert.collocationpoints(space(f), ns)

        @test c_vals[1] ≈ finitepart(stieltjes(f,RiemannDual(0.0,-im)))
        @test c_vals[2:ns[1]-1] ≈  stieltjes.(f,pts[2:ns[1]-1]⁻)
        @test c_vals[ns[1]] ≈ finitepart(stieltjes(f,RiemannDual(-1.0,-1.0)))
        @test c_vals[ns[1]+1] ≈ finitepart(stieltjes(f,RiemannDual(1.0,1.0)))
        @test c_vals[ns[1]+2:end-1] ≈  stieltjes.(f,pts[ns[1]+2:end-1]⁻)
        @test c_vals[end] ≈ finitepart(stieltjes(f,RiemannDual(0.0,-im)))
    end
end

@testset "ArraySpace" begin
    sp = ArraySpace(Legendre() ,2)
    f = Fun(x->[exp(-40(x-0.1)^2); cos(x-0.1)*exp(-40(x-0.1)^2)], sp)
    ns = ncoefficients.(Array(f))
    C = fpstieltjesmatrix(sp, ns, ns)
    C1 = fpstieltjesmatrix(sp[1], ns[1], ns[1])
    C2 = fpstieltjesmatrix(sp[2], ns[2], ns[2])

    @test C*coefficients(f) ≈ [C1*coefficients(f[1]) ; C2*coefficients(f[2])]

    sp = ArraySpace(Legendre(), 2, 2)
    f = Fun(x->[exp(-40(x-0.1)^2) sin(x-0.1)exp(-40(x-0.1)^2); cos(x-0.1)*exp(-40(x-0.1)^2) airyai(x-0.1)*exp(-40(x-0.1)^2)], sp)
    ns = ncoefficients.(Array(f))
    C = fpstieltjesmatrix(sp, ns, ns)
    C11 = fpstieltjesmatrix(sp[1,1], ns[1,1], ns[1,1])
    C21 = fpstieltjesmatrix(sp[2,1], ns[2,1], ns[2,1])
    C12 = fpstieltjesmatrix(sp[1,2], ns[1,2], ns[1,2])
    C22 = fpstieltjesmatrix(sp[2,2], ns[2,2], ns[2,2])

    @test C[1:ns[1], 1:4:end] ≈ C11
    @test norm(C[1:ns[1], 2:4:end]) ≤ 100eps()
    @test norm(C[1:ns[1], 3:4:end]) ≤ 100eps()
    @test norm(C[1:ns[1], 4:4:end]) ≤ 100eps()

    @test C[ns[1]+1:ns[1]+ns[2], 2:4:end] ≈ C21

    @test C*coefficients(f) ≈ [C11*coefficients(f[1,1]) ; C21*coefficients(f[2,1]) ; C12*coefficients(f[1,2]) ; C22*coefficients(f[2,2])]

    sp = ArraySpace(Legendre(-1 .. 0) ∪ Legendre(0 .. 1), 2)
    f = Fun(x->[exp(-40(x-0.1)^2); cos(x-0.1)*exp(-40(x-0.1)^2)], sp)
    ns = ncoefficients.(Array(f))
    C = fpstieltjesmatrix(sp, ns, ns)

    C1 = fpstieltjesmatrix(sp[1], ns[1], ns[1])
    C2 = fpstieltjesmatrix(sp[2], ns[2], ns[2])

    @test C*coefficients(f) ≈ [C1*coefficients(f[1]) ; C2*coefficients(f[2])]

    sp = ArraySpace(Legendre(-1 .. 0) ∪ Legendre(0 .. 1), 2, 2)
    f = Fun(x->[exp(-40(x-0.1)^2) sin(x-0.1)exp(-40(x-0.1)^2); cos(x-0.1)*exp(-40(x-0.1)^2) airyai(x-0.1)*exp(-40(x-0.1)^2)], sp)
    ns = ncoefficients.(Array(f))
    C = fpstieltjesmatrix(sp, ns, ns)

    C11 = fpstieltjesmatrix(sp[1,1], ns[1,1], ns[1,1])
    C21 = fpstieltjesmatrix(sp[2,1], ns[2,1], ns[2,1])
    C12 = fpstieltjesmatrix(sp[1,2], ns[1,2], ns[1,2])
    C22 = fpstieltjesmatrix(sp[2,2], ns[2,2], ns[2,2])

    @test C*coefficients(f) ≈ [C11*coefficients(f[1,1]) ; C21*coefficients(f[2,1]) ; C12*coefficients(f[1,2]) ; C22*coefficients(f[2,2])]
end

@testset "rhsolve" begin
    sp = Legendre()
    g = 1-0.3Fun(x->exp(-40x^2), sp)
    n = 2ncoefficients(g)
    g_v = RiemannHilbert.collocationvalues(g-1, n)
    u = Fun(sp, rhmatrix(g,n) \ g_v)

    @testset "-1 .. 1" begin
        n = 2ncoefficients(g)
        g = pad(g,n)
        C₋ = fpcauchymatrix(sp, n, n)
        pts = RiemannHilbert.collocationpoints(sp, n)
        @test C₋[2:end-1,:]*coefficients(g) ≈ cauchy.(g, pts[2:end-1]⁻)
        g_v = RiemannHilbert.collocationvalues(g-1, n)
        @test g_v ≈ g.(pts).-1
        G = Diagonal(g_v)

        @test G*g_v ≈ (g.(pts) .- 1).^2

        E = RiemannHilbert.evaluationmatrix(sp, pts, length(pts))
        @test E*coefficients(g) ≈ g.(pts)

        @test (g.(pts).-1).*(C₋*coefficients(g)) ≈ (g.(pts).-1).*cauchy.(g, pts.-0.000000001im)
        @test G*(C₋*coefficients(g)) ≈ (g.(pts).-1).*cauchy.(g, pts.-0.000000001im)

        L = E - G*C₋
        @test L*coefficients(g) ≈ g.(pts) - (g.(pts).-1).*cauchy.(g, pts.-0.000000001im)

        @test L == rhmatrix(g,n)

        φ = z -> 1 + cauchy(u,z)
        @test φ(0.1⁺)  ≈ g(0.1)φ(0.1⁻)

        @test 1+cauchy(u)(0.1+0.2im) ≈ φ(0.1+0.2im)
        @test (1+cauchy(u) )(0.1+0.2im) ≈ φ(0.1+0.2im)
        @test (1+cauchy(u) )(0.1⁻) ≈ φ(0.1⁻)
    end

    @testset "-1 .. 0 and 0 .. 1" begin
        sp = Legendre(-1 .. 0) ∪ Legendre(0 .. 1)
        u_1 = u
        u_ex = Fun(x->u_1(x), sp)
        g_1 = 1-0.3Fun(x->exp(-40x^2), Legendre())

        g = 1-0.3Fun(x->exp(-40x^2), sp)

        n = 2ncoefficients(g)
        g = pad(g,n)
        u_ex = pad(u_ex,n)
        @test (1+cauchy(u_ex,0.1⁺)) ≈ g(0.1)*(1+cauchy(u_ex,0.1⁻))

        C₋ = fpcauchymatrix(sp, n, n)
        pts = RiemannHilbert.collocationpoints(sp, n)

        @test transpose(C₋[1,:])*coefficients(u_ex) ≈ cauchy(u_ex,pts[1]-eps()im)
        @test transpose(C₋[5,:])*coefficients(u_ex) ≈ cauchy(u_ex,pts[5]-eps()im)
        @test transpose(C₋[end,:])*coefficients(u_ex) ≈ cauchy(u_ex, 0.0-eps()im)
        @test transpose(C₋[end-1,:])*coefficients(u_ex) ≈ cauchy(u_ex,pts[end-1]-eps()im)
        @test transpose(C₋[n÷2,:])*coefficients(u_ex) ≈ cauchy(u_ex,-1.0-eps())
        @test transpose(C₋[(n÷2)+1,:])*coefficients(u_ex) ≈ cauchy(u_ex,1.0+eps())

        @test C₋*coefficients(u_ex) ≈ cauchy.(u_ex, pts.-eps()im)
        g_v = RiemannHilbert.collocationvalues(g-1, n)
        @test g_v ≈ g.(pts).-1
        G = Diagonal(g_v)

        @test G*g_v ≈ (g.(pts) .- 1).^2

        g1 = component(g,2)
        pts1 = RiemannHilbert.collocationpoints(component(sp,2), n÷2)
        E1=RiemannHilbert.evaluationmatrix(component(sp,1), length(pts1))

        @test E1*g1.coefficients ≈ g1.(pts1)

        E = RiemannHilbert.evaluationmatrix(sp, n)
        @test E*coefficients(g) ≈ g.(pts)

        @test (g.(pts).-1).*(C₋*coefficients(g)) ≈ (g.(pts).-1).*cauchy.(g, pts.-0.000000001im)
        @test G*(C₋*coefficients(g)) ≈ (g.(pts).-1).*cauchy.(g, pts.-0.000000001im)

        L = E - G*C₋
        @test L*coefficients(g) ≈ g.(pts) - (g.(pts).-1).*cauchy.(g, pts.-0.000000001im)

        @test L == rhmatrix(g,n)

        u = Fun(sp, rhmatrix(g,n) \ g_v)
        φ = z -> 1 + cauchy(u,z)
        @test φ(0.1⁺)  ≈ g(0.1)φ(0.1⁻)

        @test 1+cauchy(u)(0.1+0.2im) ≈ φ(0.1+0.2im)
        @test (1+cauchy(u) )(0.1+0.2im) ≈ φ(0.1+0.2im)
        @test (1+cauchy(u) )(0.1⁻) ≈ φ(0.1⁻)

        @time φ = rhsolve(g, 2ncoefficients(g))
        @test φ(0.1⁺)  ≈ g(0.1)φ(0.1⁻)
    end

    @testset "0..-1 and 0..1" begin
        sp = Legendre(Segment(0 , -1)) ∪ Legendre(0 .. 1)
        g = Fun(x->x ≥ 0 ? 1-0.3exp(-40x^2) : inv(1-0.3exp(-40x^2)), sp)
        u_1 = u
        u_ex = Fun(x->sign(x)u_1(x), sp)

        n = 2ncoefficients(g)
        g = pad(g,n)
        u_ex = pad(u_ex,n)
        @test (1+cauchy(u_ex,0.1⁺)) ≈ g(0.1)*(1+cauchy(u_ex,0.1⁻))

        C₋ = fpcauchymatrix(sp, n, n)
        pts = RiemannHilbert.collocationpoints(sp, n)

        @test transpose(C₋[1,:])*coefficients(u_ex) ≈ cauchy(u_ex,pts[1]+eps()im)
        @test transpose(C₋[5,:])*coefficients(u_ex) ≈ cauchy(u_ex,pts[5]+eps()im)
        @test transpose(C₋[n÷2,:])*coefficients(u_ex) ≈ cauchy(u_ex,0.0+eps()im)
        @test transpose(C₋[(n÷2)+1,:])*coefficients(u_ex) ≈ cauchy(u_ex,1.0+eps())
        @test transpose(C₋[end,:])*coefficients(u_ex) ≈ cauchy(u_ex, pts[end]-eps()im)
        @test transpose(C₋[end,:])*coefficients(u_ex) ≈ cauchy(u_ex, 0.0-eps()im)
        @test transpose(C₋[end-1,:])*coefficients(u_ex) ≈ cauchy(u_ex,pts[end-1]-eps()im)

        @test C₋*coefficients(u_ex) ≈ [cauchy.(u_ex, pts[1:n÷2].+eps()im); cauchy.(u_ex, pts[(n÷2)+1:end].-eps()im)]
        g_v = RiemannHilbert.collocationvalues(g-1, n)

        g1 = component(g,1)
        pts1 = RiemannHilbert.collocationpoints(component(sp,1), n÷2)
        E1=RiemannHilbert.evaluationmatrix(component(sp,1), length(pts1))

        @test E1*g1.coefficients ≈ g1.(pts1)

        g2 = component(g,2)
        pts2 = RiemannHilbert.collocationpoints(component(sp,2), n÷2)

        @test g_v ≈ [g1.(pts1).-1; g2.(pts2).-1]
        G = Diagonal(g_v)

        @test g1.(pts1) ≈ g.(pts1)
        @test_broken g2.(pts2) ≈ g.(pts2) # evaluation at 0
        @test_broken G*g_v ≈ (g.(pts) .- 1).^2

        E = RiemannHilbert.evaluationmatrix(sp, n)
        @test E[1:end-1,:]*coefficients(g) ≈ g.(pts[1:end-1])
        @test transpose(E[end,:])*coefficients(g) ≈ component(g,2)(0.0)

        @test (g2.(pts2).-1).*(C₋[(n÷2)+1:end,:]*coefficients(u_ex)) ≈
            (g2.(pts2).-1).*cauchy.(u_ex, pts2.-eps()im)
        @test (g.(pts1).-1).*(C₋[1:n÷2,:]*coefficients(u_ex)) ≈
            (g.(pts1).-1).*cauchy.(u_ex, pts1.+0.000000001im)
        @test (g.(pts2).-1).*(C₋[(n÷2)+1:end,:]*coefficients(u_ex)) ≈
            (g.(pts2).-1).*cauchy.(u_ex, pts2.-0.000000001im)

        u_ex1,u_ex2 = components(u_ex)
        @test (C₋*coefficients(u_ex))[ncoefficients(u_ex1)+1:end] ≈ cauchy.(u_ex, pts2.-eps()im)
        @test G*C₋*coefficients(u_ex) ≈ [(g1.(pts1).-1).*cauchy.(u_ex, pts1.+eps()im);
                                         (g2.(pts2).-1).*cauchy.(u_ex, pts2.-eps()im)]

        L = E - G*C₋
        @test L*coefficients(u_ex) ≈ [u_ex1.(pts1) - (g1.(pts1).-1).*cauchy.(u_ex, pts1.+eps()im);
                                        u_ex2.(pts2) - (g2.(pts2).-1).*cauchy.(u_ex, pts2.-eps()im)]

        @test L == rhmatrix(g,n)

        u = Fun(sp, rhmatrix(g,n) \ g_v)
        φ1 = z -> 1 + cauchy(u,z)
        @test φ1(0.1⁺)  ≈ g(0.1)φ1(0.1⁻)
        @test φ1(-0.1⁺)  ≈ g(-0.1)φ1(-0.1⁻)

        @test 1+cauchy(u)(0.1+0.2im) ≈ φ1(0.1+0.2im)
        @test (1+cauchy(u) )(0.1+0.2im) ≈ φ1(0.1+0.2im)
        @test (1+cauchy(u) )(0.1⁻) ≈ φ1(0.1⁻)
        @test (1+cauchy(u) )(0.1⁺) ≈ φ1(0.1⁺)
        @test (1+cauchy(u) )(-0.1⁻) ≈ φ1(-0.1⁻)
        @test (1+cauchy(u) )(-0.1⁺) ≈ φ1(-0.1⁺)

        @time φ = rhsolve(g, 2ncoefficients(g))
        @test φ(0.1⁺)  ≈ g(0.1)φ(0.1⁻)
    end

    @testset "Chebyshev g" begin
        G = Fun(x -> 1 + exp(-30x^2), -1..1)
        Φ = rhsolve(G, 200)
        @test Φ(0.1⁺) ≈ G(0.1)Φ(0.1⁻)
    end
end

@testset "Matrix rhsolve" begin
    sp = ArraySpace(Legendre(), 2)
    f = Fun(Fun(x->[cos(x);sin(x)], Chebyshev()), ArraySpace(Legendre(), 2))
    G = Fun(Fun(x->[1 exp(-40x^2); 0.1exp(-40x^2) 1], Chebyshev()), ArraySpace(Legendre(), 2, 2))

    n = 2ncoefficients(G)
    E = RiemannHilbert.evaluationmatrix(sp, n)
    pts = RiemannHilbert.collocationpoints(sp, n÷2)

    @test E*coefficients(pad(f,n)) ≈ [f[1].(pts); f[2].(pts)]

    M = RiemannHilbert.multiplicationmatrix(G-I, n)
    @test M*(E*coefficients(pad(f,n))) ≈ mapreduce(f -> f.(pts), vcat, (G-I)*f)

    L = E - M*fpcauchymatrix(sp, n,n)
    u1 = L \ mapreduce(f -> f.(pts), vcat, (G-I)[:,1])
    u2 = L \ mapreduce(f -> f.(pts), vcat, (G-I)[:,2])
    z=2+I
        Φ = z -> I + [cauchy(Fun(sp, u1),z) cauchy(Fun(sp, u2),z)]

    @test Φ(0.1+eps()im) ≈ G(0.1)*Φ(0.1-eps()im)

    @test RiemannHilbert.rhmatrix(G, n) == L

    @test RiemannHilbert.collocationvalues(G-I, n) ≈
        [RiemannHilbert.collocationvalues((G-I)[:,1], n) RiemannHilbert.collocationvalues((G-I)[:,2], n)]

    Φ = rhsolve(G,n)
    @test G(0.1)*Φ(0.1⁻) ≈ Φ(0.1⁺)
end

@testset "Matrix rhsolve two interval" begin
    sp = ArraySpace(Legendre(-1 .. 0) ∪ Legendre(0 .. 1), 2)
    G = Fun(x->[1 exp(-40x^2); 0.1exp(-40x^2) 1], ArraySpace(sp[1], 2, 2))

    n=2ncoefficients(G)
    Φ = rhsolve(G,n)

    @test G(0.1)*Φ(0.1⁻) ≈ Φ(0.1⁺)
    U = RiemannHilbert.rh_sie_solve(G,n)

    U1 = U[:,1]
    g = G
    sp = space(g)[:,1]
    C₋ = fpcauchymatrix(sp, n, n)
    G = RiemannHilbert.multiplicationmatrix(g-I, n)
    E = RiemannHilbert.evaluationmatrix(sp, n)

    pts = RiemannHilbert.collocationpoints(sp, n÷2)
    @test last(E*coefficients(U1)) ≈ U1(0.00000000001)[2]
    @test last(C₋*coefficients(U1)) ≈ cauchy(U1[2], 0.00000000001-eps()*im)
    @test cauchy(U1[1], 0.00000000001-eps()*im) ≈ (C₋*coefficients(U1))[n÷2]
    @test last(G*C₋*coefficients(U1)) ≈ ((g(0.00000000001)-I)*cauchy(U1, 0.00000000001-eps()*im))[2]

    @test RiemannHilbert.collocationvalues((g-I)[:,1], n)[end] ≈ 0.1

    sp = ArraySpace(Legendre(Segment(0 , -1)) ∪ Legendre(0 .. 1), 2)
    G = g = Fun(x-> x ≥ 0 ? [1 exp(-40x^2); 0.1exp(-40x^2) 1] : inv([1 exp(-40x^2); 0.1exp(-40x^2) 1]), ArraySpace(sp[1], 2, 2))
    n=2ncoefficients(G)
    Ũ1 = Fun(x-> x ≥ 0 ? U1(x) : -U1(x), sp, n÷2)
    @test cauchy(Ũ1, 0.1⁻) ≈ cauchy(Ũ1, 0.1-eps()im)
    @test cauchy(Ũ1, 0.1+eps()im) ≈ cauchy(U1, 0.1+eps()im)
    @test cauchy(Ũ1, -0.1+eps()im) ≈ cauchy(U1, -0.1+eps()im)

    L = rhmatrix(G, n)
    rhs = RiemannHilbert.collocationvalues((g-I)[:,1], n)
    @test rhs[end] ≈ 0.1

    Φ = rhsolve(G,n)

    g = G
    sp = space(g)[:,1]
    C₋ = fpcauchymatrix(sp, n, n)
    G = RiemannHilbert.multiplicationmatrix(g-I, n)
    E = RiemannHilbert.evaluationmatrix(sp, n)

    pts = RiemannHilbert.collocationpoints(sp, n÷2)

    @test last(E*coefficients(Ũ1)) ≈ Ũ1(0.00000000001)[2]
    @test last(C₋*coefficients(Ũ1)) ≈ cauchy(Ũ1[2], 0.00000000001-eps()*im)
    @test cauchy(Ũ1[1], 0.00000000001-eps()*im) ≈ (C₋*coefficients(Ũ1))[n÷2]
    @test G[[n÷2,n],[n÷2,n]] ≈ (g(0.00000000001)-I)

    @test last(G*C₋*coefficients(Ũ1)) ≈ transpose((g(0.00000000001)-I)[2,:])*cauchy(Ũ1, 0.00000000001-eps()*im)

    @test transpose(L[end,:])*coefficients(Ũ1) ≈ Ũ1(0.00000000001)[2] - transpose((g(0.00000000001)-I)[2,:])*cauchy(Ũ1, 0.00000000001-eps()*im)

    @test U1(0.00000000001)[2] - transpose((g(0.00000000001)-I)[2,:])*cauchy(U1, 0.00000000001-eps()*im) ≈ 0.1
    @test g(0.1)*Φ(0.1⁻) ≈ Φ(0.1⁺)
end

@testset "4 rays" begin
    Γ = Segment(0, 2.5exp(im*π/6)) ∪ Segment(0, 2.5exp(5im*π/6)) ∪
            Segment(0, 2.5exp(-5im*π/6)) ∪ Segment(0, 2.5exp(-im*π/6))
    sp = ArraySpace(PiecewiseSpace(Legendre.(components(Γ))), 2,2)

    s₁ = im
    s₃ = -im

    G = Fun( z -> if angle(z) ≈ π/6
                        [1 0; s₁*exp(8im/3*z^3) 1]
                    elseif angle(z) ≈ 5π/6
                        [1 0; s₃*exp(8im/3*z^3) 1]
                    elseif angle(z) ≈ -π/6
                        [1 -s₃*exp(-8im/3*z^3); 0 1]
                    elseif angle(z) ≈ -5π/6
                        [1 -s₁*exp(-8im/3*z^3); 0 1]
                    end
                        , sp)

    Φ = transpose(rhsolve(transpose(G), 2*4*100))

    s=exp(im*π/6)
    @test Φ((s)⁻)*G(s) ≈ Φ((s)⁺)
    @test map(g->first(components(g)[1]), G)*map(g->first(components(g)[2]), G)*
        map(g->first(components(g)[3]), G)*map(g->first(components(g)[4]), G) ≈ Matrix(I,2,2)

    for _=1:5
        @test RiemannHilbert.evaluationmatrix(sp[:,1], 32) ≈ RiemannHilbert.evaluationmatrix(sp[:,1], 32)
        @test RiemannHilbert.rhmatrix(G, 32) ≈ RiemannHilbert.rhmatrix(G, 32)
    end

    @test RiemannHilbert.rhmatrix(transpose(G), 900) ≈ RiemannHilbert.rhmatrix(transpose(G), 900)
    U = RiemannHilbert.rh_sie_solve(transpose(G), 2*4*100)
    @test -0.36706155154807807 ≈ sum(U[1,2])/(-π*im)

    G̃ = Fun( z -> if angle(z) ≈ π/6
                        [1 0; s₁*exp(8im/3*z^3) 1]
                    elseif angle(z) ≈ 5π/6
                        [1 0; s₃*exp(8im/3*z^3) 1]
                    elseif angle(z) ≈ -π/6
                        [1 -s₃*exp(-8im/3*z^3); 0 1]
                    elseif angle(z) ≈ -5π/6
                        [1 -s₁*exp(-8im/3*z^3); 0 1]
                    end
                        , Γ)

    Ũ = RiemannHilbert.rh_sie_solve(transpose(G̃), 2*4*100)

    @test rhmatrix(transpose(G), 2*4*100) ≈ rhmatrix(transpose(G̃), 2*4*100)
    @test -0.36706155154807807 ≈ sum(Ũ[1,2])/(-π*im)

    V = istieltjes(Φ)
    x = Fun(domain(V))

    @test 10.0I + sum.(Array(V)) + stieltjes(x*V, 10.0) ≈ 10.0Φ(10.0)

    # check rhmatrix relationship for below
    U = V*(-2π*im)
    U1 = U[1,:]
    n= ncoefficients(U1)
    L = rhmatrix(transpose(G),n)
    vals = collocationvalues(transpose(G)-I, n)
    @test L*coefficients(U1) ≈ vals[:,1]
end

@testset "6 rays" begin
    @testset "HM on 6 rays" begin
        s₁,s₂,s₃ = -im,0,im
        @assert s₁ - s₂ + s₃ + s₁*s₂*s₃ ≈ 0

        # construct true solution using 4 rays
        Γ = Segment(0, 2.5exp(im*π/6)) ∪ Segment(0, 2.5exp(5im*π/6)) ∪
        Segment(0, 2.5exp(-5im*π/6)) ∪ Segment(0, 2.5exp(-im*π/6))
        sp = ArraySpace(PiecewiseSpace(Legendre.(components(Γ))), 2,2)

        G = Fun( z -> if angle(z) ≈ π/6
                            [1 0; s₁*exp(8im/3*z^3) 1]
                        elseif angle(z) ≈ 5π/6
                            [1 0; s₃*exp(8im/3*z^3) 1]
                        elseif angle(z) ≈ -π/6
                            [1 -s₃*exp(-8im/3*z^3); 0 1]
                        elseif angle(z) ≈ -5π/6
                            [1 -s₁*exp(-8im/3*z^3); 0 1]
                        end
                            , sp)

        Φ = transpose(rhsolve(transpose(G), 2*6*100))

        # debug
        V4 = istieltjes(Φ)
        x = 0.0

        @test stieltjes(V4,1+im)+I == Φ(1+im)

        Γ = Segment(0, 2.5exp(im*π/6))   ∪
        Segment(0, 2.5exp(im*π/2))       ∪
        Segment(0, 2.5exp(5im*π/6))      ∪
        Segment(0, 2.5exp(-5im*π/6))     ∪
        Segment(0, 2.5exp(-im*π/2))      ∪
        Segment(0, 2.5exp(-im*π/6));

        G = Fun( z -> if angle(z) ≈ π/6
                        [1                             0;
                        s₁*exp(8im/3*z^3+2im*x*z)     1]
                    elseif angle(z) ≈ π/2
                        [1                 s₂*exp(-8im/3*z^3-2im*x*z);
                        0                 1]
                    elseif angle(z) ≈ 5π/6
                        [1                             0;
                        s₃*exp(8im/3*z^3+2im*x*z)     1]
                    elseif angle(z) ≈ -π/6
                        [1                -s₃*exp(-8im/3*z^3-2im*x*z);
                        0                 1]
                    elseif angle(z) ≈ -π/2
                        [1                             0;
                        -s₂*exp(8im/3*z^3+2im*x*z)    1]
                    elseif angle(z) ≈ -5π/6
                        [1                -s₁*exp(-8im/3*z^3-2im*x*z);
                        0                 1]
                    end
                        , Γ);
        sp = ArraySpace(PiecewiseSpace(Legendre.(components(Γ))), 2,2)
        V = Fun(V4, sp)

        @test stieltjes(V,1+im)+I ≈ Φ(1+im)

        U = V*(-2π*im)

        U1 = U[1,:]
        @test cauchy(U1,1+im)+[1,0] ≈ Φ(1+im)[1,:]
        @test abs(sum(first.(components(U1[1])))) ≤ 100eps()
        @test abs(sum(first.(components(U1[2])))) ≤ 100eps()

        U11 = U1[1]
        n = ncoefficients(U11)
        C₋ = fpcauchymatrix(space(U11), n, n)
        pts = collocationpoints(space(U11), n)
        c = C₋*coefficients(U11)

        @test c[1] ≈ finitepart(cauchy(U11,orientedrightendpoint(component(Γ,1))))
        @test c[2] ≈ cauchy(U11,pts[2])
        @test c[150] ≈ finitepart(cauchy(U11,orientedleftendpoint(component(Γ,1))⁻))
        @test c[151] ≈ finitepart(cauchy(U11,orientedrightendpoint(component(Γ,2))))

        n = ncoefficients(U1)
        L = rhmatrix(transpose(G),n)
        vals = collocationvalues(transpose(G)-I, n)
        pts = collocationpoints(space(U1),n)

        c = L*coefficients(U1)
        @test c ≈ vals[:,1]

        @test coefficients(U1) ≈ L \ vals[:,1]
    end
end

include("test_nls.jl")


text:
# http://www.cs.miami.edu/home/burt/learning/Csc524.092/notes/ip_example.html
const dns_packet = [
    0x00, 0x18, 0x8b, 0x75, 0x1d, 0xe0, 0x00, 0x1f, 0xf3, 0xd8, 0x47, 0xab, 0x08, 0x00, 0x45, 0x00,
    0x00, 0x44, 0xad, 0x0b, 0x00, 0x00, 0x40, 0x11, 0x72, 0x72, 0xac, 0x14, 0x02, 0xfd, 0xac, 0x14,
    0x00, 0x06, 0xe5, 0x87, 0x00, 0x35, 0x00, 0x30, 0x5b, 0x6d, 0xab, 0xc9, 0x01, 0x00, 0x00, 0x01,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x09, 0x6d, 0x63, 0x63, 0x6c, 0x65, 0x6c, 0x6c, 0x61, 0x6e,
    0x02, 0x63, 0x73, 0x05, 0x6d, 0x69, 0x61, 0x6d, 0x69, 0x03, 0x65, 0x64, 0x75, 0x00, 0x00, 0x01,
    0x00, 0x01
]

ext: jl

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols


text:
using Interact

function handle_msg(w::InputWidget, msg)
    if msg.content["data"]["method"] == "backbone"
        Interact.recv(w, msg.content["data"]["sync_data"]["value"])
    end
end

function handle_msg(w::Button, msg)
    try
        if msg.content["data"]["method"] == "custom" &&
            msg.content["data"]["content"]["event"] == "click"
            # click event occured
            push!(w.input, nothing)
        end
    catch e
        warn(string("Couldn't handle Button message ", e))
    end
end

function handle_msg{view}(w::Options{view}, msg)
    try
        if msg.content["data"]["method"] == "backbone"
            key = msg.content["data"]["sync_data"]["value_name"]
            if haskey(w.options, key)
                Interact.recv(w, w.options[key])
            end
        end
    catch e
        warn(string("Couldn't handle ", view, " message ", e))
    end
end


text:
using StochasticDiffEq, Test

SRIW1_tab = constructSRIW1()
@test minimum(checkSRIOrder(SRIW1_tab))

SRA1_tab = constructSRA1()
@test minimum(checkSRAOrder(SRA1_tab))

DRI1_tab = constructDRI1()
@test minimum(checkRIOrder(DRI1_tab))

RI1_tab = constructRI1()
@test minimum(checkRIOrder(RI1_tab))

RI3_tab = constructRI3()
@test minimum(checkRIOrder(RI3_tab))

RI5_tab = constructRI5()
@test minimum(checkRIOrder(RI5_tab))

RI6_tab = constructRI6()
@test minimum(checkRIOrder(RI6_tab))

RDI1WM_tab = constructRDI1WM()
@test minimum(checkRIOrder(RDI1WM_tab, ps=1))

RDI2WM_tab = constructRDI2WM()
@test minimum(checkRIOrder(RDI2WM_tab))

RDI3WM_tab = constructRDI3WM()
@test minimum(checkRIOrder(RDI3WM_tab))

RDI4WM_tab = constructRDI4WM()
@test minimum(checkRIOrder(RDI4WM_tab))

RS1_tab = constructRS1()
@test minimum(checkRSOrder(RS1_tab))

RS2_tab = constructRS2()
@test minimum(checkRSOrder(RS2_tab))

# Komori's NON
NON_tab = constructNON()
@test minimum(checkNONOrder(NON_tab))

NON2_tab = constructNON2()
@test minimum(checkNONOrder(NON2_tab))


text:
# This is only here to define a namespace
module JuliaTutorial

###############################################################################
# BASICS OF TYPES AND PACKAGES
###############################################################################
# # This is how every packages that you will use work udner the hood.

# # 1/ Use the function typeof() to get the type of anything you like

# # 2/ Investigate the type hierarchy:

# 1 isa Real
# 1 isa Integer
# Integer <: Real

# # 3/ Multiple dispatch: Define a function to add 2 numbers

# function myaddition(a, b)
#     return a + b
# end

# myaddition(1, 2)

# function myaddition(a, b::AbstractFloat)
#     return a / b
# end

# myaddition(1, 2.)

# # 4/ Defining a new type and extending our function

struct Cocktail
    ingredients::Tuple
end

# sexonthebeach = Cocktail(("vodka", "schnaps", "orange juice", "cranberry juice"))
# tipunch = Cocktail(("rhum", "sugar", "lime juice"))

# specialised function
function myaddition(a::Cocktail, b::Cocktail)
    return Cocktail((a.ingredients..., b.ingredients...))
end

# olivierssecret = myaddition(sexonthebeach, tipunch)

# olivierssecret isa Cocktail
# # Make your great function `myaddition` available to the public

# c = (1, 1, 2, 5, 6)
# function toto(a, b, args...;d=1)
#     return d
# end

# toto(c...;d=10)

# to make available outside the module (e.g. available to the test function)
export myaddition, Cocktail

# 5/ Bonus: move the methods into `functions.jl` and include the file in this module

# Let's move to the ML part --> `machinelearning.jl`

end


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, LinearAlgebra, Plots, NonlinearSolve
using ProgressMeter: @showprogress

begin
    x0 = 0
    x1 = 1
    ncell = 100
    dx = (x1 - x0) / ncell
    deg = 2 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.05
    dt = cfl * dx / 2
    t = 0.0
end
ps = FRPSpace1D(x0, x1, ncell, deg)
VInv = inv(Array(ps.V))
ℓ = FR.basis_norm(ps.deg)

u = zeros(ncell, nsp, 3)
for i = 1:ncell, ppp1 = 1:nsp
    if ps.x[i] <= 0.5
        prim = [1.0, 0.0, 0.5]
    else
        prim = [0.125, 0.0, 0.625]
    end

    u[i, ppp1, :] .= prim_conserve(prim, γ)
end

function dudt!(du, u, p, t)
    du .= 0.0
    nx, nsp, J, ll, lr, lpdm, dgl, dgr, γ = p

    for i = 1:nx, j = 1:nsp
        prim = conserve_prim(u[i, j, :], γ)
        if prim[1] <= 0 || prim[end] <= 0
            @show i, j
        end
    end

    ncell = size(u, 1)
    nsp = size(u, 2)

    f = zeros(ncell, nsp, 3)
    for i = 1:ncell, j = 1:nsp
        f[i, j, :] .= euler_flux(u[i, j, :], γ)[1] ./ J[i]
    end

    u_face = zeros(ncell, 3, 2)
    f_face = zeros(ncell, 3, 2)
    for i = 1:ncell
        for j = 1:3
            # right face of element i
            u_face[i, j, 1] = dot(u[i, :, j], lr)
            f_face[i, j, 1] = dot(f[i, :, j], lr)

            # left face of element i
            u_face[i, j, 2] = dot(u[i, :, j], ll)
            f_face[i, j, 2] = dot(f[i, :, j], ll)
        end
    end

    f_interaction = zeros(nx + 1, 3)
    for i = 2:nx
        fw = @view f_interaction[i, :]
        flux_hll!(fw, u_face[i-1, :, 1], u_face[i, :, 2], γ, 1.0)
    end

    rhs1 = zeros(ncell, nsp, 3)
    for i = 1:ncell, ppp1 = 1:nsp, k = 1:3
        rhs1[i, ppp1, k] = dot(f[i, :, k], lpdm[ppp1, :])
    end

    idx = 2:ncell-1 # ending points are Dirichlet
    for i in idx, ppp1 = 1:nsp, k = 1:3
        du[i, ppp1, k] =
            -(
                rhs1[i, ppp1, k] +
                (f_interaction[i, k] / J[i] - f_face[i, k, 2]) * dgl[ppp1] +
                (f_interaction[i+1, k] / J[i] - f_face[i, k, 1]) * dgr[ppp1]
            )
    end
end

tspan = (0.0, 0.15)
p = (ps.nx, ps.deg + 1, ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, γ)
prob = ODEProblem(dudt!, u, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    for i = 1:ps.nx
        ũ = VInv * itg.u[i, :, 1]
        su = ũ[end]^2 / sum(ũ.^2)
        isShock = shock_detector(log10(su), ps.deg)

        if isShock
            for s = 1:size(itg.u, 3)
                λ = 5 * dt * sqrt(su) #/ (ps.xpg[i, 2] - ps.xpg[i, 1])
                #λ = 5 * dt * max(abs(dot(itg.u[i, :, 1], ps.dl[1, :])), abs(dot(itg.u[i, :, 1], ps.dl[2, :]))) / ps.J[i] #/ itg.u[i, 1, 1]
                û = VInv * itg.u[i, :, s]
                #FR.modal_filter!(û, λ; filter = :l2)
                #FR.modal_filter!(û, 1.5e-2; filter = :l2)
                #FR.modal_filter!(û, 2e-2; filter = :l2opt)
                #FR.modal_filter!(û, 1e-2, 1e-3; filter = :l2)
                FR.modal_filter!(û, ℓ; filter = :lasso)

                itg.u[i, :, s] .= ps.V * û
            end
        end

        tmp = @view itg.u[i, :, :]
        positive_limiter(tmp, γ, ps.wp ./ 2, ps.ll, ps.lr)
    end

    step!(itg)
end

begin
    x = zeros(ncell * nsp)
    sol = zeros(ncell * nsp, 3)
    for i = 1:ncell
        idx0 = (i - 1) * nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]

            sol[idx, :] .= conserve_prim(itg.u[i, j, :], γ)
            sol[idx, 3] = 1 / sol[idx, 3]
        end
    end
    plot(x, sol[:, :])
end

#plot(x, sol[:, 2])

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
# Model
include_model("double_integrator")

function fd(model::DoubleIntegrator{Discrete, FixedTime}, x⁺, x, u, w, h, t)
    x⁺ - [x[1] + x[2] + w[1]; x[2] + (w[2] + 1.0) * u[1]]
end

function fd(model::DoubleIntegrator{Discrete, FixedTime}, x, u, w, h, t)
    [x[1] + x[2] + w[1]; x[2] + (w[2] + 1.0) * u[1]]
end

# Horizon
T = 101

# Time step
h = 1.0

# Initial and final states
x1 = [1.0; 0.0]
xT = [0.0; 0.0]

xl, xu = state_bounds(model, T, x1 = x1, xT = xT)
_ul = -0.01
_uu = 0.01
ul, uu = control_bounds(model, T, _ul, _uu)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(ones(model.n)) : Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(100.0 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T],
		[zeros(model.m) for t = 1:T])

# Problem
prob = trajectory_optimization_problem(model,
			   obj,
			   T,
               xl = xl,
               xu = xu,
			   ul = ul,
			   uu = uu)

# Initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation for states
u0 = [ones(model.m) for t = 1:T-1]

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
@time z, info = solve(prob, copy(z0))

# Visualize
using Plots
x, u = unpack(z, prob)
plot(hcat(x...)', width = 2.0)
plot(hcat(u...)', width = 2.0, linetype = :steppost)

# Model-predictive control

function run_mpc()
	T_mpc = 25
	x_ref = [x...,[x[end] for i = 1:T_mpc]...]
	u_ref = [u...,[u[end] for i = 1:T_mpc]...]
	Q_mpc = Diagonal(10.0 * ones(model.n))
	R_mpc = Diagonal(1.0 * ones(model.m))
	K, P = tvlqr(model, x, u, h, [Q_mpc for t = 1:T], [R_mpc for t = 1:T-1])
	P_mpc = [P..., [P[end] for i = 1:T_mpc]...]

	shift = 1

	x_hist = [x_ref[shift] + 1.0e-1 * randn(model.n)]
	u_hist = []

	for t = 1:T-1
		println("T = $t")
		xl_mpc, xu_mpc = state_bounds(model, T_mpc, x1 = x_hist[end])
		ul_mpc, uu_mpc = control_bounds(model, T_mpc, _ul, _uu)
		# Objective
		obj_mpc = quadratic_tracking_objective(
		        [t < T_mpc ? Q_mpc : P_mpc[T_mpc + shift - 1] for t = 1:T_mpc],
		        [R_mpc for t = 1:T_mpc-1],
		        [x_ref[shift + t - 1] for t = 1:T_mpc],
				[u_ref[shift + t - 1] for t = 1:T_mpc-1])

		# Problem
		prob_mpc = trajectory_optimization_problem(model,
					   obj_mpc,
					   T_mpc,
		               xl = xl_mpc,
		               xu = xu_mpc,
					   ul = ul_mpc,
					   uu = uu_mpc)

		# Pack trajectories into vector
		z0_mpc = pack(x_ref[shift .+ (1:T_mpc)], u_ref[shift .+ (1:T_mpc-1)], prob_mpc)

		# Solve
		@time z_mpc, info = solve(prob_mpc, copy(z0_mpc))
		x_mpc, u_mpc = unpack(z_mpc, prob_mpc)

		# u_ctrl = -1.0 *  K[t] * (x_hist[end] - x_ref[t])
		# u_ctrl = u[t]
		u_ctrl = u_mpc[1]

		u_ctrl = min.(max.(_ul, u_ctrl[1]), _uu)

		push!(x_hist,
			fd(model, x_hist[end],
				u_ctrl,
				[0.0; 0.5] .+ 0.001 * randn(model.d), h, nothing))
		push!(u_hist, u_ctrl)
		shift += 1
	end

	return x_hist, u_hist
end

x_hist_mpc, u_hist_mpc = run_mpc()

plot(hcat(x...)',
	label = ["x (ref)" "v (ref)"], color = :black, width = 2.0)
plot!(hcat(x_hist_mpc...)',
	label = ["x (mpc)" "v (mpc)"], color = :magenta, width = 1.0)

plot(hcat(u...)',
	label = "u (ref)", color = :black, width = 2.0, linetype = :steppost)
plot!(hcat(u_hist_mpc...)',
	label = "u (mpc)", color = :magenta, width = 1.0, linetype = :steppost)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Plots
using SparseArrays

# Solver options
N = 51
integration = :rk4
opts = SolverOptions()
opts.verbose = false
opts.square_root = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-2
opts.penalty_max = 1e4

# 3 sphere obstacles + unit quaternion
model, obj_obs = TrajectoryOptimization.Dynamics.quadrotor_3obs

solver = Solver(model,obj_obs,integration=integration,N=N,opts=opts)
U = 0.5*9.81/4.0*ones(solver.model.m, solver.N-1)
results, stats = solve(solver,U)

λ_update_default!(results,solver)
update_constraints!(results,solver)
@assert max_violation(results) < opts.constraint_tolerance

J_prev = cost(solver,results)
c_max_prev = stats["c_max"][end]
J_prev = cost(solver,results)

p,pI,pE = get_num_constraints(solver)
p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

## Newton 2 ###############
results_new = copy(results)

# newton_solve!(results_new,solver)
newton_results = NewtonResults(solver)
newton_active_set!(newton_results,results_new,solver)
update_newton_results!(newton_results,results_new,solver)
newton_step!(results_new,newton_results,solver,1.0,100.0)
max_violation(results_new)

# get batch problem sizes
Nz = newton_results.Nz
Np = newton_results.Np
Nx = newton_results.Nx

# assemble KKT matrix,vector
_idx1 = Array(1:Nz)
_idx2 = Array((1:Np) .+ Nz)[newton_results.active_set]
_idx3 = Array((1:Nx) .+ (Nz + Np))
_idx4 = Array((1:Np) .+ (Nz + Np + Nx))[newton_results.active_set_ineq]

_idx5 = [_idx1;_idx2;_idx3;_idx4]
eigvals(Array(newton_results.A)[_idx5,_idx5]+1000I)

cond(Array(newton_results.A)[_idx5,_idx5])
cond(Array(newton_results.A)[_idx5,_idx5]+1000I)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module PendulumODE2

ODE-model of a mass point attached via a rod to a revolute joint with 2 states.
"""
module PendulumODE2

import ..ModiaMath

mutable struct Model <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # Parameters
    L::Float64
    m::Float64
    d::Float64
    g::Float64

    function Model(;L=1.0, m=1.0, d=0.1, g=9.81, phi0_deg=90)
        @assert(L > 0.0)
        @assert(m > 0.0)
        @assert(d >= 0.0)
        simulationState = ModiaMath.SimulationState("PendulumODE2", getModelResidues!, [phi0_deg * pi / 180; 0.0], getVariableName; nw=2)
        new(simulationState, L, m, d, g)
    end
end

getVariableName(model, vcat, vindex) = ModiaMath.getVariableName(model, vcat, vindex;
                                                         xNames=["phi"    , "w"],
                                                         derxNames=["der_phi", "der_w"],
                                                         wNames=["rx", "ry"])

function getModelResidues!(m::Model, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    phi  = _x[1]
    w    = _x[2]
    derw = _derx[2]
    _r[1] = w    - _derx[1]
    _r[2] = m.m * m.L^2 * derw + (m.m * m.g * m.L * sin(phi) + m.d * w)

    if ModiaMath.isStoreResult(m.simulationState)
        rx =  m.L * sin(phi)
        ry = -m.L * cos(phi)
        _w[1] = rx
        _w[2] = ry
    end
    return nothing
end

end


text:
using DataFrames
n = 10^3                                #  this is slow!!!
d = DataFrame(x = randn(n), y = 1:n)

w = Widget()
setWindowTitle(w, "View a data frame")

lyt = VBoxLayout(w)
setLayout(w, lyt)

## This is really slow for large data sets, but works fine at this scale (2000 items)
## for large data sets the DataFrameModel can be used.

m = StandardItemModel(w)
## currently no set_items method to populate from a data frame, so we
## do it at the Qt level
for row in 1:nrow(d), col in 1:ncol(d)
    Qt[:QStandardItem](string(d[row, col])) | item -> m[:setItem](row-1, col-1, item)
end

view = TableView(w)
setModel(view, m)
addWidget(lyt, view)

## Some adjustments for non-exported functions
view[:setAlternatingRowColors](true)
header = view[:horizontalHeader]()
header[:setStretchLastSection](true)

raise(w)


text:
using FactCheck, Images, Colors, FixedPointNumbers
using Compat

facts("Overlay") do
    gray = linspace(0.0, 1.0, 5)
    context("One") do
        ovr = Images.Overlay((2gray, 2gray), (RGB(1, 0, 0), RGB(0, 0, 1)), (Clamp{Float64}(), Clamp{Float64}()))
        @fact ovr[1] --> RGB(0, 0, 0)
        @fact ovr[2] --> RGB{U8}(0.5, 0, 0.5)
        @fact ovr[3] --> ovr[4]
        @fact ovr[4] --> ovr[5]
        @fact ovr[5] --> exactly(RGB(1, 0, 1))
        @fact eltype(ovr) --> RGB{U8}
        @fact length(ovr) --> 5
        @fact size(ovr) --> (5,)
        @fact size(ovr, 1) --> 5
        @fact size(ovr, 2) --> 1
        @fact nchannels(ovr) --> 2
        @fact raw(ovr) --> [0x00 0x80 0xff 0xff 0xff;
                            0x00 0x00 0x00 0x00 0x00;
                            0x00 0x80 0xff 0xff 0xff]
        @fact separate(ovr) --> UFixed8[0   0   0;
                                        0.5 0   0.5;
                                        1   0   1;
                                        1   0   1;
                                        1   0   1]
        iob = IOBuffer()
        show(iob, ovr)  # exercise only
    end

    context("Two") do
        ovr = Images.Overlay((gray, 0*gray), (RGB{UFixed8}(1, 0, 1), RGB{UFixed8}(0, 1, 0)), ((0, 1), (0, 1)))
        @fact eltype(ovr) --> RGB{UFixed8}
        s = similar(ovr)
        @fact typeof(s) --> Vector{RGB{UFixed8}}
        @fact length(s) --> 5
        s = similar(ovr, RGB{Float32})
        @fact isa(s, Vector{RGB{Float32}}) --> true
        @fact length(s) --> 5
        s = similar(ovr, RGB{Float32}, (3, 2))
        @fact isa(s, Matrix{RGB{Float32}}) --> true
        @fact size(s) --> (3, 2)
        buf = Images.uint32color(ovr)
        gray8 = round(UInt8, 255*gray)
        nogreen = [convert(UInt32, g)<<16 | convert(UInt32, g) for g in gray8]
        @fact buf --> nogreen
    end

    context("Three") do
        ovr = Images.Overlay((gray, 0*gray), [RGB(1, 0, 1), RGB(0, 1, 0)], ([0, 1], [0, 1]))
        @fact_throws ErrorException Images.Overlay((gray, 0*gray), (RGB(1, 0, 1), RGB(0, 1, 0)), ((0,), (0, 1)))
        @fact_throws ErrorException ovr[1] = RGB(0.2, 0.4, 0.6)
    end

    context("Four") do
        img1 = Images.Image(gray, Dict{Compat.ASCIIString, Any}([("colorspace", "Gray"), ("spatialorder",["x"])]))
        ovr = Images.OverlayImage((2gray, img1), (RGB{Float32}(1, 0, 1), RGB{Float32}(0, 1, 0)), ((0, 1),(0, 1)))
        @fact isa(ovr, Images.Image) --> true
        @fact haskey(ovr, "colorspace") --> false
        @fact Images.colorspace(ovr) --> "RGB"
        @fact ovr[2] --> RGB{Float32}(0.5, 0.25, 0.5)
        a = rand(Float32, 3, 2)
        b = rand(Float32, 3, 2)
        ovr = Images.OverlayImage((a, b), (RGB{Float32}(1, 0, 1), RGB{Float32}(0, 1, 0)), ((0, 1), (0, 1)))
        @fact isa(ovr, Images.Image) --> true
        @fact abs(ovr[1, 2] - RGB{Float32}(a[1, 2], b[1, 2], a[1, 2])) --> roughly(0, atol=1e-5)
    end

    context("permutation") do
        L1 = convert(Image{Gray}, rand(Float32, 10,10))
        L2 = convert(Image{Gray}, rand(Float32, 10,10))
        L3 = convert(Image{Gray}, rand(Float32, 10,10))
        overlay = OverlayImage((L1, L2, L3), (colorant"red", colorant"blue", colorant"green"), ((0,1),(0,1),(0,1)))
        permutedims(overlay, [2,1])
        permutedims(data(overlay), [2,1])
    end
end


text:
@testset "Receiver" begin

    receiver = Receiver(
        4e6Hz,
        intermediate_frequency = 0.0Hz,
        gain_phase_mism_crosstalk = SMatrix{4,4,ComplexF64}(I),
        attitude = RotXYZ(0.0, 0.0, 0.0),
        n0 = 1/Hz
    )
    next_receiver = GNSSSimulator.propagate(receiver, 1000, Random.GLOBAL_RNG)
    @test get_gain_phase_mism_crosstalk(next_receiver) == SMatrix{4,4}(I)
    @test get_attitude(next_receiver) == RotXYZ(0.0, 0.0, 0.0)
    @test get_noise_density(next_receiver) == 1/Hz
    @test get_sampling_frequency(next_receiver) == 4e6Hz
    @test get_intermediate_frequency(next_receiver) == 0.0Hz
    @test get_gain_phase_mism_crosstalk(next_receiver) == SMatrix{4,4}(I)
    @test @inferred(get_noise_std(next_receiver)) == sqrt(1/Hz * 4e6Hz)

    receiver = Receiver(
        4e6Hz,
        intermediate_frequency = 0.0Hz,
        gain_phase_mism_crosstalk = SMatrix{4,4,ComplexF64}(I),
        attitude = RotXYZ(0.0, 0.0, 0.0),
        n0 = 1/Hz
    )
    @test get_gain_phase_mism_crosstalk(receiver) == SMatrix{4,4}(I)
    @test get_attitude(receiver) == RotXYZ(0.0, 0.0, 0.0)
    @test get_noise_density(receiver) == 1/Hz
    @test get_sampling_frequency(receiver) == 4e6Hz
    @test get_intermediate_frequency(receiver) == 0.0Hz
    @test get_gain_phase_mism_crosstalk(receiver) == SMatrix{4,4}(I)
end


text:
"""
    System{TF, TV<:AbstractVector{TF}, TM<:AbstractMatrix{TF}, TTF<:AbstractVector{TF}}

Contains the system state, residual vector, and jacobian matrices as well as
pointers to be able to access their contents.  Also contains additional storage
needed for time domain simulations.

# Fields:
 - `static`: Flag indicating whether system matrices are only valid for static analyses
 - `x`: State vector
 - `r`: Residual vector
 - `K`: System jacobian matrix with respect to the state variables
 - `M`: System jacobian matrix with respect to the time derivative of the state variables
 - `irow_pt`: Row index of first equilibrium equation for each point
 - `irow_beam`: Row index of first equation for just this beam element
 - `irow_beam1`: Row index of first equation for the left side of each beam
 - `irow_beam2`: Row index of first equation for the right side of each beam
 - `icol_pt`: Row/Column index of first state variable for each point
 - `icol_beam`: Row/Column index of first state variable for each beam element
 - `current_step`: Current time step
 - `udot_init`: `2/dt*u + udot` for each beam element from the previous time step
 - `θdot_init`: `2/dt*θ + θdot` for each beam element from the previous time step
 - `CtCabPdot_init`: `2/dt*C'*Cab*P + C'*Cab*Pdot` for each beam element from the previous time step
 - `CtCabHdot_init`: `2/dt*C'*Cab*H + C'*Cab*Hdot` for each beam element from the previous time step
"""
struct System{TF, TV<:AbstractVector{TF}, TM<:AbstractMatrix{TF}}
    static::Bool
    x::TV
    r::TV
    K::TM
    M::TM
    irow_pt::Vector{Int}
    irow_beam::Vector{Int}
    irow_beam1::Vector{Int}
    irow_beam2::Vector{Int}
    icol_pt::Vector{Int}
    icol_beam::Vector{Int}
    current_step::Ref{Int}
    udot_init::Vector{SVector{3,TF}}
    θdot_init::Vector{SVector{3,TF}}
    CtCabPdot_init::Vector{SVector{3,TF}}
    CtCabHdot_init::Vector{SVector{3,TF}}
end
Base.eltype(::System{TF, TV, TM}) where {TF, TV, TM} = TF

"""
    System([TF=eltype(assembly),] assembly, points, static)

Initialize an object of type `System` which stores the system state, residual vector,
current time function values,and jacobian matrices as well as pointers to
be able to access their contents.

# Arguments:
 - `TF:` (optional) Used to specify floating point type used by resulting `System` object
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `points`: Point indices which should be preserved in the system of equations.
        All points with prescribed conditions should be included.
 - `static`: Flag indicating whether system matrices will be used for static simulations
"""
System(assembly, points, static) = System(eltype(assembly), assembly, points, static)

function System(TF::Type{<:AbstractFloat}, assembly, points, static)

    # get number of beams
    nbeam = length(assembly.elements)

    # get the number of beams connected to each point.
    n_connections = point_connections(assembly)

    # get the size and pointers into the system matrices
    n, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam =
        system_indices(assembly, points, n_connections, static)

    # initialize system matrices
    x = zeros(TF, n)
    r = zeros(TF, n)
    K = spzeros(TF, n, n)
    M = spzeros(TF, n, n)

    # initialize pointer to value of current time step
    current_step = Ref(0)

    # initialize storage for time domain simulations
    udot_init = [@SVector zeros(TF, 3) for i = 1:nbeam]
    θdot_init = [@SVector zeros(TF, 3) for i = 1:nbeam]
    CtCabPdot_init = [@SVector zeros(TF, 3) for i = 1:nbeam]
    CtCabHdot_init = [@SVector zeros(TF, 3) for i = 1:nbeam]

    return System(static, x, r, K, M, irow_pt, irow_beam, irow_beam1, irow_beam2,
        icol_pt, icol_beam, current_step, udot_init, θdot_init, CtCabPdot_init,
        CtCabHdot_init)
end

"""
   reset_state!(system)

Reset the state variables in `system` (stored in `system.x`) to zero.
"""
function reset_state!(system)
    system.x .= 0
    return system
end

"""
    point_connections(assembly)

Count the number of beams connected to each point
"""
function point_connections(assembly)

    npt = length(assembly.points)

    n_connections = Array{Int, 1}(undef, npt)

    for ipt = 1:npt
        n_connections[ipt] = count(x -> x == ipt, assembly.start) + count(x -> x == ipt, assembly.stop)
    end

    return n_connections
end

"""
    system_indices(assembly, points, n_connections, static)

Solve for the row indices of the first equilibrium or compatability equations for
each point and side of each beam element.  Also solve for the row/column index of
each point and beam state variable.

Note that this function includes the following logic which reduces the size of
the system of equations where possible (without sacrificing any accuracy):

If only two beams meet at a point, the 6 unknowns associated with that point as
well as the 6 compatability equations are eliminated from the system, except if
specified in the array `points`.  Points for which unknowns have been
eliminated are assigned a column index of -1.  Beams for which the compatability
equations have been eliminated are also assigned an index of -1

# Arguments:
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `points`: Point indices which should be preserved in the system matrices
 - `n_connections`: Number of connections to each point
 - `static`: flag indicating whether analysis is static

# Return Arguments:
 - `n`: total number of equations/unknowns in the system
 - `irow_pt`: Row index of first equilibrium equation for each point
 - `irow_beam`: Row index of first equation for just this beam element
 - `irow_beam1`: Row index of first equation for the left side of each beam
 - `irow_beam2`: Row index of first equation for the right side of each beam
 - `icol_pt`: Column index of first state variable for each point
 - `icol_beam`: Column index of first state variable for each beam element
"""
function system_indices(assembly, points, n_connections, static)

    npt = length(assembly.points)
    nbeam = length(assembly.elements)

    assigned = fill(false, npt)

    irow_pt = Array{Int, 1}(undef, npt)
    irow_beam = Array{Int, 1}(undef, nbeam)
    irow_beam1 = Array{Int, 1}(undef, nbeam)
    irow_beam2 = Array{Int, 1}(undef, nbeam)

    icol_pt = Array{Int, 1}(undef, npt)
    icol_beam = Array{Int, 1}(undef, nbeam)

    irow = 1
    icol = 1
    for ibeam = 1:nbeam
        ipt = assembly.start[ibeam]
        if !assigned[ipt]
            assigned[ipt] = true
            # 6 equilibrium equations + 6 compatability equations
            irow_pt[ipt] = irow
            irow_beam1[ibeam] = irow
            irow += 12
            # add unknowns for each point
            if n_connections[ipt] == 2 && !(ipt in points)
                # no additional unknowns when only two beams meet at a point
                icol_pt[ipt] = -1
            else
                # 6 additional unknowns for each point
                icol_pt[ipt] = icol
                icol += 6
            end
        else
            # add compatability equations
            if n_connections[ipt] == 2 && !(ipt in points)
                # no additional equations when only two beams meet at a point
                irow_beam1[ibeam] = -1
            else
                # 6 additional compatibility equations
                irow_beam1[ibeam] = irow
                irow += 6
            end
        end

        # 12 unknowns for each element
        icol_beam[ibeam] = icol
        icol += 12

        if static
            irow_beam[ibeam] = -1
        else
            # 6 linear and angular momentum residual equations for each element
            irow_beam[ibeam] = irow
            irow += 6
            # 6 additional unknowns for each member for unsteady analyses
            icol += 6
        end

        ipt = assembly.stop[ibeam]
        if !assigned[ipt]
            assigned[ipt] = true
            # 6 equilibrium equations + 6 compatability equations
            irow_pt[ipt] = irow
            irow_beam2[ibeam] = irow
            irow += 12
            # add unknowns for each point
            if n_connections[ipt] == 2 && !(ipt in points)
                # no additional unknowns when only two beams meet at a point
                icol_pt[ipt] = -1
            else
                # 6 additional unknowns for each point
                icol_pt[ipt] = icol
                icol += 6
            end
        else
            if n_connections[ipt] == 2 && !(ipt in points)
                # no additional compatability equations when only two beams meet at a point
                irow_beam2[ibeam] = -1
            else
                # 6 additional compatibility equations
                irow_beam2[ibeam] = irow
                irow += 6
            end
        end
    end

    # number of unknowns/equations
    n = irow - 1

    return n, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam
end

"""
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam)
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam, x0, v0, ω0)
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam, x0, v0, ω0,
        u, θ, udot, θdot)
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam, x0, v0, ω0,
        udot, θdot_init, CtCabPdot, CtCabHdot, dt)

Populate the residual vector `resid` with the results of the residual equations
for the system.

There are four implementations corresponding to the following analysis types:
 - Static
 - Dynamic - Steady State
 - Dynamic - Initial Step (for initializing time domain simulations)
 - Dynamic - Time Marching

 See "GEBT: A general-purpose nonlinear analysis tool for composite beams" by
 Wenbin Yu and "Geometrically nonlinear analysis of composite beams using
 Wiener-Milenković parameters" by Qi Wang and Wenbin Yu.

# Arguments
 - `resid`: System residual vector
 - `x`: Current state variables of the system
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `prescribed_conditions`: Dictionary of prescribed conditions at all time steps
 - `distributed_loads`: Dictionary of distributed loads at all time steps
 - `istep`: Current time step
 - `irow_pt`: Row index of first equilibrium equation for each point
 - `irow_beam`: Row index of first equation for just this beam element
 - `irow_beam1`: Row index of first equation for the left side of each beam
 - `irow_beam2`: Row index of first equation for the right side of each beam
 - `icol_pt`: Column index of first state variable for each point
 - `icol_beam`: Column index of first state variable for each beam element

# Additional Arguments for Dynamic Analyses
 - `x0`: Global frame origin (for the current time step)
 - `v0`: Global frame linear velocity (for the current time step)
 - `ω0`: Global frame angular velocity (for the current time step)

# Additional Arguments for Initial Step Analyses
 - `u`: deflection variables for each beam element
 - `θ`: rotation variables for each beam element
 - `udot`: time derivative of u for each beam element
 - `θdot`: time derivative of θ for each beam element

# Additional Arguments for Time Marching Analyses
 - `udot`: `-2/dt*u - udot` for each beam element from the previous time step
 - `θdot_init`: `-2/dt*θ - θdot` for each beam element from the previous time step
 - `CtCabPdot`: `-2/dt*C'*Cab*P - C'*Cab*Pdot` for each beam element from the previous time step
 - `CtCabHdot`: `-2/dt*C'*Cab*H - C'*Cab*Hdot` for each beam element from the previous time step
 - `dt`: time step size
"""
system_residual!

# static
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        # get pointers for element
        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        resid = element_residual!(resid, x, ibeam, assembly.elements[ibeam],
             distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        resid = point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    return resid
end

# dynamic
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
    x0, v0, ω0)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        # get pointers for element
        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        resid = element_residual!(resid, x, ibeam, assembly.elements[ibeam],
             distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2, x0, v0, ω0)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        resid = point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    return resid
end

# initial step
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
    x0, v0, ω0, u, θ, udot, θdot)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        # get pointers for element
        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        resid = element_residual!(resid, x, ibeam, assembly.elements[ibeam],
             distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2, x0, v0, ω0, u, θ, udot, θdot)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        resid = point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    return resid
end

# time marching
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
    x0, v0, ω0, udot, θdot_init, CtCabPdot, CtCabHdot, dt)

    nbeam = length(assembly.elements)
    npoint = length(assembly.points)

    # add contributions to residual equations from the beam elements
    for ibeam = 1:nbeam

        # get pointers for element
        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        resid = element_residual!(resid, x, ibeam, assembly.elements[ibeam],
             distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2, x0, v0, ω0, udot, θdot_init, CtCabPdot,
            CtCabHdot, dt)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        resid = point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    return resid
end

"""
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam)
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
        x0, v0, ω0)
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
        x0, v0, ω0, u, θ, udot, θdot)
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
        x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_init, dt)

Populate the jacobian matrix `jacob` with the jacobian of the residual vector
with respect to the state variables.

There are four implementations corresponding to the following analysis types:
 - Static
 - Dynamic - Steady State
 - Dynamic - Initial Step (for initializing time domain simulations)
 - Dynamic - Time Marching

 See "GEBT: A general-purpose nonlinear analysis tool for composite beams" by
 Wenbin Yu and "Geometrically nonlinear analysis of composite beams using
 Wiener-Milenković parameters" by Qi Wang and Wenbin Yu.

# Arguments
 - `jacob`: Jacobian matrix
 - `x`: Vector containing current state variables of the system
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `prescribed_conditions`: Dictionary of prescribed conditions at all time steps
 - `distributed_loads`: Dictionary of distributed loads at all time steps
 - `istep`: Current time step
 - `irow_pt`: Row index of first equilibrium equation for each point
 - `irow_beam`: Row index of first equation for just this beam element
 - `irow_beam1`: Row index of first equation for the left side of each beam
 - `irow_beam2`: Row index of first equation for the right side of each beam
 - `icol_pt`: Column index of first state variable for each point
 - `icol_beam`: Column index of first state variable for each beam element

# Additional Arguments for Dynamic Analyses
 - `x0`: Global frame origin (for the current time step)
 - `v0`: Global frame linear velocity (for the current time step)
 - `ω0`: Global frame angular velocity (for the current time step)

# Additional Arguments for Initial Step Analyses
 - `u`: deflection variables for each beam element
 - `θ`: rotation variables for each beam element
 - `udot`: time derivative of u for each beam element
 - `θdot`: time derivative of θ for each beam element

# Additional Arguments for Time Marching Analyses
 - `udot`: `-2/dt*u - udot` for each beam element from the previous time step
 - `θdot_init`: `-2/dt*θ - θdot` for each beam element from the previous time step
 - `CtCabPdot`: `-2/dt*C'*Cab*P - C'*Cab*Pdot` for each beam element from the previous time step
 - `CtCabHdot`: `-2/dt*C'*Cab*H - C'*Cab*Hdot` for each beam element from the previous time step
 - `dt`: time step size
"""
system_jacobian!

# static
function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        jacob = element_jacobian!(jacob, x, ibeam, assembly.elements[ibeam],
            distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        jacob = point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)

    end

    return jacob
end

# dynamic
function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
    x0, v0, ω0)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        jacob = element_jacobian!(jacob, x, ibeam, assembly.elements[ibeam],
            distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2, x0, v0, ω0)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        jacob = point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    return jacob
end

# initial step
function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
    x0, v0, ω0, u, θ, udot, θdot)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        jacob = element_jacobian!(jacob, x, ibeam, assembly.elements[ibeam],
            distributed_loads, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2, x0, v0, ω0, u, θ, udot, θdot)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        jacob = point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    return jacob
end

# time marching
function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    istep, irow_pt, irow_beam, irow_beam1, irow_beam2, icol_pt, icol_beam,
    x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_init, dt)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ibeam = 1:nbeam

        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        jacob = element_jacobian!(jacob, x, ibeam, assembly.elements[ibeam],
            distributed_loads, istep, icol, irow_b, irow_b1,
            irow_p1, irow_b2, irow_p2, x0, v0, ω0, udot_init, θdot_init,
            CtCabPdot_init, CtCabHdot_init, dt)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_pt[ipoint] <= 0
            continue
        end

        icol = icol_pt[ipoint]
        irow_p = irow_pt[ipoint]

        jacob = point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            istep, icol, irow_p, irow_beam1, irow_beam2)
    end

    # zero out near-zero values ( < eps() )
    jacob = droptol!(jacob, eps(eltype(jacob)))

    return jacob
end

"""
    system_mass_matrix!(jacob, x, assembly, irow_pt, irow_beam, irow_beam1,
        irow_beam2, icol_pt, icol_beam)

Populate the system "mass matrix", the jacobian of the residual vector
with respect to the time derivatives of the state variables.

See "GEBT: A general-purpose nonlinear analysis tool for composite beams" by
Wenbin Yu and "Geometrically nonlinear analysis of composite beams using
Wiener-Milenković parameters" by Qi Wang and Wenbin Yu.

# Arguments
 - `jacob`: Jacobian matrix
 - `x`: Vector containing current state variables of the system
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `irow_pt`: Row index of first equilibrium equation for each point
 - `irow_beam`: Row index of first equation for just this beam element
 - `irow_beam1`: Row index of first equation for the left side of each beam
 - `irow_beam2`: Row index of first equation for the right side of each beam
 - `icol_pt`: Column index of first state variable for each point
 - `icol_beam`: Column index of first state variable for each beam element
"""
function system_mass_matrix!(jacob, x, assembly, irow_pt, irow_beam, irow_beam1,
    irow_beam2, icol_pt, icol_beam)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # add contributions to "mass matrix" from the beam elements
    for ibeam = 1:nbeam

        icol = icol_beam[ibeam]
        irow_b = irow_beam[ibeam]
        irow_b1 = irow_beam1[ibeam]
        irow_p1 = irow_pt[assembly.start[ibeam]]
        irow_b2 = irow_beam2[ibeam]
        irow_p2 = irow_pt[assembly.stop[ibeam]]

        element_mass_matrix!(jacob, x, assembly.elements[ibeam], icol, irow_b,
            irow_b1, irow_p1, irow_b2, irow_p2)
    end

    # no contributions to "mass matrix" from point state variables

    # zero out near-zero values ( < eps() )
    jacob = droptol!(jacob, eps(eltype(jacob)))

    return jacob
end


text:
using SimpleDiffEq, Test, StaticArrays

# dX_t = 2u dt + dW_t
f(u,p,t) = 2u
g(u,p,t) = 1
u0 = 0.5
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u0,tspan)

sol = solve(prob,SimpleEM(),dt=0.25)

@test sol.t == collect(0:0.25:1.0)
@test length(sol.u) == 5
@test typeof(sol) <: DESolution

u0 = @SVector [0.1,0.2]
prob = SDEProblem(f,g,u0,tspan)
sol = solve(prob,SimpleEM(),dt=0.25)

@test typeof(sol.u) <: Vector{SVector{2,Float64}}

f(du,u,p,t) = du .= 2.0 * u
g(du,u,p,t) = du .= 1
u0 = 0.5ones(4)
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u0,tspan)
sol = solve(prob,SimpleEM(),dt=0.25)


text:
function f(x)
    s = Threads.Atomic{Float64}(0.0)
    n = length(x)
    @Threads.threads for k in 1:Threads.nthreads()
        y = 0.0
        @inbounds @simd for i in getrange(n)
            y += x[i]
        end
        Threads.atomic_add!(s, y)
    end
    s[]
end

function g(x)
    n = length(x)
    s = 0.0
    @inbounds @simd for i in 1:n
        s += x[i]
    end
    s
end

z = rand(10^8)
@testset "summation: getrange" begin
    println("--------")
    println("threaded")
    f(z)
    @time _f = f(z)
    println("unthreaded")
    g(z)
    @time _g = g(z)
    @test _f ≈ _g
end


text:
# tracks.jl
## https://developer.spotify.com/documentation/web-api/reference/tracks/

## https://developer.spotify.com/documentation/web-api/reference/tracks/get-audio-analysis/
@doc """
# Get Audio Analysis for a Track
**Summary**: Get a detailed audio analysis for a single track identified by it's unique Spotify ID.

`track_id` _Required_: The Spotify ID for the track.

[Reference](https://developer.spotify.com/documentation/web-api/reference/tracks/get-audio-analysis/)
""" ->
function tracks_get_audio_analysis(track_id)
    return spotify_request("audio-analysis/$track_id")
end

## https://developer.spotify.com/documentation/web-api/reference/tracks/get-audio-features/
@doc """
# Get Audio Features for a Track
**Summary**: Get audio feature information for a single track identified by it's unique Spotify ID

`track_id` _Required_: The Spotify ID for the track. Up to 50 tracks can be passed seperated with a comma. No whitespace.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/tracks/get-audio-features/)
""" ->
function tracks_get_audio_features(track_id)
    return spotify_request("audio-features/$track_id")
end

## https://developer.spotify.com/documentation/web-api/reference/tracks/get-track/
@doc """
# Get a Track
**Summary**: Get a spotify catalog information for a single track identified by it's unique Spotify ID

`track_id` _Required_: The Spotify ID for the track. Up to 50 tracks can be passed seperated with a comma. No whitespace.\n
`market` _Optional_: An ISO 3166=1 alpha-2 country code

[Reference](https://developer.spotify.com/documentation/web-api/reference/tracks/get-track/)
"""
function tracks_get(track_id, market="US")
    return spotify_request("tracks/$track_id?market=$market")
end


text:
function bar(arr)
    if "aaa" in arr && "zzz" in arr
        return "hello"
    else
        return "goodbye"
    end
end


text:
using BenchmarkTools
using CSV
using DataFrames
using HCIDatasets
using PSFModels
using PythonCall

@info "Starting PSF evaluation benchmark"

ap = pyimport("astropy")
apm = pyimport("astropy.modeling.models")
apf = pyimport("astropy.modeling.fitting")

@info "Using astropy version $(ap.__version__)"

julia_models = [
    gaussian(x=0, y=0, fwhm=(5.3, 4.7)),
    airydisk(x=0, y=0, fwhm=12.7),
    moffat(x=0, y=0, fwhm=8.1)
]
astropy_models = [
    apm.Gaussian2D(x_stddev=5.3, y_stddev=4.7),
    apm.AiryDisk2D(radius=12.7),
    apm.Moffat2D(gamma=8.1, alpha=-1)
]

names = ["Gaussian", "AiryDisk", "Moffat"]
jl_times = []
py_times = []

for (jl_mod, py_mod, name) in zip(julia_models, astropy_models, names)
    jl_time = @belapsed $jl_mod(0.5, 0.7)
    push!(jl_times, jl_time)

    py_time = @belapsed $py_mod(0.5, 0.7)
    push!(py_times, py_time)

    @info "$name" PSFModels=jl_time astropy=py_time
end

filename = joinpath(@__DIR__, "evaluation_results.csv")
DataFrame(name=names, psfmodels=jl_times, astropy=py_times) |> CSV.write(filename)

@info "Results saved to $filename"

@info "Finished PSF evaluation benchmark"

################################################################################
################################################################################
################################################################################

@info "Starting PSF fitting benchmark"

psf = HCIDatasets.BetaPictoris[:psf]
jl_times = []
py_times = []

P0 = (x=20, y=20, fwhm=(5, 5), amp=0.1)
t_jl = @belapsed PSFModels.fit(gaussian, $P0, $psf)
cartinds = CartesianIndices(psf)
axx = map(idx -> idx.I[1], cartinds)
axy = map(idx -> idx.I[2], cartinds)
t_py = @belapsed begin
    g0 = $apm.Gaussian2D(x_stddev=5/2.355, y_stddev=5/2.355, x_mean=20, y_mean=20, amplitude=0.1)
    fitter = $apf.LevMarLSQFitter()
    fitter(g0, $axx, $axy, $psf)
end
push!(jl_times, t_jl)
push!(py_times, t_py)
@info "Gaussian" PSFModels=t_jl astropy=t_py

P0 = (x=20, y=20, fwhm=5, amp=0.1)
t_jl = @belapsed PSFModels.fit(airydisk, $P0, $psf)
t_py = @belapsed begin
    a0 = $apm.AiryDisk2D(radius=5/0.973, x_0=20, y_0=20, amplitude=0.1)
    fitter = $apf.LevMarLSQFitter()
    fitter(a0, $axx, $axy, $psf)
end
push!(jl_times, t_jl)
push!(py_times, t_py)
@info "AiryDisk" PSFModels=t_jl astropy=t_py

P0 = (x=20, y=20, fwhm=5, amp=0.1)
t_jl = @belapsed PSFModels.fit(moffat, $P0, $psf)
t_py = @belapsed begin
    m0 = $apm.Moffat2D(gamma=5, x_0=20, y_0=20, amplitude=0.1)
    fitter = $apf.LevMarLSQFitter()
    fitter(m0, $axx, $axy, $psf)
end
push!(jl_times, t_jl)
push!(py_times, t_py)
@info "Moffat" PSFModels=t_jl astropy=t_py

filename = joinpath(@__DIR__, "fitting_results.csv")
DataFrame(name=names, psfmodels=jl_times, astropy=py_times) |> CSV.write(filename)

@info "Results saved to $filename"

@info "Finished PSF fitting benchmark"

################################################################################
################################################################################
################################################################################

nothing