text:
using LoopVectorization

function RMP3{T}(refwfn::Fermi.HartreeFock.RHF,alg::Fermi.MollerPlesset.DF) where T <: AbstractFloat
    Fermi.MollerPlesset.print_header()
    ints = refwfn.ints
    nocc    = refwfn.ndocc
    nvir    = refwfn.nvir
    eps = refwfn.eps
    δ2 = 0.0
    @output "\tComputing MP3 with DF algorithm\n\n"
    Fermi.Integrals.aux_ri!(ints)
    @output "\tBasis: {}\n" ints.bname["primary"]
    @output "\tDF basis: {}\n\n" ints.bname["aux"]
    ttotal = @elapsed begin
    @output "\tComputing and transforming integrals ... "
    t = @elapsed begin
        Bov = ints["BOV"]
        Boo = ints["BOO"]
        Bvv = ints["BVV"]
        Bvo = ints["BVO"]
        #Δ   = [(1/(eps[i] + eps[j] - eps[a] - eps[b])) for i=1:nocc,j=1:nocc,a=nocc+1:nocc+nvir,b=nocc+1:nocc+nvir]
    end
    naux    = size(ints["BOV"],1)
    @assert Bov ≈ permutedims(Bvo,(1,3,2))
    @output " done in {:>5.2f} s\n" t
    Bi = zeros(T,size(Bov[:,1,:]))
    Bj = zeros(T,size(Bi))
    BAB = zeros(T,nvir,nvir)
    for i in 1:nocc
        for j in 1:nocc
            Bi[:,:] = Bov[:,i,:]
            Bj[:,:] = Bov[:,j,:]
            @tensoropt begin
                BAB[a,b] = Bi[Q,a]*Bj[Q,b]
            end
            BBA = transpose(BAB)
            for b in 1:nvir
                for a in 1:nvir
                    iajb = BAB[a,b]
                    ibja = BBA[a,b]
                    δ2 += iajb*(2*iajb - ibja)/(eps[i] + eps[j] - eps[a+nocc] - eps[b+nocc])
                end
            end
        end
    end
    @output "\n\tDF-MP2 energy is {:>14.10f}\n" δ2

    δ3 = 0.0
    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]
    v = 1+drop_occ:nvir
    o = 1:nocc-drop_vir

    @output "\tEvaluating 12 diagrams ..."
    #there are 12 of these bad bois
    δ3 = zero(T)
    # Diagram 1
    δ3_1 = zeros(T,Threads.nthreads())
    @output " 1 "
    ruvvs = [zeros(T,nvir,nvir) for i=1:Threads.nthreads()]
    oorus = [zeros(T,nocc,nocc) for i=1:Threads.nthreads()]
    tsoos = [zeros(T,nocc,nocc) for i=1:Threads.nthreads()]

    Threads.@threads for r in v
        for u in v
            @views B1 = Bvv[:,:,r]
            @views B2 = Bvv[:,u,:]
            ruvv = ruvvs[Threads.threadid()]
            @tensor ruvv[t,s] = B1[Q,t]*B2[Q,s]
            @views B1  = Bov[:,:,r]
            @views B2  = Bov[:,:,u]
            ooru = oorus[Threads.threadid()]
            @tensor ooru[a,b] = B1[Q,a]*B2[Q,b]
            Δ1 = [(1/(eps[i] + eps[j] - eps[r+nocc] - eps[u+nocc])) for i in o,j in o]
            for t in v,s in v
                @views B1 = Bov[:,:,t]
                @views B2 = Bov[:,:,s]
                tsoo = tsoos[Threads.threadid()]
                @tensor tsoo[a,b] = B1[Q,a]*B2[Q,b]
                ruts = ruvv[t,s]
                Δ2 = [(1/(eps[i] + eps[j] - eps[t+nocc] - eps[s+nocc])) for i in o,j in o]
                #@tensor δ3 += 2*ooru[a,b]*ruts*tsoo[a,b]*Δ[a,b]
                _δ3 = 0.0
                @avx for b in o, a in o
                    _δ3 += ooru[a,b]*tsoo[a,b]*Δ1[a,b]*Δ2[a,b]
                end
                δ3_1[Threads.threadid()] += ruts*_δ3
            end
        end
    end
    δ3_1 = 2.0*sum(δ3_1)

    # Diagram 2
    δ3_2 = zeros(T,Threads.nthreads())
    @output "2 "
    Threads.@threads for a in o
        for d in o
            @views B1 = Bov[:,a,:]
            @views B2 = Bov[:,d,:]
            @tensor advv[r,s] := B1[Q,r]*B2[Q,s]
            @views B1 = Boo[:,a,:]
            @views B2 = Boo[:,d,:]
            @tensor ooad[c,b] := B1[Q,c]*B2[Q,b]
            Δ1 = [(1/(eps[a]+eps[d]-eps[r+nocc]-eps[s+nocc])) for r in v, s in v]
            for c in o, b in o
                cbad = ooad[c,b]
                @views B1 = Bov[:,c,:]
                @views B2 = Bov[:,b,:]
                @tensor vvcb[r,s] := B1[Q,r]*B2[Q,s]
                Δ2 = [(1/(eps[c]+eps[b]-eps[r+nocc]-eps[s+nocc])) for r in v, s in v]
                _δ3 = 0.0
                @avx for s in v, r in v
                    _δ3 += advv[r,s]*vvcb[r,s]*Δ1[r,s]*Δ2[r,s]
                end
                δ3_2[Threads.threadid()] += cbad*_δ3
            end
        end
    end
    δ3_2 = 2.0*sum(δ3_2)

    # Diagram 3
    @output "3 "
    δ3_3 = zeros(T,Threads.nthreads())
    Threads.@threads for r in v
        for c in o
            @views B1 = Bov[:,:,r]
            @views B2 = Bov[:,c,:]
            @tensor ocrv[a,t] := B1[Q,a]*B2[Q,t]
            @views B1 = Bvv[:,r,:]
            @views B2 = Boo[:,:,c]
            @tensor rovc[t,a] := B1[Q,a]*B2[Q,t]
            Δ1 = [(1/(eps[a]+eps[c]-eps[r+nocc]-eps[t+nocc])) for a in o, t in v]
            for b in o, s in v
                rbsc = rovc[b,s]
                @views B1 = Bov[:,:,s]
                @views B2 = Bov[:,b,:]
                @tensor svob[t,a] := B1[Q,a]*B2[Q,t]
                Δ2 = [(1/(eps[a]+eps[b]-eps[s+nocc]-eps[t+nocc])) for a in o, t in v]
                _δ3 = 0.0
                @avx for t in v, a in o
                    _δ3 += ocrv[a,t]*svob[t,a]*Δ1[a,t]*Δ2[a,t]
                end
                δ3_3[Threads.threadid()] += rbsc*_δ3
            end
        end
    end
    δ3_3 = -4.0*sum(δ3_3)

    # Diagram 4
    # r b
    #   a s
    #      c t
    #           +
    @output "4 "
    δ3_4 = zeros(T,Threads.nthreads())
    Threads.@threads for r in v
        for b in o
            @views B1 = Bov[:,b,r]
            @views B2 = Bov[:,:,:]
            @tensor borv[c,t] := B1[Q]*B2[Q,c,t]
            @views B1 = Bvv[:,r,:]
            @views B2 = Boo[:,:,b]
            @tensor rovb[a,s] := B1[Q,s]*B2[Q,a]
            Δ1 = [(1/(eps[b]+eps[c]-eps[r+nocc]-eps[t+nocc])) for c in o, t in v]
            for a in o, s in v
                rasb = rovb[a,s]
                @views B1 = Bov[:,a,s]
                @views B2 = Bov[:,:,:]
                @tensor svao[t,c] := B1[Q]*B2[Q,c,t]
                Δ2 = [(1/(eps[a]+eps[c]-eps[s+nocc]-eps[t+nocc])) for c in o, t in v]
                _δ3 = 0.0
                @avx for t in v, c in o
                    _δ3 += borv[c,t]*svao[t,c]*Δ1[c,t]*Δ2[c,t]
                end
                δ3_4[Threads.threadid()] += rasb*_δ3
            end
        end
    end
    δ3_4 = -4.0*sum(δ3_4)

    # Diagram 5
    # c t ->
    #   b s ->
    #     a r ->
    #           +
    @output "5 "
    δ3_5 = zeros(T,Threads.nthreads())
    Threads.@threads for c in o
        for t in v
            @views B1 = Bov[:,:,:]
            @views B2 = Bov[:,c,t]
            @tensor ocvt[a,r] := B1[Q,a,r]*B2[Q]
            @views B1 = Bov[:,:,:]
            @views B2 = Bov[:,c,t]
            @tensor otvc[b,s] := B1[Q,b,s]*B2[Q]
            Δ1 = [(1/(eps[a]+eps[c]-eps[r+nocc]-eps[t+nocc])) for a in o, r in v]
            for b in o, s in v
                btsc = otvc[b,s]
                @views B1 = Bov[:,b,s]
                @views B2 = Bov[:,:,:]
                @tensor vsob[a,r] := B1[Q]*B2[Q,a,r]
                Δ2 = [(1/(eps[a]+eps[b]-eps[r+nocc]-eps[s+nocc])) for a in o, r in v]
                _δ3 = 0.0
                @avx for r in v,a in o
                    _δ3 += ocvt[a,r]*vsob[a,r]*Δ1[a,r]*Δ2[a,r]
                end
                δ3_5[Threads.threadid()] += btsc*_δ3
            end
        end
    end
    δ3_5 = 8.0*sum(δ3_5)

    # Diagram 6
    # c t -> o v
    #   a s -> o v
    #     b r -> o v
    @output "6 "
    δ3_6 = zeros(T,Threads.nthreads())
    Threads.@threads for c in o
        for t in v
            @views B1 = Bov[:,c,:]
            @views B2 = Bov[:,:,t]
            @tensor covt[b,r] := B1[Q,r]*B2[Q,b]
            @views B1 = Bov[:,:,:]
            @views B2 = Bov[:,c,t]
            @tensor otvc[a,s] := B1[Q,a,s]*B2[Q]
            Δ1 = [(1/(eps[c]+eps[b]-eps[r+nocc]-eps[t+nocc])) for b in o, r in v]
            for a in o, s in v
                atsc = otvc[a,s]
                @views B1 = Bov[:,a,:]
                @views B2 = Bov[:,:,s]
                @tensor vsao[r,b] := B1[Q,r]*B2[Q,b]
                Δ2 = [(1/(eps[a]+eps[b]-eps[r+nocc]-eps[s+nocc])) for b in o, r in v]
                _δ3 = 0.0
                @avx for r in v, b in o
                    _δ3 += covt[b,r]*vsao[r,b]*Δ1[b,r]*Δ2[b,r]
                end
                δ3_6[Threads.threadid()] += atsc*_δ3
            end
        end
    end
    δ3_6 = 2.0*sum(δ3_6)

    # Diagram 7
    # a c -> o o
    #   d b -> o o
    #     r s -> v v
    @output "7 "
    δ3_7 = zeros(T,Threads.nthreads())
    Threads.@threads for a in o
        for c in o
            @views B1 = Bov[:,a,:]
            @views B2 = Bov[:,c,:]
            @tensor acrs[r,s] := B1[Q,r]*B2[Q,s]
            @views B1 = Boo[:,:,a]
            @views B2 = Boo[:,:,c]
            @tensor dbac[d,b] := B1[Q,d]*B2[Q,b]
            Δ1 = [(1/(eps[a]+eps[c]-eps[r+nocc]-eps[s+nocc])) for r in v, s in v]
            for d in o, b in o
                _dbac = dbac[d,b]
                @views B1 = Bov[:,d,:]
                @views B2 = Bov[:,b,:]
                @tensor srdb[s,r] := B1[Q,s]*B2[Q,r]
                Δ2 = [(1/(eps[d]+eps[b]-eps[r+nocc]-eps[s+nocc])) for r in v, s in v]
                _δ3 = 0.0
                @avx for s in v, r in v
                    _δ3 += acrs[r,s]*srdb[s,r]*Δ1[r,s]*Δ2[r,s]
                end
                δ3_7[Threads.threadid()] += _dbac*_δ3
            end
        end
    end
    δ3_7 = -1.0*sum(δ3_7)

    # Diagram 8
    # t r -> v v
    #   u s -> v v
    #     a b -> o o
    @output "8 "
    δ3_8 = zeros(T,Threads.nthreads())
    abrts = [zeros(T,nocc,nocc) for i=1:Threads.nthreads()]
    truss = [zeros(T,nvir,nvir) for i=1:Threads.nthreads()]
    usabs = [zeros(T,nocc,nocc) for i=1:Threads.nthreads()]
    Threads.@threads for t in v
        for r in v
            @views B1 = Bov[:,:,r]
            @views B2 = Bov[:,:,t]
            abrt = abrts[Threads.threadid()]
            @tensor abrt[a,b] = B1[Q,a]*B2[Q,b]
            @views B1 = Bvv[:,t,:]
            @views B2 = Bvv[:,r,:]
            trus = truss[Threads.threadid()]
            @tensor trus[u,s] = B1[Q,u]*B2[Q,s]
            Δ1 = [(1/(eps[a]+eps[b]-eps[t+nocc]-eps[r+nocc])) for a in o, b in o]
            for u in v, s in v
                _trus = trus[u,s]
                @views B1 = Bov[:,:,u]
                @views B2 = Bov[:,:,s]
                usab = usabs[Threads.threadid()]
                @tensor usab[a,b] = B1[Q,a]*B2[Q,b]
                Δ2 = [(1/(eps[a]+eps[b]-eps[u+nocc]-eps[s+nocc])) for a in o, b in o]
                _δ3 = 0.0
                @avx for b in o, a in o
                    _δ3 += abrt[a,b]*usab[a,b]*Δ1[a,b]*Δ2[a,b]
                end
                δ3_8[Threads.threadid()] += _trus*_δ3
            end
        end
    end
    δ3_8 = -1.0*sum(δ3_8)

    # Diagram 9
    # b r -> o v
    #   s a -> v o
    #     c t -> o v
    @output "9 "
    δ3_9 = zeros(T,Threads.nthreads())
    Threads.@threads for b in o
        for r in v
            @views B1 = Bov[:,b,r]
            @views B2 = Bov[:,:,:]
            @tensor bcrt[c,t] := B1[Q]*B2[Q,c,t]
            @views B1 = Boo[:,:,b]
            @views B2 = Bvv[:,r,:]
            @tensor arbs[a,s] := B1[Q,a]*B2[Q,s]
            Δ1 = [(1/(eps[c]+eps[b]-eps[r+nocc]-eps[t+nocc])) for c in o, t in v]
            for s in v, a in o
                _arbs = arbs[a,s]
                @views B1 = Bov[:,a,:]
                @views B2 = Bov[:,:,s]
                @tensor tsac[c,t] := B1[Q,t]*B2[Q,c]
                Δ2 = [(1/(eps[a]+eps[c]-eps[s+nocc]-eps[t+nocc])) for c in o, t in v]
                _δ3 = 0.0
                @avx for t in v, c in o
                    _δ3 += bcrt[c,t]*tsac[c,t]*Δ1[c,t]*Δ2[c,t]
                end
                δ3_9[Threads.threadid()] += _arbs*_δ3
            end
        end
    end
    δ3_9 = 2.0*sum(δ3_9)

    # Diagram 10
    # b r -> o v
    #   a s -> o v
    #     c t -> o v
    @output "10 "
    δ3_10 = zeros(T,Threads.nthreads())
    Threads.@threads for b in o
        for r in v
            @views B1 = Bov[:,:,r]
            @views B2 = Bov[:,b,:]
            @tensor cbrt[c,t] := B1[Q,c]*B2[Q,t]
            @views B1 = Bvv[:,r,:]
            @views B2 = Boo[:,:,b]
            @tensor rasb[a,s] := B1[Q,s]*B2[Q,a]
            Δ1 = [(1/(eps[c]+eps[b]-eps[r+nocc]-eps[t+nocc])) for c in o, t in v]
            for a in o, s in v
                _rasb = rasb[a,s]
                @views B1 = Bov[:,a,s]
                @views B2 = Bov[:,:,:]
                @tensor stac[c,t] := B1[Q]*B2[Q,c,t]
                Δ2 = [(1/(eps[a]+eps[c]-eps[s+nocc]-eps[t+nocc])) for c in o, t in v]
                _δ3 = 0.0
                @avx for t in v, c in o
                    _δ3 += cbrt[c,t]*stac[c,t]*Δ1[c,t]*Δ2[c,t]
                end
                δ3_10[Threads.threadid()] += _rasb*_δ3
            end
        end
    end
    δ3_10 = 2.0*sum(δ3_10)

    # Diagram 11
    # a s -> o v
    #   c t -> o v
    #     r b -> v o
    @output "11 "
    δ3_11 = zeros(T,Threads.nthreads())
    Threads.@threads for a in o
        for s in v
            @views B1 = Bov[:,a,:]
            @views B2 = Bov[:,:,s]
            @tensor abrs[b,r] := B1[Q,r]*B2[Q,b]
            @views B1 = Bov[:,a,s]
            @views B2 = Bov[:,:,:]
            @tensor scat[c,t] := B1[Q]*B2[Q,c,t]
            Δ1 = [(1/(eps[a]+eps[b]-eps[r+nocc]-eps[s+nocc])) for  b in o, r in v]
            for c in o, t in v
                _scat = scat[c,t]
                @views B1 = Bov[:,:,:]
                @views B2 = Bov[:,c,t]
                @tensor rtbc[b,r] := B1[Q,b,r]*B2[Q]
                Δ2 = [(1/(eps[c]+eps[b]-eps[r+nocc]-eps[t+nocc])) for b in o, r in v]
                _δ3 = 0.0
                @avx for r in v, b in o
                    _δ3 += abrs[b,r]*rtbc[b,r]*Δ1[b,r]*Δ2[b,r]
                end
                δ3_11[Threads.threadid()] += _scat*_δ3
            end
        end
    end
    δ3_11 = -4.0*sum(δ3_11)

    # Diagram 12
    # c t -> o v
    #   a s -> o v
    #     r b -> v o
    @output "12 "
    δ3_12 = zeros(T,Threads.nthreads())
    Threads.@threads for c in o
        for t in v
            @views B1 = Bov[:,:,:]
            @views B2 = Bov[:,c,t]
            @tensor bcrt[b,r] := B1[Q,b,r]*B2[Q]
            @views B1 = Bov[:,:,:]
            @views B2 = Bov[:,c,t]
            @tensor atsc[a,s] := B1[Q,a,s]*B2[Q]
            Δ1 = [(1/(eps[b]+eps[c]-eps[t+nocc]-eps[r+nocc])) for b in o, r in v]
            for a in o, s in v
                _atsc = atsc[a,s]
                @views B1 = Bov[:,a,:]
                @views B2 = Bov[:,:,s]
                @tensor rsab[r,b] := B1[Q,r]*B2[Q,b]
                Δ2 = [(1/(eps[a]+eps[b]-eps[r+nocc]-eps[s+nocc])) for b in o, r in v]
                _δ3 = 0.0
                @avx for r in v, b in o
                    _δ3 += bcrt[b,r]*rsab[r,b]*Δ1[b,r]*Δ2[b,r]
                end
                δ3_12[Threads.threadid()] += _atsc*_δ3
            end
        end
    end
    δ3_12 = -4.0*sum(δ3_12)
    end

    δ3 += δ3_1 + δ3_2 + δ3_3 + δ3_4 + δ3_5 + δ3_6 + δ3_7 + δ3_8 + δ3_9 + δ3_10 + δ3_11 + δ3_12
    @output "\n\tDF-MP3 energy is {:>14.10f}\n" δ3

    @output "\tDF-MP3 done in {:>5.2f} s\n" ttotal
    #@output repeat("-",80)*"\n"
    RMP3{T}(δ2 + δ3,Fermi.MemTensor{T}(zeros(T,0,0,0,0)))
end


text:
# tests for the State and UAVState types and related functions.

# using Base.Test

# push!(LOAD_PATH, ".")
# include("includes.jl")

# State tests
function test_state_constructor()
    uav_state = UAVState(collect(1.:8)...)
    s = State([uav_state])
    @test s.num_uavs == 1
end

function test_to_polar()
    uav_state_own = UAVState(ones(Float64, 8)...)
    uav_state_int = UAVState(ones(Float64, 8)...)
    s = State([uav_state_own; uav_state_int])
    s.responding_states = [1,1]
    cvec = to_polar(s)
    @test cvec == [0,0,0,sqrt(2),sqrt(2),1,1]

    uav_state_1 = UAVState(0., 10., 0., 10., 0., 0., 0., 0.)
    uav_state_2 = UAVState(100., 10., 100., 10., .2, 0., 0., 0.)
    state = State([uav_state_1, uav_state_2])
    state.responding_states = [1,1]
    cvec = to_polar(state)
    r = sqrt(100^2 + 100^2)
    theta = atan2(100, 100)
    heading = norm_angle(.2 - 2 * pi)
    vown = sqrt(200)
    vint = sqrt(200)
    resown = 1
    resint = 1
    polar_state = [r, theta, heading, vown, vint, resown, resint]
    @test cvec == polar_state

end

function test_state_equality()
    uav_state_1 = UAVState(1.,2.,3.,4.,5.,6.,7.,8.)
    uav_state_2 = UAVState(1.,2.,3.,4.,5.,6.,7.,8.)
    uav_state_3 = UAVState(0.,0.,0.,0.,0.,0.,0.,0.)
    state_1 = State([uav_state_1, uav_state_2])
    state_1.responding_states = [1,1]
    state_2 = State([uav_state_1, uav_state_2])
    state_2.responding_states = [1,1]
    @test state_1 == state_2
    state_3 = State([uav_state_1, uav_state_3])
    state_3.responding_states = [1,1]
    @test state_1 != state_3
end

# UAVState tests
function test_convert_to_vector()
    uav_state = UAVState(1.,2.,3.,4.,5.,6.,7.,8.)
    v = convert(Vector{Float64}, uav_state)
    @test v == collect(1.:8)
end

function test_convert_from_vector()
    uav_state = UAVState(collect(1.:8)...)
    @test uav_state.x == 1.
    @test uav_state.xdot == 2.
    @test uav_state.y == 3.
    @test uav_state.ydot == 4.
    @test uav_state.psi == 5.
    @test uav_state.psidot == 6.
    @test uav_state.phi == 7.
    @test uav_state.phidot == 8.
end

function test_get_speed()
    uav_state_vector = collect(1.:8)
    @test get_speed(uav_state_vector) == sqrt(2^2 + 4^2)
end

function test_get_horizontal_distance()
    uav_state_own = UAVState(ones(Float64, 8)...)
    uav_state_int = UAVState(zeros(Float64, 8)...)
    s = State([uav_state_own; uav_state_int])
    d = get_horizontal_distance(s)
    @test d == sqrt(2)
end

function test_uav_state_equality()
    uav_state_1 = UAVState(1.,2.,3.,4.,5.,6.,7.,8.)
    uav_state_2 = UAVState(1.,2.,3.,4.,5.,6.,7.,8.)
    @test uav_state_1 == uav_state_2
    uav_state_3 = UAVState(0.,0.,0.,0.,0.,0.,0.,0.)
    @test uav_state_1 != uav_state_3
end
tic()
# State tests
test_state_constructor()
# commented tests out of order until decide on polar formulation
test_to_polar()
test_state_equality()

# UAVState tests
test_convert_to_vector()
test_convert_from_vector()
test_get_speed()
test_get_horizontal_distance()
test_uav_state_equality()
toc()

ext: jl


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Base.Test
u0=rand(300,20).*ones(300,20)/2
prob = prob_ode_2Dlinear_notinplace
prob2 = prob_ode_2Dlinear

sol =solve(prob,Euler(),dt=1//2^(4),save_timeseries=false)
sol =solve(prob2,Euler(),dt=1//2^(4),save_timeseries=false)

@time sol =solve(prob,Euler(),dt=1//2^(6),save_timeseries=false)
@time sol2 =solve(prob2,Euler(),dt=1//2^(6),save_timeseries=false)

alloc1 = @allocated sol =solve(prob,Euler(),dt=1//2^(6),save_timeseries=false)
alloc2 = @allocated sol2 =solve(prob2,Euler(),dt=1//2^(6),save_timeseries=false)

alloc1 = @allocated sol =solve(prob,Euler(),dt=1//2^(6),save_timeseries=false)
alloc2 = @allocated sol2 =solve(prob2,Euler(),dt=1//2^(6),save_timeseries=false)

@test alloc2 <= alloc1

sol = solve(prob_ode_linear,Euler(),dt=1//2^(6),save_timeseries=true)
sol2 = solve(prob_ode_linear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_timeseries=true)

sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(6),save_timeseries=true)
sol2 = solve(prob_ode_large2Dlinear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_timeseries=true)

sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(6),save_timeseries=true)
alloc1 = @allocated sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(8),save_timeseries=true)
alloc2 = @allocated sol2 = solve(prob_ode_large2Dlinear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_timeseries=true)

@test alloc2 <= alloc1

sol =solve(prob,Tsit5())


text:
### add modules
using Interpolations, ColorVectorSpace, ImageFiltering, Images
using Images
using TestImages

### setup
function imWarp(img, dx, dy)
    itp = interpolate(img, BSpline(Linear()))
    inds = indices_spatial(img)
    rng = extrema.(inds)
    imw = similar(img, eltype(itp))
    for I in CartesianIndices(inds)
        dxi, dyi = dx[I], dy[I]
        y, x = clamp(I[1]+dyi, rng[1]...), clamp(I[2]+dxi, rng[2]...)
        imw[I] = itp(y, x)
    end
    return imw
end

# get image
n = 20
mgrid = ones(512, 512);
@views mgrid[1:n:end, :] .= 0;
@views mgrid[:, 1:n:end] .= 0;
img = mgrid
# fig = PyPlot.figure(dpi = 300, figsize = (5, 5))
# PyPlot.imshow(mgrid, cmap = :gray); gcf()

#rand flow
flow = gen_tiled_flow(size(mgrid), 50);

#center flow:
ix = size(img)[1]
iy = size(img)[2]
ixhalf = Int64(ix/2)
iyhalf = Int64(iy/2)
flow = Array{ComplexF64}(undef, ix, iy)
flow[1:ixhalf, 1:iyhalf] .= -20 .+ -20 .* im;
flow[1:ixhalf, iyhalf:iy] .= 20 .+ -20 .* im;
flow[ixhalf:ix, 1:iyhalf] .= -20 .+ 20 .* im;
flow[ixhalf:ix, iyhalf:iy] .= 20 .+ 20 .* im;

dx, dy = real(flow), imag(flow);
imgw = imWarp(mgrid, dx, dy);

fig = PyPlot.figure(dpi = 300, figsize = (5, 5))
PyPlot.imshow(imgw, cmap = :gray);
ax = gca()
ax.get_ylim()
ax.set_ylim(0, size(imgw)[1])

flow[1:ixhalf, 1:iyhalf] .= 0;

showflow(flow, 20, fig=fig);

flow[1,1]
gcf()

showflow(flow, 20);
gcf()

u = gen_tiled_flow((51, 123), 10)
showflow(u)

1
LAP_julia.visualise.imgshowflow(imgw, flow)

LAP_julia.visualise.warp_imgshowflow(img, flow .* (-1))


text:
using AxisArrays
using AxisKeys
using CSV
using Combinatorics
using DataFrames
using Dates
using Distances
using Documenter
using HypothesisTests
using LinearAlgebra
using Random
using Statistics
using StatsBase
using TableOperations
using Tables
using Test

using Impute
using Impute:
    Impute,
    Imputor,
    Chain,
    DropObs,
    DropVars,
    Interpolate,
    Fill,
    KNN,
    LOCF,
    NOCB,
    Replace,
    SRS,
    DeclareMissings,
    Substitute,
    WeightedSubstitute,
    SVD,
    Filter,
    Threshold,
    WeightedThreshold,
    ThresholdError,
    apply,
    apply!,
    impute,
    impute!,
    interp,
    run,
    threshold,
    wthreshold,
    validate

@testset "Impute" begin
    include("testutils.jl")

    include("validators.jl")
    include("declaremissings.jl")
    include("chain.jl")
    include("data.jl")
    include("deprecated.jl")
    include("filter.jl")
    include("imputors/interp.jl")
    include("imputors/knn.jl")
    include("imputors/locf.jl")
    include("imputors/nocb.jl")
    include("imputors/replace.jl")
    include("imputors/srs.jl")
    include("imputors/substitute.jl")
    include("imputors/svd.jl")
    include("utils.jl")

    # Start running doctests before we wrap up technical changes and work
    # on more documentation
    # The doctests fail on x86, so only run them on 64-bit hardware & Julia 1.6
    Sys.WORD_SIZE == 64 && v"1.6" <= VERSION < v"1.7" && doctest(Impute)
end


text:
using Dates
using LibGit2
using Pkg
using Pkg: GitTools
using UUIDs
using DelimitedFiles
using SHA
# using Logging

user_original_directory = pwd()                 # where are we right now (this is the place that we have opend our terminal)
project_root = @__DIR__                         # this is the root directory of the program.

# cd(project_root)

# Internal constants

# =========================================> change to the better version ================
include("./modules/inputs.jl")
# include(joinpath(project_root,SCREIPTS_DIR_NAME,"inputs.jl"))
include("./modules/functions.jl")
# include(joinpath(project_root,SCREIPTS_DIR_NAME,"functions.jl"))

if CLEAN_INSTALLATION
    remove_dir!(PREFIX)
end

mkpath(PREFIX)

handle_dublicate_file(joinpath(BUILD_INFORMATION_PATH,"BUILD_LOG.log"))
julia_version_compatibility_check()
if LOGGING_TO_FILE
logprintln("== In case of crash, build log is available at: $(joinpath(BUILD_INFORMATION_PATH,"BUILD_LOG.log")).")
end

Saving_the_build_information(BUILD_INFORMATION_PATH,:start)
if SILENT
    logprintln("== Silent installation mode.")
end
logprintln("WARNING: DO NOT type your password into this program.")
@warn("DO NOT type your password into this program.")

# ====================== PHASE 1 =======================
# reading configuration.toml
# downloading files provided in configuration.toml
# saving the metadata

logprintln("== downloading Project.toml and Manifest.toml files...")
special_project_files = handle_files(configuration["include"],"projects",SP_PROJECTS_PATH)
@debug("DEBUG: Project Addresses: ", special_project_files)
special_manifest_files = handle_files(configuration["include"],"manifests",SP_MANIFESTS_PATH)
@debug("DEBUG: Manifest Addresses: ", special_manifest_files)
keep_online_files = handle_files(configuration["include"],"keep_online",KEEP_ONLINE_PATH)
remove_from_repo_files = handle_files(configuration["include"],"remove_from_repo",REMOVE_FROM_REPO_PATH)

logprintln("    >> Saving metadata...")

if !isempty(special_manifest_files) save_data(special_manifest_files,META_DATA_PATH,"SPMANIFESTS_PATH.DATA") end
if !isempty(special_project_files)  save_data(special_project_files,META_DATA_PATH,"SPPROJECTS_PATH.DATA")   end
if !isempty(keep_online_files)  save_data(keep_online_files,META_DATA_PATH,"KEEPONLINE_PATH.DATA")   end
if !isempty(remove_from_repo_files)  save_data(remove_from_repo_files,META_DATA_PATH,"REMOVEFROMREPO_PATH.DATA")   end
# logprintln("<> all done. successfully downloaded all Project.toml and Manifest.toml files")
logprintln("<> Done.")

logprintln("== Cloning registeries...")
cloned_registeries_path = cloning_registeries(configuration["include"],"registries",SP_REG_PATH)
@debug("DEBUG: cloned_registeries_path: ",cloned_registeries_path)
logprintln("    >> Saving metadata...")
save_data(cloned_registeries_path,META_DATA_PATH,"CLONED_REGISTERIES_PATH.DATA")
logprintln("<> Done.")

# ======================== PHASE TWO ==========================
# Reading metadata
# Preparing a download list

spr_path, spm_path, spp_path, ko_path, rm_path = read_meta_data_phase2()

# I found more conviniet to work with a simpler dictionary structure.
# reconstruction of packages metadata.
all_packages_names_dict = Dict{String,String}()
all_packages_addresses_dict = Dict{String,String}()
all_packages_urls_dict = Dict{String,String}()
download_list = String[]
download_list_urls_dict = Dict{String,String}()

logprintln("== Preparing metadata...")
fetch_packages_names_and_addresses!(all_packages_names_dict, all_packages_addresses_dict, spr_path)
packages_metadata = Dict("names" => all_packages_names_dict, "addresses" => all_packages_addresses_dict)
get_urls!(all_packages_urls_dict,packages_metadata["addresses"])
if rm_path != nothing
    remove_from_repo_packages = read_remove_from_repo_packages(rm_path,packages_metadata)
else
    remove_from_repo_packages = nothing
end

if configuration["options"]["download_all_packages"]
    logprintln("== All repo mode ...")
    download_list = collect(keys(packages_metadata["names"]))
else
    logprintln("== Partial repo mode ...")
    if spp_path != nothing
        read_packages_in_projects!(download_list,spp_path)
    end
    if spm_path != nothing
        read_packages_in_manifests!(download_list,spm_path)
    end
    remove_stdlib_packages!(download_list)
end

logprintln("    >> Preparing download list ...")
if remove_from_repo_packages != nothing
    remove!(download_list,remove_from_repo_packages)
end
if ko_path != nothing
    remove_keep_online_packages!(download_list,ko_path)
end
make_download_list!(download_list_urls_dict,download_list,all_packages_urls_dict)

# ============== Phase 3 [Downloading] ==========

logprintln("== Downloading packages...")
broken_packages_dict , downloaded_packages_path_dict, downloaded_packages_rel_path_dict = make_repo(download_list_urls_dict,packages_metadata)

handel_brocken_packages_metadata(broken_packages_dict,packages_metadata)

# =========== Phase 4 [Making new repo] =========

if MAKE_LOCAL_ACCESS_REG
    logprintln("== Making new registry for local access ...")
    type = "local"
    reg_root = joinpath(NEW_REG_PATH, LOCAL_REG_DIR)
    make_reg(packages_metadata,
                download_list_urls_dict,
                broken_packages_dict,
                reg_root,
                type,
                downloaded_packages_path_dict,
                downloaded_packages_rel_path_dict)
end

if INSTALL_ON_NETWORK
    logprintln("== Making new registry for public access ...")
    type = "public_access"
    reg_root = joinpath(NEW_REG_PATH, PUBLIC_REG_DIR)
        make_reg(packages_metadata,
        download_list_urls_dict,
        broken_packages_dict,
        reg_root,
        type,
        downloaded_packages_path_dict,
        downloaded_packages_rel_path_dict)
end

logprintln("""
==========================================
               Finished
==========================================
""")
Saving_the_build_information(BUILD_INFORMATION_PATH,:finish)


text:
using JuMP, Gurobi

mutable struct Node
    Level::Int
    Model::JuMP.Model
    Zbound::Float64
    Xrelax::Array{Float64}
    Status::Symbol
end

mutable struct Best
    Zstar::Float64
    Zbound::Float64
    Xstar::Array{Float64}
    Intsols::Int
    Visited::Int
end

function IndeciseVariable(node::Node)
    sizeX = length(node.Xrelax)
    nonInteger = Array{Float64}(sizeX)
    for i=1:sizeX
        if node.Model.colCat[i] == :Bin
            nonInteger[i] = abs(node.Xrelax[i] - 0.5)
        else
            nonInteger[i] = 2 #Variable is not binary
        end
    end
    mostInDoubt = indmin(nonInteger) #get index of the most fractionary variable
    if nonInteger[mostInDoubt] == 0.5  #all variables are 0 or 1
        return false
    end
    return mostInDoubt
end

function Bound(node::Node, best::Best, sense::Symbol)
    errortolerance = 1e-4 #if 0.01% away from optimal -> current answer is optimal
    if node.Status != :Optimal #Bound by Infeasabilty
        return false
    end
    if sense == :Max #Bound by Limit
        if node.Zbound < best.Zbound
            return false
        end
    else
        if node.Zbound > best.Zbound
            return false
        end
    end
    if abs(node.Zbound - best.Zstar) <= errortolerance #Bound by Optimality
        return false
    end
end

function Branch(node::Node)
    indeciseVariable = IndeciseVariable(node)
    if indeciseVariable == false #if all binary variables are integer -> no branching
        return false, false
    end

    leftChild = deepcopy(node.Model)
    # Forces the most indecise variable to be 0
    leftChild.colUpper[indeciseVariable] = 0
    leftChild.colLower[indeciseVariable] = 0

    rightChild = deepcopy(node.Model)
    # Forces the most indecise variable to be 0
    rightChild.colUpper[indeciseVariable] = 1
    rightChild.colLower[indeciseVariable] = 1

    return leftChild, rightChild
end

function InitializeHeadNode(model::JuMP.Model, sense::Symbol)
            #Level,    Model,   Zbound,             Xrelax,                 Status
    node = Node(0, deepcopy(model), Inf, Array{Float64}(length(model.colUpper)), :None)
    if sense == :Min
        node.Zbound = -Inf
    end
    return node
end

function InitializeBest(model::JuMP.Model)
    sense = getobjectivesense(model)
               #Zstar, Zbound        Xstar,                    Intsols    Visited
    best = Best( -Inf, Inf, Array{Float64}(length(model.colUpper)), 0 ,0) #initialize Best as a Max problem
    if sense == :Min
        best.Zstar = Inf
        best.Zbound = -Inf
    end
    return best
end

function SolveRelax(model::JuMP.Model, solver::MathProgBase.AbstractMathProgSolver)
    setsolver(model, solver)
    status = solve(model, relaxation = true) #solve model relaxation
    return status, getobjectivevalue(model), model.colVal
end

function UpdateBest(best::Best, node::Node, sense::Symbol, binaryVariables::Array{Int})
    if sense == :Max
        if (node.Zbound > best.Zstar && all(isinteger, node.Xrelax[binaryVariables])) #if the bound is better and if all binary variables are integer
            best.Zstar = node.Zbound
            best.Zbound = node.Zbound
            best.Xstar = node.Xrelax
            best.Intsols = best.Intsols + 1
        end
        if (node.Zbound < best.Zbound && node.Zbound > best.Zstar)
            best.Zbound = node.Zbound
        end
    else
        if (node.Zbound < best.Zstar && all(isinteger, node.Xrelax[binaryVariables]))
            best.Zstar = node.Zbound
            best.Zbound = node.Zbound
            best.Xstar = node.Xrelax
            best.Intsols = best.Intsols + 1
        end
        if (node.Zbound > best.Zbound && node.Zbound < best.Zstar)
            best.Zbound = node.Zbound
        end
    end
end

function BinaryVariables(model::JuMP.Model)
    binaryVariableIndexes = Array{Int}(0)
    for i=1:length(model.colCat)
        if model.colCat[i] == :Bin
            push!(binaryVariableIndexes,i) #gets index of all binary variables
        end
    end
    return binaryVariableIndexes
end

function SolveMIP(model::JuMP.Model)
    tic()
    t1 = time_ns()
    t2 = 0.01
    iter = 0
    maxiter = 1000 #max iterations
    tmax = 180 #max time in nanoseconds
    level = 0
    solver = GurobiSolver(OutputFlag=0)
    sense = getobjectivesense(model)
    best = InitializeBest(model)
    nodes = Array{Node}(1)
    binaryVariableIndexes = BinaryVariables(model) #gets index of all binary variables
    nodes[1] = InitializeHeadNode(model, sense) #root node
    nodes[1].Status, nodes[1].Zbound, nodes[1].Xrelax = SolveRelax(nodes[1].Model,solver)
    model.ext[:status] = nodes[1].Status
    UpdateBest(best, nodes[1], sense, binaryVariableIndexes) #updates best answer
    best.Visited = best.Visited + 1

    while (!isempty(nodes) && iter <= maxiter && t2 <= tmax)
        level = level + 1
        leftChild, rightChild = Branch(nodes[end])
        pop!(nodes)

        if leftChild != false
            statusLeftChild, boundLeftChild, xRelaxLeftChild = SolveRelax(leftChild, solver)
            best.Visited = best.Visited + 1
            nodeLeftChild = Node(level, leftChild, boundLeftChild, xRelaxLeftChild, statusLeftChild)
            UpdateBest(best, nodeLeftChild, sense, binaryVariableIndexes)
            if Bound(nodeLeftChild, best, sense) != false
                push!(nodes, nodeLeftChild)
            end
        end

        if rightChild != false
            statusRightChild, boundRightChild, xRelaxRightChild = SolveRelax(rightChild, solver)
            best.Visited = best.Visited + 1
            nodeRightChild = Node(level, rightChild, boundRightChild, xRelaxRightChild, statusRightChild)
            UpdateBest(best,nodeRightChild, sense, binaryVariableIndexes)
            if Bound(nodeRightChild, best, sense) != false
                push!(nodes, nodeRightChild)
            end
        end
        iter = iter + 1
        if iter%50 == 0
            t2 = (time_ns() - t1)/1e9
        end
    end

    if ((iter >= maxiter || t2 >= tmax) && all(isinteger, best.Xstar[binaryVariableIndexes])) #have found not optimal integer solutions
        model.ext[:status] = :SubOptimal
    elseif ((iter >= maxiter || t2 >= tmax) && ~(all(isinteger, best.Xstar[binaryVariableIndexes]))) #have not found any solutions
        model.ext[:status] = :NotSolved
    end

    model.ext[:nodes] = best.Visited
    model.ext[:intsols] = best.Intsols
    model.colVal = best.Xstar
    model.objVal = best.Zstar
    model.objBound = best.Zbound
    time = toc()
    model.ext[:time] = time
end


text:
module ConformingFESpacesTests

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces

# testing compute_conforming_cell_dofs

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

order = 2
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]
conf = GradConformity()

face_labeling = get_face_labeling(model)
dirichlet_tags = ["tag_1","tag_6"]

cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
  reffes,conf,grid_topology, face_labeling, dirichlet_tags)

r = [
  [-1,1,4,5,14,15,16,17,35],[1,2,5,6,18,19,17,20,36],[2,3,6,7,21,22,20,23,37],
  [4,5,8,9,15,24,25,26,38],[5,6,9,10,19,27,26,28,39],[6,7,10,11,22,29,28,30,40],
  [8,9,12,-2,24,-4,31,32,41],[9,10,-2,-3,27,-5,32,33,42],[10,11,-3,13,29,-6,33,34,43]]
test_array(cell_dofs,r)
@test nfree == 43
@test ndiri == 6
@test dirichlet_dof_tag == [1, 2, 2, 2, 2, 2]
@test dirichlet_cells == [1, 7, 8, 9]

order = 1
reffes = [LagrangianRefFE(VectorValue{2,Float64},p,order) for p in polytopes]

dirichlet_components = [(true,true), (false,true)]

cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
  reffes,conf,grid_topology, face_labeling, dirichlet_tags, dirichlet_components)

r = [
  [-1,1,7,9,-2,2,8,10],[1,3,9,11,2,4,10,12],[3,5,11,13,4,6,12,14],
  [7,9,15,17,8,10,16,18],[9,11,17,19,10,12,18,20],[11,13,19,21,12,14,20,22],
  [15,17,23,25,16,18,24,-3],[17,19,25,26,18,20,-3,-4],[19,21,26,27,20,22,-4,28]]

test_array(cell_dofs,r)
@test nfree==28
@test ndiri==4
@test dirichlet_dof_tag == [1, 1, 2, 2,]
@test dirichlet_cells == [1, 7, 8, 9]

order = 3
reffes = [LagrangianRefFE(VectorValue{2,Float64},p,order) for p in polytopes]

dirichlet_components = [(true,true), (false,true)]

cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
  reffes, conf, grid_topology, face_labeling, dirichlet_tags, dirichlet_components)

V = GradConformingFESpace(reffes,model,dirichlet_tags)
test_single_field_fe_space(V)

matvecdata = ([],[],[])
matdata = ([],[],[])
vecdata = ([],[])
test_single_field_fe_space(V,matvecdata,matdata,vecdata)

V = GradConformingFESpace(reffes,model,dirichlet_tags,dirichlet_components)
test_single_field_fe_space(V)

V = ConformingFESpace(reffes,model,get_face_labeling(model),dirichlet_tags,dirichlet_components)
test_single_field_fe_space(V)

end  # module


text:
import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp

machine = (function ()
    # Primitives
    label      = re"[^\(\);,:]+"
    len        = re":[^\(\);,:]+"
    cladestart = re"\("
    cladestop  = re"\)"
    sep        = re","
    finish     = re";"

    # Composites
    metadata   = (label * len) | label | len
    nodechange = cladestart | cladestop | sep
    element    = nodechange * re.opt(metadata)
    tree       = re.opt(metadata) * re.rep(element) * finish

    # Actions
    label.actions[:enter] = [:mark]
    label.actions[:exit]  = [:label]
    len.actions[:enter] = [:mark]
    len.actions[:exit] = [:len]
    cladestart.actions[:enter] = [:push]
    cladestop.actions[:enter] = [:pop]
    sep.actions[:enter] = [:pop, :push]
    finish.actions[:enter] = [:finish]

    # Finally, compile the final FASTA pattern into a state machine.
    return Automa.compile(tree)
end)()

mutable struct Node
    label::Union{Nothing, String}
    len::Union{Nothing, Float64}
    parent::Union{Nothing, Node}
end

function Base.show(io::IO, node::Node)
    parent = node.parent === nothing ? nothing : "Node(\"$(node.parent.label)\")"
    print(io, "Node(\"$(node.label)\", $(node.len), $parent)")
end

@noinline function popstack!(stack::Vector{Node})
    length(stack) > 1 || error("Too many closing parenthesis in file")
    pop!(stack)
    return nothing
end

actions = Dict(
    :finish => :(finished = true),
    :mark => :(mark = p),
    :label => :(last(stack).label = String(data[mark:p-1])),
    :len => :(last(stack).len = Base.parse(Float64, String(data[mark+1:p-1]))),
    :pop => :(popstack!(stack)),
    :push => quote
        push!(stack, Node(nothing, nothing, last(stack)))
        push!(nodes, last(stack))
    end
)

context = Automa.CodeGenContext()
@eval function parse(data::Union{String,Vector{UInt8}})
    mark = 0
    nodes = [Node(nothing, nothing, nothing)]
    stack = copy(nodes)
    finished = false

    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine, actions))

    if (cs != 0) & !finished
        error("failed to parse on byte ", p)
    end

    # Final stack should only have the root node.
    if length(stack) != 1
        error("Not enough closing parentheses in file")
    end

    return nodes
end

max_stars_repo_name: jakobnissen/play

max_issues_repo_name: jakobnissen/play

max_forks_repo_name: jakobnissen/play


text:
module LogCompose

using Pkg, Logging

include("compose.jl")
include("connectors.jl")

end # module LogCompose


text:
using FinSight
using Dates
using Test
@testset "(x)irr and (x)npv functions" begin
# Examples from
# https://support.microsoft.com/en-us/office/npv-function-8672cb67-2576-4d07-b67b-ac28acf2a568
    @test isapprox(npv(0.1, [-10000, 3000, 4200, 6800]), 1188.44, rtol = 1e-3)
    @test isapprox(npv(0.08, [8000, 9200, 10000, 12000, 14500]) + (-40000), 1922.06, rtol=1e-3)

    let values = [-70e3, 12e3, 15e3, 18e3, 21e3, 26e3]
        @test isapprox(irr(values[1:end-1]), -0.021, atol=1e-3)
        @test isapprox(npv(irr(values), values), 0.0, atol=1e-6)
        v2 = @view values[1:5]
        @test isapprox(npv(irr(v2), v2), 0.0, atol=1e-6)
    end

    # From https://support.microsoft.com/en-us/office/xnpv-function-1b42bbf6-370f-4532-a0eb-d67c16b664b7
    let (dates, values) =
        ([DateTime(2008, 1, 1), DateTime(2008, 3, 1), DateTime(2008, 10, 30), DateTime(2009,2,15), DateTime(2009,4, 1)],
         [-10e3, 2750, 4250, 3250, 2750])
        @test isapprox(xnpv(0.09, values, dates), 2086.65, atol=1e-2)
        @test isapprox(xnpv(xirr(values, dates), values, dates), 0.0, atol=1e-6)
    end
end

max_stars_repo_name: davidtaylor-fastlanebi/FinSight

max_issues_repo_name: davidtaylor-fastlanebi/FinSight

max_forks_repo_name: davidtaylor-fastlanebi/FinSight


text:
#######################################################################################################################################################################################################
#
# Changes to this type
# General
#     2022-Jan-18: add abstract colimitation type
#     2022-Jan-25: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Hierachy of `AbstractColimit`
- [`MinimumColimit`](@ref)
- [`QuadraticColimit`](@ref)
- [`SerialColimit`](@ref)
"""
abstract type AbstractColimit{FT<:AbstractFloat} end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-14: add minimum colimitation struct
#     2022-Jan-25: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Empty structure to indicate minimum colimitation.

---
# Examples
```julia
col = MinimumColimit{Float64}();
```
"""
struct MinimumColimit{FT<:AbstractFloat} <:AbstractColimit{FT} end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-14: add quadratic colimitation struct
#     2022-Jan-25: add documentation
#     2022-Mar-01: add documentation to field
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure to indicate quadratic colimitation.

# Fields

$(TYPEDFIELDS)

---
# Examples
```julia
col = QuadraticColimit{Float64}(0.98);
```
"""
mutable struct QuadraticColimit{FT<:AbstractFloat} <: AbstractColimit{FT}
    "Curvature factor"
    CURVATURE::FT
end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Feb-28: add serial colimitation struct
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Empty structure to indicate serial colimitation.

---
# Examples
```julia
col = SerialColimit{Float64}();
```
"""
struct SerialColimit{FT<:AbstractFloat} <:AbstractColimit{FT} end


text:
module WordTokenizers

using HTML_Entities
using StrTables
using Unicode
using DataDeps

abstract type PretrainedTokenizer end

export poormans_tokenize, punctuation_space_tokenize,
       penn_tokenize, improved_penn_tokenize, nltk_word_tokenize,
       tweet_tokenize,
       tokenize,
       rulebased_split_sentences,
       split_sentences,
       set_tokenizer, set_sentence_splitter,
       rev_tokenize, rev_detokenize,
       toktok_tokenize
export ALBERT_V1, ALBERT_V2, load, tokenizer, sentence_from_tokens, ids_from_tokens
export PretrainedTokenizer, tokenizer_files
include("words/fast.jl")

include("words/simple.jl")
include("words/nltk_word.jl")
include("words/reversible_tokenize.jl")
include("words/sedbased.jl")
include("words/tweet_tokenizer.jl")
include("sentences/sentence_splitting.jl")
include("words/TokTok.jl")

include("set_method_api.jl")
include("split_api.jl")

include("statistical/unigram.jl")

const pretrained = Dict{DataType, Vector{String}}()
function tokenizer_files(::Type{T}) where T<:PretrainedTokenizer
    get!(pretrained,T) do
        String[]
    end
end

function __init__()
    include(joinpath(@__DIR__, "statistical/Vocab_DataDeps.jl"))
end
end # module


text:
module PostProcess

include("../Types.jl")

export run

function run(input:: Channel, time_threshold)
    current = nothing
    return function (output:: Channel)
        for dp in input
            if current === nothing
                current = dp
            else
                # If the time difference exceeds the threshold, we have a confirmed step
                if (dp.time - current.time) > time_threshold
                    current = dp
                    put!(output, dp)
                else
                    # Keep the point with the largest magnitude.
                    if dp.magnitude > current.magnitude
                        current = dp
                    end
                end
            end
        end
    end
end

end

max_stars_repo_name: eliabiraschi/pedometer

max_issues_repo_name: eliabiraschi/pedometer

max_forks_repo_name: eliabiraschi/pedometer


text:
using PtFEM

data = Dict(
  # Beam(ndim, nst, nxe, nip, direction, finite_element(nod, nodof), axisymmet
  :struc_el => Beam(2, 1, 4, 1, :x, Line(2, 2), false),
  :properties => [4.0e4; 2.0e4],
  :etype => [1, 1, 2, 2],
  :x_coords => [0.0, 2.5, 5.0, 8.0, 10.0],
  :support => [
      (1, [0 1]),
      (4, [0 1])
    ],
  :loaded_nodes => [
      (2, [-20.0 0.0]),
      (3, [-6.0 -3.0]),
      (4, [-8.8 2.2]),
      (5, [-1.2 0.5333])
    ],
  :fixed_freedoms => [
      (1, 2, -0.001),
      (3, 1, -0.005)
    ],
  :penalty => 1e19,
  :eq_nodal_forces_and_moments => [
    (3, [-6.0 -3.0 -6.0 3.0]),
    (4, [-2.8 -0.8 -1.2  0.5333])
  ]
)

data |> display
println()

@time fem, dis_df, fm_df = p43(data)
println()

display(dis_df)
println()
display(fm_df)


text:
function PulseDiff(; name, Val=1, dt=0.1)
    @named d = DigitalPin()
    @variables val(t)
    D = ModelingToolkit.Difference(t; dt=dt)

    eqs = [
        D(val) ~ Val
        val ~ d.val
    ]

    ODESystem(eqs, t, [val], [], systems=[d], defaults=Dict(Val=>0), name=name)
end

function Set(; name)
    @named d = DigitalPin()

    eqs = [
        d.val ~ 1
    ]
    ODESystem(eqs, t, [], [], systems=[d],  name=name)
end

function Reset(; name)
    @named d = DigitalPin()

    eqs = [
        d.val ~ 0
    ]
    ODESystem(eqs, t, [], [], systems=[d], name=name)
end

function Pulse(; name, duty_cycle=0.5, T=1.0)
    @named d = DigitalPin()

    eqs = [
        d.val ~ IfElse.ifelse(t%T > duty_cycle*T, 1, 0)
    ]
    ODESystem(eqs, t, [], [], systems=[d], name=name)
end


text:
module  TestTypesDatasets

using PopGen
using DataFrames
using PooledArrays
using Test

cats = @nancycats;
sharks = @gulfsharks;

@testset "Dataset types" begin
    @test typeof(cats) == PopData
    @test typeof(cats.meta) == DataFrame
    @test typeof(cats.loci) == DataFrame
    @test typeof(sharks) == PopData
    @test typeof(sharks.meta) == DataFrame
    @test typeof(sharks.loci) == DataFrame
end

@testset "Dataset dimensions" begin
    @test size(cats.meta) == (237,5)
    @test size(cats.loci) == (2133,4)
    @test size(sharks.meta) == (212,5)
    @test size(sharks.loci) == (468308,4)
end

@testset "Dataset column names" begin
    meta_names_sorted = ["latitude","longitude","name","ploidy","population"]
    loci_names_sorted = ["genotype","locus","name","population"]
    @test sort(names(cats.meta)) == meta_names_sorted
    @test sort(names(cats.loci)) == loci_names_sorted
    @test sort(names(sharks.meta)) == meta_names_sorted
    @test sort(names(sharks.loci)) == loci_names_sorted
end

@testset "Nancycats column types" begin
    @test typeof(cats.meta.name) == Vector{String}
    @test typeof(cats.meta.population) == Vector{String}
    @test typeof(cats.meta.ploidy) == Vector{Int8}
    @test typeof(cats.meta.latitude) ==  Vector{Union{Missing, Float32}}
    @test typeof(cats.meta.longitude) == Vector{Union{Missing, Float32}}
    @test typeof(cats.loci.name) <: PooledArray
    @test eltype(cats.loci.name) == String
    @test typeof(cats.loci.population) <: PooledArray
    @test eltype(cats.loci.population) == String
    @test typeof(cats.loci.locus) <: PooledArray
    @test eltype(cats.loci.locus) == String
    @test typeof(cats.loci.genotype) <: GenoArray
    @test eltype(cats.loci.genotype) <: Union{Missing, Genotype}
end

@testset "Gulfsharks column types" begin
    @test typeof(sharks.meta.name) == Vector{String}
    @test typeof(sharks.meta.population) == Vector{String}
    @test typeof(sharks.meta.ploidy) == Vector{Int8}
    @test typeof(sharks.meta.latitude) ==  Vector{Float64}
    @test typeof(sharks.meta.longitude) == Vector{Float64}
    @test typeof(sharks.loci.name) <: PooledArray
    @test eltype(sharks.loci.name) == String
    @test typeof(sharks.loci.population) <: PooledArray
    @test eltype(sharks.loci.population) == String
    @test typeof(sharks.loci.locus) <: PooledArray
    @test eltype(sharks.loci.locus) == String
    @test typeof(sharks.loci.genotype) <: GenoArray
    @test eltype(sharks.loci.genotype) <: Union{Missing, Genotype}
end

end  # module


text:
###################################
##    Cairn Client Functionality ##
###################################
module Client
import Cairn
using Logging
@Logging.configure(level=DEBUG)
import Rock
import SimpleFileServer
import SHA
immutable t
    id::Int64
    host::ASCIIString
    port::Int64
    t(host::ASCIIString, port::Int64) = new(view(host, port, Cairn.RPC.GetID()), host, port)
end
make(host, port) = t(host, port)
function command(c::t, comm::Cairn.RPC.Command)
    r = Rock.Client.make(c.host, c.port)
    Rock.Client.command(r, comm)
end
function view(host, port, v::Cairn.RPC.View)
    r = Rock.Client.make(host, port)
    Rock.Client.command(r, v)
end
function view(c::Cairn.Client.t, v::Cairn.RPC.View)
    r = Rock.Client.make(c.host, c.port)
    Rock.Client.command(r, v)
end
function delete(cli::Cairn.Client.t, name::ASCIIString)
    # First, un track the file, then delete it for real
    resp = Cairn.Client.command(cli, Cairn.RPC.DeleteRequest(name))::DeleteResponse
    if !isnull(resp.err)
        throw(resp.err)
    end
    # for (host, port) = nodes
    #     c = SimpleFileServer.Client.make(host, port)
    #     SimpleFileServer.Client.delete(c, name)
    # end
end
function command_result(c::t, command::Cairn.RPC.Command)
    arg  = Cairn.RPC.CacheResultOfAction(c.id,command)
    @error arg
    Cairn.Client.view(c, Cairn.RPC.PrepareCache(arg))
    l = Cairn.Client.command(c, arg)
    @error l
    fetch_cache = Cairn.RPC.FetchCache(arg, true)
    while true
        cache_result = view(c, fetch_cache)
        if cache_result.state == Cairn.RPC.Ready
            return cache_result.val
        elseif cache_result.state == Cairn.RPC.Gone
            error("Nothing in Cache")
        elseif cache_result.state == Cairn.RPC.Loading
            wait(Timer(1))
        end
    end
end
function download(c::Cairn.Client.t, key_span::Cairn.RPC.KeySpan,
                  to_dir::ASCIIString, rename::Vector{ASCIIString} = [],
                  stream::Nullable{Channel{ASCIIString}} = Nullable())
    files = command_result(c, Cairn.RPC.GetKeySpanRequest(key_span)).files
    @sync for r::Cairn.RPC.GetKeySpanSubResponse = files
        if isnull(r.err)
            @async for (hash, sfs_client) = zip(r.hashes, map(x -> SimpleFileServer.Client.make(x...), r.replicas))
                try
                    # Download the chunk, appending it to the file
                    name = joinpath(to_dir, length(rename) > 0 ? pop!(rename): r.key )
                    SimpleFileServer.Client.download(sfs_client, hash, name, true)
                    if !isnull(stream)
                        put!(stream,(r.key, name))
                    end

                catch err
                    @debug err
                    # TODO Better story here?
                    error("Could not locate file")
                end
            end

        else
            throw(r.err.value)
        end
    end
end
function upload(c::Cairn.Client.t, name::ASCIIString, m::Vector{UInt8}, replication::Int64)
    r = command_result(c, Cairn.RPC.CreateRequest(name, replication))::Cairn.RPC.CreateResponse
    if isnull(r.err)

        name = SHA.sha256(m)
        file = SimpleFileServer.File(name,name, length(m))
        n = RemoteRef()
        put!(n, 0)
        @sync for client = map(a ->SimpleFileServer.Client.make(a...), r.replicas)
            @async begin
                try
                    SimpleFileServer.Client.upload(client, file, m)

                    put!(n, take!(n)+1)
                catch err
                    @error err
                    # One Retry
                    try
                        wait(Timer(3))
                        SimpleFileServer.Client.upload(client, file, m)
                        put!(n, take!(n)+1)
                    catch err
                        @error err
                    end

                end
            end
        end
        v = take!(n)
        if v != replication
            @warn "Could not get the requested number of servers, got $n"
        end
        v
    else
        throw(r.err.value)
    end
end
function find_connection(cs::Vector{Cairn.Client.t})
    Base.Random.shuffle!(cs)
    cs[1]
    # TODO; Ping, find slave, etc
end

end


text:
abstract type AbstractStateChanges end

struct StateChanges{A<:AbstractArray,B<:AbstractArray} <: AbstractStateChanges
    to_change::A
    new_values::B
end

function StateChanges(state::State)
    return StateChanges(Int[], eltype(state)[])
end

function Base.resize!(c::StateChanges, n)
    resize!(c.to_change,  n)
    resize!(c.new_values, n)
end

zero!(c::StateChanges) = resize!(c, 0)

Base.eltype(c::StateChanges{A,B}) where {A,B} = (eltype(A), eltype(B))
Base.length(c::StateChanges) = length(c.to_change)
Base.size(c::StateChanges) = size(c.to_change)

function Base.:(==)(a::StateChanges, b::StateChanges)
    sk_a = sortperm(a.to_change)
    sk_b = sortperm(b.to_change)

    return ( a.to_change[sk_a]  == b.to_change[sk_b] &&
             a.new_values[sk_a] == b.new_values[sk_b] )
end

function Base.push!(c::StateChanges, (to_change, new_values))
    push!(c.to_change, to_change)
    push!(c.new_values, new_values)
    return c
end

function Base.append!(c::StateChanges, (to_change, new_values))
    append!(c.to_change, to_change)
    append!(c.new_values, new_values)
    return c
end

Base.getindex(c::StateChanges, i) = (c.to_change[i], c.new_values[i])

function Base.iterate(iter::StateChanges, state=1)
    state > length(iter) && return nothing
    return (iter[state], state + 1)
end

function Base.show(io::IO, ch::StateChanges)
    print(io, "StateChanges{$(eltype(ch))} : [")
    content = ""
    for (id, val)=ch
        content *= "$id->$val, "
    end
    print(io, content[1:end-2]*"]")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Birch_Standard"
# Birch's version is `(git describe --long || echo) | sed -E 's/v([0-9]+)-([0-9]+)-g[0-9a-f]+/\1.\2/'`
# (see https://github.com/JuliaPackaging/Yggdrasil/pull/2156#discussion_r528224861)
version = v"1.104"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/lawmurray/Birch.git", "25154eac433bdeaf55dffda2432c666659286350"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
# add libraries for building birch on host
apk add boost-dev yaml-dev

# build and install driver on host system
cd ${WORKSPACE}/srcdir/Birch/driver/
mkdir -p ${WORKSPACE}/srcdir/build_host
autoreconf -vi
CXX=${HOSTCXX} CPPFLAGS=-I/usr/include LDFLAGS=-L/usr/lib ./configure --prefix=${WORKSPACE}/srcdir/build_host --build=${MACHTYPE} --host=${MACHTYPE}
make clean
make -j${nproc}
make install

# build and install the standard library
cd ${WORKSPACE}/srcdir/Birch/libraries/Standard/
export BIRCH_PREFIX=${WORKSPACE}/srcdir/build_host
export PATH=${WORKSPACE}/srcdir/build_host/bin:$PATH
# musl c library does not define stdin, stdout, and stderr macros as objects
if [[ "${target}" == *-musl* ]]; then
  atomic_patch -p3 ${WORKSPACE}/srcdir/stdio.patch
fi
birch bootstrap
CPPFLAGS="-I${prefix}/include" ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-debug --enable-release
make -j${nproc}
make install

# install license
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# Native support for Windows is not yet provided
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libbirch-standard", :libbirch_standard),
    LibraryProduct("libbirch-standard-debug", :libbirch_standard_debug),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("boost_jll"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency("Eigen_jll"),
    Dependency(PackageSpec(name = "LibBirch_jll", version = version)),
    Dependency("LibYAML_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               preferred_gcc_version=v"5")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
include("./kernel_system.jl")

function test_data(d, p, m, capacity)
    U = rand(d, p)
    Fu = rand(p, p)
    X = data_array(rand(d, capacity), m)
    Fx = data_array(rand(p, capacity), m)
    y = data_array(rand(capacity), p + m)

    return U, Fu, X, Fx, y
end

function test_const(U, Fu, X, Fx, capacity)
    d, p = size(U)
    m = size(X, 2)
    k = FixedParam(Gaussian(), 0.01)
    K = Matrix{Float64}(undef, m + p, m + p)
    correlation!(K, k, [U X])
    M = [zeros(p, p) [Fu Fx]; [Fu Fx]' K]
    sd = system_data!(U, Fu, X, Fx)
    Ku = Matrix{Float64}(undef, p, p)
    Kux = data_array(Float64, capacity, p)
    Kx = incremental_cholesky(Float64, capacity)
    kd = kernel_data!(k, Ku, Kux, Kx, sd.U)
    compute_next!(kd, sd)
    ks = KernelSystem(sd, kd)

    return ks, M
end

function test_solve(ks, M, y)
    cu, cx = solve!(similar(y), ks, y)
    d = fit_tail!(similar(cu), ks, y[1:length(cu)], cu, cx)
    c = M \ [zeros(length(cu)); y]
    display([d; cu; cx])
    display(c)
    return nothing
end

function test_system()
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@doc raw"""

"""
module LotkaVolterra4d

    using GeometricEquations
    using Parameters

    export hamiltonian, ϑ, ϑ₁, ϑ₂, ω

    export lotka_volterra_4d_ode,
           lotka_volterra_4d_pode, lotka_volterra_4d_pdae,
           lotka_volterra_4d_iode, lotka_volterra_4d_idae,
           lotka_volterra_4d_lode, lotka_volterra_4d_ldae,
           lotka_volterra_4d_dg

    Δt = 0.1
    nt = 1000

    const q₀ = [2.0, 1.0, 1.0, 1.0]
    const p  = (a₁=1.0, a₂=1.0, a₃=1.0, a₄=1.0, b₁=-1.0, b₂=-2.0, b₃=-1.0, b₄=-1.0)

    const reference_solution = [1.6390462434739954, 1.3764800055785835, 0.37903204434372284, 1.4399236281802124]

    # const q₀ = [2.0, 1.0, 2.0, 1.0]
    # const p  = (a₁=1.0, a₂=1.0, a₃=1.0, a₄=1.0, b₁=-1.0, b₂=-2.0, b₃=-1.0, b₄=-2.0)

    # const q₀ = [2.0, 1.0, 2.0, 1.0]
    # const p  = (a₁=1.0, a₂=1.0, a₃=1.0, a₄=1.0, b₁=-1.0, b₂=-4.0, b₃=-2.0, b₄=-3.0)

    ϑ₁(t, q) = 0.5 * ( + log(q[2]) - log(q[3]) + log(q[4]) ) / q[1]
    ϑ₂(t, q) = 0.5 * ( - log(q[1]) + log(q[3]) - log(q[4]) ) / q[2]
    ϑ₃(t, q) = 0.5 * ( + log(q[1]) - log(q[2]) + log(q[4]) ) / q[3]
    ϑ₄(t, q) = 0.5 * ( - log(q[1]) + log(q[2]) - log(q[3]) ) / q[4]

    # ϑ₁(t, q) = ( log(q[2]) + log(q[4]) ) / q[1]
    # ϑ₂(t, q) = ( log(q[3]) ) / q[2]
    # ϑ₃(t, q) = ( log(q[1]) + log(q[4]) ) / q[3]
    # ϑ₄(t, q) = ( log(q[2]) ) / q[4]

    # ϑ₁(t, q) = ( + log(q[2]) - log(q[3]) + log(q[4]) ) / q[1] / 2 + q[2] + q[3] + q[4]
    # ϑ₂(t, q) = ( - log(q[1]) + log(q[3]) - log(q[4]) ) / q[2] / 2 + q[1] + q[3] + q[4]
    # ϑ₃(t, q) = ( + log(q[1]) - log(q[2]) + log(q[4]) ) / q[3] / 2 + q[1] + q[2] + q[4]
    # ϑ₄(t, q) = ( - log(q[1]) + log(q[2]) - log(q[3]) ) / q[4] / 2 + q[1] + q[2] + q[3]

    # ϑ₁(t, q) = ( + log(q[2]) - log(q[3]) + log(q[4]) ) / q[1]
    # ϑ₂(t, q) = ( + log(q[3]) - log(q[4]) ) / q[2]
    # ϑ₃(t, q) = log(q[4]) / q[3]
    # ϑ₄(t, q) = zero(eltype(q))

    function v₁(t, q, params)
        @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄ = params
        q[1] * (+ a₂*q[2] + a₃*q[3] + a₄*q[4] + b₂ + b₃ + b₄)
    end

    function v₂(t, q, params)
        @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄ = params
        q[2] * (- a₁*q[1] + a₃*q[3] + a₄*q[4] - b₁ + b₃ + b₄)
    end

    function v₃(t, q, params)
        @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄ = params
        q[3] * (- a₁*q[1] - a₂*q[2] + a₄*q[4] - b₁ - b₂ + b₄)
    end

    function v₄(t, q, params)
        @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄ = params
        q[4] * (- a₁*q[1] - a₂*q[2] - a₃*q[3] - b₁ - b₂ - b₃)
    end

    dϑ₁dx₁(t, q) = ( - log(q[2]) + log(q[3]) - log(q[4]) ) / q[1]^2 / 2
    dϑ₁dx₂(t, q) = + 1 / (q[1] * q[2]) / 2
    dϑ₁dx₃(t, q) = - 1 / (q[1] * q[3]) / 2
    dϑ₁dx₄(t, q) = + 1 / (q[1] * q[4]) / 2

    dϑ₂dx₁(t, q) = - 1 / (q[2] * q[1]) / 2
    dϑ₂dx₂(t, q) = ( + log(q[1]) - log(q[3]) + log(q[4]) ) / q[2]^2 / 2
    dϑ₂dx₃(t, q) = + 1 / (q[2] * q[3]) / 2
    dϑ₂dx₄(t, q) = - 1 / (q[2] * q[4]) / 2

    dϑ₃dx₁(t, q) = + 1 / (q[3] * q[1]) / 2
    dϑ₃dx₂(t, q) = - 1 / (q[3] * q[2]) / 2
    dϑ₃dx₃(t, q) = ( - log(q[1]) + log(q[2]) - log(q[4]) ) / q[3]^2 / 2
    dϑ₃dx₄(t, q) = + 1 / (q[3] * q[4]) / 2

    dϑ₄dx₁(t, q) = - 1 / (q[4] * q[1]) / 2
    dϑ₄dx₂(t, q) = + 1 / (q[4] * q[2]) / 2
    dϑ₄dx₃(t, q) = - 1 / (q[4] * q[3]) / 2
    dϑ₄dx₄(t, q) = ( + log(q[1]) - log(q[2]) + log(q[3]) ) / q[4]^2 / 2

    # dϑ₁dx₁(t, q) = - ( log(q[2]) + log(q[4]) ) / q[1]^2
    # dϑ₁dx₂(t, q) = 1 / (q[1] * q[2])
    # dϑ₁dx₃(t, q) = zero(eltype(q))
    # dϑ₁dx₄(t, q) = 1 / (q[1] * q[4])

    # dϑ₂dx₁(t, q) = zero(eltype(q))
    # dϑ₂dx₂(t, q) = - ( log(q[3]) ) / q[2]^2
    # dϑ₂dx₃(t, q) = 1 / (q[2] * q[3])
    # dϑ₂dx₄(t, q) = zero(eltype(q))

    # dϑ₃dx₁(t, q) = + 1 / (q[3] * q[1])
    # dϑ₃dx₂(t, q) = zero(eltype(q))
    # dϑ₃dx₃(t, q) = - ( log(q[1]) + log(q[4]) ) / q[3]^2
    # dϑ₃dx₄(t, q) = + 1 / (q[3] * q[4])

    # dϑ₄dx₁(t, q) = zero(eltype(q))
    # dϑ₄dx₂(t, q) = 1 / (q[4] * q[2])
    # dϑ₄dx₃(t, q) = zero(eltype(q))
    # dϑ₄dx₄(t, q) = - ( log(q[2]) ) / q[4]^2

    # dϑ₁dx₁(t, q) = ( - log(q[2]) + log(q[3]) - log(q[4]) ) / q[1]^2 / 2
    # dϑ₁dx₂(t, q) = 1 + 1 / (q[1] * q[2]) / 2
    # dϑ₁dx₃(t, q) = 1 - 1 / (q[1] * q[3]) / 2
    # dϑ₁dx₄(t, q) = 1 + 1 / (q[1] * q[4]) / 2

    # dϑ₂dx₁(t, q) = 1 - 1 / (q[2] * q[1]) / 2
    # dϑ₂dx₂(t, q) = ( + log(q[1]) - log(q[3]) + log(q[4]) ) / q[2]^2 / 2
    # dϑ₂dx₃(t, q) = 1 + 1 / (q[2] * q[3]) / 2
    # dϑ₂dx₄(t, q) = 1 - 1 / (q[2] * q[4]) / 2

    # dϑ₃dx₁(t, q) = 1 + 1 / (q[3] * q[1]) / 2
    # dϑ₃dx₂(t, q) = 1 - 1 / (q[3] * q[2]) / 2
    # dϑ₃dx₃(t, q) = ( - log(q[1]) + log(q[2]) - log(q[4]) ) / q[3]^2 / 2
    # dϑ₃dx₄(t, q) = 1 + 1 / (q[3] * q[4]) / 2

    # dϑ₄dx₁(t, q) = 1 - 1 / (q[4] * q[1]) / 2
    # dϑ₄dx₂(t, q) = 1 + 1 / (q[4] * q[2]) / 2
    # dϑ₄dx₃(t, q) = 1 - 1 / (q[4] * q[3]) / 2
    # dϑ₄dx₄(t, q) = ( + log(q[1]) - log(q[2]) + log(q[3]) ) / q[4]^2 / 2

    # dϑ₁dx₁(t, q) = ( - log(q[2]) + log(q[3]) - log(q[4]) ) / q[1]^2
    # dϑ₁dx₂(t, q) = + 1 / (q[1] * q[2])
    # dϑ₁dx₃(t, q) = - 1 / (q[1] * q[3])
    # dϑ₁dx₄(t, q) = + 1 / (q[1] * q[4])

    # dϑ₂dx₁(t, q) = zero(eltype(q))
    # dϑ₂dx₂(t, q) = ( - log(q[3]) + log(q[4]) ) / q[2]^2
    # dϑ₂dx₃(t, q) = + 1 / (q[2] * q[3])
    # dϑ₂dx₄(t, q) = - 1 / (q[2] * q[4])

    # dϑ₃dx₁(t, q) = zero(eltype(q))
    # dϑ₃dx₂(t, q) = zero(eltype(q))
    # dϑ₃dx₃(t, q) = ( - log(q[4]) ) / q[3]^2
    # dϑ₃dx₄(t, q) = + 1 / (q[3] * q[4])

    # dϑ₄dx₁(t, q) = zero(eltype(q))
    # dϑ₄dx₂(t, q) = zero(eltype(q))
    # dϑ₄dx₃(t, q) = zero(eltype(q))
    # dϑ₄dx₄(t, q) = zero(eltype(q))

    function ϑ(t::Number, q::AbstractVector, Θ::AbstractVector)
        Θ[1] = ϑ₁(t,q)
        Θ[2] = ϑ₂(t,q)
        Θ[3] = ϑ₃(t,q)
        Θ[4] = ϑ₄(t,q)
        nothing
    end

    function ϑ(t::Number, q::AbstractVector)
        Θ = zero(q)
        ϑ(t, q, Θ)
        return Θ
    end

    function ϑ(t::Number, q::AbstractVector, k::Int)
        if k == 1
            ϑ₁(t, q)
        elseif k == 2
            ϑ₂(t, q)
        elseif k == 3
            ϑ₃(t, q)
        elseif k == 4
            ϑ₄(t, q)
        else
            throw(BoundsError(ϑ,k))
        end
    end

    function ω(t, q, Ω)
        Ω[1,1] = 0
        Ω[1,2] = dϑ₁dx₂(t,q) - dϑ₂dx₁(t,q)
        Ω[1,3] = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
        Ω[1,4] = dϑ₁dx₄(t,q) - dϑ₄dx₁(t,q)

        Ω[2,1] = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)
        Ω[2,2] = 0
        Ω[2,3] = dϑ₂dx₃(t,q) - dϑ₃dx₂(t,q)
        Ω[2,4] = dϑ₂dx₄(t,q) - dϑ₄dx₂(t,q)

        Ω[3,1] = dϑ₃dx₁(t,q) - dϑ₁dx₃(t,q)
        Ω[3,2] = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
        Ω[3,3] = 0
        Ω[3,4] = dϑ₃dx₄(t,q) - dϑ₄dx₃(t,q)

        Ω[4,1] = dϑ₄dx₁(t,q) - dϑ₁dx₄(t,q)
        Ω[4,2] = dϑ₄dx₂(t,q) - dϑ₂dx₄(t,q)
        Ω[4,3] = dϑ₄dx₃(t,q) - dϑ₃dx₄(t,q)
        Ω[4,4] = 0

        nothing
    end

    f₁(t, q, v) = dϑ₁dx₁(t,q) * v[1] + dϑ₂dx₁(t,q) * v[2] + dϑ₃dx₁(t,q) * v[3] + dϑ₄dx₁(t,q) * v[4]
    f₂(t, q, v) = dϑ₁dx₂(t,q) * v[1] + dϑ₂dx₂(t,q) * v[2] + dϑ₃dx₂(t,q) * v[3] + dϑ₄dx₂(t,q) * v[4]
    f₃(t, q, v) = dϑ₁dx₃(t,q) * v[1] + dϑ₂dx₃(t,q) * v[2] + dϑ₃dx₃(t,q) * v[3] + dϑ₄dx₃(t,q) * v[4]
    f₄(t, q, v) = dϑ₁dx₄(t,q) * v[1] + dϑ₂dx₄(t,q) * v[2] + dϑ₃dx₄(t,q) * v[3] + dϑ₄dx₄(t,q) * v[4]

    g₁(t, q, v) = dϑ₁dx₁(t,q) * v[1] + dϑ₁dx₂(t,q) * v[2] + dϑ₁dx₃(t,q) * v[3] + dϑ₁dx₄(t,q) * v[4]
    g₂(t, q, v) = dϑ₂dx₁(t,q) * v[1] + dϑ₂dx₂(t,q) * v[2] + dϑ₂dx₃(t,q) * v[3] + dϑ₂dx₄(t,q) * v[4]
    g₃(t, q, v) = dϑ₃dx₁(t,q) * v[1] + dϑ₃dx₂(t,q) * v[2] + dϑ₃dx₃(t,q) * v[3] + dϑ₃dx₄(t,q) * v[4]
    g₄(t, q, v) = dϑ₄dx₁(t,q) * v[1] + dϑ₄dx₂(t,q) * v[2] + dϑ₄dx₃(t,q) * v[3] + dϑ₄dx₄(t,q) * v[4]

    function hamiltonian(t, q, params)
        @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄ = params
        a₁*q[1] + a₂*q[2] + a₃*q[3] + a₄*q[4] + b₁*log(q[1]) + b₂*log(q[2]) + b₃*log(q[3]) + b₄*log(q[4])
    end

    hamiltonian_iode(t, q, v, params) = hamiltonian(t, q, params)

    hamiltonian_pode(t, q, p, params) = hamiltonian(t, q, params)

    function dHd₁(t, q, params)
        @unpack a₁, b₁ = params
        a₁ + b₁ / q[1]
    end

    function dHd₂(t, q, params)
        @unpack a₂, b₂ = params
        a₂ + b₂ / q[2]
    end

    function dHd₃(t, q, params)
        @unpack a₃, b₃ = params
        a₃ + b₃ / q[3]
    end

    function dHd₄(t, q, params)
        @unpack a₄, b₄ = params
        a₄ + b₄ / q[4]
    end

    function lotka_volterra_4d_dH(t, q, dH, params)
        dH[1] = dHd₁(t, q, params)
        dH[2] = dHd₂(t, q, params)
        dH[3] = dHd₃(t, q, params)
        dH[4] = dHd₄(t, q, params)
        nothing
    end

    lotka_volterra_4d_ϑ(t, q, Θ, params) = ϑ(t, q, Θ)
    lotka_volterra_4d_ϑ(t, q, v, Θ, params) = ϑ(t, q, Θ)
    lotka_volterra_4d_ω(t, q, Ω, params) = ω(t, q, Ω)

    function lotka_volterra_4d_v(t, q, v, params)
        v[1] = v₁(t, q, params)
        v[2] = v₂(t, q, params)
        v[3] = v₃(t, q, params)
        v[4] = v₄(t, q, params)
        nothing
    end

    function lotka_volterra_4d_v(t, q, p, v, params)
        lotka_volterra_4d_v(t, q, v, params)
    end

    function lotka_volterra_4d_v_ham(t, q, p, v, params)
        v .= 0
        nothing
    end

    function lotka_volterra_4d_f(t::Real, q::Vector, v::Vector, f::Vector, params)
        f[1] = f₁(t,q,v) - dHd₁(t, q, params)
        f[2] = f₂(t,q,v) - dHd₂(t, q, params)
        f[3] = f₃(t,q,v) - dHd₃(t, q, params)
        f[4] = f₄(t,q,v) - dHd₄(t, q, params)
        nothing
    end

    function lotka_volterra_4d_f_ham(t::Real, q::Vector, f::Vector, params)
        f[1] = - dHd₁(t, q, params)
        f[2] = - dHd₂(t, q, params)
        f[3] = - dHd₃(t, q, params)
        f[4] = - dHd₄(t, q, params)
        nothing
    end

    function lotka_volterra_4d_f_ham(t::Real, q::Vector, v::Vector, f::Vector, params)
        lotka_volterra_4d_f_ham(t, q, f, params)
    end

    function lotka_volterra_4d_g(t::Real, q::Vector, v::Vector, g::Vector, params)
        g[1] = f₁(t,q,v)
        g[2] = f₂(t,q,v)
        g[3] = f₃(t,q,v)
        g[4] = f₄(t,q,v)
        nothing
    end

    function lotka_volterra_4d_g(t::Real, q::Vector, p::Vector, v::Vector, g::Vector, params)
        lotka_volterra_4d_g(t, q, v, g, params)
    end

    function lotka_volterra_4d_g̅(t::Real, q::Vector, v::Vector, g::Vector, params)
        g[1] = g₁(t,q,v)
        g[2] = g₂(t,q,v)
        g[3] = g₃(t,q,v)
        g[4] = g₄(t,q,v)
        nothing
    end

    function lotka_volterra_4d_g̅(t::Real, q::Vector, p::Vector, v::Vector, g::Vector, params)
        lotka_volterra_4d_g̅(t, q, v, g, params)
    end

    function lotka_volterra_4d_u(t, q, v, u, params)
        u .= v
        nothing
    end

    function lotka_volterra_4d_u(t, q, p, v, u, params)
        lotka_volterra_4d_u(t, q, v, u, params)
    end

    function lotka_volterra_4d_ϕ(t, q, p, ϕ, params)
        ϕ[1] = p[1] - ϑ₁(t,q)
        ϕ[2] = p[2] - ϑ₂(t,q)
        ϕ[3] = p[3] - ϑ₃(t,q)
        ϕ[4] = p[4] - ϑ₄(t,q)
        nothing
    end

    function lotka_volterra_4d_ψ(t, q, p, v, f, ψ, params)
        ψ[1] = f[1] - g₁(t,q,v)
        ψ[2] = f[2] - g₂(t,q,v)
        ψ[3] = f[3] - g₃(t,q,v)
        ψ[4] = f[4] - g₄(t,q,v)
        nothing
    end

    function lotka_volterra_4d_ode(q₀=q₀, params=p)
        ODE(lotka_volterra_4d_v, q₀; parameters=params, invariants=(h=hamiltonian,))
    end

    function lotka_volterra_4d_pode(q₀=q₀, p₀=ϑ(0, q₀), params=p)
        PODE(lotka_volterra_4d_v, lotka_volterra_4d_f,
             q₀, p₀; parameters=params, invariants=(h=hamiltonian_pode,))
    end

    function lotka_volterra_4d_iode(q₀=q₀, p₀=ϑ(0, q₀), params=p)
        IODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
             lotka_volterra_4d_g, q₀, p₀;
             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v)
    end

    function lotka_volterra_4d_lode(q₀=q₀, p₀=ϑ(0, q₀), params=p)
        LODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
             lotka_volterra_4d_g, q₀, p₀;
             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v,
             Ω=lotka_volterra_4d_ω, ∇H=lotka_volterra_4d_dH)
    end

    function lotka_volterra_4d_idae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
        IDAE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
             lotka_volterra_4d_u, lotka_volterra_4d_g,
             lotka_volterra_4d_ϕ, q₀, p₀, λ₀;
             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v)
    end

    function lotka_volterra_4d_pdae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
        PDAE(lotka_volterra_4d_v_ham, lotka_volterra_4d_f_ham,
             lotka_volterra_4d_u, lotka_volterra_4d_g,
             lotka_volterra_4d_ϕ, q₀, p₀, λ₀;
             v̄=lotka_volterra_4d_v, f̄=lotka_volterra_4d_f,
             parameters=params, invariants=(h=hamiltonian_pode,))
    end

    function lotka_volterra_4d_ldae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
        LDAE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f_ham,
             lotka_volterra_4d_g, lotka_volterra_4d_g̅,
             lotka_volterra_4d_ϕ, lotka_volterra_4d_ψ,
             q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
             v̄=lotka_volterra_4d_v, f̄=lotka_volterra_4d_f,)
    end

    function lotka_volterra_4d_dg(q₀=q₀, p₀=ϑ(0, q₀), params=p)
        IODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
             lotka_volterra_4d_g, q₀, p₀;
             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v)
    end

end