text:
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using Test

@testset "Tests" begin
    @testset "ppf" begin
        @test isapprox(ttt.ppf(ttt.N01,0.3),-0.52440044)
        @test isapprox(ttt.ppf(ttt.Gaussian(2.,3.),0.3),0.42679866)
    end
    @testset "cdf" begin
        @test isapprox(ttt.cdf(ttt.N01,0.3),0.617911409)
        @test isapprox(ttt.cdf(ttt.Gaussian(2.,3.),0.3),0.28547031)
    end
    @testset "pdf" begin
        @test isapprox(ttt.pdf(ttt.N01,0.3),0.38138781)
        @test isapprox(ttt.pdf(ttt.Gaussian(2.,3.),0.3),0.11325579)
    end
    @testset "compute_margin" begin
        @test isapprox(ttt.compute_margin(0.25,sqrt(2)*25.0/6),1.8776005988)
        @test isapprox(ttt.compute_margin(0.25,sqrt(3)*25.0/6),2.29958170)
        @test isapprox(ttt.compute_margin(0.0,sqrt(3)*25.0/6),2.7134875810435737e-07)
        @test isapprox(ttt.compute_margin(1.0,sqrt(3)*25.0/6),Inf)
    end
    @testset "approx trunc" begin
        res = ttt.approx(ttt.N01,0.0,false)
        @test isapprox(res,ttt.Gaussian(0.79788453,0.602810306),1e-5)
        margin = 1.8776005988
        res = ttt.approx(ttt.Gaussian(0.,sqrt(2)*(25/6) ),margin,false)
        @test isapprox(res, ttt.Gaussian(5.958, 3.226), 1e-3)
        res = ttt.approx(ttt.Gaussian(0.,sqrt(2)*(25/6) ),margin,true)
        @test isapprox(res, ttt.Gaussian(0.,1.0767055),1e-6)
        res = ttt.approx(ttt.Gaussian(12.,sqrt(2)*(25/6) ),margin,true)
        @test isapprox(res, ttt.Gaussian(0.3900995,1.0343979),1e-6)
    res.sigma
    end
    @testset "Gaussian" begin
        N, M = ttt.Gaussian(25.0, 25.0/3), ttt.Gaussian(0.0, 1.0)
        @test isapprox(M/N, ttt.Gaussian(-0.365, 1.007),1e-3)
        @test isapprox(N*M, ttt.Gaussian(0.355, 0.993),1e-3)
        @test isapprox(N+M, ttt.Gaussian(25.00, 8.393),1e-3)
        M = ttt.Gaussian(1.0, 1.0)
        @test isapprox(N-M, ttt.Gaussian(24.00, 8.393),1e-3)
    end
    @testset "1vs1" begin
        ta = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb], [0.,1.], 0.0)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1], ttt.Gaussian(20.794779,7.194481), 1e-4)
        @test isapprox(post[2][1], ttt.Gaussian(29.205220,7.194481), 1e-4)

        g = ttt.Game([[ttt.Player(ttt.Gaussian(29.,1.),25.0/6)] ,[ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6)]], [0.,1.])
        post = ttt.posteriors(g)
        @test isapprox(post[1][1], ttt.Gaussian(28.896,0.996), 1e-3)
        @test isapprox(post[2][1], ttt.Gaussian(32.189,6.062), 1e-3)

        ta = [ttt.Player(ttt.Gaussian(1.139,0.531),1.0,0.2125)]
        tb = [ttt.Player(ttt.Gaussian(15.568,0.51),1.0,0.2125)]
        g = ttt.Game([ta,tb], [0.,1.], 0.0)
        @test isapprox(g.likelihoods[1][1].sigma, Inf)
        @test isapprox(g.likelihoods[1][1].mu, 0.0)
        @test isapprox(g.likelihoods[2][1].sigma, Inf)

    end
    @testset "1vs1vs1" begin
        teams = [[ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
                ,[ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
                ,[ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]]

        g = ttt.Game(teams, [1.,2.,0.])
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(25.000000,6.238469796),1e-5)
        @test isapprox(post[2][1],ttt.Gaussian(31.3113582213,6.69881865) ,1e-5)

        g = ttt.Game(teams)
        post = ttt.posteriors(g)
        @test isapprox(post[2][1],ttt.Gaussian(25.000000,6.238469796),1e-5)
        @test isapprox(post[1][1],ttt.Gaussian(31.3113582213,6.69881865) ,1e-5)

        g = ttt.Game(teams, [1.,2.,0.], 0.5)
        post = ttt.posteriors(g)

        @test isapprox(post[1][1],ttt.Gaussian(25.000,6.093),1e-3)
        @test isapprox(post[2][1],ttt.Gaussian(33.379,6.484),1e-3)
        @test isapprox(post[3][1],ttt.Gaussian(16.621,6.484),1e-3)
    end
    @testset "1vs1 draw" begin
        ta = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb], [0.,0.], 0.25)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(25.000,6.469481),1e-4)
        @test isapprox(post[2][1],ttt.Gaussian(25.000,6.469481),1e-4)

        ta = [ttt.Player(ttt.Gaussian(25.,3.),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(29.,2.),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb], [0.,0.], 0.25)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(25.736,2.709956),1e-4)
        @test isapprox(post[2][1],ttt.Gaussian(28.67289,1.916471),1e-4)
    end
    @testset "1vs1vs1 draw" begin
        ta = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        tc = [ttt.Player(ttt.Gaussian(25.0,25.0/3),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb,tc], [0.,0.,0.],0.25)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(25.000,5.729),1e-3)
        @test isapprox(post[2][1],ttt.Gaussian(25.000,5.707),1e-3)

        ta = [ttt.Player(ttt.Gaussian(25.,3.),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(25.,3.),25.0/6,25.0/300)]
        tc = [ttt.Player(ttt.Gaussian(29.,2.),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb,tc], [0.,0.,0.],0.25)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(25.489,2.638),1e-3)
        @test isapprox(post[2][1],ttt.Gaussian(25.511,2.629),1e-3)
        @test isapprox(post[3][1],ttt.Gaussian(28.556,1.886),1e-3)
    end
    @testset "NvsN Draw" begin
        ta = [ttt.Player(ttt.Gaussian(15.,1.),25.0/6,25.0/300)
             ,ttt.Player(ttt.Gaussian(15.,1.),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(30.,2.),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb], [0.,0.], 0.25)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(15.000,0.9916),1e-3)
        @test isapprox(post[1][2],ttt.Gaussian(15.000,0.9916),1e-3)
        @test isapprox(post[2][1],ttt.Gaussian(30.000,1.9320),1e-3)
    end
    @testset "NvsNvsN mixt" begin
        ta = [ttt.Player(ttt.Gaussian(12.,3.),25.0/6,25.0/300)
             ,ttt.Player(ttt.Gaussian(18.,3.),25.0/6,25.0/300)]
        tb = [ttt.Player(ttt.Gaussian(30.,3.),25.0/6,25.0/300)]
        tc = [ttt.Player(ttt.Gaussian(14.,3.),25.0/6,25.0/300)
             ,ttt.Player(ttt.Gaussian(16.,3.),25.0/6,25.0/300)]
        g = ttt.Game([ta,tb, tc], [1.,0.,0.], 0.25)
        post = ttt.posteriors(g)
        @test isapprox(post[1][1],ttt.Gaussian(13.051,2.864),1e-3)
        @test isapprox(post[1][2],ttt.Gaussian(19.051,2.864),1e-3)
        @test isapprox(post[2][1],ttt.Gaussian(29.292,2.764),1e-3)
        @test isapprox(post[3][1],ttt.Gaussian(13.658,2.813),1e-3)
        @test isapprox(post[3][2],ttt.Gaussian(15.658,2.813),1e-3)
    end
    @testset "Game evidence" begin
        @testset "1vs1" begin
            ta = [ttt.Player(ttt.Gaussian(27.,2.),1.0)]
            tb = [ttt.Player(ttt.Gaussian(25.,3.0),1.0)]
            d = ta[1].prior.mu - tb[1].prior.mu
            v = sqrt(2+ta[1].prior.sigma^2 + tb[1].prior.sigma^2)
            evidence = 1-ttt.cdf(ttt.Gaussian(d,v),0.0)
            g = ttt.Game([ta,tb], [1.,0.])
            @test isapprox(g.evidence , evidence)

            ta = [ttt.Player(ttt.Gaussian(27.,2.),1.0),ttt.Player(ttt.Gaussian(28.,1.),1.0) ]
            tb = [ttt.Player(ttt.Gaussian(25.,3.0),1.0),ttt.Player(ttt.Gaussian(26.,1.0),1.0)]
            d = (ta[1].prior.mu + ta[2].prior.mu) - (tb[1].prior.mu + tb[2].prior.mu)
            v = sqrt(4+ta[1].prior.sigma^2 + ta[2].prior.sigma^2+ tb[1].prior.sigma^2 + tb[2].prior.sigma^2 )
            evidence = 1-ttt.cdf(ttt.Gaussian(d,v),0.0)
            g = ttt.Game([ta,tb], [1.,0.])
            @test isapprox(g.evidence , evidence)

            ta = [ttt.Player(ttt.Gaussian(25.,1e-7),25.0/6,25.0/300)]
            tb = [ttt.Player(ttt.Gaussian(25.,1e-7),25.0/6,25.0/300)]
            g = ttt.Game([ta,tb], [0.,0.], 0.25)
            @test isapprox(g.evidence,0.25)
            g = ttt.Game([ta,tb], [1.,0.], 0.25)
            @test isapprox(g.evidence,0.375)
        end
        @testset "1vs1vs1 margin 0" begin
            ta = [ttt.Player(ttt.Gaussian(25.,1e-7),25.0/6,25.0/300)]
            tb = [ttt.Player(ttt.Gaussian(25.,1e-7),25.0/6,25.0/300)]
            tc = [ttt.Player(ttt.Gaussian(25.,1e-7),25.0/6,25.0/300)]

            g_abc = ttt.Game([ta,tb,tc], [3.,2.,1.], 0.)
            g_acb = ttt.Game([ta,tb,tc], [3.,1.,2.], 0.)
            g_bac = ttt.Game([ta,tb,tc], [2.,3.,1.], 0.)
            g_bca = ttt.Game([ta,tb,tc], [1.,3.,2.], 0.)
            g_cab = ttt.Game([ta,tb,tc], [2.,1.,3.], 0.)
            g_cba = ttt.Game([ta,tb,tc], [1.,2.,3.], 0.)

            d1 = ttt.performance(g_abc,1)-ttt.performance(g_abc,2)

            proba = 0
            proba += g_abc.evidence
            proba += g_acb.evidence
            proba += g_bac.evidence
            proba += g_bca.evidence
            proba += g_cab.evidence
            proba += g_cba.evidence
            println("Corregir la evidencia multiequipos para que sume 1")
            @test  isapprox(proba, 1.49999991)
        end
    end
    @testset "Forget" begin
        gamma = 0.15*25.0/3
        g = ttt.Gaussian(25.,1e-7)
        @test isapprox(ttt.forget(g,gamma,5).sigma,sqrt(5*gamma^2))
        @test isapprox(ttt.forget(g,gamma,1).sigma,sqrt(1*gamma^2))
    end
    @testset "One event each" begin
        agents = Dict{String,ttt.Agent}()
        for k in ["a", "b", "c", "d", "e", "f"]
            agents[k] = ttt.Agent(ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 25.0/300 ) , ttt.Ninf, ttt.minInt64)
        end
        composition = [ [["a"],["b"]], [["c"],["d"]] , [["e"],["f"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        b = ttt.Batch(composition = composition, results = results, agents = agents)
        @test isapprox(ttt.posterior(b,"a"),ttt.Gaussian(29.205,7.194),1e-3)
        @test isapprox(ttt.posterior(b,"b"),ttt.Gaussian(20.795,7.194),1e-3)
        @test isapprox(ttt.posterior(b,"d"),ttt.Gaussian(29.205,7.194),1e-3)
        @test isapprox(ttt.posterior(b,"c"),ttt.Gaussian(20.795,7.194),1e-3)
        @test isapprox(ttt.posterior(b,"e"),ttt.Gaussian(29.205,7.194),1e-3)
        @test isapprox(ttt.posterior(b,"f"),ttt.Gaussian(20.795,7.194),1e-3)
        iter = ttt.convergence(b)
        @test iter == 1
    end
    @testset "Same strength" begin
        agents= Dict{String,ttt.Agent}()
        for k in ["a", "b", "c", "d", "e", "f"]
            agents[k] = ttt.Agent(ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 25.0/300 ) , ttt.Ninf, ttt.minInt64)
        end
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        b = ttt.Batch(composition = composition, results = results, agents = agents)
        @test isapprox(ttt.posterior(b,"a"),ttt.Gaussian(24.96097,6.29954),1e-3)
        @test isapprox(ttt.posterior(b,"b"),ttt.Gaussian(27.09559,6.01033),1e-3)
        @test isapprox(ttt.posterior(b,"c"),ttt.Gaussian(24.88968,5.86631),1e-3)
        iter = ttt.convergence(b)
        @test isapprox(ttt.posterior(b,"a"),ttt.Gaussian(25.000,5.419),1e-3)
        @test isapprox(ttt.posterior(b,"b"),ttt.Gaussian(25.000,5.419),1e-3)
        @test isapprox(ttt.posterior(b,"c"),ttt.Gaussian(25.000,5.419),1e-3)
    end
    @testset "Add events into a batch" begin
        agents= Dict{String,ttt.Agent}()
        for k in ["a", "b", "c", "d", "e", "f"]
            agents[k] = ttt.Agent(ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 25.0/300 ) , ttt.Ninf, ttt.minInt64)
        end
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        b = ttt.Batch(composition = composition, results = results, agents = agents)
        iter = ttt.convergence(b)
        @test isapprox(ttt.posterior(b,"a"),ttt.Gaussian(25.000,5.419),1e-3)
        @test isapprox(ttt.posterior(b,"b"),ttt.Gaussian(25.000,5.419),1e-3)
        @test isapprox(ttt.posterior(b,"c"),ttt.Gaussian(25.000,5.419),1e-3)
        ttt.add_events(b,composition,results)
        @test length(b.events) == 6
        iter = ttt.convergence(b,1e-6,20)
        @test iter == 20
        @test isapprox(ttt.posterior(b,"a"),ttt.Gaussian(25.000,3.88),1e-3)
        @test isapprox(ttt.posterior(b,"b"),ttt.Gaussian(25.000,3.88),1e-3)
        @test isapprox(ttt.posterior(b,"c"),ttt.Gaussian(25.000,3.88),1e-3)
        end
    @testset "TrueSkill initalization" begin
        composition = [ [["aa"],["b"]], [["aa"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        priors = Dict{String,ttt.Player}()
        for k in ["aa", "b", "c"]
            priors[k] = ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 0.15*25.0/3)
        end

        h = ttt.History(composition, results, [1,2,3], priors)

        @test isapprox(ttt.posterior(h.batches[1],"aa"),ttt.Gaussian(29.205,7.19448),1e-3)

        observed = h.batches[2].skills["aa"].forward.sigma
        gamma = 0.15*25.0/3
        expected = sqrt((gamma*1)^2 + ttt.posterior(h.batches[1],"aa").sigma^2)
        @test isapprox(observed, expected)

        observed = ttt.posterior(h.batches[2],"aa")

        g = ttt.Game(ttt.within_priors(h.batches[2], 1),[0.,1.])
        expected = ttt.posteriors(g)[1][1]
        @test isapprox(observed, expected, 1e-7)
    end
    @testset "One-batch history" begin
        composition = [ [["aj"],["bj"]],[["bj"],["cj"]], [["cj"],["aj"]] ]
        results = [[1.,0.],[1.,0.],[1.,0.]]
        times = [1,1,1]
        priors = Dict{String,ttt.Player}()
        for k in ["aj", "bj", "cj"]
            priors[k] = ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 0.15*25.0/3)
        end

        h1 = ttt.History(composition,results, times, priors)
        # TrueSkill
        @test isapprox(ttt.posterior(h1.batches[1],"aj"),ttt.Gaussian(22.904,6.010),1e-3)
        @test isapprox(ttt.posterior(h1.batches[1],"cj"),ttt.Gaussian(25.110,5.866),1e-3)
        # TTT
        step , i = ttt.convergence(h1)
        @test isapprox(ttt.posterior(h1.batches[1],"aj"),ttt.Gaussian(25.000,5.419),1e-3)
        @test isapprox(ttt.posterior(h1.batches[1],"cj"),ttt.Gaussian(25.000,5.419),1e-3)

        priors = Dict{String,ttt.Player}()
        for k in ["aj", "bj", "cj"]
            priors[k] = ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 25.0/300)
        end
        h2 = ttt.History(composition,results, [1,2,3], priors  )
        # TrueSkill
        @test isapprox(ttt.posterior(h2.batches[3],"aj"),ttt.Gaussian(22.904,6.011),1e-3)
        @test isapprox(ttt.posterior(h2.batches[3],"cj"),ttt.Gaussian(25.111,5.867),1e-3)
        # TTT
        step , i = ttt.convergence(h2)
        @test isapprox(ttt.posterior(h2.batches[3],"aj"),ttt.Gaussian(24.999,5.420),1e-3)
        @test isapprox(ttt.posterior(h2.batches[3],"cj"),ttt.Gaussian(25.001,5.420),1e-3)
    end
    @testset "TrueSkill Through Time" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        priors = Dict{String,ttt.Player}()
        for k in ["a", "b", "c"]
            priors[k] = ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 25.0/300)
        end

        h = ttt.History(composition, results, Int64[], priors)
        step , iter = ttt.convergence(h)
        @test (h.batches[3].skills["b"].elapsed == 1) & (h.batches[3].skills["c"].elapsed == 1)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian(25.0002673,5.41938162),1e-5)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(24.999465,5.419425831),1e-5)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian(25.00053219,5.419696790),1e-5)
        #ttt.learning_curves(h)
    end
    @testset "Env TTT" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]

        h = ttt.History(composition=composition, results=results, mu=25.,sigma=25.0/3, beta=25.0/6, gamma=25.0/300)
        step , iter = ttt.convergence(h)
        @test (h.batches[3].skills["b"].elapsed == 1) & (h.batches[3].skills["c"].elapsed == 1)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian(25.0002673,5.41938162),1e-5)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(24.999465,5.419425831),1e-5)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian(25.00053219,5.419696790),1e-5)

    end
    @testset "Env 0 TTT" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]

        @time h = ttt.History(composition, results, mu=0.0,sigma=6.0, beta=1.0, gamma=0.05)
        @time step , iter = ttt.convergence(h, iterations=100)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian( 0.001,2.395),1e-3)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(-0.001,2.396),1e-3)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian( 0.001,2.396),1e-3)

        composition = [ [["a"],["b"]], [["c"],["a"]] , [["b"],["c"]] ]
        @time h = ttt.History(composition=composition, mu=0.0,sigma=6.0, beta=1.0, gamma=0.05)
        @time step , iter = ttt.convergence(h, iterations=100)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian( 0.001,2.395),1e-3)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(-0.001,2.396),1e-3)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian( 0.001,2.396),1e-3)

        @time h = ttt.History(composition, mu=0.0,sigma=6.0, beta=1.0, gamma=0.05)
        @time step , iter = ttt.convergence(h, iterations=100)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian( 0.001,2.395),1e-3)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(-0.001,2.396),1e-3)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian( 0.001,2.396),1e-3)
    end
    @testset "Env 0 TTT" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]

        @time h = ttt.History(composition=composition, results=results, times = [0, 10, 20], mu=0.0,sigma=6.0, beta=1.0, gamma=0.05)
        @time step , iter = ttt.convergence(h, iterations=100)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian( 0.005,2.404),1e-3)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(-0.016,2.404),1e-3)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian( 0.017,2.406),1e-3)
    end
    @testset "Teams" begin
        composition = [ [["a","b"],["c","d"]], [["e","f"] , ["b","c"]], [["a","d"], ["e","f"]]  ]
        results = [[1.,0.],[0.,1.],[1.,0.]]

        h = ttt.History(composition=composition, results=results, mu=0.0,sigma=6.0, beta=1.0, gamma=0.0)
        trueskill_log_evidence = ttt.log_evidence(h)
        trueskill_log_evidence_online =  ttt.log_evidence(h, online=true)
        @test isapprox(trueskill_log_evidence , trueskill_log_evidence_online ,rtol=1e-5)

        @test isapprox(ttt.posterior(h.batches[1],"b").mu, -1*ttt.posterior(h.batches[1],"c").mu,rtol=1e-5)

        evidence_second_event = exp(ttt.log_evidence(h, agents = ["b"]))*2
        @test isapprox(0.5,evidence_second_event,rtol=1e-5)
        evidence_third_event = exp(ttt.log_evidence(h, agents = ["a"]))*2
        @test isapprox(0.669885,evidence_third_event ,rtol=1e-5)

        step , iter = ttt.convergence(h)

        loocv_hat = exp(ttt.log_evidence(h))
        p_d_m_hat = exp(ttt.log_evidence(h, online=true))

        @test isapprox( ttt.posterior(h.batches[1],"a"), ttt.posterior(h.batches[1],"b"), 1e-3)
        @test isapprox( ttt.posterior(h.batches[1],"c"), ttt.posterior(h.batches[1],"d"), 1e-3)
        @test isapprox( ttt.posterior(h.batches[2],"e"), ttt.posterior(h.batches[2],"f"), 1e-3)

        @test isapprox( ttt.posterior(h.batches[1],"a"), ttt.Gaussian(mu=4.085 ,sigma=5.107), 1e-3 )
        @test isapprox( ttt.posterior(h.batches[1],"c"), ttt.Gaussian(mu=-0.533 ,sigma=5.107), 1e-3)
        @test isapprox( ttt.posterior(h.batches[3],"e"), ttt.Gaussian(mu=-3.552 ,sigma=5.155), 1e-3 )
    end
    @testset "sigma-beta 0" begin
        composition = [ [["a","a_b","b"],["c","c_d","d"]]
                 , [["e","e_f","f"],["b","b_c","c"]]
                 , [["a","a_d","d"],["e","e_f","f"]]  ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        priors = Dict{String,ttt.Player}()
        for k in ["a_b", "c_d", "e_f", "b_c", "a_d", "e_f"]
            priors[k] = ttt.Player(prior=ttt.Gaussian(mu=0.0, sigma=1e-7), beta=0.0, gamma=0.2)
        end

        @time h = ttt.History(composition=composition, results=results, priors=priors, mu=0.0,sigma=6.0, beta=1.0, gamma=0.0)
        @time step , iter = ttt.convergence(h)
        @test isapprox( ttt.posterior(h.batches[1],"a_b"), ttt.Gaussian(mu=0.0 ,sigma=0.0), 1e-4 )
        @test isapprox( ttt.posterior(h.batches[3],"e_f"), ttt.Gaussian(mu=-0.002 ,sigma=0.2), 1e-4)
    end
    @testset "Memory Size" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]

        h = ttt.History(composition=composition, results=results, times = [0, 10, 20], mu=0.0,sigma=6.0, beta=1.0, gamma=0.05)
        @test Base.summarysize(h) < 3700
        @test Base.summarysize(h.batches) - Base.summarysize(h.agents) < 2900
        @test Base.summarysize(h.agents) < 700
        @test Base.summarysize(h.batches[2]) - Base.summarysize(h.agents)  < 946

        @test Base.summarysize(h) < Base.summarysize(composition)  * 6.5

        @test Base.summarysize(h.batches[1].skills) == 618
        @test Base.summarysize(h.batches[1].events) == 314
    end
    @testset "Learning curves" begin
        composition = [ [["aj"],["bj"]],[["bj"],["cj"]], [["cj"],["aj"]] ]
        results = [[1.,0.],[1.,0.],[1.,0.]]
        priors = Dict{String,ttt.Player}()
        for k in ["aj", "bj", "cj"]
            priors[k] = ttt.Player(ttt.Gaussian(25., 25.0/3), 25.0/6, 25.0/300)
        end
        h = ttt.History(composition, results, [5,6,7], priors)
        ttt.convergence(h)
        lc = ttt.learning_curves(h)

        @test lc["aj"][1][1] == 5
        @test lc["aj"][end][1] == 7
        @test isapprox(lc["aj"][end][2],ttt.Gaussian(24.999,5.420),1e-3)
        @test isapprox(lc["cj"][end][2],ttt.Gaussian(25.001,5.420),13-3)
    end
    @testset "Add events into History 1" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]

        h = ttt.History(composition=composition, results=results, mu=0.0, sigma=2.0, beta=1.0, gamma=0.0)
        step , iter = ttt.convergence(h)
        @test (h.batches[3].skills["b"].elapsed == 1) & (h.batches[3].skills["c"].elapsed == 1)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian(0.0,1.301),1e-3)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(0.0,1.301),1e-3)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian(0.0,1.301),1e-3)
        ttt.add_events(h, composition, results)

        @test length(h.batches) == 6
        @test [ttt.get_composition(b.events) for b in h.batches] == [ [[["a"], ["b"]]], [[["a"], ["c"]]], [[["b"], ["c"]]]
        , [[["a"], ["b"]]], [[["a"], ["c"]]], [[["b"], ["c"]]] ]

        step , iter = ttt.convergence(h)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian(0.0,0.931),1e-3)
        @test isapprox(ttt.posterior(h.batches[4],"a"),ttt.Gaussian(0.0,0.931),1e-3)
        @test isapprox(ttt.posterior(h.batches[4],"b"),ttt.Gaussian(0.0,0.931),1e-3)
        @test isapprox(ttt.posterior(h.batches[6],"b"),ttt.Gaussian(0.0,0.931),1e-3)

        composition = [ [["a"],["b"]], [["c"],["a"]] , [["b"],["c"]] ]

        h = ttt.History(composition=composition, mu=0.0, sigma=2.0, beta=1.0, gamma=0.0)
        step , iter = ttt.convergence(h)
        @test (h.batches[3].skills["b"].elapsed == 1) & (h.batches[3].skills["c"].elapsed == 1)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian(0.0,1.301),1e-3)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.Gaussian(0.0,1.301),1e-3)
        @test isapprox(ttt.posterior(h.batches[3],"b"),ttt.Gaussian(0.0,1.301),1e-3)
        ttt.add_events(h, composition)

        @test length(h.batches) == 6
        @test [ttt.get_composition(b.events) for b in h.batches] == [ [[["a"], ["b"]]], [[["c"], ["a"]]], [[["b"], ["c"]]]
        , [[["a"], ["b"]]], [[["c"], ["a"]]], [[["b"], ["c"]]] ]

        step , iter = ttt.convergence(h)
        @test isapprox(ttt.posterior(h.batches[1],"a"),ttt.Gaussian(0.0,0.931),1e-3)
        @test isapprox(ttt.posterior(h.batches[4],"a"),ttt.Gaussian(0.0,0.931),1e-3)
        @test isapprox(ttt.posterior(h.batches[4],"b"),ttt.Gaussian(0.0,0.931),1e-3)
        @test isapprox(ttt.posterior(h.batches[6],"b"),ttt.Gaussian(0.0,0.931),1e-3)
    end
    @testset "Add events into History 2" begin
        composition = [ [["a"],["b"]], [["a"],["c"]] , [["b"],["c"]] ]
        results = [[1.,0.],[0.,1.],[1.,0.]]
        times = [0,10,20]
        h = ttt.History(composition=composition, results=results, times=times, mu=0.0, sigma=2.0, beta=1.0, gamma=0.0)
        step , iter = ttt.convergence(h)
        times = [15,10,0]
        ttt.add_events(h, composition, results, times)
        @test length(h.batches) == 4
        @test [length(b) for b in h.batches] == [2,2,1,1]
        @test [ttt.get_composition(b.events) for b in h.batches] == [
        [[["a"], ["b"]], [["b"], ["c"]]] ,
        [[["a"], ["c"]], [["a"], ["c"]]] ,
        [[["a"], ["b"]]] ,
        [[["b"], ["c"]]]
        ]
        @test [ttt.get_results(b.events) for b in h.batches] == [ [[1., 0.], [1., 0.]], [[0., 1.], [0., 1.]], [[1., 0.]], [[1., 0.]] ]
        @test (h.batches[1].skills["c"].elapsed == 0) & (h.batches[end].skills["c"].elapsed == 10)
        @test (h.batches[1].skills["a"].elapsed == 0) & (h.batches[3].skills["a"].elapsed == 5)
        @test (h.batches[1].skills["b"].elapsed == 0) & (h.batches[end].skills["b"].elapsed == 5)

        step , iter = ttt.convergence(h)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.posterior(h.batches[end],"b"),1e-4)
        @test isapprox(ttt.posterior(h.batches[1],"c"),ttt.posterior(h.batches[end],"c"),1e-4)
        @test isapprox(ttt.posterior(h.batches[1],"c"),ttt.posterior(h.batches[1],"b"),1e-4)

        # # # # # # # # #

        composition = [ [["a"],["b"]], [["c"],["a"]] , [["b"],["c"]] ]
        times = [0,10,20]
        h = ttt.History(composition=composition, times=times, mu=0.0, sigma=2.0, beta=1.0, gamma=0.0)
        step , iter = ttt.convergence(h)
        times = [15,10,0]
        ttt.add_events(h, composition, times=times)
        @test length(h.batches) == 4
        @test [length(b) for b in h.batches] == [2,2,1,1]
        @test [ttt.get_composition(b.events) for b in h.batches] == [
        [[["a"], ["b"]], [["b"], ["c"]]] ,
        [[["c"], ["a"]], [["c"], ["a"]]] ,
        [[["a"], ["b"]]] ,
        [[["b"], ["c"]]]
        ]
        @test [ttt.get_results(b.events) for b in h.batches] == [ [[1., 0.], [1., 0.]], [[1., 0.], [1., 0.]], [[1., 0.]], [[1., 0.]] ]
        @test (h.batches[1].skills["c"].elapsed == 0) & (h.batches[end].skills["c"].elapsed == 10)
        @test (h.batches[1].skills["a"].elapsed == 0) & (h.batches[3].skills["a"].elapsed == 5)
        @test (h.batches[1].skills["b"].elapsed == 0) & (h.batches[end].skills["b"].elapsed == 5)

        step , iter = ttt.convergence(h)
        @test isapprox(ttt.posterior(h.batches[1],"b"),ttt.posterior(h.batches[end],"b"),1e-4)
        @test isapprox(ttt.posterior(h.batches[1],"c"),ttt.posterior(h.batches[end],"c"),1e-4)
        @test isapprox(ttt.posterior(h.batches[1],"c"),ttt.posterior(h.batches[1],"b"),1e-4)

    end
    @testset "Add events into History 2" begin
        composition = [[["a"],["b"]],[["b"],["a"]]]
        h = ttt.History(composition)
        p_d_m_2 = exp(ttt.log_evidence(h))*2
        @test isapprox(p_d_m_2, 0.17650911, rtol=1e-5)
        @test isapprox(p_d_m_2, exp(ttt.log_evidence(h, online=true))*2, rtol=1e-5)
        @test isapprox(p_d_m_2, exp(ttt.log_evidence(h, online=true, agents = ["a"]))*2, rtol=1e-5)
        @test isapprox(p_d_m_2, exp(ttt.log_evidence(h, online=false, agents = ["a"]))*2, rtol=1e-5)

        ttt.convergence(h,iterations=11)
        loocv_approx_2 = sqrt(exp(ttt.log_evidence(h)))
        @test isapprox(loocv_approx_2, 0.001976774, rtol=1e-5)

        p_d_m_approx_2 = exp(ttt.log_evidence(h, online=true))*2
        @test loocv_approx_2-p_d_m_approx_2 < 1e-4
        @test isapprox(loocv_approx_2, exp(ttt.log_evidence(h, online=true, agents= ["b"]))*2, atol=1e-5)

    end
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
#----------------------------------------------------------------------------------*
#*
#*-------------------------  Fundamental Variables ----------------------------------*
#*
#*----------------------------------------------------------------------------------
# Constants
# Pressure
# Temperature
# Time
# Size related
# Eletric
export potency
typealias Danapotency DanaRealParametric
type _potency
	function _potency(_::Dict{Symbol,Any})
		fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
		fields[:Brief]="Potency"
		fields[:Default]=10
		fields[:Lower]=-1e3
		fields[:Upper]=1e3
		fields[:finalUnit]="kW"
		drive!(fields,_)
		new(fields)
	end
	value::Dict{Symbol,Any}
end
typealias potency DanaRealParametric{_potency}

ext: jl


text:
struct SinkhornDivergence{A<:Sinkhorn,B<:Sinkhorn,C<:Sinkhorn}
    algμν::A
    algμμ::B
    algνν::C
end

"""
    sinkhorn_divergence(
        μ::AbstractVecOrMat,
        ν::AbstractVecOrMat,
        C,
        ε,
        alg::SinkhornDivergence=SinkhornDivergence(
            SinkhornGibbs(), SymmetricSinkhornGibbs(), SymmetricSinkhornGibbs()
        );
        kwargs...,
    )

Compute the Sinkhorn Divergence between finite discrete
measures `μ` and `ν` with respect to a common cost matrix `C`,
entropic regularization parameter `ε` and algorithm `alg`.

In the default case where `regularization = false`, the Sinkhorn Divergence is that of [^GPC18] and is computed as
```math
\\operatorname{S}_{ε}(μ,ν) := \\operatorname{W}_{ε}(μ,ν)
- \\frac{1}{2}(\\operatorname{W}_{ε}(μ,μ) + \\operatorname{W}_{ε}(ν,ν)),
```
and ``\\operatorname{W}_{ε}`` is defined as
```math
\\operatorname{W}_{ε}(μ, ν) = \\langle C, γ^\\star \\rangle,
```
where ``γ^\\star`` is the entropy-regularised transport plan between `μ` and `ν`.
For `regularization = true`, the Sinkhorn Divergence is that of [^FeydyP19] and is computed as above
where ``\\operatorname{W}_{ε}`` is replaced by ``\\operatorname{OT}_{ε}``, the entropy-regularised optimal transport
cost with regulariser penalty.

The default algorithm for computing the term ``\\operatorname{W}_{ε}(μ, ν)`` is the `SinkhornGibbs` algorithm.
For the terms ``\\operatorname{W}_{ε}(μ, μ)`` and ``\\operatorname{W}_{ε}(ν, ν)``, the symmetric fixed point iteration of [^FeydyP19] is used.
Alternatively, a pre-computed optimal transport `plan` between `μ` and `ν` may be provided.

[^GPC18]: Aude Genevay, Gabriel Peyré, Marco Cuturi, Learning Generative Models with Sinkhorn Divergences, Proceedings of the Twenty-First International Conference on Artficial Intelligence and Statistics, (AISTATS) 21, 2018
[^FeydyP19]: Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2681–2690. PMLR, 2019.
See also: [`sinkhorn2`](@ref)
"""
function sinkhorn_divergence(
    μ::AbstractVecOrMat,
    ν::AbstractVecOrMat,
    C,
    ε,
    alg::SinkhornDivergence=SinkhornDivergence(
        SinkhornGibbs(), SymmetricSinkhornGibbs(), SymmetricSinkhornGibbs()
    );
    kwargs...,
)
    return sinkhorn_divergence(μ, ν, C, C, C, ε, alg; kwargs...)
end

"""
    sinkhorn_divergence(
        μ,
        ν,
        Cμν,
        Cμ,
        Cν,
        ε,
        alg::SinkhornDivergence=SinkhornDivergence(
            SinkhornGibbs(), SymmetricSinkhornGibbs(), SymmetricSinkhornGibbs()
        );
        kwargs...,
    )

Compute the Sinkhorn Divergence between finite discrete
measures `μ` and `ν` with respect to the precomputed cost matrices `Cμν`,
`Cμμ` and `Cνν`, entropic regularization parameter `ε` and algorithm `alg`.

A pre-computed optimal transport `plan` between `μ` and `ν` may be provided.

See also: [`sinkhorn2`](@ref), [`sinkhorn_divergence`](@ref)
"""
function sinkhorn_divergence(
    μ,
    ν,
    Cμν,
    Cμ,
    Cν,
    ε,
    alg::SinkhornDivergence=SinkhornDivergence(
        SinkhornGibbs(), SymmetricSinkhornGibbs(), SymmetricSinkhornGibbs()
    );
    kwargs...,
)
    OTμν = sinkhorn2(μ, ν, Cμν, ε, alg.algμν; kwargs...)
    OTμμ = sinkhorn2(μ, Cμ, ε, alg.algμμ; kwargs...)
    OTνν = sinkhorn2(ν, Cν, ε, alg.algνν; kwargs...)
    return max.(0, OTμν .- (OTμμ .+ OTνν) / 2)
end


text:
module IPCTests

using Test

using IPC

@testset "Basic Functions       " begin
    pid = IPC.getpid()
    @test pid.value == getpid()
    @test isa(string(pid), String)
    @test isa(string(IPC.getppid()), String)
    @test isa(string(IPC.getuid()), String)
    @test isa(string(IPC.geteuid()), String)
end

@testset "File Functions        " begin
    path = "/tmp/test-$(getpid())"
    data = rand(10)
    open(path, "w") do io
        write(io, data)
    end
    f = open(IPC.FileDescriptor, path, "r")
    @test fd(f) ≥ 0
    @test filesize(f) == sizeof(data)
    @test position(f) == 0
    @test seekend(f) == position(f) == sizeof(data)
    @test seekstart(f) == position(f) == 0
    pos = (sizeof(data)>>1)
    @test seek(f, pos) == position(f) == pos
    close(f)
    @test fd(f) == -1
end

@testset "Time Functions        " begin
    # compile and warmup
    time()
    float(gettimeofday())
    float(clock_gettime(CLOCK_MONOTONIC))
    float(clock_gettime(CLOCK_REALTIME))
    float(clock_getres(CLOCK_MONOTONIC))
    float(clock_getres(CLOCK_REALTIME))

    # compare times
    ms = 0.001
    @test abs(time() - float(gettimeofday())) ≤ 1ms
    @test abs(time() - float(clock_gettime(CLOCK_REALTIME))) ≤ 1ms
    @test float(clock_getres(CLOCK_MONOTONIC)) ≤ 1ms
    @test float(clock_getres(CLOCK_REALTIME)) ≤ 1ms
    @test float(nanosleep(0.01)) == 0
    s = 10ms
    t0 = float(clock_gettime(CLOCK_MONOTONIC))
    nanosleep(s)
    t1 = float(clock_gettime(CLOCK_MONOTONIC))
    @test abs(t1 - (t0 + s)) ≤ 4ms
end

@testset "Signals               " begin
    # Sanity check:
    @test isbitstype(IPC._typeof_sigset)
    @test sizeof(IPC._typeof_sigset) == IPC._sizeof_sigset
    @test isbitstype(IPC._typeof_siginfo)
    @test sizeof(IPC._typeof_siginfo) == IPC._sizeof_siginfo

    # SigSet:
    set = SigSet()
    sigrtmin = (isdefined(IPC, :SIGRTMIN) ? IPC.SIGRTMIN : 1)
    sigrtmax = (isdefined(IPC, :SIGRTMAX) ? IPC.SIGRTMAX : 8*sizeof(set))
    @test all(set[i] == false for i in sigrtmin:sigrtmax)
    fill!(set, true)
    @test all(set[i] for i in sigrtmin:sigrtmax)
    i = rand(sigrtmin:sigrtmax)
    set[i] = false
    @test !set[i]
    set[i] = true
    @test set[i]
    set = sigpending()
    @test all(set[i] == false for i in 1:sigrtmax)

    # SigInfo:
    si = SigInfo()

    # SigAction:
    sa = SigAction()
end

@testset "BSD System V Keys     " begin
    path = "/tmp"
    key1 = IPC.Key(path, '1')
    key2 = IPC.Key(key1.value)
    @test key1 == key2
    shmds = ShmInfo()
    semds = IPC.SemInfo()
end

# Skip semaphore tests for Apple.
@static if !Sys.isapple()

    @testset "Named Semaphores      " begin
        begin
            name = "/sem-$(getpid())"
            rm(Semaphore, name)
            @test_throws SystemError Semaphore(name)
            sem1 = Semaphore(name, 0)
            sem2 = Semaphore(name)
            @test sem1[] == 0
            @test sem2[] == 0
            post(sem1)
            @test sem2[] == 1
            post(sem1)
            @test sem2[] == 2
            wait(sem2)
            @test sem2[] == 1
            @test trywait(sem2) == true
            @test trywait(sem2) == false
            @test_throws TimeoutError timedwait(sem2, 0.1)
        end
        GC.gc() # call garbage collector to exercise the finalizers
    end

    @testset "Anonymous Semaphores  " begin
        begin
            buf = DynamicMemory(sizeof(Semaphore))
            sem1 = Semaphore(buf, 0)
            sem2 = Semaphore(buf)
            @test sem1[] == 0
            @test sem2[] == 0
            post(sem1)
            @test sem2[] == 1
            post(sem1)
            @test sem2[] == 2
            wait(sem2)
            @test sem2[] == 1
            @test trywait(sem2) == true
            @test trywait(sem2) == false
            @test_throws TimeoutError timedwait(sem2, 0.1)
        end
        GC.gc() # call garbage collector to exercise the finalizers
    end

end

@testset "Wrapped Arrays        " begin
    begin
        T = Float32
        dims = (5,6)
        buf = DynamicMemory(sizeof(T)*prod(dims))
        A = WrappedArray(buf, T, dims)
        B = WrappedArray(buf) # indexable byte buffer
        C = WrappedArray(buf, T) # indexable byte buffer
        D = WrappedArray(buf, b -> (T, dims, 0)) # all parameters provided by a function
        A[:] = 1:length(A) # fill A before the copy
        E = copy(A)
        n = prod(dims)
        @test ndims(A) == ndims(D) == ndims(E) == length(dims)
        @test size(A) == size(D) == size(E) == dims
        @test all(size(A,i) == size(D,i) == dims[i] for i in 1:length(dims))
        @test eltype(A) == eltype(C) == eltype(D) == T
        @test Base.elsize(A) == Base.sizeof(T)
        @test length(A) == length(C) == div(length(B), sizeof(T)) == length(D) == n
        @test sizeof(A) == sizeof(B) == sizeof(C) == sizeof(D) == sizeof(buf)
        @test pointer(A) == pointer(B) == pointer(C) == pointer(D) == pointer(buf)
        @test isa(A.arr, Array{T,length(dims)})
        @test A[1] == 1 && A[end] == prod(dims)
        @test all(A[i] == i for i in 1:n)
        @test all(C[i] == i for i in 1:n)
        @test all(D[i] == i for i in 1:n)
        @test all(E[i] == i for i in 1:n)
        B[:] .= 0
        @test all(A[i] == 0 for i in 1:n)
        C[:] = randn(T, n)
        flag = true
        for i in eachindex(A, D)
            if A[i] != D[i]
                flag = false
            end
        end
        @test flag
        B[:] .= 0
        A[2,3] = 23
        A[3,2] = 32
        @test D[2,3] == 23 && D[3,2] == 32
        # Test copy back.
        copyto!(A, E)
        @test all(A[i] == E[i] for i in 1:n)
    end
    GC.gc() # call garbage collector to exercise the finalizers
end

@testset "Shared Memory (BSD)   " begin
    begin
        T = Int
        n = 10
        len = n*sizeof(Int)
        A = SharedMemory(IPC.PRIVATE, len)
        id = shmid(A)
        @test isa(id, ShmId)
        @test id == ShmId(A)
        @test "$id" == "ShmId($(id.value))"
        @test id.value == convert(Int64, id)
        @test id.value == convert(Int32, id)
        B = SharedMemory(id; readonly=false)
        C = SharedMemory(id; readonly=true)
        @test shmid(A) == shmid(B) == shmid(C)
        @test sizeof(A) == sizeof(B) == sizeof(C) == len
        Aptr = convert(Ptr{T},pointer(A))
        Bptr = convert(Ptr{T},pointer(B))
        Cptr = convert(Ptr{T},pointer(C))
        for i in 1:n
            unsafe_store!(Aptr,i,i)
        end
        @test all(unsafe_load(Bptr, i) == i for i in 1:n)
        @test all(unsafe_load(Cptr, i) == i for i in 1:n)
        for i in 1:n
            unsafe_store!(Bptr,-i,i)
        end
        @test all(unsafe_load(Aptr, i) == -i for i in 1:n)
        @test all(unsafe_load(Cptr, i) == -i for i in 1:n)
        @test_throws ReadOnlyMemoryError unsafe_store!(Cptr,42)
        @test ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t),
                    Aptr, 0, len) == Aptr
        @test all(unsafe_load(Cptr, i) == 0 for i in 1:n)
    end
    GC.gc() # call garbage collector to exercise the finalizers
end

@testset "Shared Memory (POSIX) " begin
    begin
        T = Int
        n = 10
        len = n*sizeof(Int)
        name = "/shm-$(getpid())"
        rm(SharedMemory, name)
        @test_throws SystemError SharedMemory(name)
        A = SharedMemory(name, len)
        id = shmid(A)
        @test isa(id, String)
        @test id == name
        B = SharedMemory(id; readonly=false)
        C = SharedMemory(id; readonly=true)
        @test shmid(A) == shmid(B) == shmid(C) == name
        @test sizeof(A) == sizeof(B) == sizeof(C) == len
        Aptr = convert(Ptr{T},pointer(A))
        Bptr = convert(Ptr{T},pointer(B))
        Cptr = convert(Ptr{T},pointer(C))
        for i in 1:n
            unsafe_store!(Aptr,i,i)
        end
        @test all(unsafe_load(Bptr, i) == i for i in 1:n)
        @test all(unsafe_load(Cptr, i) == i for i in 1:n)
        for i in 1:n
            unsafe_store!(Bptr,-i,i)
        end
        @test all(unsafe_load(Aptr, i) == -i for i in 1:n)
        @test all(unsafe_load(Cptr, i) == -i for i in 1:n)
        @test_throws ReadOnlyMemoryError unsafe_store!(Cptr,42)
        @test ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t),
                    Aptr, 0, len) == Aptr
        @test all(unsafe_load(Cptr, i) == 0 for i in 1:n)
    end
    GC.gc() # call garbage collector to exercise the finalizers
end

@testset "Wrapped Shared Arrays " begin
    begin
        T = Float32
        dims = (3,4,5)
        for key in (IPC.PRIVATE, "/wrsharr-$(getpid())")
            if isa(key, String)
                try; shmrm(key); catch err; end
            end
            A = WrappedArray(key, T, dims)
            id = shmid(A)
            if isa(id, ShmId)
                info = ShmInfo(id)
                @test info.segsz ≥ sizeof(A) + 64
            end
            B = WrappedArray(id; readonly=false)
            C = WrappedArray(id; readonly=true)
            n = length(A)
            @test shmid(A) == shmid(B) == shmid(C) == id
            @test sizeof(A) == sizeof(B) == sizeof(C) == n*sizeof(T)
            @test eltype(A) == eltype(B) == eltype(C) == T
            @test size(A) == size(B) == size(C) == dims
            @test length(A) == length(B) == length(C) == prod(dims)
            @test all(size(A,i) == size(B,i) == size(C,i) == dims[i] for i in 1:length(dims))
            A[:] = 1:n
            @test first(A) == 1
            @test last(A) == n
            @test A[end] == n
            @test all(B[i] == i for i in 1:n)
            @test all(C[i] == i for i in 1:n)
            B[:] = -(1:n)
            @test extrema(A[:] + (1:n)) == (0, 0)
            @test all(C[i] == -i for i in 1:n)
            @test_throws ReadOnlyMemoryError C[end] = 42
            @test ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t),
                        A, 0, sizeof(A)) == pointer(A)
            @test extrema(C) == (0, 0)
        end
    end
    GC.gc() # call garbage collector to exercise the finalizers
end

end # module


text:
module FindClosest

using Distances

export findclosest

const BRUTEFORCE_THRESHOLD = 33

"""
    findclosest(pts[, dist = Euclidean()]) -> (d, (i, j))

Find the distance `d` and indices `(i, j)` of the two closest points `pᵢ` and `pⱼ` in `pts`,
according to distance `dist(pᵢ, pⱼ)` (see the Distance.jl package). If `length(pts) < 2`,
return `nothing`.

The search is performed using a single-threaded, multidimensional, divide-and-conquer
algorithm, with runtime complexity `O(n log(n))`.

# Examples
```
julia> using StaticArrays; pts = rand(SVector{3,Float64}, 10^3);

julia> findclosest(pts)
(0.007292429825826281, (75, 535))
```
"""
function findclosest(pts, dist = Euclidean())
    scratch = scratch_arrays(pts)
    pts´ = last(scratch)
    copy_itr!(pts´, enumerate(pts))

    len = length(pts´)
    len < 2 && return nothing
    len <= BRUTEFORCE_THRESHOLD && return bruteforce(pts´, dist)

    dim = length(scratch)
    sort_dim!(pts´, dim)
    T = typeof(dist_inds(pts´, (firstindex(pts´), firstindex(pts´)), dist))
    return findclosest_sorted(pts´, dist, T, scratch)
end

function findclosest_sorted(pts, dist, ::Type{T}, scratch::NTuple{D})::T where {D,T}
    len = length(pts)
    if len <= BRUTEFORCE_THRESHOLD
        d = bruteforce(pts, dist)
    elseif D == 1
        d0 = dist_inds(pts, (firstindex(pts), lastindex(pts)), dist)
        d = findclosest_sorted_1D(pts, dist, d0)
    else
        half = len ÷ 2
        ptsL = view(pts, 1:half)
        ptsR = view(pts, half + 1:len)
        dL = findclosest_sorted(ptsL, dist, T, scratch)
        dR = findclosest_sorted(ptsR, dist, T, scratch)
        dLR = min_dinds(dL, dR)
        δ = first(dLR)
        inter_range = inter_range_dim(ptsL, ptsR, δ, D, half, len)
        if length(inter_range) < 2
            d = dLR
        else
            scratch´ = Base.front(scratch)
            ptsLR = copy_view!(last(scratch´), pts, inter_range)
            sort_dim!(ptsLR, D - 1)
            if D - 1 == 1
                dinter = findclosest_sorted_1D(ptsLR, dist, dLR)
            else
                dinter = findclosest_sorted(ptsLR, dist, T, scratch´)
            end
            d = min_dinds(dLR, dinter)
        end
    end
    return d
end

function findclosest_sorted_1D(pts, dist, dLR)
    d = dLR
    ifirst, ilast = firstindex(pts), lastindex(pts)
    for i in ifirst:ilast-1
        j = i + 1
        while j <= ilast && first(last(pts[j])) - first(last(pts[i])) < first(d)
            d = min_dinds(d, dist_inds(pts, (i, j), dist))
            j += 1
        end
    end
    return d
end

function bruteforce(ps, dist = Euclidean())
    ifirst, ilast = firstindex(ps), lastindex(ps)
    d = dist_inds(ps, (ifirst, ilast), dist)
    for i in ifirst:ilast-1, j in i+1:ilast
        d = min_dinds(d, dist_inds(ps, (i, j), dist))
    end
    return d
end

## Auxiliary functions

function scratch_arrays(pts)
    p = first(pts)
    T = typeof(p)
    D = length(p)
    scratch = ntuple(i -> Vector{Tuple{Int,T}}(undef, 0), D)
    return scratch
end

function dist_inds(ps::AbstractVector{<:Tuple}, (i, j), dist)
    (i1, p1), (i2, p2) = ps[i], ps[j]
    return dist(p1, p2), sortedtuple(i1, i2)
end

dist_inds(ps, (i, j), dist) = dist(ps[i], ps[j]), sortedtuple(i, j)

min_dinds((d, is), (d´, is´)) = ifelse(d < d´, (d, is), (d´, is´))

sortedtuple(i, j) = ifelse(i < j, (i, j), (j, i))

sort_dim!(pts, dim) = sort!(pts; alg = QuickSort, by = p -> last(p)[dim])

center_dim((i, p), (i´, p´), dim) = (p[dim] + p´[dim])/2

function inter_range_dim(ptsL, ptsR, δ, D, half, len)
    center = center_dim(last(ptsL), first(ptsR), D)
    rmax = findfirst(p -> last(p)[D] > center + δ, ptsR)
    rmin = findlast( p -> last(p)[D] < center - δ, ptsL)
    inter_range =
        (rmin === nothing ? 1 : rmin + 1):(rmax === nothing ? len : half + rmax - 1)
    return inter_range
end

function copy_itr!(dst, src)
    resize!(dst, 0)
    @inbounds for s in src
        push!(dst, s)
    end
    return dst
end

function copy_view!(dst, src, inter_range)
    len = length(inter_range)
    resize!(dst, max(length(dst), len))
    dstview = view(dst, 1:len)
    copy!(dstview, view(src, inter_range))
    return dstview
end

end # module


text:
function value_of(::Type{Transformer}, arg0::JString)
    return jcall(Transformer, "valueOf", Transformer, (JString,), arg0)
end

function values(::Type{Transformer})
    return jcall(Transformer, "values", Vector{Transformer}, ())
end


text:
module util

export Vocab, setdefault!, Report, submittedEarlier, calculateDayTimestamp, DocFreq, updateDocFreq!,isduplicate, get_basic_args, aggregatePositionByFunction, sigmoid, set_doc_len!, generate_ngrams, add_ngram_to_vocab!, vocabLength, aggregatePositionByFuncKeepDim

using ArgParse
using Dates
using JSON
using DataStructures

sigmoid(x)=1.0/(1.0+exp(-x))

struct Vocab
    vocab::Dict{String, UInt32}
    ukn_set::Set{UInt32}
end

function vocabLength(vocab::Vocab)::Int64
    return length(vocab.vocab)
end

function setdefault!(vocab::Vocab, functionName::String, isUkn=false)
    tkn_id = get!(vocab.vocab, functionName, vocabLength(vocab) + 1)

    if isUkn
        push!(vocab.ukn_set, tkn_id)
    end

    return tkn_id
end

mutable struct  ReportStats
    doc_len::Int64
    doc_lens::Vector{Int64}
    weight_sum::Float64
    sorted_idxs::Vector{Vector{UInt16}}
end

mutable struct Report
    id::Int64
    dup_id::Int64
    creation_date::Int64
    creation_day_ts::UInt64
    stacks::Vector{Vector{UInt32}}
    bow::Vector{Pair{UInt32,Float64}}
    positional_bow::Vector{Pair{UInt32, Vector{Int16}}}

    v_bow::Vector{Vector{Pair{UInt32,Float64}}}
    v_positional_bow::Vector{Vector{Pair{UInt32, Vector{Int16}}}}

    stats::ReportStats
end

@inline function Report(id, dup_id, creation_date, creation_day_ts, stacks)
    return Report(id, dup_id, creation_date, creation_day_ts, stacks, Vector{Pair{UInt32,Float64}}(), Vector{Pair{UInt32, Vector{Int16}}}(), Vector{Vector{Pair{UInt32,Float64}}}(), Vector{Vector{Pair{UInt32, Vector{Int16}}}}(), ReportStats(-999999999, Vector{Int64}(), -Inf, Vector{Vector{UInt16}}()))
end

isduplicate(report::Report) = report.id != report.dup_id

submittedEarlier(r1::Report, r2::Report) = r1.creation_date < r2.creation_date || (r1.creation_date == r2.creation_date && r1.id < r2.id)

function calculateDayTimestamp(date)
    return UInt64(floor(date / (24 * 60 * 60)))
end

function set_doc_len!(report::Report)
    doc_len = 0

    for stack in report.stacks
        doc_len += length(stack)
    end

    report.stats.doc_len = doc_len

    return doc_len
end

mutable struct DocFreq
    doc_freq::Vector{Float64}
    nDocs::Float64
    freq_by_stacks::Bool
    static_df_ukn::Bool
    vocab::Vocab

    cache::Vector{Float64}
    previous_n_docs::Float64
end

DocFreq(doc_freq::Vector{Float64}, nDocs::Float64, freq_by_stacks::Bool, static_df_ukn::Bool, vocab::Vocab) = DocFreq(doc_freq, nDocs, freq_by_stacks, static_df_ukn, vocab, Vector{Float64}(), 0)

function updateDocFreq!(docfreqObj::DocFreq, reports)
    for report in reports
        if docfreqObj.freq_by_stacks
            for stack in report.stacks
                for funcId in Set(stack)
                    docfreqObj.doc_freq[funcId] += 1
                end
                docfreqObj.nDocs += 1
            end
        else
            function_set = Set(( funcId for stack in report.stacks for funcId in stack))
            docfreqObj.nDocs += 1

            for funcId in function_set
                docfreqObj.doc_freq[funcId] += 1
            end
        end
    end

    if docfreqObj.static_df_ukn
        for funcId in docfreqObj.vocab.ukn_set
            docfreqObj.doc_freq[funcId]= docfreqObj.nDocs - 1
        end
    end
end

function aggregatePositionByFunction(report::Report)::Vector{Pair{UInt32, Vector{Int16}}}
    bow = Dict{UInt32, Vector{Int16}}()

    for stack in report.stacks
        for (pos, func_id) in enumerate(stack)
            if haskey(bow, func_id)
                push!(bow[func_id], pos)
            else
                bow[func_id] = [pos]
            end
        end
    end

    for positions in values(bow)
        sort!(positions)
    end

    return collect(bow)
end

function aggregatePositionByFuncKeepDim(report::Report)::Vector{Vector{Pair{UInt32, Vector{Int16}}}}
    v_bow = Vector{Vector{Pair{UInt32, Vector{Int16}}}}(undef, length(report.stacks))

    for (idx, stack) in enumerate(report.stacks)
        bow = Dict{UInt32, Vector{Int16}}()
        for (pos, func_id) in enumerate(stack)
            if haskey(bow, func_id)
                push!(bow[func_id], pos)
            else
                bow[func_id] = [pos]
            end
        end

        for positions in values(bow)
            sort!(positions)
        end

        v_bow[idx] = sort!(collect(bow))
    end

    return v_bow
end

function add_ngram_to_vocab!(vocab::Vocab, reportid2report::Dict{Int64, Report}, ngram)
    for report in values(reportid2report)
        for ngram_tkn in generate_ngrams(report.stacks, ngram)
            setdefault!(vocab, ngram_tkn)
        end
    end
end

function generate_ngrams(stacks::Vector{Vector{UInt32}}, ngram)
    ngram_vec = Vector{String}()

    for stack in stacks
        for n=2:ngram
            v = Vector{UInt32}(undef, n)

            for start_idx=1:(length(stack) - n + 1)
                p = 1
                for idx=start_idx:(start_idx + n - 1)
                    v[p] = stack[idx]
                    p+=1
                end

                push!(ngram_vec, join(v, ","))
            end

        end
    end

    return ngram_vec
end

function generate_ngrams(stack::Vector{UInt32}, ngram)
    ngram_vec = Vector{String}()

    for n=2:ngram
        v = Vector{UInt32}(undef, n)

        for start_idx=1:(length(stack) - n + 1)
            p = 1
            for idx=start_idx:(start_idx + n - 1)
                v[p] = stack[idx]
                p+=1
            end

            push!(ngram_vec, join(v, ","))
        end

    end

    return ngram_vec
end

function get_basic_args()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--window", "-w"
            help = "Time window"
            arg_type = Int
            default = -1
        "--max_depth"
            help = ""
            arg_type = Int
            default = 300
        "--filter_recursion"
            help = "Options: none, modani and brodie"
            default = "none"
        "--keep_ukn"
            action = :store_true
        "--aggregate"
            help = "Options: max, avg_query, avg_cand, avg_short, avg_long avg_query_cand"
            default = "max"
        "--static_df_ukn"
            action = :store_true
        "--filter_func"
            help = "Options: none, top_k, top_k_trim"
            default = "none"
        "--filter_func_k"
            help = "Percentage of the k-top functions (values between 0 and 100)"
            arg_type = Float32
            default = Float32(-1.0)
        "--freq_by_stacks"
            action = :store_true
        "--upd_doc_freq"
            action = :store_true
        "--smoothing"
            action = :store_true
        "--interesting_stacks"
            action = :store_true
            help = "Following https://bazaar.launchpad.net/~bgo-maintainers/bugzilla-traceparser/3.4/view/head:/lib/TraceParser/Trace.pm,
                we select the interesting stacks that might contain the problem. This option only works for Gnome."
    end

    s.add_help = false

    return s
end

end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
using CPIMC
using CPIMC.UniformElectronGas
import CPIMC.UniformElectronGas: w

import CPIMC: orbs

@testset "λ(N, rs, d) and rs(N::Int, λ::Float64, d) are consistent" begin
    N = 18
    r = 0.56
    @test rs(N, λ(N, r, 3), 3) ≈ r
    @test rs(N, λ(N, r, 2), 2) ≈ r
end

@testset " consistence of β(Θ::Float64, rs::Float64, ξ, d) with λ(N, rs, d) and rs(N::Int, λ::Float64, d) " begin
    N = 18
    r = 0.56
    Θ = 0.123
    ξ = 0.254
    α = (4 / (9π))^(1/3)
    d = 3
    @test β(Θ, N, ξ, d) ≈ (α * r)^2 * 16 / ( (1 + ξ)^(2/d) * (2π)^4 * λ(N, r, d)^2 * Θ )
end

e = CEnsemble(2,5.68089,7)

S = sphere_with_same_spin(PlaneWave((0,0,0)),dk=1)
a = PlaneWave((-2,0,0))
b = PlaneWave((3,0,0))
c = PlaneWave((0,0,0))
d = PlaneWave((1,0,0))

sd = Kinks( ImgTime(0.5) => T4(a,b,c,d),
            ImgTime(0.6) => T4(c,d,a,b))

conf = Configuration(sphere(PlaneWave((0,0,0),Up),dk=1),sd)
conf_pol = Configuration(sphere_with_same_spin(PlaneWave((0,0,0),Up),dk=1),sd)

@testset "ΔW_diag for polarized occupation" begin

    mod = UEG()

    a = PlaneWave((-2,0,0))
    b = PlaneWave((3,0,0))
    c = PlaneWave((0,0,0))
    d = PlaneWave((1,0,0))

    ## occupation to be changed
    occ = Set([a,b,c,d])

    ### Test 2-particle excitation
    # choose creator orbitals
    i = PlaneWave((1,1,1))
    j = PlaneWave((0,-1,-1))
    # choose annihilator orbitals
    k = c
    l = d
    @assert iszero( i.vec + j.vec - k.vec - l.vec ) " momentum not conserved for this excitation "
    @assert (i.spin == k.spin) & (j.spin == l.spin) " spin is not conserved for this excitation "

    # this functions sums all combinations of orbitals o ∈ occ with annihilated orbitals k (o ≠ k) and l (o ≠ l)
    # and substracts the sum of all combinations of orbitals o ∈ ( occ ̸ {k,l} ) ∪ {i,j} of the new occupation created with orbitals i (o ≠ i) and j (o ≠ j)
    @test ΔW_diag(mod, i, j, k, l, occ) ≈ ( w(mod, a,k,k,a) + w(mod, b,k,k,b) + w(mod, a,l,l,a) + w(mod, b,l,l,b) + w(mod, k,l,l,k)
                                      - w(mod, a,i,i,a) - w(mod, b,i,i,b) - w(mod, a,j,j,a) - w(mod, b,j,j,b) - w(mod, i,j,j,i) )

    ### Test 1-particle excitation
    i = PlaneWave((-3,0,0))
    j = b
    # this functions sums all combinations of orbitals o ∈ occ with the annihilated orbital j (o ≠ j)
    # and substracts the sum of all combinations of orbitals o ∈ ( occ ̸ {j} ) ∪ {i} of the new occupation created with orbital i (o ≠ i)
    @test ΔW_diag(mod, i, j, occ) ≈ ( w(mod, a,j,j,a) + w(mod, c,j,j,c) + w(mod, d,j,j,d)
                                - w(mod, a,i,i,a) - w(mod, c,i,i,c) - w(mod, d,i,i,d) )
end

@testset "ΔW_diag for unpolarized occupation" begin

    mod = UEG()

    a = PlaneWave((-2,0,0),Up)
    b = PlaneWave((3,0,0),Down)
    c = PlaneWave((0,0,0),Up)
    d = PlaneWave((1,0,0),Down)

    ## occupation to be changed
    occ = Set([a,b,c,d])

    ### Test 2-particle excitation
    # choose creator orbitals
    i = PlaneWave((1,1,1),Up)
    j = PlaneWave((0,-1,-1),Down)
    # choose annihilator orbitals
    k = c
    l = d
    @assert iszero( i.vec + j.vec - k.vec - l.vec ) " momentum not conserved for this excitation "
    @assert (i.spin == k.spin) & (j.spin == l.spin) " spin is not conserved for this excitation "

    # this functions sums all combinations of orbitals o ∈ occ with annihilated orbitals k (o ≠ k) and l (o ≠ l)
    # and substracts the sum of all combinations of orbitals o ∈ ( occ ̸ {k,l} ) ∪ {i,j} of the new occupation created with orbitals i (o ≠ i) and j (o ≠ j)
    @test ΔW_diag(mod, i, j, k, l, occ) ≈ ( w(mod, a,k,k,a) + w(mod, b,k,k,b) + w(mod, a,l,l,a) + w(mod, b,l,l,b) + w(mod, k,l,l,k)
                                      - w(mod, a,i,i,a) - w(mod, b,i,i,b) - w(mod, a,j,j,a) - w(mod, b,j,j,b) - w(mod, i,j,j,i) )

    ## test for orbitals where a momentum vector of one creator orbital is equal to the momentum vector of one orbital in the occupation
    # choose creator orbitals
    i = PlaneWave((-2,0,0),Down)
    j = PlaneWave((5,0,0),Up)
    # choose annihilator orbitals
    k = b
    l = c
    @assert iszero( i.vec + j.vec - k.vec - l.vec ) " momentum not conserved for this excitation "
    @assert (i.spin == k.spin) & (j.spin == l.spin) " spin is not conserved for this excitation "

    # the contribution w(a,i,i,a) does not occur here because a.vec == i.vec
    @test ΔW_diag(mod, i, j, k, l, occ) ≈ ( w(mod, a,k,k,a) + w(mod, d,k,k,d) + w(mod, a,l,l,a) + w(mod, d,l,l,d) + w(mod, k,l,l,k)
                                      # - w(a,i,i,a) does not occur here because a.vec == i.vec
                                       - w(mod, d,i,i,d) - w(mod, a,j,j,a) - w(mod, d,j,j,d) - w(mod, i,j,j,i) )

    ### Test 1-particle excitation
    i = PlaneWave((-3,0,0),Down)
    j = b
    # this functions sums all combinations of orbitals o ∈ occ with the annihilated orbital j (o ≠ j)
    # and substracts the sum of all combinations of orbitals o ∈ ( occ ̸ {j} ) ∪ {i} of the new occupation created with orbital i (o ≠ i)
    @test ΔW_diag(mod, i, j, occ) ≈ ( w(mod, a,j,j,a) + w(mod, c,j,j,c) + w(mod, d,j,j,d)
                                - w(mod, a,i,i,a) - w(mod, c,i,i,c) - w(mod, d,i,i,d) )
end


text:
using PEPOKit
using TensorKit
using KrylovKit
using LinearAlgebra
using Test
using TestSetExtensions

@testset ExtendedTestSet "All tests" begin
    @includetests ARGS
end

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
include("polyagamma/polyagamma_mcmc.jl")
include("mmle/mmle.jl")

function calibrate_item!(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{Response};
    method = "mmle",
    kwargs...
    )
    if method == "pg"
        return calibrate_item_pg!(item, examinees, responses; kwargs...)
    elseif method =="mmle"
        return calibrate_item_mmle!(item, examinees, responses; kwargs...)
    end
end


text:
using DiffEqFlux, Flux, GalacticOptim, OrdinaryDiffEq, Test

n = 2  # number of ODEs
tspan = (0.0, 1.0)

d = 5  # number of data pairs
x = rand(n, 5)
y = rand(n, 5)

cb = function (p,l)
  @show l
  false
end

using Random
Random.seed!(100)

NN = Chain(Dense(n, 10n, tanh),
           Dense(10n, n))

@info "ROCK4"
nODE = NeuralODE(NN, tspan, ROCK4(), reltol=1e-4, saveat=[tspan[end]])

loss_function(θ) = Flux.mse(y, nODE(x, θ))
l1 = loss_function(nODE.p)

res = DiffEqFlux.sciml_train(loss_function, nODE.p, NewtonTrustRegion(), GalacticOptim.AutoZygote(), maxiters = 200, cb=cb)
@test 5loss_function(res.minimizer) < l1
res = DiffEqFlux.sciml_train(loss_function, nODE.p, Optim.KrylovTrustRegion(), GalacticOptim.AutoZygote(), maxiters = 200, cb=cb)
@test 5loss_function(res.minimizer) < l1

NN = FastChain(FastDense(n, 10n, tanh),
               FastDense(10n, n))

@info "ROCK2"
nODE = NeuralODE(NN, tspan, ROCK2(), reltol=1e-4, saveat=[tspan[end]])

loss_function(θ) = Flux.mse(y, nODE(x, θ))
l1 = loss_function(nODE.p)
optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_function(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, nODE.p,)

res = GalacticOptim.solve(optprob, NewtonTrustRegion(), maxiters = 200, cb=cb)
@test 5loss_function(res.minimizer) < l1
res = GalacticOptim.solve(optprob, Optim.KrylovTrustRegion(), maxiters = 200, cb=cb, allow_f_increases = true)
@test 5loss_function(res.minimizer) < l1


text:
 using Agents: vertices
 using MetaGraphs: add_vertex!, has_edge
 using Agents, AgentsPlots, Plots, LightGraphs, MetaGraphs, GraphPlot, GraphRecipes, NetworkLayout
 using Tables, DataFrames
 using CSV, JSON, Serialization, DelimitedFiles
 using BritishNationalGrid
 using ZipFile, Shapefile
 using Random
 using Match
 using LinearAlgebra
 using StatsBase
 using Distributions
 using StatsPlots
 using SparseArrays
 using Laplacians
 using DataStructures
 using RollingFunctions
 using BenchmarkTools
 using Statistics
 using Logging,LoggingExtras, LoggingFacilities
 using Dates,TimeZones
 using PyCall
#Core structures
 include("core/core_structs.jl")
#graph-related
 include("core/graph_functions.jl")
#events
 include("eve/artificial_events.jl")
#various util functions
 include("utils/util_functions.jl")
#logging
 include("utils/logging_functions.jl")
#Plotting functions
 include("utils/plotting_functions.jl")
# time-to-event functions
 include("utils/tte_functions.jl")
#running functions
 include("utils/running_functions.jl")

#maintenance model
 include("phy/maintenance_model.jl")

 include("ntw/of_switch.jl")
 include("ntw/of_control.jl")
 include("ctl/agent_control.jl")

 include("phy/physical_model.jl")

 include("ntw/network_model.jl")
 include("phy/geo_model.jl")

#Main Functions
 include("model/netManFunctions.jl")

#Agents.jl function implementation for this model
 include("model/netManModel.jl")

#Basic queries
#include("model/queries_basic.jl")
#Multiple queries
 include("ctl/queries_multiple.jl")

export load_run_configs, single_run_with_logging

# GurobiPy Optimisation config
script_dir = pwd() * "/src/pyopt/" #@__DIR__
pushfirst!(PyVector(pyimport("sys")."path"), script_dir)
opt_init = pyimport("optimisation_initialisation")
opt_run = pyimport("maintenance_planning")
np = pyimport("numpy")

# load sim config file
base_cfgs = load_base_cfgs("data/configs/configs.csv")
#first row of the file
bcfg = base_cfgs[1]

#obtain random services that provide this coverage of the network
# coverage = 0.95
# ntw_topo = GraphModel(bcfg.ntw_topo_n)
# ntw_services = get_end_points(bcfg.seed,get_graph(bcfg.seed,bcfg.size,ntw_topo;k=bcfg.k,adj_m_csv=bcfg.ntw_csv_adj_matrix),coverage)

#small scale scenario
# remember to set ntw topo and other params in config.csv file. ntw_topo_n: 4, size = 16, mnt_policy: 2
ntw_services = [(1,7),(4,1),(5,14),(12,8)]
init_sne_params = (ids=[2,6,9,13,10,15,16],ruls=[65,78,84,86,93,49,90])

#large scale scenario
# services = [(23,80),(18,48),(56,92),(64,31),(48,70)]
# init_sne_params = (ids = [36,6,3,5,12,57,2,19,1,25,8], ruls = [34,42,48,50,53,73,82,88,90,94,97])

cfg = config(bcfg,ntw_services,init_sne_params)

single_run_with_logging(cfg)


text:
@doc """
Generate a closure over the roll-call data set that will compute the
unnormalized negative log-posterior at any set of parameters, `θ`.
""" ->
function make_nlp(
    D::Int,
    σ_x::Float64,
    σ_y::Float64,
    n_legislators::Int,
    n_bills::Int,
    legislators::Vector{Int},
    bills::Vector{Int},
    votes::Vector{Float64},
)
    function nlp(θ::Vector)
        n = length(legislators)

        lp = 0.0

        # Log-posterior contribution from the likelihood.
        for i in 1:n
            l = legislators[i]
            b = bills[i]
            v = votes[i]
            z = 0.0
            for d in 1:D
                x_l_d = θ[D * (l - 1) + d]
                x_b_d = θ[D * n_legislators + D * (b - 1) + d]
                z += x_l_d * x_b_d
            end
            y_b = θ[D * n_legislators + D * n_bills + b]
            z += y_b
            p = logistic(z)
            lp += log(v * p + (1 - v) * (1 - p))
        end

        # Log-posterior contribution from the prior.
        for i in 1:(D * n_legislators + D * n_bills)
            lp += logpdf(Normal(0, σ_x), θ[i])
        end
        for i in (D * n_legislators + D * n_bills + 1):(D * n_legislators + D * n_bills + n_bills)
            lp += logpdf(Normal(0, σ_y), θ[i])
        end

        -lp
    end

    nlp
end

@doc """
Generate a closure over the roll-call data set that will compute the gradient
of the unnormalized negative log-posterior at any set of parameters, `θ`. The
gradient will be stored in the second argument, `gr`, which will be mutated.
""" ->
function make_nlp_gr(
    D::Int,
    σ_x::Float64,
    σ_y::Float64,
    n_legislators::Int,
    n_bills::Int,
    legislators::Vector{Int},
    bills::Vector{Int},
    votes::Vector{Float64},
)
    function nlp_gr!(θ::Vector, gr::Vector)
        n = length(legislators)

        fill!(gr, 0.0)

        # Log-posterior contribution from the likelihood.
        for i in 1:n
            l = legislators[i]
            b = bills[i]
            v = votes[i]
            z = 0.0
            for d in 1:D
                x_l_d = θ[D * (l - 1) + d]
                x_b_d = θ[D * n_legislators + D * (b - 1) + d]
                z += x_l_d * x_b_d
            end
            y_b = θ[D * n_legislators + D * n_bills + b]
            z += y_b
            p = logistic(z)
            for d in 1:D
                x_l_d = θ[D * (l - 1) + d]
                x_b_d = θ[D * n_legislators + D * (b - 1) + d]
                gr[D * (l - 1) + d] += (p - v) * x_b_d
                gr[D * n_legislators + D * (b - 1) + d] += (p - v) * x_l_d
            end
            gr[D * n_legislators + D * n_bills + b] += (p - v)
        end

        # Log-posterior contribution from the prior.
        for i in 1:(D * n_legislators + D * n_bills)
            gr[i] += θ[i] / σ_x^2
        end
        for i in (D * n_legislators + D * n_bills + 1):(D * n_legislators + D * n_bills + n_bills)
            gr[i] += θ[i] / σ_y^2
        end

        nothing
    end

    nlp_gr!
end

@doc """
Fit an ideal points model to roll call data from an ORD file. Defaults to
fitting latent space positions in 1 dimension; this can changed by increasing
the value of the second argument. To ensure that left-wing candidates take on
negative valued latent space positions, all legislators from a specific party
idea are not in the non-negative orthant. For long-running model fits, it is
possible to enable a trace of the optimization process to ensure that the
model fit is improving with each optimization pass.
""" ->
function ideal_points(
    ord_file::ORDFile,
    D::Integer = 1,
    σ_x::Real = 1.0,
    σ_y::Real = 25.0,
    party_id::Integer = 100,
    show_progress::Bool = false,
)
    parties = RollCallDataIO.parties(ord_file)
    roll_call = convert(SparseRollCall, RollCallDataIO.roll_call(ord_file))
    n_legislators, n_bills = roll_call.n_legislators, roll_call.n_bills

    nlp = make_nlp(
        D,
        σ_x,
        σ_y,
        n_legislators,
        n_bills,
        roll_call.legislators,
        roll_call.bills,
        roll_call.votes,
    )

    nlp_gr! = make_nlp_gr(
        D,
        σ_x,
        σ_y,
        n_legislators,
        n_bills,
        roll_call.legislators,
        roll_call.bills,
        roll_call.votes,
    )

    # Force members of party with ID == party_id to one part of the latent
    # space to ensure consistent placement across runs. The default party_id
    # will always place US Democrats to the left for consistency with
    # left-wing terminology.
    Θ₀ = fill(0.0, D * n_legislators + D * n_bills + n_bills)
    for l in 1:n_legislators
        for d in 1:D
            Θ₀[D * (l - 1) + d] = ifelse(
                parties[l] == party_id,
                -1.0 + d / l,
                1.0 + d / l,
            )
        end
    end

    fit = Optim.optimize(
        nlp,
        nlp_gr!,
        Θ₀,
        method = :l_bfgs,
        show_trace = show_progress,
    )

    legislator_slopes = fit.minimum[1:(D * n_legislators)]
    bill_slopes = fit.minimum[
        (D * n_legislators + 1):(D * n_legislators + D * n_bills)
    ]
    bill_intercepts = fit.minimum[
        (D * n_legislators + D * n_bills + 1):(D * n_legislators + D * n_bills + n_bills)
    ]

    (
        legislator_slopes,
        bill_slopes,
        bill_intercepts,
        fit,
    )
end


text:
module DomainSets

using StaticArrays
using LinearAlgebra, Statistics
import LinearAlgebra: cross, ×, pinv

using IntervalSets
using CompositeTypes, CompositeTypes.Display, CompositeTypes.Indexing

# deprecations in v0.5
@deprecate IntersectionDomain IntersectDomain
@deprecate DifferenceDomain SetdiffDomain
@deprecate FlexibleUnitCube DynamicUnitCube
@deprecate FlexibleUnitSphere DynamicUnitSphere
@deprecate FlexibleUnitBall DynamicUnitBall
@deprecate Composition ComposedMap

@deprecate element component
@deprecate elements components
@deprecate numelements ncomponents

################################
## Exhaustive list of imports
################################

# Generated functions
import Base: @ncall

# Operator symbols
import Base: *, +, -, /, \, ^,
    |, &,
    ∪, ∩,
    ==, isapprox,
    ∘,
    # Set operations
    intersect, union, setdiff, in, isempty, minimum, maximum,
    issubset,
    # Arrays
    eltype, hash, isreal,
    # Types, promotions and conversions
    convert, promote,
    # Display
    show

# IntervalSets
import IntervalSets: (..), endpoints, Domain, AbstractInterval, TypedEndpointsInterval,
                        leftendpoint, rightendpoint, isleftopen, isleftclosed,
                        isrightopen, isrightclosed, isopenset, isclosedset,
                        infimum, supremum
export ..

import CompositeTypes: component, components

################################
## Exhaustive list of exports
################################

## Utils

# from util/common.jl
export prectype, numtype,
    convert_numtype, promote_numtype,
    convert_prectype, promote_prectype,
    iscomposite, component, components, ncomponents

## Maps

# from maps/map.jl
export AbstractMap, Map, TypedMap,
    applymap,
    domaintype, codomaintype,
    inverse, leftinverse, rightinverse,
    mapsize, jacobian, jacdet, diffvolume
# from maps/composite.jl
export ComposedMap, composedmap, ∘
# from maps/product.jl
export ProductMap, productmap
# from maps/basic.jl
export IdentityMap,
    StaticIdentityMap, VectorIdentityMap,
    ZeroMap, UnityMap, ConstantMap,
    isconstant, constant
# from maps/affine.jl
export AffineMap, Translation, LinearMap,
    matrix, vector,
    islinear, isaffine

## Generic domains

# from generic/domain.jl
export Domain, EuclideanDomain, VectorDomain,
    dimension,
    approx_in,
    isopenset, isclosedset, iscompact,
    boundary, ∂,
    boundingbox,
    interior, closure,
    volume,
    point_in_domain,
    normal, tangents, distance_to,
    canonicaldomain, mapto_canonical, mapfrom_canonical, hascanonicaldomain,
    mapto,
    parameterdomain, parameterization, hasparameterization,
    mapfrom_parameterdomain, mapto_parameterdomain

# from generic/lazy.jl
export superdomain

# from generic/productdomain.jl
export ProductDomain, productdomain,
    VcatDomain, VectorProductDomain, TupleProductDomain

# from generic/mapped.jl
export MappedDomain,
    map_domain,
    mapped_domain,
    forward_map,
    inverse_map

# from generic/setoperations.jl
export UnionDomain, uniondomain,
    IntersectDomain, intersectdomain,
    SetdiffDomain, setdiffdomain

# from applications/rotations.jl
export rotate

export infimum, supremum

## Specific domains

# from domains/trivial.jl
export EmptySpace, FullSpace, TypeDomain,
    emptyspace, fullspace, typedomain
# from domains/numbers.jl
export Integers, RealNumbers, Rationals, ComplexNumbers,
    ℕ, ℤ, ℚ, ℝ, ℂ, ℝ1, ℝ2, ℝ3, ℝ4
# from domains/interval.jl
export AbstractInterval, Interval, UnitInterval, ChebyshevInterval,
    OpenInterval, ClosedInterval,
    leftendpoint, rightendpoint, isleftopen, isrightopen,
    HalfLine, NegativeHalfLine
# from domains/simplex.jl
export UnitSimplex,
    StaticUnitSimplex, DynamicUnitSimplex,
    EuclideanUnitSimplex, VectorUnitSimplex,
    corners
# from domains/point.jl
export Point
# from domains/ball.jl
export Ball, UnitBall,
    center, radius,
    StaticUnitBall, DynamicUnitBall,
    EuclideanUnitBall, VectorUnitBall,
    Sphere, UnitSphere,
    StaticUnitSphere, DynamicUnitSphere,
    VectorUnitSphere, EuclideanUnitSphere,
    Disk, UnitDisk, VectorUnitDisk,
    UnitCircle, VectorUnitCircle,
    ComplexUnitCircle, ComplexUnitDisk,
    ellipse, ellipse_shape, cylinder
# from domains/cube.jl
export UnitCube,
    StaticUnitCube, DynamicUnitCube,
    EuclideanUnitCube, VectorUnitCube,
    UnitSquare, UnitCube,
    Rectangle
# from domain/levelset.jl
export LevelSet, ZeroSet,
    SublevelSet, SubzeroSet,
    SuperlevelSet, SuperzeroSet,
    pseudolevel
# from domain/indicator.jl
export IndicatorFunction

## Applications
# from applications/rotation.jl
export rotation_map,
    CartToPolarMap, PolarToCartMap

include("util/common.jl")

include("maps/map.jl")
include("maps/lazy.jl")
include("maps/inverse.jl")
include("maps/jacobian.jl")
include("maps/composite.jl")
include("maps/product.jl")
include("maps/isomorphism.jl")
include("maps/basic.jl")
include("maps/affine.jl")
include("maps/arithmetics.jl")

include("generic/domain.jl")
include("generic/canonical.jl")
include("generic/lazy.jl")
include("generic/productdomain.jl")
include("generic/setoperations.jl")
include("generic/mapped.jl")
include("generic/broadcast.jl")

include("domains/trivial.jl")
include("domains/numbers.jl")
include("domains/levelset.jl")
include("domains/point.jl")
include("domains/interval.jl")
include("domains/simplex.jl")
include("domains/ball.jl")
include("domains/cube.jl")
include("domains/indicator.jl")
include("domains/boundingbox.jl")

include("applications/coordinates.jl")
include("applications/rotation.jl")

end # module


text:
function get_half_track(obj::HalfTrackSampler)
    return jcall(obj, "getHalfTrack", List, ())
end

function handle_step(obj::HalfTrackSampler, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "handleStep", void, (SpacecraftState, jboolean), arg0, arg1)
end

function init(obj::OrekitFixedStepHandler, arg0::SpacecraftState, arg1::AbsoluteDate, arg2::jdouble)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate, jdouble), arg0, arg1, arg2)
end


text:
module MacrosTests

end # module


text:
module PTMCMC
using AbstractMCMC
using Distributions
using Requires

import Random

# Abstract type for PT-style samplers

end


text:
module mnist

using Statistics

# big endian (network order) !!! must use ntoh(x)

export MnistData,
    getImages,
    getMeanImageLabel

struct MnistData
    images::Array{UInt8,3}
    labels::Array{UInt8,1}
    #
    function MnistData(filenameImg, filenameLabel)
        imagesT = readidxImagefile(filenameImg)
        labelsT = readidxLabelsfile(filenameLabel)
        new(imagesT,labelsT)
    end
end

function readidxImagefile(filename::String)
    idfile = open(filename)
    # read 4 bytes (magic)
    val = read(idfile, UInt32)
    magic = Int32(ntoh(val))
    val = read(idfile, UInt32)
    nbitem = Int32(ntoh(val))
    val = read(idfile, UInt32)
    nbrow = Int32(ntoh(val))
    val = read(idfile, UInt32)
    nbcol = Int32(ntoh(val))
    images = Array{UInt8,3}(undef, nbrow,nbcol, nbitem)
    # then read matrix of bytes by row
    for k in 1:nbitem
        for i in 1:nbrow
            row = Vector{UInt8}(undef, nbcol)
            unsafe_read(idfile, pointer(row), nbcol*sizeof(UInt8))
            images[i,:,k] = row
        end
    end
    images
end

function readidxLabelsfile(filenameLabel)
   idfile = open(filenameLabel)
    # read 4 bytes (magic)
    val = read(idfile, UInt32)
    magic = Int32(ntoh(val))
    val = read(idfile, UInt32)
    nbitem = Int32(ntoh(val))
    #
    labels = Vector{UInt8}(undef, nbitem)
    unsafe_read(idfile, pointer(labels), sizeof(UInt8)*nbitem)
    labels
end

function getImages(mnist::MnistData , label::UInt8)
    idx=findall(x-> x==label, mnist.labels)
    return mnist.images[:,:,idx]
end

function getMeanImageLabel(mnist::MnistData,  label::UInt8)
    idx=findall(x-> x==label, mnist.labels)
    meanIL=zeros(Float64, mnist.images[:,:,1])
    meanIL = mean(mnist.images[:,:,idx],3)
    meanIL[:,:,1]
end

end


text:
export
    DiophantineEquationQuadraticXY,
    diophantine_equation_quadratic_xy

# Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
struct DiophantineEquationQuadraticXY{T<:Integer}
    cx²::T # A
    cxy::T # B
    cy²::T # C
    cx::T  # D
    cy::T  # E
    c0::T  # F
end

diophantine_equation_quadratic_xy(;cx²::T=0, cxy::T=0, cy²::T=0, cx::T=0, cy::T=0, c0::T=0) where {T<:Integer} = DiophantineEquationQuadraticXY(cx², cxy, cy², cx, cy, c0)

Base.show(io::IO, eq::DiophantineEquationQuadraticXY) = print(io, "$(eq.cx²)x² + $(eq.cxy)xy + $(eq.cy²)y² + $(eq.cx)x + $(eq.cy)y + $(eq.c0) = 0")

evaluate(eq::DiophantineEquationQuadraticXY{T}, sol::DiophantineSolutionXY{T}) where {T<:Integer} = evaluate(eq, sol.x, sol.y)
evaluate(eq::DiophantineEquationQuadraticXY{T}, x::T, y::T) where {T<:Integer} = eq.cx² * x^2 + eq.cxy * x * y + eq.cy² * y^2 + eq.cx * x + eq.cy * y + eq.c0

solve(eq::DiophantineEquationQuadraticXY{T}) where {T<:Integer} = begin
    cx², cxy, cy², cx, cy, c0 = eq.cx², eq.cxy, eq.cy², eq.cx, eq.cy, eq.c0

    cx² == cxy == cy² == 0 && return AbstractDiophantineSolutions{DiophantineSolutionXY{T}}[solve(diophantine_equation_linear_xy(cx=cx, cy=cy, c0=c0))]
    cx² == cy² == 0 && return solve_simplehyperbolic(eq)

    discr = cxy^2 - 4*cx²*cy²
    discr < 0 && return solve_elliptical(eq)
    discr == 0 && return solve_parabolic(eq)
    solve_hyperbolic(eq)
end

solve_simplehyperbolic(eq::DiophantineEquationQuadraticXY{T}) where {T<:Integer} = begin
    cxy, cx, cy, c0 = eq.cxy, eq.cx, eq.cy, eq.c0

    solutions = AbstractDiophantineSolutions{DiophantineSolutionXY{T}}[]

    d = cx * cy - cxy * c0
    if d == 0
        cy % cxy == 0 && push!(solutions, diophantine_onex_anyy(-(cy ÷ cxy)))
        cx % cxy == 0 && push!(solutions, diophantine_anyx_oney(-(cx ÷ cxy)))
    else
        xytuples = Tuple{T,T}[]
        for f in factors(abs(d), true)
            x, r = divrem(f - cy, cxy)
            r == 0 || continue
            y, r = divrem(d ÷ f - cx, cxy)
            r == 0 || continue
            push!(xytuples, (x, y))
        end
        push!(solutions, isempty(xytuples) ? diophantine_nonex_noney(T) : diophantine_solutions(xytuples))
    end

    solutions
end

solve_elliptical(eq::DiophantineEquationQuadraticXY{T}) where {T<:Integer} = begin
    cx², cxy, cy², cx, cy, c0 = eq.cx², eq.cxy, eq.cy², eq.cx, eq.cy, eq.c0

    # f(t) = (cxy^2 - 4*cx²*cy²)*t^2 + 2(cxy*cy - 2*cy²*cx)*t + (cy^2 - 4*cy²*c0)
    fa = (cxy^2 - 4*cx²*cy²)
    fb = (cxy*cy - 2*cy²*cx)
    fc = (cy^2 - 4*cy²*c0)
    f(t) = fa*t^2 + 2*fb*t + fc
    D = fb^2 - fa*fc

    D < 0 && return AbstractDiophantineSolutions{DiophantineSolutionXY{T}}[diophantine_nonex_noney(T)]

    z = map(x -> Base.trunc(T, x), Float64[(-fb-sqrt(D)) / fa, (-fb+sqrt(D)) / fa])
    xytuples = Tuple{T,T}[]

    sort!(z)
    for x in z[1]:z[2]
        u = f(x)
        v = isqrt(u)
        v*v == u || continue

        y, r = divrem(-(cxy*x + cy) + v, 2*cy²)
        r == 0 && push!(xytuples, (x, y))

        v != 0 || continue
        y, r = divrem(-(cxy*x + cy) - v, 2*cy²)
        r == 0 && push!(xytuples, (x, y))
    end

    solution = isempty(xytuples) ? diophantine_nonex_noney(T) : diophantine_solutions(xytuples)
    AbstractDiophantineSolutions{DiophantineSolutionXY{T}}[solution]
end

solve_parabolic(eq::DiophantineEquationQuadraticXY{T}) where {T<:Integer} = begin
    cx², cxy, cy², cx, cy, c0 = eq.cx², eq.cxy, eq.cy², eq.cx, eq.cy, eq.c0

    solutions = AbstractDiophantineSolutions{DiophantineSolutionXY{T}}[]

    g = copysign(gcd(cx², cy²), cx²)
    a, b, c = cx² ÷ g, cxy ÷ g, cy² ÷ g
    ra, rc = isqrt(a), copysign(isqrt(c), b)

    fa = ra * g
    fb = cx
    fc = ra * c0
    f(t) = fa * t^2 + fb * t + fc

    d = rc * cx - ra * cy
    if d == 0
        D = fb^2 - 4*fa*fc
        if D < 0
            push!(solutions, diophantine_nonex_noney(T))
        else
            u = map(x -> Base.trunc(T, x), Float64[(-fb-sqrt(D)) / (2*fa), (-fb+sqrt(D)) / (2*fa)])

            for v in unique(u)
                f(v) == 0 && push!(solutions, solve(diophantine_equation_linear_xy(cx=ra, cy=rc, c0=-v)))
            end
        end
    else
        for u in 0:(abs(d) - 1)
            f(u) % d == 0 || continue
            push!(solutions, diophantine_quadraticx_quadraticy(
                rc * g * (-d), -(cy + 2 * rc * g * u), -((rc * g * u^2 + cy * u + rc * c0) ÷ d),
                ra * g *   d,    cx + 2 * ra * g * u,    (ra * g * u^2 + cx * u + ra * c0) ÷ d
            ))
        end
    end

    solutions
end

solve_hyperbolic(eq::DiophantineEquationQuadraticXY{T}) where {T<:Integer} = begin
    (eq.cx == eq.cy == 0) ?
        solve_hyperbolic_homogeneous(eq) :
        solve_hyperbolic_general_quadratic(eq)
end

solve_hyperbolic_homogeneous(eq::DiophantineEquationQuadraticXY{T}) where {T<:Integer} = begin
    cx², cxy, cy², c0 = eq.cx², eq.cxy, eq.cy², eq.c0

    solutions = AbstractDiophantineSolutions{DiophantineSolutionXY{T}}[]

    discr = cxy^2 - 4*cx²*cy²
    k = isqrt(discr)

    if c0 == 0
        push!(solutions, diophantine_solutions((0, 0)))
        if k*k == discr
            push!(solutions, solve(diophantine_equation_linear_xy(cx=2*cx², cy=cxy+k, c0=0)))
            push!(solutions, solve(diophantine_equation_linear_xy(cx=2*cx², cy=cxy-k, c0=0)))
        end
        return solutions
    end

    if k*k == discr
        xytuples = Tuple{T,T}[]
        for f in factors(abs(-4*cx²*c0), true)
            y, r = divrem(f + (4*cx²*c0) ÷ f, 2*k)
            r == 0 || continue
            x, r = divrem(f - (cxy + k)*y, 2*cx²)
            r == 0 || continue
            push!(xytuples, (x, y))
        end
        push!(solutions, isempty(xytuples) ? diophantine_nonex_noney(T) : diophantine_solutions(xytuples))
        return solutions
    end

    # WIP
end


text:
α = abs.(rand(Uniform(-1.0, 1.0), (ns))) .+ 0.5

n_nodes = ns;
nn_W = Chain(x -> x,
              Dense(ns, n_nodes, gelu),
              Dense(n_nodes, n_nodes, gelu),
              Dense(n_nodes, n_nodes, gelu),
              Dense(n_nodes, n_nodes, gelu),
              Dense(n_nodes, ns))
pW, re = Flux.destructure(nn_W);

p = vcat(α, pW)

function cellbox!(du, u, p, t)
    Wu = re(@view(p[ns+1:end]))(u)
    pα = @view(p[1:ns])
    # du = Wu - μ - u
    @. du = pα * tanh(Wu - μ) - u
end

tspan = (0, tfinal);
ts = 0:tspan[2]/ntotal:tspan[2];
ts = ts[2:end];
u0 = zeros(ns);
prob = ODEProblem(cellbox!, u0, tspan, saveat=ts);

predict_neuralode(u0, p, 1);
@show loss_neuralode(p, 1);

cbi = function (p, i_exp)
    ode_data = ode_data_list[i_exp, :, :]
    _ts = 0:ts[end]/nplot:ts[end]
    pred = predict_neuralode(u0, p, i_exp, nothing, false)
    # gold = predict_neuralode(u0, p_gold, i_exp, nothing, false)
    l_plt = []
    for i in 1:minimum([10, ns])
        plt = Plots.scatter(ts, ode_data[i,:], label="Data");
        plot!(plt, _ts, pred[i,:], label="Prediction");
        # plot!(plt, _ts, gold[i,:], label="Ground truth");
        ylabel!(plt, "x$i")
        xlabel!(plt, "Time")
        push!(l_plt, plt)
    end
    plt_all = plot(l_plt..., legend=false, size=(1000, 1000));
    png(plt_all, string(fig_path, "/conditions/i_exp_", i_exp))
    return false
end

function show_network(p)
    nothing
end

function loss_network(p)
     return 1.0, 1.0
end

cbp = function (p, iter)
    nothing
end

cbi(p, 1)


text:
# Equation references to
#
# - Bethe, H., & Salpeter, E. (1977). Quantum mechanics of one-and
#   two-electron atoms. New York: Plenum Pub. Corp.

@doc raw"""
    non_relativistic_energy(n[, Z=1])

Returns the non-relativistic energy of a hydrogen orbital with
principal quantum number ``n`` (and an effective nuclear charge ``Z``)
according to

```math
\begin{equation}
\tag{BS2.11}
E_n =
-\frac{Z^2}{2n^2}
\end{equation}
```

"""
non_relativistic_energy(n, Z=1) = -Z^2/2n^2

"""
    non_relativistic_energy(o)

Return the non-relativistic energy of the orbital `o`.

# Examples

```jldoctest
julia> non_relativistic_energy(o"1s")
-0.5

julia> non_relativistic_energy(o"2p")
-0.125

julia> non_relativistic_energy(ro"2p")
-0.125

julia> non_relativistic_energy(ro"2p-")
-0.125
```
"""
non_relativistic_energy(o::AbstractOrbital) = non_relativistic_energy(o.n)
non_relativistic_energy(o::SpinOrbital) = non_relativistic_energy(o.orb)

@doc raw"""
    relativistic_energy(n, κ[, Z=1])

Returns the relativistic energy of a hydrogen orbital with principal
quantum number ``n`` and "angular momentum" quantum number ``κ`` (and
an effective nuclear charge ``Z``) according to

```math
\begin{equation}
\tag{BS14.29*}
E_{nκ} =
c^2\left[
-1+
\frac{1}{\sqrt{1 +
\left(
\frac{Zα}{n - \abs{κ} + \sqrt{κ^2-Z^2α^2}}
\right)^2}}
\right]
\end{equation}
```

where the energy has been shifted by ``-c^2`` compared to
``\mathrm{(BS14.29)}``.

"""
function relativistic_energy(n, κ, Z=1)
    r = Z*α/(n-abs(κ)+√(κ^2-Z^2*α^2))
    -c^2*(1-inv(√(1+r^2)))
end

"""
    relativistic_energy(o::RelativisticOrbital)

Return the relativistic energy for the orbital `o`.

# Example

```jldoctest
julia> relativistic_energy(ro"1s")
-0.5000066565961734

julia> relativistic_energy(ro"2p-")
-0.12500208019062972

julia> relativistic_energy(ro"2p")
-0.1250004160284539
```
"""
relativistic_energy(o::RelativisticOrbital) = relativistic_energy(o.n, o.κ)
relativistic_energy(o::SpinOrbital) = relativistic_energy(o.orb)

"""
    orbital_energy(o)

Return the energy of orbital `o`, automatically choosing between
[`non_relativistic_energy`](@ref) and [`relativistic_energy`](@ref),
depending on the type of `o`.

# Examples

```jldoctest
julia> orbital_energy(o"1s")
-0.5

julia> orbital_energy(ro"1s")
-0.5000066565961734

julia> orbital_energy(o"2p")
-0.125

julia> orbital_energy(ro"2p")
-0.1250004160284539
```

"""
orbital_energy(o::Orbital) = non_relativistic_energy(o)
orbital_energy(o::RelativisticOrbital) = relativistic_energy(o)
orbital_energy(o::SpinOrbital) = orbital_energy(o.orb)

ground_state( ::Orbital) = o"1s"
ground_state( ::RelativisticOrbital) = ro"1s"
ground_state(o::SpinOrbital) = ground_state(o.orb)

"""
    atomic_hamiltonian(orbitals[; Efun=orbital_energy, neutral_zero=true])

Return the diagonal matrix whose elements corresponds to the energies
of each of the `orbitals`, by default [`orbital_energy`](@ref) is used
to determine the energy, which dispatches on the orbital type to
return the appropriate energy (non-relativistic or relativistic), and
all energies are shifted such that `1s` has zero energy, if
`neutral_zero==true`.

# Examples

```jldoctest
julia> atomic_hamiltonian(os"1[s] 2[s-p]")
3×3 LinearAlgebra.Diagonal{Float64,Array{Float64,1}}:
 0.0   ⋅      ⋅
  ⋅   0.375   ⋅
  ⋅    ⋅     0.375

julia> atomic_hamiltonian(ros"1[s] 2[s-p]")
4×4 LinearAlgebra.Diagonal{Float64,Array{Float64,1}}:
 0.0   ⋅         ⋅         ⋅
  ⋅   0.375005   ⋅         ⋅
  ⋅    ⋅        0.375005   ⋅
  ⋅    ⋅         ⋅        0.375006
```
"""
function atomic_hamiltonian(orbitals; Efun::Function=orbital_energy,
                            neutral_zero=true)
    E = zeros(length(orbitals))

    for (i,o) in enumerate(orbitals)
        E[i] = Efun(o)
    end

    if neutral_zero && !isempty(orbitals)
        E .-= Efun(ground_state(first(orbitals)))
    end

    Diagonal(E)
end

export non_relativistic_energy, relativistic_energy, orbital_energy, atomic_hamiltonian


text:
# ReverseDiff pullbacks for tracing AD

# using ChainRulesCore
# using ReverseDiff
# ChainRulesCore.@scalar_rule(ptransport_inner_coefficient(θ),
#              ptransport_inner_coefficient(θ)^2*sin(θ))
#
#
# ReverseDiff.@grad
# f(x::ReverseDiff.TrackedVector) = ReverseDiff.track(f, x)
# ReverseDiff.@grad function f(x)
#     xv = ReverseDiff.value(x)
#     return dot(xv, xv), Δ -> (Δ * 2 * xv,)
# end

# DiffCache to reduce cache craziness

# DiffCache - how do?
# immutable DiffCache{T<:AbstractArray, S<:AbstractArray}
#     du::T
#     dual_du::S
# end
#
#
# function DiffCache{chunk_size}(T, size, ::Type{Val{chunk_size}})
#     DiffCache(zeros(T, size...), zeros(ForwardDiff.Dual{nothing,T,chunk_size}, size...))
# end
#
# DiffCache(u::AbstractArray) = DiffCache(eltype(u),size(u),Val{ForwardDiff.pickchunksize(length(u))})
# DiffCache(u::AbstractArray, size) = DiffCache(eltype(u),size,Val{ForwardDiff.pickchunksize(length(u))})
# get_tmp{T<:ForwardDiff.Dual}(dc::DiffCache, ::Type{T}) = dc.dual_du
# get_tmp(dc::DiffCache, T) = dc.du

# function fused_update_strains_energy(r::basic_rod, Δr::rod_update, ϵr::rod_strains)
#
#     # Unpack
#     x0, d0 = r
#     Δx, Δθ = Δr
#     l0, κ0, τ0 = ϵr
#
#     # Tangents (Should maybe save these?)
#     t0 = x0[:,2:end] - x0[:,1:end-1]
#     t = t0 + Δx[:,2:end] - Δx[:,1:end-1]
#     t0 = t0./norm3(t0)
#     l = norm3(t)
#     t = t./l
#
#     # Update
#     d1 = ptransport(d0, t0, t)
#     rotate_orthogonal_unit!(d1, t, Δθ)
#     d2 = cross3(t, d1)
#
#     # Curvatures
#     k = 2*cross3(t[:,1:end-1], t[:,2:end])./(1 + dot3(t[:,1:end-1], t[:,2:end]))
#     κ = Matrix{Float64}(undef, 2, length(l)-1)
#     κ[1,:] = dot3(k, (d2[:,1:end-1] + d2[:,2:end])/2)
#     κ[2,:] = -dot3(k, (d1[:,1:end-1] + d1[:,2:end])/2)
#
#     # Twists (could move in part of ptransport to reduce intermediates)
#     τ = (τ0 + Δθ[2:end] - Δθ[1:end-1]
#             + spherical_excess(t0[:,1:end-1], t0[:,2:end], t[:,1:end-1])
#             + spherical_excess(t[:,1:end-1], t0[:,2:end], t[:,2:end]))
#
#     rod_strains(l, κ, τ)
# end


text:
"""
Mixed-fidelity unsteady aerodynamics simulation engine.

    # AUTHORSHIP
      * Author    : Eduardo J. Alvarez
      * Email     : Edo.AlvarezR@gmail.com
      * Created   : Oct 2019
      * License   : MIT
"""
module FLOWUnsteady

# ------------ GENERIC MODULES -------------------------------------------------
import Dierckx
import CSV
import DataFrames
import JLD
import Dates
using PyPlot
using LinearAlgebra: cross, I

# ------------ FLOW CODES ------------------------------------------------------
# FLOWVLM https://github.com/byuflowlab/FLOWVLM
import FLOWVLM
const vlm = FLOWVLM

# FLOWVPM https://github.com/byuflowlab/FLOWVPM.jl
try                     # Load FLOWVPM if available
    import FLOWVPM
catch e                 # Otherwise load a dummy version of FLOWVPM
    @warn("FLOWVPM module not found. Using dummy module instead.")
    include("FLOWUnsteady_dummy_FLOWVPM.jl")
end
const vpm = FLOWVPM

# GeometricTools https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
const gt = GeometricTools

import FLOWNoise
const noise = FLOWNoise

# BPM https://github.com/byuflowlab/BPM.jl
import BPM

# ------------ GLOBAL VARIABLES ------------------------------------------------
const module_path = splitdir(@__FILE__)[1]                # Path to this module
const def_data_path = joinpath(module_path, "../data/")   # Default path to data folder

# ------------ HEADERS ---------------------------------------------------------
# Load modules
for module_name in ["vehicle", "vehicle_vlm",
                    "maneuver", "rotor",
                    "simulation_types", "simulation", "utils",
                    "processing", "monitors",
                    "noise_wopwop", "noise_bpm"]
    include("FLOWUnsteady_"*module_name*".jl")
end

end # END OF MODULE

max_stars_repo_name: rachealerhard/FLOWUnsteady

max_issues_repo_name: rachealerhard/FLOWUnsteady

max_forks_repo_name: rachealerhard/FLOWUnsteady


text:
"""
    GridOperators

Provides local finite-difference
operator functions.
"""
module GridOperators

export ∇_z, Δ_z, adv_upwind, adv_upwind_conservative, upwind_conservative
using ..Grids

"""
    ∇_z(f, grid::Grid)

Computes the local derivative of field ``f``:
``∇f``
"""
∇_z(f, grid::Grid) = (f[3] - f[1])/(2*grid.dz)

"""
    Δ_z(f, grid::Grid)

Computes the local Laplacian of field ``f``:
``∇ • (∇f)``
"""
Δ_z(f, grid::Grid) = (f[3] - 2*f[2] + f[1])*grid.dzi2

"""
    Δ_z(f, grid::Grid, K)

Computes the local Laplacian of field ``f``
with a variable coefficient:
``∇ • (K ∇f)``
"""
function Δ_z(f, grid::Grid, K)
  K∇f_left  = (K[1]+K[2])/2*(f[2]-f[1])*grid.dzi
  K∇f_right = (K[2]+K[3])/2*(f[3]-f[2])*grid.dzi
  return (K∇f_right - K∇f_left)*grid.dzi
end

"""
    adv_upwind(ϕ, w, grid::Grid)

Local upwind advection operator ``w • ∇ϕ``. This
operator is stable but numerically diffusive.
"""
function adv_upwind(ϕ, w, grid::Grid)
  if w[2] > 0
    return w[ 2 ]*(ϕ[2] - ϕ[1]) * grid.dzi
  else
    return w[ 2 ]*(ϕ[3] - ϕ[2]) * grid.dzi
  end
end

"""
    adv_upwind_conservative(ϕ, w, grid::Grid)

Local upwind advection operator ``∇ • (wϕ)``. This
operator is stable but numerically diffusive.
"""
function adv_upwind_conservative(ϕ, w, grid::Grid)
  w_dual = [(w[1]+w[2])/2,
            (w[2]+w[3])/2]
  g_dual = [w_dual[1] > 0 ? ϕ[1] : ϕ[2],
            w_dual[2] > 0 ? ϕ[2] : ϕ[3]]
  return (w_dual[2]*g_dual[2] - w_dual[1]*g_dual[1]) * grid.dzi
end

"""
    upwind_conservative(ϕ, w, grid::Grid)

Local upwind gradient operator ``∇ϕ``. This
operator is stable but numerically diffusive.
"""
function upwind_conservative(ϕ, w, grid::Grid)
  w_dual = [(w[1]+w[2])/2,
            (w[2]+w[3])/2]
  g_dual = [w_dual[1] > 0 ? ϕ[1] : ϕ[2],
            w_dual[2] > 0 ? ϕ[2] : ϕ[3]]
  return (g_dual[2] - g_dual[1]) * grid.dzi
end

end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


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

"""
    module SimpleStateEvents

Model of a mass controlled by a two-point controller leading to state events.

# Modelica model

    f = if s > 0 then 0 else fMax;
    v = der(s);
    m*der(v) + d*v + k*s = f
"""
module SimpleStateEvents

import ModiaMath

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

    # Parameters
    m::Float64
    k::Float64
    d::Float64
    fMax::Float64

    function Model(;m=1.0, k=1.0, d=0.1, fMax=1.5, s0=2.0)
        @assert(m > 0.0)
        @assert(k > 0.0)
        @assert(d >= 0.0)
        @assert(fMax > 0.0)
        simulationState = ModiaMath.SimulationState("SimpleStateEvents", getModelResidues!, [s0;0.0], getVariableName;
                                nz=1, nw=1)
        new(simulationState, m, k, d, fMax)
    end
end

getVariableName(model, vcat, vindex) = ModiaMath.getVariableName(model, vcat, vindex;
                                                         xNames=["s", "v"],
                                                         wNames=["f"])

function getModelResidues!(m::Model, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    sim = m.simulationState

    s = _x[1]
    v = _x[2]

    sPos = ModiaMath.positive!(sim, 1, s, "s")
    f    = sPos ? 0.0 : m.fMax
    derv = (f - m.d * v - m.k * s) / m.m

    _r[1] = v    - _derx[1]
    _r[2] = derv - _derx[2]

    if ModiaMath.isStoreResult(m.simulationState)
        _w[1] = f
    end
    return nothing
end

end