text:
using Test, LinearAlgebra

@testset "/src/Degeneracy.jl" begin

using BellScenario: Degeneracy

@testset "symmetry_groups()" begin
    @testset "2x2 permutation group" begin
        transforms = [[1 0;0 1],[0 1;1 0]]
        elements = [[0 0],[1 0],[0 1],[1 1]]

        sym_groups = Degeneracy.symmetry_groups(transforms,elements)

        @test issetequal(elements, collect(Iterators.flatten(sym_groups)))
        @test length(sym_groups) == 3

        @test sym_groups[1] == [[0 0]]
        @test sym_groups[2] == [[1 0],[0 1]]
        @test sym_groups[3] == [[1 1]]
    end

    @testset "3x3 permutation group" begin
        transforms = Math.permutation_maps(3)
        elements = [
            [0 0 0],[0 0 1],[0 1 0],[0 1 1],
            [1 0 0],[1 0 1],[1 1 0],[1 1 1]
        ]

        sym_groups = Degeneracy.symmetry_groups(transforms,elements)

        @test issetequal(elements, collect(Iterators.flatten(sym_groups)))
        @test length(sym_groups) == 4

        @test sym_groups[1] == [[0 0 0]]
        @test sym_groups[2] == [[0 0 1],[0 1 0],[1 0 0]]
        @test sym_groups[3] == [[0 1 1],[1 0 1],[1 1 0]]
        @test sym_groups[4] == [[1 1 1]]
    end

    @testset "4x4 permutation group" begin
        transforms = Math.permutation_maps(4)
        elements = [
            [0 0 0 0],[0 0 0 1],[0 0 1 0],[0 0 1 1],
            [0 1 0 0],[0 1 0 1],[0 1 1 0],[0 1 1 1],
            [1 0 0 0],[1 0 0 1],[1 0 1 0],[1 0 1 1],
            [1 1 0 0],[1 1 0 1],[1 1 1 0],[1 1 1 1]
        ]

        sym_groups = Degeneracy.symmetry_groups(transforms,elements)

        @test issetequal(elements, collect(Iterators.flatten(sym_groups)))
        @test length(sym_groups) == 5

        @test sym_groups[1] == [[0 0 0 0]]
        @test sym_groups[2] == [[0 0 0 1],[0 0 1 0],[0 1 0 0],[1 0 0 0]]
        @test sym_groups[3] == [[0 0 1 1],[0 1 0 1],[0 1 1 0],[1 0 0 1],[1 0 1 0],[1 1 0 0]]
        @test sym_groups[4] == [[0 1 1 1],[1 0 1 1],[1 1 0 1],[1 1 1 0]]
        @test sym_groups[5] == [[1 1 1 1]]
    end
end

@testset "merge_symmetry_groups()" begin
    group1_sets = [[[1,1],[2,2]],[[3,3],[4,4],[5,5]],[[6,6]],[[7,7]]]
    group2_sets = [[[1,1],[3,3]],[[6,6],[7,7]],[[2,2],[4,4],[5,5]]]

    merged_groups = Degeneracy.merge_symmetry_groups(group1_sets,group2_sets)
    @test length(merged_groups) == 2
    @test all( # returned group sets do not overlap
        (x) -> x === nothing,
        indexin(merged_groups[1],merged_groups[2])
    )
    @test issetequal(merged_groups[1], [[1,1],[2,2],[3,3],[4,4],[5,5]])
    @test issetequal(merged_groups[2], [[6,6],[7,7]])
end

@testset "canonical_facets()" begin
    α_expt = (2,2)
    β_expt = (2,2)

    facets = [
        [0 0 0 0 0 -1 0 0 0],[0 0 0 0 0 0 -1 0 0],[0 0 0 0 0 0 0 -1 0],
        [0 0 0 0 0 0 0 0 -1],[0 -1 0 0 0 0 1 0 0],[0 -1 0 0 0 1 0 0 0],
        [0 0 -1 0 0 0 0 0 1],[0 0 -1 0 0 0 0 1 0],[0 0 0 -1 0 0 0 1 0],
        [0 0 0 -1 0 1 0 0 0],[0 0 0 0 -1 0 0 0 1],[0 0 0 0 -1 0 1 0 0],
        [0 -1 0 -1 0 1 1 1 -1],[0 -1 0 0 -1 1 1 -1 1],[0 0 -1 -1 0 1 -1 1 1],
        [0 0 -1 0 -1 -1 1 1 1],[-1 0 1 0 1 0 0 0 -1],[-1 0 1 1 0 0 0 -1 0],
        [-1 1 0 0 1 0 -1 0 0],[-1 1 0 1 0 -1 0 0 0],[-1 0 1 0 1 1 -1 -1 -1],
        [-1 0 1 1 0 -1 1 -1 -1],[-1 1 0 0 1 -1 -1 1 -1],[-1 1 0 1 0 -1 -1 -1 1]
    ]

    canonical_groups = Degeneracy.canonical_facets(α_expt, β_expt, facets, "no-signaling")

    @test issetequal(facets, collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 2
    @test issetequal(canonical_groups[1], [
        [0 0 0 0 0 -1 0 0 0],[0 0 0 0 0 0 -1 0 0],[0 0 0 0 0 0 0 -1 0],[0 0 0 0 0 0 0 0 -1],
        [0 -1 0 0 0 0 1 0 0],[0 -1 0 0 0 1 0 0 0],[0 0 -1 0 0 0 0 0 1],[0 0 -1 0 0 0 0 1 0],
        [0 0 0 -1 0 0 0 1 0],[0 0 0 -1 0 1 0 0 0],[0 0 0 0 -1 0 0 0 1],[0 0 0 0 -1 0 1 0 0],
        [-1 0 1 0 1 0 0 0 -1],[-1 0 1 1 0 0 0 -1 0],[-1 1 0 0 1 0 -1 0 0],[-1 1 0 1 0 -1 0 0 0]
    ])
    @test issetequal(canonical_groups[2], [
        [0 -1 0 -1 0 1 1 1 -1],[0 -1 0 0 -1 1 1 -1 1],[0 0 -1 -1 0 1 -1 1 1],[0 0 -1 0 -1 -1 1 1 1],
        [-1 0 1 0 1 1 -1 -1 -1],[-1 0 1 1 0 -1 1 -1 -1],[-1 1 0 0 1 -1 -1 1 -1],[-1 1 0 1 0 -1 -1 -1 1]
    ])
end

@testset "canonical_vertices()" begin
    @testset "chsh vertices" begin
        α1 = (2,2)
        β1 = (2,2)

        vertices = [
            [1;0;0;0;0;0;0;0;0],[1;0;0;0;1;0;0;0;0],[1;0;0;1;0;0;0;0;0],[1;0;0;1;1;0;0;0;0],
            [1;0;1;0;0;0;0;0;0],[1;0;1;0;1;0;0;0;1],[1;0;1;1;0;0;0;1;0],[1;0;1;1;1;0;0;1;1],
            [1;1;0;0;0;0;0;0;0],[1;1;0;0;1;0;1;0;0],[1;1;0;1;0;1;0;0;0],[1;1;0;1;1;1;1;0;0],
            [1;1;1;0;0;0;0;0;0],[1;1;1;0;1;0;1;0;1],[1;1;1;1;0;1;0;1;0],[1;1;1;1;1;1;1;1;1]
        ]

        canonical_vertices = Degeneracy.canonical_vertices(α1, β1, vertices, "no-signaling")

        @test length(canonical_vertices) == 1
        @test length(canonical_vertices[1]) == 16
        @test issetequal(canonical_vertices[1], vertices)
    end

    @testset "prepare and measure vertices" begin
        α = (3,1)
        β = (1,3)

        vertices = [
            [1 1 0 0 0 0 0]',[1 1 1 0 0 0 0]',[1 1 0 1 0 0 0]',[1 1 1 1 0 0 0]',[1 1 1 0 0 0 1]',[1 1 0 1 0 1 0]',
            [1 1 0 0 0 1 1]',[1 0 0 0 1 0 0]',[1 0 0 0 1 1 0]',[1 0 0 0 1 0 1]',[1 0 1 1 1 0 0]',[1 0 1 0 1 0 1]',
            [1 0 0 1 1 1 0]',[1 0 0 0 1 1 1]',[1 0 0 0 0 0 0]',[1 0 1 0 0 0 0]',[1 0 0 1 0 0 0]',[1 0 1 1 0 0 0]',
            [1 0 0 0 0 1 1]',[1 0 0 0 0 1 0]',[1 0 0 0 0 0 1]'
        ]

        canonical_vertices = Degeneracy.canonical_vertices(α, β, vertices, "normalized")

        @test length(canonical_vertices) == 2
        @test length(canonical_vertices[1]) == 3
        @test length(canonical_vertices[2]) == 18
        @test issetequal(canonical_vertices[1], [[1 1 1 1 0 0 0]',[1 0 0 0 0 0 0]',[1 0 0 0 1 1 1]'])
        @test issetequal( canonical_vertices[2],[
            [1 1 0 0 0 0 0]',[1 1 1 0 0 0 0]',[1 1 0 1 0 0 0]',[1 1 1 0 0 0 1]',[1 1 0 1 0 1 0]',
            [1 1 0 0 0 1 1]',[1 0 0 0 1 0 0]',[1 0 0 0 1 1 0]',[1 0 0 0 1 0 1]',[1 0 1 1 1 0 0]',
            [1 0 0 1 1 1 0]',[1 0 1 0 0 0 0]',[1 0 0 1 0 0 0]',[1 0 1 1 0 0 0]',[1 0 1 0 1 0 1]',
            [1 0 0 0 0 1 1]',[1 0 0 0 0 1 0]',[1 0 0 0 0 0 1]'
        ])
    end
end

@testset "bipartite_input_relabels()" begin

    @testset "no-signaling subspace" begin
        id = diagm(0 => fill(1, 9))
        flip_alice_input = [
            1 0 0 0 0 0 0 0 0;
            0 0 1 0 0 0 0 0 0;
            0 1 0 0 0 0 0 0 0;
            0 0 0 1 0 0 0 0 0;
            0 0 0 0 1 0 0 0 0;
            0 0 0 0 0 0 0 1 0;
            0 0 0 0 0 0 0 0 1;
            0 0 0 0 0 1 0 0 0;
            0 0 0 0 0 0 1 0 0;
        ];
        flip_bob_input = [
            1 0 0 0 0 0 0 0 0;
            0 1 0 0 0 0 0 0 0;
            0 0 1 0 0 0 0 0 0;
            0 0 0 0 1 0 0 0 0;
            0 0 0 1 0 0 0 0 0;
            0 0 0 0 0 0 1 0 0;
            0 0 0 0 0 1 0 0 0;
            0 0 0 0 0 0 0 0 1;
            0 0 0 0 0 0 0 1 0;
        ];
        flip_ab_input = [
            1 0 0 0 0 0 0 0 0;
            0 0 1 0 0 0 0 0 0;
            0 1 0 0 0 0 0 0 0;
            0 0 0 0 1 0 0 0 0;
            0 0 0 1 0 0 0 0 0;
            0 0 0 0 0 0 0 0 1;
            0 0 0 0 0 0 0 1 0;
            0 0 0 0 0 0 1 0 0;
            0 0 0 0 0 1 0 0 0;
        ]

        @test issetequal(
            Degeneracy.bipartite_input_relabels((2,2), (2,2), "no-signaling"),
            [id, flip_alice_input, flip_bob_input, flip_ab_input]
        )
    end

    test_cases = [
        ("no-signaling", 8),
        ("fixed-direction", 10),
        ("normalized", 12),
        ("generalized", 16)
    ]
    @testset "basic tests for subspace: $(case[1])" for case in test_cases
        maps = Degeneracy.bipartite_input_relabels((2,2),(2,2),case[1])

        dim = case[2]

        @test size(maps) == size(unique(maps))

        for i in 1:size(maps)[1]
            map = maps[i]

            @test size(map) == (dim+1,dim+1)
            @test map[1,1] == 1
            @test map[2:end,1] == zeros(Int64, dim)
            @test map[1,2:end] == zeros(Int64, dim)
        end
    end
end

@testset "bipartite_output_relabels()" begin
    @testset "dichotomic case" begin
        @testset "no-signaling" begin
            (α_relabels, β_relabels) = Degeneracy.bipartite_output_relabels((2,2),(2,2),"no-signaling")

            @test size(α_relabels) == (4,)
            @test size(β_relabels) == (4,)

            @test α_relabels[1] == diagm(0=>fill(1, 9))
            @test β_relabels[1] == diagm(0=>fill(1, 9))

            @test α_relabels[2] == [
                1 0 0 0 0 0  0 0 0;
                1 -1 0 0 0 0 0 0 0;
                0 0 1 0 0 0 0 0 0;
                0 0 0 1 0 0 0 0 0;
                0 0 0 0 1 0 0 0 0;
                0 0 0 1 0 -1 0 0 0;
                0 0 0 0 1 0 -1 0 0;
                0 0 0 0 0 0 0 1 0;
                0 0 0 0 0 0 0 0 1;
            ]
            @test α_relabels[3] == [
                 1 0 0 0 0 0 0 0 0;
                 0 1 0 0 0 0 0 0 0;
                 1 0 -1 0 0 0 0 0 0;
                 0 0 0 1 0 0 0 0 0;
                 0 0 0 0 1 0 0 0 0;
                 0 0 0 0 0 1 0 0 0;
                 0 0 0 0 0 0 1 0 0;
                 0 0 0 1 0 0 0 -1 0;
                 0 0 0 0 1 0 0 0 -1;
             ]
             @test β_relabels[2] == [
                1 0 0 0 0 0 0 0 0;
                0 1 0 0 0 0 0 0 0;
                0 0 1 0 0 0 0 0 0;
                1 0 0 -1 0 0 0 0 0;
                0 0 0 0 1 0 0 0 0;
                0 1 0 0 0 -1 0 0 0;
                0 0 0 0 0 0 1 0 0;
                0 0 1 0 0 0 0 -1 0;
                0 0 0 0 0 0 0 0 1;
            ]
            @test β_relabels[3] == [
                1 0 0 0 0 0 0 0 0;
                0 1 0 0 0 0 0 0 0;
                0 0 1 0 0 0 0 0 0;
                0 0 0 1 0 0 0 0 0;
                1 0 0 0 -1 0 0 0 0;
                0 0 0 0 0 1 0 0 0;
                0 1 0 0 0 0 -1 0 0;
                0 0 0 0 0 0 0 1 0;
                0 0 1 0 0 0 0 0 -1;
            ]
            @test β_relabels[4] == β_relabels[2]*β_relabels[3]
        end
    end

    @testset "relabel set contains the inverse of each element" begin
        test_cases = [
            ((2,2),(2,2)),
            ((3,2),(2,2)),
            ((2,3),(2,3)),
            ((3,1),(1,3))
        ]

        @testset "generalized $(case[1]), $(case[2])" for case in test_cases
            (α_relabels, β_relabels) = Degeneracy.bipartite_output_relabels(case[1],case[2],"generalized")

            dim = Behavior.dimension(case[1],case[2],"generalized")

            # in generalized coordinates, output relabelings inverted through
            # their own transpose.
            for relabel in α_relabels
                parity = det(relabel)
                @test (parity == 1) | (parity == -1)
                @test relabel*relabel' == diagm( 0 => fill(1,dim))
            end

            for relabel in β_relabels
                parity = det(relabel)
                @test (parity == 1) | (parity == -1)
                @test relabel*relabel' == diagm( 0 => fill(1,dim))
            end
        end

        @testset "no-signaling $(case[1]), $(case[2])" for case in test_cases
            (α_relabels, β_relabels) = Degeneracy.bipartite_output_relabels(case[1],case[2],"no-signaling")

            dim = Behavior.dimension(case[1],case[2],"no-signaling")
            Id = diagm(0 => fill(1,dim))

            foreach( relabel_test->begin
                matches = filter(relabel_inv->relabel_test*relabel_inv == Id, α_relabels)
                @test (1,) == size(matches)
            end, α_relabels)

            foreach( relabel_test->begin
                matches = filter(relabel_inv->relabel_test*relabel_inv == Id, β_relabels)
                @test (1,) == size(matches)
            end, β_relabels)
        end

        @testset "fixed-direction $(case[1]), $(case[2])" for case in test_cases
            (α_relabels, β_relabels) = Degeneracy.bipartite_output_relabels(case[1],case[2],"fixed-direction")

            dim = Behavior.dimension(case[1],case[2],"fixed-direction")
            Id = diagm(0 => fill(1,dim))

            # check that there is exactly 1 inverting map for every map in the
            # relabeling set
            foreach( relabel_test->begin
                matches = filter(relabel_inv->relabel_test*relabel_inv == Id, α_relabels)
                @test (1,) == size(matches)
            end, α_relabels)

            foreach( relabel_test->begin
                matches = filter(relabel_inv->relabel_test*relabel_inv == Id, β_relabels)
                @test (1,) == size(matches)
            end, β_relabels)
        end
    end
end

end


text:
# =============================================================================
# Cut chocolate bar
#      x1   x2   x3    n
#  ___  ___  ___  ___
# |___||___||___||___|  y1
# |___||___||___||___|  y2
# |___||___||___||___|  y3
# |___||___||___||___|  ml
#
# Luis Lopez Diaz - llopeza@uninorte.edu.co
# =============================================================================
using JuMP, GLPK, Gurobi;
# Input data
m, n = 6, 4;                # Chocolate bar dimension
x = [2 1 3];                # Cost of vertical cuts
y = [1 4 4 1 2];            # Cost of horizontal cuts
nx = length(x);             # Number of horizontal cuts
ny = length(y);             # Number of vertical cuts
nX = 1:nx;                  # Sets for vertical cuts
nY = 1:ny;                  # Sets for horizontal cuts
ml = Model(with_optimizer(Gurobi.Optimizer));
@variable(ml, 1 ≤ X[i in nX] ≤ m,Int);
@variable(ml, 1 ≤ Y[i in nY] ≤ n, Int);
@variable(ml, f);
@constraint(ml, sumX_l, nx .<= ones(1,nx)*X)
@constraint(ml, sumY_l, ny .<= ones(1,ny)*Y)
@constraint(ml, sumX_u, ones(1,nx)*X .<= nx*(ny+1))
@constraint(ml, sumY_u, ones(1,ny)*Y .<= (nx+1)*ny)
@constraint(ml, cuts, ones(1,ny)*Y + ones(1,nx)*X .== m*n-1)
@constraint(ml, con, f .== x*X + y*Y);
@objective(ml, Min, f)

optimize!(ml)
println(value.(X))
println(value.(Y))
println(value(f))

ext: jl


text:
# test to see if the likelihood is correctly calculated
# and if the networks are correctly estimated
# Claudia August 2015

# -------------------5taxon tree------------------

PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")

df=DataFrame(t1=["6","6","10","6","6"],
             t2=["7","7","7","10","7"],
             t3=["4","10","4","4","4"],
             t4=["8","8","8","8","10"],
             obsCF12=[0.2729102510259939, 0.3967750546426937, 0.30161247267865315, 0.24693940689390592, 0.2729102510259939],
             obsCF13=[0.45417949794801216, 0.30161247267865315, 0.30161247267865315, 0.5061211862121882, 0.45417949794801216],
             obsCF14=[0.2729102510259939, 0.30161247267865315, 0.3967750546426937, 0.24693940689390592, 0.2729102510259939])
d = readTableCF(df)
@test_throws ErrorException PhyloNetworks.writeExpCF(d)
@test writeTableCF(d) == rename(df, [:obsCF12 => :CF12_34, :obsCF13 => :CF13_24, :obsCF14 => :CF14_23])
@test tipLabels(d) ==  ["4","6","7","8","10"]
@test_logs PhyloNetworks.descData(d, devnull)

df[!,:ngenes] = [10,10,10,10,20]
allowmissing!(df, :ngenes)
d = readTableCF(df)
df[1,:ngenes] = missing; d.quartet[1].ngenes = -1.0
newdf = writeTableCF(d)
@test newdf[!,1:7] == rename(df, [:obsCF12 => :CF12_34, :obsCF13 => :CF13_24, :obsCF14 => :CF14_23])[!,1:7]
@test ismissing(newdf[1,:ngenes])
@test newdf[2:end,:ngenes] == df[2:end,:ngenes]

# starting tree:
tree = "((6,4),(7,8),10);"
currT = readTopologyLevel1(tree);
#printEdges(currT)

@testset "correct pseudo likelihood and snaq" begin
@testset "lik on tree" begin
extractQuartet!(currT,d)
calculateExpCFAll!(d)
tmp = (@test_logs PhyloNetworks.writeExpCF(d))
for i in [5,7] for j in 2:5 @test tmp[j,i] ≈ 0.12262648039048077; end end
for j in 2:5 @test tmp[j,6] ≈ 0.7547470392190385; end
lik = logPseudoLik(d)
@test lik ≈ 193.7812623319291
#estTree = optTopRun1!(currT,d,0,5454) # issue with printCounts, TravisCI?
#@test estTree.loglik ≈ 0.0 atol=1e-8
#println("passed optTopRun1! on tree")
end

# ------------------5taxon network 1 hybridization: Case H-----------------
# starting topology: Case G
global tree = "((((6:0.1,4:1.5)1:0.2,(7)11#H1)5:0.1,(11#H1,8)),10:0.1);" # Case G
global currT = readTopologyLevel1(tree);
# real network: Case H
global df=DataFrame(t1=["6","6","10","6","6"],t2=["7","7","7","10","7"],t3=["4","10","4","4","4"],t4=["8","8","8","8","10"],CF1234=[0.13002257237728915, 0.36936019721747243, 0.34692592933269173, 0.12051951084152591, 0.11095702789935982], CF1324=[0.7399548552454217, 0.28371387344983595, 0.28371387344983595, 0.7589609783169482, 0.7780859442012804],CF1423=[0.13002257237728915, 0.34692592933269173, 0.36936019721747243, 0.12051951084152591, 0.11095702789935982])
global d = readTableCF(df)

@testset "lik of network" begin
extractQuartet!(currT,d)
calculateExpCFAll!(d)
lik = logPseudoLik(d)
@test lik ≈ 50.17161079450669
end

@testset "network estimation h=1" begin
estNet = optTopRun1!(currT, 0.01,75, d,1, 1e-5,1e-6,1e-3,1e-4,
                     false,true,Int[], 5454, stdout,false,0.3)
# topology, likAbs,Nfail, data,hmax, fRel,fAbs,xRel,xAbs,
# verbose,closeN,numMoves, seed, logfile,writelog,probST,sout)
@test estNet.loglik < 0.00217
end

@testset "snaq! in serial and in parallel" begin
  global tree = readTopology("((((6:0.1,4:1.5),9)1:0.1,8),10:0.1);")
  @test_throws ErrorException snaq!(tree, d) # some taxa are in quartets, not in tree
  originalstdout = stdout
  redirect_stdout(devnull)
  global net = readTopology("((((6:0.1,4:1.5)1:0.2,((7,60))11#H1)5:0.1,(11#H1,8)),10:0.1);")
  @test_logs (:warn, r"^these taxa will be deleted") snaq!(net, d, # taxon "60" in net: not in quartets
    hmax=1, runs=1, Nfail=1, seed=1234, ftolRel=1e-2,ftolAbs=1e-2,xtolAbs=1e-2,xtolRel=1e-2)
  global n1 = snaq!(currT, d, hmax=1, runs=1, Nfail=1, seed=1234,
             ftolRel=1e-2,ftolAbs=1e-2,xtolAbs=1e-2,xtolRel=1e-2,
             verbose=true)
  addprocs(1)
  @everywhere using PhyloNetworks
  global n2 = snaq!(currT, d, hmax=1, runs=2, Nfail=1, seed=1234,
             ftolRel=1e-2,ftolAbs=1e-2,xtolAbs=1e-2,xtolRel=1e-2)
  redirect_stdout(originalstdout)
  rmprocs(workers())
  @test writeTopology(n1, round=true)==writeTopology(n2, round=true)
  @test n1.loglik == n2.loglik
  n3 = readSnaqNetwork("snaq.out")
  @test writeTopology(n3, round=true)==writeTopology(n2, round=true)
  @test n3.loglik > 0.0
  rm("snaq.out")
  rm("snaq.networks")
  rm("snaq.log") # .log and .err should be git-ignored, but still
  rm("snaq.err")
end
end


text:
using HDF5, JLD

examplefolder = dirname(@__FILE__)
datafolder = joinpath(examplefolder,"..","..","data","datasets","victoriapark")
include(joinpath(examplefolder,"VicPrkEstimator.jl"))

d = jldopen(joinpath(datafolder,"VicPrkBasic.jld"), "r") do file
  read(file, "dBasic")
end
f = jldopen(joinpath(datafolder,"VicPrkBasic.jld"), "r") do file
  read(file, "fBasic")
end
# MM = jldopen("test/data/VicPrkBasic.jld", "r") do file
#   read(file, "MMBasic")
# end
include(joinpath(datafolder,"Vic120MM.jl"))
MMr = MMwithRedirects(MM);

# isamdict = jldopen(joinpath(datafolder,"fgWithISAMrefERR01.jld"), "r") do file
#   read(file, "isamdict")
# end
# fgu = jldopen(joinpath(datafolder,"fgWithISAMrefERR01.jld"), "r") do file
#   read(file, "fgu")
# end


text:
function eccall(func_name::Symbol, libname::Symbol, rtype, args, types)
    :(ccall(($(QuoteNode(func_name)), $(QuoteNode(libname))),
        $rtype,
        $(Expr(:tuple, types...)),
        $(args...))
     )
end

function is_ptr_type_expr(@nospecialize t)
    (t === :Cstring || t === :Cwstring) && return true
    isa(t, Expr) || return false
    t = t::Expr
    t.head === :curly && t.args[1] === :Ptr
end

function efunsig(name::Symbol, args::Vector{Symbol}, types)
    x = [is_ptr_type_expr(t) ? a : Expr(:(::), a, t) for (a,t) in zip(args,types)]
    Expr(:call, name, x...)
end

function sort_cursors(cursors::Vector{CLMacroDefinition})
    cursor_names = spelling.(cursors)
    cursor_dict = Dict(zip(cursor_names, cursors))
    sort!(unique!(cursor_names))
    sorted_cursors = [cursor_dict[cursor_names[j]] for j in 1:length(cursor_names)]
    return sorted_cursors
end

function wrap_vx_fun!(ctx::AbstractContext, cursor::CLFunctionDecl)

    func_type = type(cursor)

    if kind(func_type) == CXType_FunctionNoProto
        @warn "No Prototype for $cursor - assuming no arguments"
    elseif isvariadic(func_type)
        @warn "Skipping VarArg Function $cursor"
        return ctx
    end

    func_name = isempty(ctx.force_name) ? Symbol(spelling(cursor)) : ctx.force_name
    ret_type = clang2julia(return_type(cursor)) # reformat name
    args = function_args(cursor)
    arg_types = [argtype(func_type, i) for i in 0:length(args)-1]

    # do a check to dealias/reformat args

    arg_reps = clang2julia.(arg_types)

    for (i, arg) in enumerate(arg_reps)
        if Meta.isexpr(arg, :curly) && first(arg.args) == :NTuple
            arg_reps[i] = Expr(:curly, :Ref, last(arg.args))
        end
        if Meta.isexpr(arg, :curly) && first(arg.args) == :Ptr
            arg_reps[i] = Expr(:curly, :Ref, last(arg.args))
        end
    end

    arg_count = 0
    arg_names = map(args) do x
        n = name_safe(name(x))
        s = !isempty(n) ? n : "arg"*string(arg_count+=1)
        Symbol(s)
    end

    signature = efunsig(Symbol(spelling(cursor)), arg_names, arg_reps)
    body = eccall(func_name, Symbol("libopenvx"), ret_type, arg_names, arg_reps)
    push!(ctx.api_buffer, Expr(:function, signature, Expr(:block, body)))

    return ctx
end

function wrap_vx_struct!(ctx::AbstractContext, cursor::CLStructDecl)
    # make sure a empty struct is indeed an opaque struct typedef/typealias
    # cursor = canonical(cursor)  # this won't work
    cursor = type(cursor) |> canonical |> typedecl
    cursor_name = name(cursor)
    # handle typedef anonymous struct
    idx = ctx.children_index
    if 0 < idx < length(ctx.children)
        next_cursor = ctx.children[idx+1]
        if is_typedef_anon(cursor, next_cursor)
            cursor_name = name(next_cursor)
        end
    end
    !isempty(ctx.force_name) && (cursor_name = ctx.force_name;)
    cursor_name == "" && (@warn("Skipping unnamed StructDecl: $cursor"); return ctx)

    struct_sym = symbol_safe(cursor_name)
    ismutable = get(ctx.options, "is_struct_mutable", false)
    buffer = ctx.common_buffer

    # generate struct declaration
    block = Expr(:block)
    expr = Expr(:struct, ismutable, struct_sym, block)
    deps = OrderedSet{Symbol}()
    struct_fields = children(cursor)
    for (field_idx, field_cursor) in enumerate(struct_fields)
        field_name = name(field_cursor)
        field_kind = kind(field_cursor)
        if field_kind == CXCursor_StructDecl || field_kind == CXCursor_UnionDecl
            continue
        elseif field_kind == CXCursor_FirstAttr && field_idx == 1
            continue
        elseif field_kind != CXCursor_FieldDecl || field_kind == CXCursor_TypeRef || kind(type(field_cursor)) == CXType_FunctionNoProto
            buffer[struct_sym] = ExprUnit(Clang.Poisoned())
            @warn "Skipping struct: \"$cursor\" due to unsupported field: $field_cursor"
            return ctx
        elseif isempty(field_name)
            error("Unnamed struct member in: $cursor ... cursor: $field_cursor")
        end

        if occursin("anonymous", string(clang2julia(field_cursor)))
            idx = field_idx-1
            anonymous_record = struct_fields[idx]
            while idx != 0 && kind(anonymous_record) == CXCursor_FieldDecl
                idx -= 1
                anonymous_record = struct_fields[idx]
            end
            if idx == field_idx-1
                ctx.anonymous_counter += 1
                anon_name = "ANONYMOUS$(ctx.anonymous_counter)_"*spelling(field_cursor)
                ctx.force_name = anon_name
                wrap!(ctx, anonymous_record)
                ctx.force_name = ""
                repr = symbol_safe(anon_name)
            else
                anon_name = "ANONYMOUS$(ctx.anonymous_counter)_"*spelling(struct_fields[idx+1])
                repr = symbol_safe(anon_name)
            end
        else
            repr = clang2julia(field_cursor)
        end
        push!(block.args, Expr(:(::), symbol_safe(field_name), repr))
        push!(deps, target_type(repr))
    end

    # check for a previous forward ordering
    if !(struct_sym in keys(buffer)) || buffer[struct_sym].state == :empty
        if !isempty(struct_fields)
            buffer[struct_sym] = ExprUnit(expr, deps)
        else
            # opaque struct typedef/typealias
            buffer[struct_sym] = ExprUnit(:(const $struct_sym = Cvoid), deps, state=:empty)
        end
    end

    return ctx
end


text:
@testset "Original tests" begin

function test_order(lru, keys, vals)
    for (i, (k,v)) = enumerate(lru)
        @test k == keys[i]
        @test v == vals[i]
    end
end

CACHE = LRU{Int, Int}(; maxsize = 20)
# Test insertion ordering
kvs = 1:10
for i in reverse(kvs)
    CACHE[i] = i
end
test_order(CACHE, 1:10, 1:10)

# Test reinsertion ordering
kvs = [1, 3, 5, 7, 9, 2, 4, 6, 8, 10]
for i in reverse(kvs)
    CACHE[i] = i
end
test_order(CACHE, kvs, kvs)

# Test least recently used items are evicted
resize!(CACHE; maxsize = 5)
test_order(CACHE, kvs[1:5], kvs[1:5])

resize!(CACHE; maxsize = 10)
test_order(CACHE, kvs[1:5], kvs[1:5])

kvs = 1:11
for i in reverse(kvs)
    CACHE[i] = i
end
test_order(CACHE, kvs[1:end-1], kvs[1:end-1])

# Test lookups, and that lookups reorder
kvs = 1:10
for i in reverse(kvs)
    CACHE[i] = 2*i
end
for i in kvs
    @test CACHE[i] == 2*i
end
test_order(CACHE, reverse(kvs), reverse(2*kvs))

# Test empty!
empty!(CACHE)
@test length(CACHE) == 0

# Test get! with default parameter
@test get!(CACHE, 1, 2)  == 2    # not in cache
@test get!(CACHE, 1, 2)  == 2    # in cache

# Test get! with function call
val = get!(CACHE, 2) do
          3
      end
@test val == 3
get!(CACHE, 2) do
    error("this shouldn't have been called!")
end

# Test Abstract typed cache. All we're checking for here is that the container
# is able to hold abstract types without issue. Insertion order is already
# tested above.
CACHE2 = LRU{String, Integer}(; maxsize = 5)
CACHE2["test"] = 4
@test CACHE2["test"] == 4

end


text:
export Container
abstract AbstractContainer

type Container <: AbstractContainer
  image::AbstractImage
  volumes::Dict{String, String}
  workdir::String
  name::String
  interactive::Bool
  rm::Bool
  env::Dict{Symbol, String}
end

const DEFAULT_CONTAINER = convert(Dict{Symbol, Any}, {
  :workdir => "",
  :name => "",
  :image => Image("ubuntu:14.04"),
  :volumes => Dict{String, String}(),
  :interactive => false,
  :rm => true,
  :env => Dict{Symbol, String}()
})

function Container(image::AbstractImage; kwargs...)
  # Treat volume separately
  volumes = Util.getarg(:volumes, kwargs, DEFAULT_CONTAINER)
  for (vol, paths) in filter(x -> x[1] == :volume, kwargs)
    if length(paths) != 2 || isa(paths, String)
      error("volume keyword argument should consist of (vm, container) path tuples")
    end
    volumes[string(paths[2])] = string(paths[1])
  end
  filter!(x -> x[1] != :volume, kwargs)

  # Remove special keyword arg volumes
  if length(kwargs) !=(map(x->x[1], kwargs) |> unique |> length)
    error("Keyword arguments, other than volume should be unique")
  end

  Container(
    image, volumes,
    map(x -> Util.getarg(x, kwargs, DEFAULT_CONTAINER), names(Container)[3:end])...
  )
end
Container(image::String; kwargs...) = Container(Image(image); kwargs...)
Container(; kwargs...) = Container(deepcopy(DEFAULT_CONTAINER[:image]); kwargs...)

function Container(container::Container; kwargs...)
  volumes = map(x -> x[2], filter(x -> x[1] == :volume, kwargs))
  filter!(x -> x[1] != :volume, kwargs)

  result = deepcopy(container)
  for (k, v) in kwargs
    setfield!(result, k, v)
  end
  for paths in volumes
    if length(paths) != 2 || isa(paths, String)
      error("volume keyword argument should consist of (vm, container) path tuples")
    end
    result.volumes[string(paths[2])] = string(paths[1])
  end

  result
end

function command(machine::MachineEnv, container::Container, cmd::Cmd=``)
  volargs = ""
  for (cont, host) in container.volumes
    volargs *= "-v $(chomp(abspath(host))):$(chomp(cont)):rw "
  end
  intargs = if container.interactive; "-i" else "" end
  rmargs = if container.rm "--rm=true" else "" end
  nmargs = if length(container.name) > 0 "--name=\"$(container.name)\"" else "" end

  envs = join(["-e $(string(k))=\"$(string(v))\"" for (k, v) in container.env], " ")
  if !is(cmd.env, nothing)
    envs = envs * " " * join(["-e $u" for u in cmd.env], " ")
  end

  if length(cmd.dir) > 0
    wkd = "--workdir=\"$(cmd.dir)\""
  elseif length(container.workdir) > 0
    wkd = "--workdir=\"$(container.workdir)\""
  else
    wkd = ""
  end
  image = container.image.name
  result = command(machine, "run", "$rmargs $envs $volargs $intargs $nmargs $wkd $image")
  `$result $cmd`
end

|>(machine::MachineEnv, container::AbstractContainer) = (machine, container)
|>(env::(MachineEnv, AbstractContainer), cmd::Union(Cmd, Vector{Cmd})) = (env[1], env[2], cmd)
function |> (image::AbstractImage, container::AbstractContainer)
  result = deepcopy(container)
  result.image = image
  result
end
|> (env::(MachineEnv, AbstractImage), cont::AbstractContainer) = (env[1], env[2] |> cont)
function |>(env::(MachineEnv, AbstractImage), cmd::Union(Cmd, Vector{Cmd}))
  env[1], Container(env[2]), cmd
end
|> (env::MachineEnv, cmd::Union(Cmd, Vector{Cmd})) = (env, Container(), cmd)
|> (m::MachineEnv, c::(Container, Union(Cmd, Vector{Cmd}))) = (m, c[1], c[2])

function command_impl(func::Function, vm::MachineEnv, container::Container, cmd::Cmd)
  if isa(container.image, BuildImage)
    vm |> container.image |> run
  end
  command(vm, container, cmd) |> func
end

for runner in [:run, :readchomp, :readall]
  @eval begin
    $runner(cmd::(MachineEnv, Container, Cmd)) = command_impl($runner, cmd[1], cmd[2], cmd[3])
    $runner(cmd::(MachineEnv, Container)) = command_impl($runner, cmd[1], cmd[2], ``)
    $runner(vm::MachineEnv, cont::Container, cmd::Cmd=``) = command_impl($runner, vm, cont, cmd)
  end
end

function bash(cmd::Cmd)
  result = ""
  # Adds environment variable
  if !is(cmd.env, nothing)
    result *= join(cmd.env, "\n")
  end
  # Move to workdir
  if length(cmd.dir) > 0
    result *= "cd $(cmd.dir)\n"
  end
  # Add command
  result *= "$(`$cmd`)"[2:end-1] * "\n"
  result
end
function run(vm::MachineEnv, container::Container, cmds::Vector{Cmd})
  if isa(container.image, BuildImage)
    vm |> container.image |> run
  end
  container = deepcopy(container)
  container.interactive = true
  open(command(vm, container), "w", STDOUT) do stream
    for cmd in cmds
      write(stream, bash(cmd))
    end
  end
end
run(cmd::(MachineEnv, Container, Vector{Cmd})) = run(cmd[1], cmd[2], cmd[3])


text:
module AlgebraOfGraphics

using Tables: columns, getcolumn, columnnames, AbstractColumns, Tables
using CategoricalArrays: categorical, cut, levelcode, refs, levels, CategoricalArray
using Observables: AbstractObservable, Observable, to_value
using NamedDims: NamedDimsArray, dimnames
using StructArrays: GroupPerm, refine_perm!, StructArray
using KernelDensity: kde
import GLM, Loess

using AbstractPlotting: Point2f0,
                        Attributes,
                        AbstractPlot,
                        Node,
                        automatic,
                        lift,
                        @lift,
                        RGBAf0,
                        AbstractPlotting

using AbstractPlotting.MakieLayout: LAxis,
                                    LText,
                                    LRect,
                                    GridLayout,
                                    linkxaxes!,
                                    linkyaxes!,
                                    hidexdecorations!,
                                    hideydecorations!,
                                    Top,
                                    Bottom,
                                    Left,
                                    Right,
                                    MakieLayout

export data, dims, draw, spec, style
export categorical, cut

include("algebraic_list.jl")
include("utils.jl")
include("context.jl")
include("specs.jl")
include("scales.jl")
include("analysis/analysis.jl")
include("analysis/smooth.jl")
include("analysis/density.jl")
include("legend.jl")
include("draw.jl")

end # module


text:
@testset "Fifth order" begin
    @testset "Bits $i" for i in 11:12
        x = 0:1 << i
        y_approx = fpcos.(Int32.(x), Val(i))
        y = cos.(x / 1 << i * π / 2) .* 1 << i
        @test all(abs.(y_approx .- y) .< 2.0)
    end
end


text:
using Documenter
using Unfold
using DocStringExtensions
using Plots
gr()
#unicodeplots()

makedocs(sitename="Unfold.jl",
        root = joinpath(dirname(pathof(Unfold)), "..", "docs"),
        prettyurls = get(ENV, "CI", nothing) == "true",
        pages = [
            "index.md",
            "LM Tutorial" =>"lm_tutorial.md",
            "LMM Tutorial" =>"lmm_tutorial.md",
        ])

deploydocs(; repo = "github.com/unfoldtoolbox/Unfold.jl", push_preview = true,        devbranch = "main")


text:
using SMLMData
using NearestNeighbors

"""
    initialdensity = estimatedensities(smld::SMLMData.SMLD2D,
        clusterdata::Vector{Matrix{Float32}}, params::ParamStruct)

Estimate local emitter densities for clusters in `smld` and `clusterdata`.

# Description
The initial local densities `initialdensity` around each pre-cluster present in
`smld`/`clusterdata` are estimated based on the local density of pre-clusters
throughout the entire set of data as well as some of the rate parameters
provided in `params`.
"""
function estimatedensities(smld::SMLMData.SMLD2D,
    clusterdata::Vector{Matrix{Float32}},
    params::ParamStruct)

    # Define some new parameters.
    dutycycle = params.k_on / (params.k_on + params.k_off + params.k_bleach)
    nclusters = length(clusterdata)
    if isempty(smld.nframes)
        maxframe = maximum(smld.framenum)
    else
        maxframe = smld.nframes
    end

    # If only one cluster is present, we should return an answer right away and stop.
    if nclusters == 1
        if size(clusterdata[1], 1) == 1
            initialdensity = 1
        else
            clusterarea = (maximum(clusterdata[1][:, 1]) - minimum(clusterdata[1][:, 1])) *
                          (maximum(clusterdata[1][:, 1]) - minimum(clusterdata[1][:, 1]))
            initialdensity = (1 / clusterarea) *
                             ((params.k_bleach / params.k_off) / (1 - params.p_miss)) /
                             (1 - exp(-params.k_bleach * dutycycle * (maxframe - 1)))
        end
        return initialdensity
    end

    # Determine the center of all clusters assuming each arose from the same
    # emitter.
    clustercenters = zeros(2, nclusters)
    for nn = 1:nclusters
        # Isolate some arrays to improve readability.
        x = clusterdata[nn][:, 1]
        y = clusterdata[nn][:, 2]
        σ_x = clusterdata[nn][:, 3]
        σ_y = clusterdata[nn][:, 4]

        # Compute the cluster centers based on MLE of position.
        clustercenters[1:2, nn] = [sum(x ./ σ_x .^ 2) / sum(1 ./ σ_x .^ 2)
            sum(y ./ σ_y .^ 2) / sum(1 ./ σ_y .^ 2)]
    end

    # Estimate the local cluster density based on the distance to
    # nearest-neighbors.
    kneighbors = minimum([params.nnearestclusters nclusters - 1])
    kdtree = NearestNeighbors.KDTree(clustercenters)
    _, nndist = NearestNeighbors.knn(kdtree, clustercenters, kneighbors + 1)
    nndist = getindex.(nndist, 1:(kneighbors-1))
    clusterdensity = (kneighbors + 1) ./ (pi * nndist .^ 2)

    # Estimate the density of underlying emitters based on cluster density.
    lambda1 = params.k_bleach * dutycycle
    lambda2 = (params.k_on + params.k_off + params.k_bleach) - lambda1
    initialdensity = clusterdensity *
                     (1.0 / dutycycle) * (1.0 / params.k_off) * (1.0 / (1.0 - params.p_miss)) ./
                     ((1.0 / lambda1) * (1.0 - exp(-lambda1 * (maxframe - 1.0))) -
                      (1.0 / lambda2) * (1.0 - exp(-lambda2 * (maxframe - 1.0))))

    return initialdensity
end


text:
import BangBang
import OnlineStats
using Documenter
using Transducers
# Not importing `JSON` and `Literate` here so that this file can be
# included from test without adding them as dependencies.

EXAMPLE_PAGES = [
    "Tutorial: Missing values" => "examples/tutorial_missings.md",
    "Tutorial: Parallelism" => "examples/tutorial_parallel.md",
    "Empty result handling" => "examples/empty_result_handling.md",
    "Parallel word count" => "examples/words.md",
    "Writing transducers" => "examples/transducers.md",
    "Writing reducibles" => "examples/reducibles.md",
]

transducers_rm_examples() =
    rm(joinpath(@__DIR__, "src/examples"); force=true, recursive=true)

function transducers_literate(;
        inputbase = joinpath(@__DIR__, "..", "examples"),
        outputbase = joinpath(@__DIR__, "src", "examples"),
        examples = EXAMPLE_PAGES,
        kwargs...)
    transducers_rm_examples()
    for (_, outpath) in examples
        name, = splitext(basename(outpath))
        inputfile = joinpath(inputbase, "$name.jl")
        outputdir = outputbase
        Literate.markdown(
            inputfile, outputdir;
            documenter = true,
            kwargs...)
    end
end

function transducers_makedocs(;
        examples = EXAMPLE_PAGES,
        kwargs...)
    if isempty(examples)
        # Make some dummy examples
        dummypage = joinpath("examples", "dummy.md")
        path = joinpath(@__DIR__, "src", dummypage)
        mkpath(dirname(path))
        write(
            path,
            """
            # Dummy page
            ## Empty result handling
            """
        )
        examples = dummypage
    end
    makedocs(;
        modules = [Transducers],
        pages = [
            "Home" => "index.md",
            "Parallelism" => "parallelism.md",
            "Manual" => "manual.md",
            "Interface" => "interface.md",
            "Examples" => examples,
            hide("Internals" => "internals.md"),
            hide("doctests/show_xf.md"),
            hide("doctests/show_rf.md"),
        ],
        repo = "https://github.com/tkf/Transducers.jl/blob/{commit}{path}#L{line}",
        sitename = "Transducers.jl",
        authors = "Takafumi Arakaki",
        root = @__DIR__,
        strict = true,
        kwargs...)
end

function should_push_preview(event_path = get(ENV, "GITHUB_EVENT_PATH", nothing))
    event_path === nothing && return false
    event = JSON.parsefile(event_path)
    pull_request = get(event, "pull_request", nothing)
    pull_request === nothing && return false
    labels = [x["name"] for x in pull_request["labels"]]
    # https://developer.github.com/v3/activity/events/types/#pullrequestevent
    yes = "push_preview" in labels
    if yes
        @info "Trying to push preview as label `push_preview` is specified." labels
    else
        @info "Not pushing preview as label `push_preview` is not specified." labels
    end
    return yes
end


text:
_view(::Nothing, i) = nothing
_view(A::Fill{T,2,Axes}, i) where {T,Axes} = view(A, :, 1)
_view(A::AbstractMatrix, idx) = view(A, :, idx)

aggregate(aggr::typeof(+), X) = vec(sum(X, dims=2))
aggregate(aggr::typeof(-), X) = -vec(sum(X, dims=2))
aggregate(aggr::typeof(*), X) = vec(prod(X, dims=2))
aggregate(aggr::typeof(/), X) = 1 ./ vec(prod(X, dims=2))
aggregate(aggr::typeof(max), X) = vec(maximum(X, dims=2))
aggregate(aggr::typeof(min), X) = vec(minimum(X, dims=2))
aggregate(aggr::typeof(mean), X) = vec(aggr(X, dims=2))

abstract type GraphNet <: AbstractGraphLayer end

@inline update_edge(gn::GraphNet, e, vi, vj, u) = e
@inline update_vertex(gn::GraphNet, ē, vi, u) = vi
@inline update_global(gn::GraphNet, ē, v̄, u) = u

@inline function update_batch_edge(gn::GraphNet, adj, E, V, u)
    n = size(adj, 1)
    edge_idx = edge_index_table(adj)
    mapreduce(i -> apply_batch_message(gn, i, adj[i], edge_idx, E, V, u), hcat, 1:n)
end

@inline apply_batch_message(gn::GraphNet, i, js, edge_idx, E, V, u) =
    mapreduce(j -> update_edge(gn, _view(E, edge_idx[(i,j)]), _view(V, i), _view(V, j), u), hcat, js)

@inline update_batch_vertex(gn::GraphNet, Ē, V, u) =
    mapreduce(i -> update_vertex(gn, _view(Ē, i), _view(V, i), u), hcat, 1:size(V,2))

@inline function aggregate_neighbors(gn::GraphNet, aggr, E, accu_edge)
    @assert !iszero(accu_edge) "accumulated edge must not be zero."
    cluster = generate_cluster(E, accu_edge)
    NNlib.scatter(aggr, E, cluster)
end

@inline function aggregate_neighbors(gn::GraphNet, aggr::Nothing, E, accu_edge)
    @nospecialize E accu_edge
    return nothing
end

@inline aggregate_edges(gn::GraphNet, aggr, E) = aggregate(aggr, E)

@inline function aggregate_edges(gn::GraphNet, aggr::Nothing, E)
    @nospecialize E
    return nothing
end

@inline aggregate_vertices(gn::GraphNet, aggr, V) = aggregate(aggr, V)

@inline function aggregate_vertices(gn::GraphNet, aggr::Nothing, V)
    @nospecialize V
    return nothing
end

function propagate(gn::GraphNet, fg::FeaturedGraph, naggr=nothing, eaggr=nothing, vaggr=nothing)
    E, V, u = propagate(gn, adjacency_list(fg), fg.ef, fg.nf, fg.gf, naggr, eaggr, vaggr)
    FeaturedGraph(fg, nf=V, ef=E, gf=u)
end

function propagate(gn::GraphNet, adj::AbstractVector{S}, E::R, V::Q, u::P,
                   naggr=nothing, eaggr=nothing, vaggr=nothing) where {S<:AbstractVector,R,Q,P}
    E = update_batch_edge(gn, adj, E, V, u)
    Ē = aggregate_neighbors(gn, naggr, E, accumulated_edges(adj))
    V = update_batch_vertex(gn, Ē, V, u)
    ē = aggregate_edges(gn, eaggr, E)
    v̄ = aggregate_vertices(gn, vaggr, V)
    u = update_global(gn, ē, v̄, u)
    return E, V, u
end


text:
using Plots
pyplot()

function plot_pieces(pieces...; gap=false)
    plt = plot(pieces[1], legend=false)
    offset = length(pieces[1])
    last_y = pieces[1][end]

    for piece in pieces[2:end]
        x = range(gap ? offset + 1 : offset, stop = offset + length(piece))
        y = gap ? piece : vcat(last_y, piece)
        plot!(plt, x, y)
    end

    plt
end


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-
#
# Force field transform
#
# Aaron LI
# 2015/07/14
#

using FITSIO;
#include("../julia/ndgrid.jl");

@everywhere function meshgrid(vx, vy)
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

# Calculate the forces between the specified point with respect to the image.
@everywhere function force(p0, img)
    img = copy(img);
    x0, y0 = p0;
    v0 = img[y0, x0];
    img[y0, x0] = 0.0;
    rows, cols = size(img);
    x, y = meshgrid(1:cols, 1:rows);
    x[y0, x0] = -1;
    y[y0, x0] = -1;
    f_x = v0 .* img .* (x-x0) ./ ((x-x0).^2 + (y-y0).^2).^1.5;
    f_y = v0 .* img .* (y-y0) ./ ((x-x0).^2 + (y-y0).^2).^1.5;
    #return (f_x, f_y);
    return (sum(f_x), sum(f_y));
end

# Perform the "force field transform" for the input image.
#
# Return:
#   (amplitudes, angles)
#   amplitudes: the amplitudes of the resulting forces of each pixel
#   angles: the directions of the resulting forces of each pixel,
#           in unit radian.
@everywhere function force_field_transform_serial(img, rowstart=1, rowend="end")
    rows, cols = size(img)
    if rowend == "end"
        rowend = rows
    end
    amplitudes = zeros(rows, cols)
    angles = zeros(rows, cols)
    t0 = time()
    t_p = t0 + 30  # in 30 seconds
    for y = rowstart:rowend
        for x = 1:cols
            t1 = time()
            if (t1 >= t_p)
                percent = 100*((y-rowstart)*cols + x+1) / ((rowend-rowstart+1)*cols)
                @printf("Worker #%d: progress: %.3f%%; %.1f min\n",
                        myid(), percent, (t1-t0)/60.0)
                t_p += 30  # in 30 seconds
            end
            F_x, F_y = force((x, y), img)
            #@printf("F_x, F_y = (%f, %f)\n", F_x, F_y);
            amplitudes[y, x] = sqrt(F_x^2 + F_y^2)
            angles[y, x] = atan2(F_y, F_x)
        end
    end
    t1 = time()
    @printf("Worker #%d: finished in %.1f min!\n", myid(), (t1-t0)/60.0)
    return (amplitudes, angles)
end

# parallel-capable
function force_field_transform(img)
    t0 = time()
    rows, cols = size(img)
    np = nprocs()
    amplitudes = cell(np)
    angles = cell(np)
    # split rows for each process
    rows_chunk = div(rows, np)
    rowstart = cell(np)
    rowend = cell(np)
    @sync begin
        for p = 1:np
            rowstart[p] = 1 + rows_chunk * (p-1)
            if p == np
                rowend[p] = rows
            else
                rowend[p] = rowstart[p] + rows_chunk - 1
            end
            # perform transform
            @async begin
                amplitudes[p], angles[p] = remotecall_fetch(p,
                        force_field_transform_serial,
                        img, rowstart[p], rowend[p])
            end
        end
    end
    t1 = time()
    @printf("Finished in %.1f min!\n", (t1-t0)/60.0)
    return (sum(amplitudes), sum(angles))
end

# arguments
#println(ARGS);
if length(ARGS) != 3
    println("Usage: PROG <input_fits_img> <out_fits_amplitudes> <out_fits_angles>");
    exit(1);
end

infile = ARGS[1];
outfile_ampl = ARGS[2];
outfile_angles = ARGS[3];

fits_img = FITS(infile);
img = read(fits_img[1]);
header = read_header(fits_img[1]);

# perform force field transform
ampl, angles = force_field_transform(img);

outfits_ampl = FITS(outfile_ampl, "w");
outfits_angles = FITS(outfile_angles, "w");
write(outfits_ampl, ampl; header=header);
write(outfits_angles, angles; header=header);

close(fits_img);
close(outfits_ampl);
close(outfits_angles);

#= vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=julia: =#

max_stars_repo_name: liweitianux/atoolbox

max_issues_repo_name: liweitianux/atoolbox

max_forks_repo_name: liweitianux/atoolbox


text:
# # Wrapper blocks

abstract type WrapperBlock <: AbstractBlock end

wrapped(w::WrapperBlock) = w.block
wrapped(b::Block) = b
function setwrapped(w::WrapperBlock, b)
    return Setfield.@set w.block = b
end
mockblock(w::WrapperBlock) = mockblock(wrapped(w))
checkblock(w::WrapperBlock, obs) = checkblock(wrapped(w), obs)

# If not overwritten, encodings are applied to the wrapped block
"""
    abstract type PropagateWrapper

Defines the default propagation behavior of a `WrapperBlock` when
an encoding is applied to it.

Propagation refers to what happens when an encoding is applied to
a `WrapperBlock`. If no `encode` method is defined for a wrapper block
`wrapper`, `encode` is instead called on the wrapped block.
Propagating the wrapper block means that the block resulting from
encoding the wrapped block is rewrapped in `wrapper.`.

```
wrapper = Wrapper(block)
# propagate
encodedblock(enc, wrapper) = Wrapper(encodedblock(enc, wrapped(wrapper)))

# don't propagate
encodedblock(enc, wrapper) = encodedblock(enc, wrapped(wrapper))
```

The following wrapping behaviors exist:

- `PropagateAlways`: Always propagate. This is the default behavior.
- `PropagateNever`: Never propagate
- `PropagateSameBlock`: Only propagate if the wrapped block is unchanged
    by the encoding
"""
abstract type PropagateWrapper end

struct PropagateAlways <: PropagateWrapper end
struct PropagateSameBlock <: PropagateWrapper end
struct PropagateNever <: PropagateWrapper end

propagatewrapper(::WrapperBlock) = PropagateAlways()

encodedblock(enc::Encoding, wrapper::WrapperBlock) =
    encodedblock(enc, wrapper, propagatewrapper(wrapper))

function encodedblock(enc::Encoding, wrapper::WrapperBlock, ::PropagateAlways)
    inner = encodedblock(enc, wrapped(wrapper))
    return isnothing(inner) ? nothing : setwrapped(wrapper, inner)
end

function encodedblock(enc::Encoding, wrapper::WrapperBlock, ::PropagateNever)
    return encodedblock(enc, wrapped(wrapper))
end

function encodedblock(enc::Encoding, wrapper::WrapperBlock, ::PropagateSameBlock)
    inner = encodedblock(enc, wrapped(wrapper))
    inner == wrapped(block) && return setwrapped(wrapper, inner)
    return inner
end

decodedblock(enc::Encoding, wrapper::WrapperBlock) =
    decodedblock(enc, wrapper, propagatewrapper(wrapper))

function decodedblock(enc::Encoding, wrapper::WrapperBlock, ::PropagateAlways)
    inner = decodedblock(enc, wrapped(wrapper))
    return isnothing(inner) ? nothing : setwrapped(wrapper, inner)
end

function decodedblock(enc::Encoding, wrapper::WrapperBlock, ::PropagateNever)
    return decodedblock(enc, wrapped(wrapper))
end

function decodedblock(enc::Encoding, wrapper::WrapperBlock, ::PropagateSameBlock)
    inner = decodedblock(enc, wrapped(wrapper))
    inner == wrapped(block) && return setwrapped(wrapper, inner)
    return inner
end

# Encoding and decoding, if not overwritten for specific wrapper, are fowarded
# to wrapped block.

function encode(enc::Encoding, ctx, wrapper::WrapperBlock, obs; kwargs...)
    return encode(enc, ctx, wrapped(wrapper), obs; kwargs...)
end

function decode(enc::Encoding, ctx, wrapper::WrapperBlock, obs; kwargs...)
    return decode(enc, ctx, wrapped(wrapper), obs; kwargs...)
end

# ### Training interface
#
# The default behavior for `WrapperBlock`s is to forward behavior that is not
# explicitly overwritten to the wrapped block. So the same is done here for
# the training interface.

blockbackbone(wrapper::WrapperBlock) = blockbackbone(wrapped(wrapper))
blockmodel(wrapper::WrapperBlock, out, args...) = blockmodel(wrapped(wrapper), out, args...)
blockmodel(in::Block, out::WrapperBlock, args...) = blockmodel(in, wrapped(out), args...)

blocklossfn(wrapper::WrapperBlock, out) = blocklossfn(wrapped(wrapper), out)
blocklossfn(in::Block, out::WrapperBlock) = blocklossfn(in, wrapped(out))


text:
# Filtered catamorphisms

function cata(f, xtree, selector::Selector)
    return cata(xtree) do x
        matches(selector, x) ? f(x) : x
end
end

function catafirst(f, xtree, selector::Selector)
    xtree_, state_ = catafold(xtree, false) do x, done
        (!done && matches(selector, x)) ? (f(x), true) : (x, done)
    end
    return xtree_
end

# Replace

function replace(xtree, xnode, selector::Selector)
    cata(x -> xnode, xtree, selector)
end

function replacefirst(xtree, xnode, selector::Selector)
    return catafirst(x -> xnode, xtree, selector)
end

function replacemany(xtree, xnodes, selector::Selector)
    xtree_, state_ = catafold(xtree, 1) do x, i
        matches(selector, x) ? (xnodes[i], i + 1) : (x, i)
    end
    return xtree_
end

# Filter

function Base.filter(f, xtree::XNode)
    return cata(xtree) do x
        if x isa XLeaf
            return x
        else
            return withchildren(x, collect(filter(f, children(x))))
        end
    end
end

Base.filter(xtree::XNode, sel::Selector) = filter(x -> matches(sel, x), xtree)

# Insertion

abstract type Position end

struct NthChild <: Position
    n::Int
    selector::Selector
end

FirstChild(selector) = NthChild(1, selector)

struct After <: Position
    selector::Selector
end

struct Before <: Position
    selector::Selector
end

function insert(xtree, x, pos::Position)
    cata(xtree) do child
        i = insertionindex(child, children(child), pos)
        if isnothing(i)
            return child
        else
            return withchildren(child, _insert(children(child), i, x))
        end
    end
end

function insertfirst(xtree, x, pos::Position)
    xtree_, inserted = catafold(xtree, false) do child, inserted
        inserted && return (child, inserted)
        i = insertionindex(child, children(child), pos)
        if isnothing(i)
            return child, false
        else
            return withchildren(child, _insert(children(child), i, x)), true
        end
    end
    return xtree_
end

function insertionindex(parent, children, pos::NthChild)
    matches(pos.selector, parent) ? pos.n : nothing
end

function insertionindex(parent, children, pos::Before)
    for (i, c) in enumerate(children)
        if matches(pos.selector, c)
            return i
        end
    end
    return nothing
end

function insertionindex(parent, children, pos::After)
    for (i, c) in enumerate(children)
        if matches(pos.selector, c)
            return i + 1
        end
    end
    return nothing
end

_insert(xs::AbstractVector{T}, i, x::T) where T = insert!(copy(xs), i, x)
_insert(xs::AbstractVector{<:XTree}, i, x::XTree) = vcat(xs[1:i - 1], [x], xs[i:end])


text:
export Rocchio, RocchioBagging, fit, predict, transform
import SimilaritySearch: KnnResult
using StatsBase: countmap

mutable struct Rocchio
    protos::Vector{SVEC}  # prototypes
    pops::Vector{Int} # population per class
end

const RocchioBagging = Vector{Rocchio}

function fit(::Type{Rocchio}, X::AbstractVector{A}, y::AbstractVector{I}; nclasses=0) where {I<:Integer,A<:DVEC}
    if nclasses == 0
        nclasses = unique(y) |> length
    end

    prototypes = [SVEC() for i in 1:nclasses]
    populations = zeros(Int, nclasses)
    for i in 1:length(X)
        c = y[i]
        add!(prototypes[c], X[i])
        populations[c] += 1
    end

    n = length(X[1])
    for p in prototypes
        normalize!(p)
    end
    Rocchio(prototypes, populations)
end

function predict(rocchio::Rocchio, x::SVEC)
    res = KnnResult(1)
    for i in 1:length(rocchio.protos)
        d = cosine_distance(rocchio.protos[i], x)
        push!(res, i, d)
    end

    first(res).objID
end

function transform(rocchio::Rocchio, x::SVEC)
    [dot(rocchio.protos[i], x) for i in 1:length(rocchio.protos)]
end

function broadcastable(rocchio::Rocchio)
    (rocchio,)
end

function fit(::Type{RocchioBagging}, X::AbstractVector, y::AbstractVector{Int}, nrocchios=5; nclasses=0)::RocchioBagging
    if nclasses == 0
        nclasses = unique(y) |> length
    end

    s = ceil(Int, length(X) / 2)
    L = RocchioBagging()
    for i in 1:nrocchios
        I = rand(1:length(X), s)
        rocchio = fit(Rocchio, X[I], y[I], nclasses=nclasses)
        push!(L, rocchio)
    end

    L
end

function predict(rocchio::RocchioBagging, x)
    P = [predict(r, x) for r in rocchio] |> countmap |> collect
    sort!(P, by=x->x[end])
    P[end][1]
end

function transform(rocchio::RocchioBagging, x)
    sum(transform(r, x) for r in rocchio) |> normalize!
end

function broadcastable(rocchio::RocchioBagging)
    (rocchio,)
end


text: # gradient methods


text:
import Images.imread, Images.imwrite

function imread{S<:IO}(stream::S, ::Type{Images.Dummy})
    pixels = reshape(uint8(1:12), 3, 4)
    Image(pixels, ["colorspace" => "Gray", "spatialorder" => "xy"])
end

function imwrite(img, filename::String, ::Type{Images.Dummy})
    stream = open(filename, "w")
    println(stream, "Dummy Image")
    close(stream)
end


text:
module RRSTExperiments

using JuMP, CPLEX
using Random, Distributions
using Graphs
using SparseArrays

export InputEdge,
    solve_rec_st_with_algorithm, solve_rec_st_model, solve_minmax_st, solve_rec_st_hurwicz, solve_inc_st_with_model, solve_inc_st_with_algorithm,

    parse_graph_data, generate_uncertain_costs, generate_scenario,

    calculate_cost

struct InputEdge
    i::Int64
    j::Int64
    C::Float64
    c::Float64
    d::Float64

    InputEdge(i::Int, j::Int, C::Float64) = new(i, j, C, 0.0, 0.0)
    InputEdge(i::Int, j::Int, C::Float64, c::Float64) = new(i, j, C, c, 0.0)
    InputEdge(i::Int, j::Int, C::Float64, c::Float64, d::Float64) = new(i, j, C, c, d)
end

include("./utils.jl")

include("./inc_st.jl")
include("./rec_st_combinatorial.jl")
include("./rec_st_model.jl")
include("./rob_st_model.jl")
include("./rec_st_hurwicz.jl")

end


text:
__precompile__()
module AbstractTrees

export print_tree, TreeCharSet, Leaves, PostOrderDFS, Tree,
    AnnotationNode, StatelessBFS, treemap, treemap!, PreOrderDFS,
    ShadowTree, children, Leaves

import Base: getindex, setindex!, iterate, nextind, print, show,
    eltype, IteratorSize, length, push!, pop!
using Base: SizeUnknown
using Markdown

abstract type AbstractShadowTree end

include("traits.jl")
include("implicitstacks.jl")

# This package is intended to provide an abstract interface for working.
# Though the package itself is not particularly sophisticated, it defines
# the interface that can be used by other packages to talk about trees.

# By default assume that if an object is iterable, it's iteration gives the
# children. If an object is not iterable, assume it does not have children by
# default.
function children(x)
    if Base.isiterable(typeof(x)) && !isa(x, Integer) && !isa(x, Char) && !isa(x, Task)
        return x
    else
        return ()
    end
end
has_children(x) = children(x) !== ()

# Print a single node. Override this if you want your print function to print
# part of the tree by default
printnode(io::IO, node) = show(IOContext(io, :compact => true), node)

# Special cases

# Don't consider strings or reals tree-iterable in general
children(x::AbstractString) = ()
children(x::Real) = ()

# Define this here, there isn't really a good canonical package to define this
# elsewhere
children(x::Expr) = x.args

# To support iteration over associatives, define printnode on Tuples to return
# the first element. If this doesn't work well in practice it may be better to
# create a special iterator that `children` on `Associative` returns.
# Even better, iteration over associatives should return pairs.

printnode(io::IO, kv::Pair{K,V}) where {K,V} = printnode(io,kv[1])
children(kv::Pair{K,V}) where {K,V} = (kv[2],)

printnode(io::IO, d::Dict{K,V}) where {K,V} = print(io, Dict{K,V})
printnode(io::IO, d::Vector{T}) where {T} = print(io, Vector{T})

# Node equality predicate
nodeequal(a, b) = a === b

# Utilities

# Printing
struct TreeCharSet
    mid
    terminator
    skip
    dash
end

# Default charset
TreeCharSet() = TreeCharSet('├','└','│','─')

function print_prefix(io, depth, charset, active_levels)
    for current_depth in 0:(depth-1)
        if current_depth in active_levels
            print(io,charset.skip," "^(textwidth(charset.dash)+1))
        else
            print(io," "^(textwidth(charset.skip)+textwidth(charset.dash)+1))
        end
    end
end

@doc doc"""
# Usage
Prints an ASCII formatted representation of the `tree` to the given `io` object.
By default all children will be printed up to a maximum level of 5, though this
valud can be overriden by the `maxdepth` parameter. The charset to use in
printing can be customized using the `charset` keyword argument.

# Examples
```julia
julia> print_tree(STDOUT,Dict("a"=>"b","b"=>['c','d']))
Dict{String,Any}("b"=>['c','d'],"a"=>"b")
├─ b
│  ├─ c
│  └─ d
└─ a
   └─ b

julia> print_tree(STDOUT,Dict("a"=>"b","b"=>['c','d']);
        charset = TreeCharSet('+','\\','|',"--"))
Dict{String,Any}("b"=>['c','d'],"a"=>"b")
+-- b
|   +-- c
|   \-- d
\-- a
   \-- b
```

"""
function _print_tree(printnode::Function, io::IO, tree, maxdepth = 5; depth = 0, active_levels = Int[],
    charset = TreeCharSet(), withinds = false, inds = [], from = nothing, to = nothing, roottree = tree)
    nodebuf = IOBuffer()
    isa(io, IOContext) && (nodebuf = IOContext(nodebuf, io))
    if withinds
        printnode(nodebuf, tree, inds)
    else
        tree != roottree && isa(treekind(roottree), IndexedTree) ?
            printnode(nodebuf, roottree[tree]) :
            printnode(nodebuf, tree)
    end
    str = String(take!(isa(nodebuf, IOContext) ? nodebuf.io : nodebuf))
    for (i,line) in enumerate(split(str, '\n'))
        i != 1 && print_prefix(io, depth, charset, active_levels)
        println(io, line)
    end
    c = isa(treekind(roottree), IndexedTree) ?
        childindices(roottree, tree) : children(roottree, tree)
    if c !== ()
        s = Iterators.Stateful(from === nothing ? pairs(c) : Iterators.Rest(pairs(c), from))
        while !isempty(s)
            ind, child = popfirst!(s)
            ind === to && break
            active = false
            child_active_levels = active_levels
            print_prefix(io, depth, charset, active_levels)
            if isempty(s)
                print(io, charset.terminator)
            else
                print(io, charset.mid)
                child_active_levels = push!(copy(active_levels), depth)
            end
            print(io, charset.dash, ' ')
            print_tree(printnode, io, child; depth = depth + 1,
              active_levels = child_active_levels, charset = charset, withinds=withinds,
              inds = withinds ? [inds; ind] : [], roottree = roottree)
        end
    end
end
print_tree(f::Function, io::IO, tree, args...; kwargs...) = _print_tree(f, io, tree, args...; kwargs...)
print_tree(io::IO, tree, args...; kwargs...) = print_tree(printnode, io, tree, args...; kwargs...)
print_tree(tree, args...; kwargs...) = print_tree(stdout::IO, tree, args...; kwargs...)

# Tree Indexing
struct Tree
    x::Any
end
Tree(x::Tree) = x
Tree(x::AbstractShadowTree) = x
show(io::IO, tree::Tree) = print_tree(io, tree.x)

mutable struct AnnotationNode{T}
    val::T
    children::Array{AnnotationNode{T}}
end

children(x::AnnotationNode) = x.children
printnode(io::IO, x::AnnotationNode) = print(io, x.val)

struct ShadowTree <: AbstractShadowTree
    tree::Tree
    shadow::Tree
    ShadowTree(x::Tree,y::Tree) = new(x,y)
    ShadowTree(x,y) = ShadowTree(Tree(x),Tree(y))
end
first_tree(x::ShadowTree) = x.tree
second_tree(x::ShadowTree) = x.shadow

function zip_min(c1, c2)
    n1, n2 = length(c1), length(c2)
    if n1 < n2
        c2 = Iterators.take(c2,n1)
    elseif n2 < n1
        c1 = Iterators.take(c1,n2)
    end
    zip(c1, c2)
end

make_zip(x::AbstractShadowTree) = zip_min(children(x.tree.x), children(x.shadow.x))

function children(x::AbstractShadowTree)
    map(res->typeof(x)(res[1], res[2]),make_zip(x))
end

iterate(x::AbstractShadowTree, state...) = iterate(make_zip(x), state...)

function make_annotations(cb, tree, parent, s)
    s = cb(tree, parent, s)
    AnnotationNode{Any}(s, AnnotationNode{Any}[make_annotations(cb, child, tree, s) for child in children(tree)])
end

function getindex(tree::Tree, indices)
    node = tree.x
    for idx in indices
        node = children(node)[idx]
    end
    node
end
getindex(tree::Tree, indices::T) where {T<:ImplicitNodeStack} =
    getindex(tree, indices.idx_stack.stack)

function getindexhighest(tree::Tree, indices)
    node = tree.x
    for (i,idx) in enumerate(indices)
        cs = children(node)
        if idx > length(cs)
            return (indices[1:i-1],node)
        end
        node = children(node)[idx]
    end
    (indices, node)
end

function setindex!(tree::Tree, val, indices)
    setindex!(children(getindex(tree,indices[1:end-1])),val,indices[end])
end
setindex!(tree::Tree, val, indices::T) where {T<:ImplicitNodeStack} =
    setindex!(tree, val, indices.idx_stack.stack)

function getindex(tree::AbstractShadowTree, indices)
    typeof(tree)(Tree(first_tree(tree))[indices],Tree(second_tree(tree))[indices])
end

function setindex!(tree::AbstractShadowTree, val, indices)
    setindex!(Tree(first_tree(tree)), first_tree(val), indices)
    setindex!(Tree(second_tree(tree)), second_tree(val), indices)
end

function setindex!(tree::AbstractShadowTree, val::Tuple, indices)
    setindex!(Tree(first_tree(tree)), val[1], indices)
    setindex!(Tree(second_tree(tree)), val[2], indices)
end

# Tree Iterators

abstract type TreeIterator{T} end

"""
Iterator to visit the leaves of a tree, e.g. for the tree

Any[1,Any[2,3]]
├─ 1
└─ Any[2,3]
   ├─ 2
   └─ 3

we will get [1,2,3]
"""
struct Leaves{T} <: TreeIterator{T}
    tree::T
end
IteratorSize(::Type{Leaves{T}}) where {T} = SizeUnknown()

"""
Iterator to visit the nodes of a tree, guaranteeing that children
will be visited before their parents.

e.g. for the tree

Any[1,Any[2,3]]
├─ 1
└─ Any[2,3]
   ├─ 2
   └─ 3

we will get [1,2,3,Any[2,3],Any[1,Any[2,3]]]
"""
struct PostOrderDFS <: TreeIterator{Any}
    tree::Any
    PostOrderDFS(x::Any) = new(x)
end
PostOrderDFS(tree::Tree) = PostOrderDFS(tree.x)
IteratorSize(::Type{PostOrderDFS}) = SizeUnknown()

"""
Iterator to visit the nodes of a tree, guaranteeing that parents
will be visited before their children.

Optionally takes a filter function that determines whether the iterator
should continue iterating over a node's children (if it has any) or should
consider that node a leaf.

e.g. for the tree

Any[Any[1,2],Any[3,4]]
├─ Any[1,2]
|  ├─ 1
|  └─ 2
└─ Any[3,4]
   ├─ 3
   └─ 4

we will get [Any[Any[1,2],Any[3,4]],Any[1,2],1,2,Any[3,4],3,4]

# Invalidation
Modifying the underlying tree while iterating over it, is allowed, however,
if parents and sibling links are not explicitly stored, the identify of any
parent of the last obtained node does not change (i.e. mutation is allowed,
replacing nodes is not).

"""
struct PreOrderDFS{T} <: TreeIterator{T}
    tree::T
    filter::Function
    PreOrderDFS{T}(tree,filter::Function=(args...)->true) where {T} = new{T}(tree,filter)
end
PreOrderDFS(tree::T,filter::Function=(args...)->true) where {T} = PreOrderDFS{T}(tree,filter)
PreOrderDFS(tree::Tree,filter::Function=(args...)->true) = PreOrderDFS(tree.x,filter)
IteratorSize(::Type{PreOrderDFS{T}}) where {T} = SizeUnknown()

# State depends on what kind of tree we have:
#   - Parents/Siblings are not stored:
#       - RegularTree: ImplicitNodeStack
#       - IndexedTree: ImplicitIndexStack
#   - Parents/Siblings are stored:
#       - RegularTree: Nodes
#       - IndexedTree: Indices
#
childstates(tree, state, ::IndexedTree) = childindices(tree, state)
childstates(tree, state, ::RegularTree) = children(tree, state)
parentstate(tree, state, ::IndexedTree) = parentind(tree, state)
parentstate(tree, state, ::RegularTree) = parent(tree, state)

parentstate(tree, state) = parentstate(tree, state, treekind(tree))

update_state!(old_state, cs, idx) = next(cs, idx)[1]

struct ImplicitRootState
end
getindex(x, ::ImplicitRootState) = x
getindex(x::AbstractArray, ::ImplicitRootState) = x

"""
Trees must override with method if the state of the root is not the same as the
tree itself (e.g. IndexedTrees should always override this method).
"""
rootstate(x) = ImplicitRootState()

function firststate(ti::PreOrderDFS{T}) where T
    if isa(parentlinks(ti.tree), StoredParents) &&
            isa(siblinglinks(ti.tree), SiblingLinks)
        rootstate(ti.tree)
    else
        state = ImplicitIndexStack(idxtype(ti.tree)[])
        if !isa(treekind(typeof(ti.tree)), IndexedTree)
            state = ImplicitNodeStack(nodetype(ti.tree)[], state)
        end
        state
    end
end
function firststate(ti::Union{Leaves, PostOrderDFS})
    state = firststate(PreOrderDFS(ti.tree))
    while true
        css = childstates(ti.tree, state)
        isempty(css) && break
        state = first(css)
    end
    state
end

nextind(::Base.Generator, idx) = idx + 1
relative_state(tree, parentstate, childstate::ImplicitIndexStack) =
    childstate.stack[end]
relative_state(tree, parentstate, childstate::ImplicitNodeStack) =
    relative_state(tree, parentstate, childstate.idx_stack)
function nextsibling(tree, state)
    ps = parentstate(tree, state)
    cs = childstates(tree, ps)
    isempty(cs) && return nothing
    new_state = nextind(cs, relative_state(tree, ps, state))
    iterate(cs, new_state) === nothing && return nothing
    update_state!(tree, ps, children(tree, ps), new_state)
end

function nextsibling(node, ::StoredParents, ::ImplicitSiblings, ::RegularTree)
    isroot(node) && return nothing
    p = parent(node)
    last_was_node = false
    for c in children(p)
        last_was_node && return c
        (c == node) && (last_was_node = true)
    end
    last_was_node && return nothing
    error("Tree inconsistency: node not a child of parent")
end
nextsibling(node, ::Any, ::StoredSiblings, ::Any) = error("Trees with explicit siblings must override the `prevsibling` method explicitly")
nextsibling(node) = nextsibling(node, parentlinks(node), siblinglinks(node), treekind(node))

function prevsibling(node, ::StoredParents, ::ImplicitSiblings, ::RegularTree)
    isroot(node) && return nothing
    p = parent(node)
    last_c = nothing
    for c in children(p)
        (c == node) && return last_c
        last_c = c
    end
    @show p
    @show node
    error("Tree inconsistency: node not a child of parent")
end
prevsibling(node, ::Any, ::StoredSiblings, ::Any) = error("Trees with explicit siblings must override the `prevsibling` method explicitly")
prevsibling(node) = prevsibling(node, parentlinks(node), siblinglinks(node), treekind(node))
prevsibling(tree, node) = prevsibling(node)

isroot(tree, state, ::RegularTree) = tree == state
isroot(tree, state, ::IndexedTree) = state == rootstate(tree)
isroot(tree, state) = isroot(tree, state, treekind(tree))

struct Subtree{T,S}
    tree::T
    state::S
end
children(tree::Subtree) = children(tree.tree, tree.state)
nodetype(tree::Subtree) = nodetype(tree.tree)
idxtype(tree::Subtree) = idxtype(tree.tree)

joinstate(tree, a, b) = b

if isdefined(Base, :UnionAll)
    Base.@pure function get_primary(T::DataType)
        T.name.wrapper
    end
else
    Base.@pure function get_primary(T::DataType)
        T.name.primary
    end
end

function stepstate(ti::TreeIterator, state)
    if isa(ti, PreOrderDFS) && ti.filter(getnode(ti.tree, state))
        ccs = childstates(ti.tree, state)
        !isempty(ccs) && return first(ccs)
    end
    while !isroot(ti.tree, state)
        nextstate = nextsibling(ti.tree, state)
        if nextstate !== nothing
            return joinstate(ti.tree, nextstate, firststate(
                get_primary(typeof(ti))(Subtree(ti.tree, nextstate))))
        end
        state = parentstate(ti.tree, state)
        isa(ti, PostOrderDFS) && return state
    end
    nothing
end

getnode(tree, ns) = isa(treekind(tree), IndexedTree) ? tree[ns] : ns
getnode(tree::AbstractShadowTree, ns::ImplicitNodeStack) = tree[ns.idx_stack.stack]

function iterate(ti::TreeIterator)
    state = firststate(ti)
    (getnode(ti.tree, state), state)
end
function iterate(ti::TreeIterator, state)
    state = stepstate(ti, state)
    state === nothing && return nothing
    (getnode(ti.tree, state), state)
end

"""
    Acends the tree, at each node choosing whether or not to continue.
    Note that the parent is computed before the callback is exectuted, allowing
    modification of the argument to the callback (as long as the overall tree
    structure is not altered).
"""
function ascend(select, node)
    isroot(node) && (select(node); return node)
    p = parent(node)
    while select(node) && !isroot(node)
        node = p
        p = parent(node)
    end
    node
end

"""
    Descends the tree, at each node choosing the child given by select callback
    or the current node if 0 is returned.
"""
function descend(select, tree)
    idx = select(tree)
    idx == 0 && return tree
    node = children(tree)[idx]
    while true
        idx = select(node)
        idx == 0 && return node
        node = children(node)[idx]
    end
end

"""
Iterator to visit the nodes of a tree, all nodes of a level will be visited
before their children

e.g. for the tree

Any[1,Any[2,3]]
├─ 1
└─ Any[2,3]
   ├─ 2
   └─ 3

we will get [Any[1,Any[2,3]],1,Any[2,3],2,3]

WARNING: This is O(n^2), only use this if you know you need it, as opposed to
a more standard statefull approach.
"""
struct StatelessBFS <: TreeIterator{Any}
    tree::Any
end
IteratorSize(::Type{StatelessBFS}) = SizeUnknown()

function descend_left(newinds, next_node, level)
    # Go down until we are at the correct level or a dead end
    while length(newinds) != level
        cs = children(next_node)
        if isempty(cs)
            break
        end
        push!(newinds, 1)
        next_node = first(cs)
    end
    return newinds
end

function nextind_or_deadend(tree, ind, level)
    current_level = active_level = length(ind)
    active_inds = copy(ind)
    # Go up until there is a right neighbor
    while current_level > 0
        # Check for next node at the current level
        active_inds = ind[1:current_level-1]
        parent = Tree(tree)[active_inds]
        cur_child = ind[current_level]
        ni = nextind(children(parent), cur_child)
        current_level -= 1
        if iterate(children(parent), ni) !== nothing
            newinds = [active_inds; ni]
            next_node = children(parent)[ni]
            return descend_left(newinds, next_node, level)
        end
    end
    return nothing
end

iterate(ti::StatelessBFS) = (Tree(ti.tree)[[]], [])
"""
Stateless level-order bfs iteration. The algorithm is as follows:

Go up. If there is a right neighbor, go right, then left until you reach the
same level. If you reach the root, go left until you reach the next level.
"""
function iterate(ti::StatelessBFS, ind)
    org_level = active_level = length(ind)
    newinds = ind
    while true
        newinds = nextind_or_deadend(ti.tree, newinds, active_level)
        if newinds === nothing
            active_level += 1
            if active_level > org_level + 1
                return nothing
            end
            newinds = descend_left([], ti.tree, active_level)
        end
        if length(newinds) == active_level
            break
        end
    end
    Tree(ti.tree)[newinds], newinds
end

# Mapping over trees
function treemap(f::Function, tree::PostOrderDFS)
    new_tree = Any[Union{}[]]
    current_length = 0
    for (ind, node) in pairs(tree)
        while length(new_tree) < length(ind)
            push!(new_tree, Union{}[])
        end
        thechildren = Union{}[]
        if length(ind) < length(new_tree)
            thechildren = pop!(new_tree)
        end
        if ind == []
            return f(ind, node, thechildren)
        end
        siblings = new_tree[end]
        el = f(ind, node, thechildren)
        S = typeof(el)
        T = eltype(siblings)
        if S === T || S <: T
            push!(siblings, el)
        else
            R = typejoin(T, S)
            new = similar(siblings, R)
            copy!(new,1,siblings,1,length(siblings))
            push!(new,el)
            new_tree[end] = new
        end
    end
end

function treemap!(f::Function, ti::PreOrderDFS)
    state = firststate(ti)
    while state !== nothing
        ind = state
        node = getnode(ti.tree, ind)
        new_node = f(node)
        if new_node !== node
            if isempty(ind)
                return treemap!(PreOrderDFS(new_node)) do x
                    x == new_node && return x
                    f(x)
                end
            end
            Tree(ti.tree)[ind] = new_node
        end
        state = stepstate(ti, ind)
    end
    ti.tree
end

end # module


text:
function K!(Kg, g, grid, β, ∂u∂c, f, f′, shocks)
# This function requires the container of the output value as argument Kg

    # Construct linear interpolation object
    g_func = LinearInterpolation(grid, g, extrapolation_bc=Line())

    # solve for updated consumption value
    for (i, y) in enumerate(grid)
        function h(c)
            vals = ∂u∂c.(g_func.(f(y - c) * shocks)) .* f′(y - c) .* shocks
            return ∂u∂c(c) - β * mean(vals)
        end
        Kg[i] = find_zero(h, (1e-10, y - 1e-10))
    end
    return Kg
end

# The following function does NOT require the container of the output value as argument
K(g, grid, β, ∂u∂c, f, f′, shocks) =
    K!(similar(g), g, grid, β, ∂u∂c, f, f′, shocks)

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using LatexTables
using Test

body  = reshape(["v", 1, 2, 3], 1, :)
#body  = vcat(["v" 1 2 3],  ["a" 2 1 0.5])
n_row = size(body)[1]

run_finished(s::String) = return true
compare_strings(s1::String, s2::String) = remove_spaces(s1) == remove_spaces(s2)
compare_strings_rows(s1::String, s2::String) = split_string(s1) == split_string(s2)

function compare_strings_rows_subset(s1::String, s2::String)
    s1_spl = split_string(s1)
    s2_spl = split_string(s2)
    return is_subset(s1_spl, s2_spl) && !is_subset(s2_spl, s1_spl)
end

function compile_success(file_name, s::String)

    dir_name = "output"
    file_pdf = file_name * ".pdf"
    file_tex = file_name * ".tex"

    isfile(joinpath(dir_name, file_pdf)) && rm(joinpath(dir_name, file_pdf))

    make_tex_file(file_tex, s; dir_name=dir_name, compile=true)

    return isfile(joinpath(dir_name, file_pdf))
end

result1 = "\\begin{table}[!ht] \n\\centering \n\\caption{cap} \n\\label{lab} \n\\begin{tabular}{@{} llll @{}} \n\\toprule \nv & \$1\$ & \$2\$ & \$3\$ \\\\ \n\\bottomrule \n\\end{tabular} \n\\end{table} \n"
result2 = "\\begin{table}[!ht] \n\\centering \n\\caption{cap} \n\\label{lab} \n\\begin{tabular}{@{} llll @{}} \n\\hline \nv & \$1\$ & \$2\$ & \$3\$ \\\\ \n\\hline \n\\end{tabular} \n\\end{table} \n"

function split_string(s::String)
    s = remove_spaces(s)
    s = sort(split(s, "\n"))
    return s[.!isempty.(s)]
end

function remove_spaces(s::String)
    ii = [t[1] for t in collect.(findall(" ", s))]
    jj = setdiff(1:length(s), ii)
    return s[jj]
end

function is_subset(s1::AbstractArray, s2::AbstractArray)
    all([any(x1 .== s2) for x1 in s1])
end

@testset "Introduction" begin
    @test_throws Exception error("Kdo testuje, neveri koderum.")
end

@testset "Highlighting" begin
    @test_throws Exception table_to_tex(body; highlight_max_row=true, highlight_max_col=true)
    @test_throws Exception table_to_tex(body; highlight_min_row=true, highlight_min_col=true)
    @test run_finished(table_to_tex(body; highlight_max_row=true, highlight_max_style=Color(:blue)))
    @test run_finished(table_to_tex(body; highlight_max_row=true, highlight_max_style=[Color(:blue), Style(:italic), CellColor(:red)]))
    @test run_finished(table_to_tex(body; highlight_min_row=true, highlight_min_style=Color(:blue)))
    @test run_finished(table_to_tex(body; highlight_min_row=true, highlight_min_style=[Color(:blue), Style(:italic), CellColor(:red)]))
end

@testset "Types" begin
    @test_throws Exception table_to_tex(body; col_format=["", "a"])
    @test run_finished(table_to_tex(body; col_format=["s", "3f", "2d", 'd']))
end

@testset "Leading col" begin
    @test run_finished(table_to_tex(body; leading_col=repeat(["1"], n_row), alignment=""))
    @test run_finished(table_to_tex(body; leading_col=repeat([1], n_row), alignment='l'))
    @test run_finished(table_to_tex(body; leading_col=repeat(["1"], n_row), alignment="l"))
    @test run_finished(table_to_tex(body; leading_col=repeat([1], n_row), alignment="llll"))
    @test run_finished(table_to_tex(body; leading_col=repeat(["1"], n_row), alignment="lllll"))
end

@testset "Comparison" begin
    @test compare_strings_rows(result1, table_to_tex(body; table_type=:booktabs, position="!ht", caption="cap", label="lab", centering=true, alignment="l", caption_position_top=true))
    @test compare_strings_rows(result1, table_to_tex(body; table_type=:booktabs, position="!ht", caption="cap", label="lab", centering=true, alignment="l", caption_position_top=false))
    @test compare_strings_rows(result2, table_to_tex(body; table_type=:tabular, position="!ht", caption="cap", label="lab", centering=true, alignment="l", caption_position_top=true))
    @test compare_strings_rows(result2, table_to_tex(body; table_type=:tabular, position="!ht", caption="cap", label="lab", centering=true, alignment="l", caption_position_top=false))
    @test !compare_strings_rows_subset(result1, table_to_tex(body; table_type=:booktabs, position="!ht", centering=true, alignment="l", caption_position_top=true, floating_table=false))
    @test compare_strings_rows_subset(table_to_tex(body; table_type=:booktabs, position="!ht", centering=true, alignment="l", caption_position_top=true, floating_table=false), result1)
end

@testset "Compile" begin
    @test compile_success("test12344321", table_to_tex(body; table_type=:booktabs, position="!ht", caption="cap", label="lab", centering=true, alignment="l", caption_position_top=true, highlight_max_row=true, highlight_min_row=true))
    @test compile_success("test12344322", table_to_tex(body; table_type=:booktabs, position="!ht", caption="cap", label="lab", centering=true, alignment="l", caption_position_top=true, highlight_max_col=true, highlight_min_col=true))
end


text:
module  A3

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\A3.jl"),
  informat="script",
  out_path = "$(pwd())\\A3.html",
  doctype = "md2html")
=#
if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt, ForwardDiff

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const DATE_FORMAT_STR = "m/d/y"
const SP500_NAME = "GSPC"
const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"

const DEFAULT_NUM_SIMS = 10000
const PlotContainer = Union{Plot,Gadfly.Compose.Context}

#this holds the parameters of an AR1 Process
struct AR1Process
  α::Float64
  β::Float64
  params::Vector{Float64}
  ε::UnivariateDistribution
end

#main constructor for an AR1 with sensible defaults
AR1Process(α::Float64, β::Float64; params::Vector{Float64}=Vector{Float64}(),
  dist::Type = Normal)::AR1Process =
  AR1Process(α, β, params, Normal(params...))

#constructor to get a normal AR1 Process
AR1Process(α::Float64, β::Float64, μ::Float64, σ::Float64)::AR1Process =
          AR1Process(α, β, params=[μ,σ], dist=Normal)

#this function gets a series given an AR1 process
function getAR1Series!(p::AR1Process, T::Int,
  p₀::Float64=(p.β==1.0)?rand(p.ε)+p.α:rand(Normal(p.α,(p.params[2]^2/(1-p.β^2)))),
  series::Vector{Float64} = Vector{Float64}(T+1))::Vector{Float64}

  series[1] = p₀
  series[2:(T+1)] = rand(p.ε, T)

  for i::Int ∈ 1:T
    series[1+i] += p.α + p.β*series[i]
  end

  return series
end

function getAR1ForOLS!(p::AR1Process, T::Int;
                      p₀::Float64=rand(p.ε)+p.α,
                      AR1DF::DataFrame = DataFrame([Float64, Float64], [:Y, :X], T),
                      series::Vector{Float64} = Vector{Float64}(T+1))::DataFrame

  getAR1Series!(p, T, p₀, series)

  AR1DF[:,:X] .= series[1:(end-1)]
  AR1DF[:,:Y] .= series[2:end]

  return AR1DF
end

#this holds information from the simulation
mutable struct AR1OLSSim
  AR1::AR1Process
  α::Vector{Float64}
  β::Vector{Float64}
  σ::Vector{Float64}
  #Σ::Vector{Matrix{Float64}}
end

#this function constructs and runs the OLS sim
function AR1OLSSim(p::AR1Process, T::Int, N::Int)::AR1OLSSim

  #pre-allocate information from the regression
  α::Vector{Float64} = Vector{Float64}(N)
  β::Vector{Float64} = Vector{Float64}(N)
  σ::Vector{Float64} = Vector{Float64}(N)
  #Σ::Vector{Matrix{Float64}} = [Matrix{Float64}(2,2) for i::Int ∈ 1:N]

  #use this to store the data for a given OLS (ONLY FOR SINGLE THREADED RUNS)
  #series::Vector{Float64} = Vector{Float64}(T+1)
  #AR1DF::DataFrame = DataFrame([Float64, Float64], [:Y, :X], T)

  #Use some cores
  @fastmath Threads.@threads for i::Int ∈ 1:N
    AR1DF::DataFrame = getAR1ForOLS!(p,T)

    reg::CTLM = CTLM(AR1DF,  :X, :Y, XNames = [:intercept, :X], YName=:Y,
        eliminateNulls=false)

      #get the coefficients
      α[i], β[i] = reg.β
      σ[i] = std(reg.ε)

    #getModWhiteΣ!(reg, Σ[i])
  end

  return AR1OLSSim(deepcopy(p),α, β, σ)
end

function makeSummaryTables(sims::Vector{AR1OLSSim})
  ########################make the summary data table
  summaryDecimals::Int = 4
  numCols = length(sims)
  summaryCols::Vector{String} = ["(\$\\beta=1.0, T=50\$)",
    "(\$\\beta=1.0, T=600\$)", "(\$\\beta=0.95, T=50\$)"]
  summaryRowNames::Vector{String} =
    ["\$\\beta_{true}\$", "\$Mean(\\beta_{est})\$", "\$\\beta\$ Standard Deviation",
    "\$\\alpha_{true}\$", "\$Mean(\\alpha_{est})\$", "\$\\alpha\$ Standard Deviation",
    "\$\\sigma_{true}\$", "\$Mean(\\sigma_{est})\$", "\$\\sigma\$ Standard Deviation",
      "T-Stat (\$\\beta_{est}-\\beta_{true}\$)", "\$\\beta_{true}\$ - 5\\% T-Tail",
      "\$\\beta_{true}\$ - 1\\% T-Tail", "N"]
  summaryNumRows::Int = length(summaryRowNames)

  βTrues::Vector{Float64} = ((i::Int)->(sims[i].AR1.β)).(1:numCols)
  βEstimates::Vector{Float64} = ((sim::AR1OLSSim)->mean(sim.β)).(sims)
  βStdDevs::Vector{Float64} = ((sim::AR1OLSSim)->std(sim.β)).(sims)
  βTStats::Vector{Float64} = (βEstimates .- βTrues) ./ βStdDevs
  Ns::Vector{Int} = ((sim::AR1OLSSim)->length(sim.β)).(sims)
  βT5Tails = βTrues .+ βStdDevs .* ((n::Int)->quantile(TDist(n-2), 0.05)).(Ns)
  βT1Tails = βTrues .+ βStdDevs .* ((n::Int)->quantile(TDist(n-2), 0.01)).(Ns)

  αTrues::Vector{Float64} = ((i::Int)->(sims[i].AR1.α)).(1:numCols)
  αEstimates::Vector{Float64} = ((sim::AR1OLSSim)->mean(sim.α)).(sims)
  αStdDevs::Vector{Float64} = ((sim::AR1OLSSim)->std(sim.α)).(sims)

  σTrues::Vector{Float64} = ((i::Int)->(sims[i].AR1.params[2])).(1:numCols)
  σEstimates::Vector{Float64} = ((sim::AR1OLSSim)->mean(sim.σ)).(sims)
  σStdDevs::Vector{Float64} = ((sim::AR1OLSSim)->std(sim.σ)).(sims)

  #make the content matrix
  summaryContent::Vector{Vector{String}} = Vector{Vector{String}}()

  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(βTrues)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(βEstimates)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(βStdDevs)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(αTrues)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(αEstimates)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(αStdDevs)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(σTrues)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(σEstimates)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(σStdDevs)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(βTStats)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(βT5Tails)))
  push!(summaryContent, string.(((x::Float64)->round(x,summaryDecimals)).(βT1Tails)))
  push!(summaryContent, string.(Ns))

  tableSummary::String = texTable( "Summary of Results",
    """See Tex File""", #caption
    [summaryCols], #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    summaryRowNames, #descRowNames
    summaryContent, #descContent
    Vector{String}(),
    columnSepPt = -5)
end

#problem specific script
function simulationScript(; N::Int=10000,
  T::Int = 50,
  α::Float64 = 0.0,
  β::Float64 = 1.0,
  μ::Float64 = 0.0,
  σ::Float64 = 0.2,
  staticsPoints::Int = 50,
  discountForStatics::Float64 = 0.1)::AR1OLSSim

  #main specification
  mainProblem::AR1Process = AR1Process(α, β, μ, σ)
  mainSim::AR1OLSSim = AR1OLSSim(mainProblem, T, N)
  println("""α mean: $(mean(mainSim.α))
            β mean: $(mean(mainSim.β)))""")

  #comparitive statics specifications
  σStatics::Vector{AR1Process} =  AR1Process.(α, β, μ, [i / staticsPoints
    for i::Int ∈ 1:staticsPoints]::Vector{Float64})
  αStatics::Vector{AR1Process} =  AR1Process.([i / staticsPoints - 0.5
    for i::Int ∈ 1:staticsPoints]::Vector{Float64}, β, μ, σ)
  simσStatics::Vector{AR1OLSSim} = AR1OLSSim.(σStatics, T, Int(floor(N*discountForStatics)))
  simαStatics::Vector{AR1OLSSim} = AR1OLSSim.(αStatics, T, Int(floor(N*discountForStatics)))

  Gadfly.push_theme(:default)
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  push!(plotNames, "Statics_Beta$(β)_N$(N)_T$(T)")
  push!(plots, vstack(
    plot(x=((ar1::AR1Process)->ar1.α).(αStatics),
      y=((sim::AR1OLSSim)->mean(sim.β)).(simαStatics), Geom.point,
      Guide.ylabel("β"), Guide.xlabel("α"),
      Guide.title("Estimates of β vs. α (True β=1.0, N=$(Int(floor(discountForStatics*N))))")),
    plot(x=((ar1::AR1Process)->ar1.params[2]).(σStatics),
      y=((sim::AR1OLSSim)->mean(sim.β)).(simσStatics), Geom.point,
      Guide.ylabel("β"), Guide.xlabel("σ"),
      Guide.title("Estimates of β vs. σ (True β=1.0, N=$(Int(floor(discountForStatics*N))))"))))

      push!(plotNames, "Statics_AlphaSigma$(β)_N$(N)_T$(T)")
      push!(plots, vstack(
        plot(x=((ar1::AR1Process)->ar1.α).(αStatics),
          y=((sim::AR1OLSSim)->mean(sim.σ)).(simαStatics), Geom.point,
          Guide.ylabel("σ"), Guide.xlabel("α"),
          Guide.title("Estimates of σ vs. α (True σ=0.2, N=$(Int(floor(discountForStatics*N))))")),
        plot(x=((ar1::AR1Process)->ar1.params[2]).(σStatics),
          y=((sim::AR1OLSSim)->mean(sim.α)).(simσStatics), Geom.point,
          Guide.ylabel("α"), Guide.xlabel("σ"),
          Guide.title("Estimates of α vs. σ (True α=0.0, N=$(Int(floor(discountForStatics*N))))"))))

    for i ∈ 1:length(plots)
      draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 7.5inch, 9inch),plots[i])
    end

  return mainSim
end

function A3Script()::Void

  sims::Vector{AR1OLSSim} = Vector{AR1OLSSim}()

  push!(sims, simulationScript())
  push!(sims, simulationScript(T=600))
  push!(sims, simulationScript(β=0.95))

  writeTables2File([makeSummaryTables(sims)],
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "SummaryTable.tex")

  return nothing
end

#enable to run
@time A3Script()

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# print information over iterations of nndescent
using NearestNeighborDescent
using LightGraphs
using Distances
using MLDatasets
using Statistics

include("../benchmark/benchutils.jl")

function diameter_stats(graph)
    diameters = [knn_diameter(graph, v) for v in vertices(graph)]
    min_diam, med_diam, max_diam = minimum(diameters), median(diameters), maximum(diameters)
    mean_diam = mean(diameters)
    return min_diam, mean_diam, med_diam, max_diam
end

num_new_edges(graph) = sum(KNNGraphs.flag(e) for e in edges(graph))

function do_iter(graph, data, metric, sample_rate, true_indices)
    c = NearestNeighborDescent.local_join!(graph, data, metric; sample_rate=sample_rate)
    min_diam, mean_diam, med_diam, max_diam = diameter_stats(graph)
    nne = num_new_edges(graph)
    indices = knn_matrices(graph)[1]
    rec = recall(indices, true_indices)
    return c, (min_diam, mean_diam, med_diam, max_diam), nne, rec
end

function join_with_history(data, k, metric, sample_rate, true_indices)
    graph = HeapKNNGraph(data, k, metric)
    indices = knn_matrices(graph)[1]
    count_history = []
    diam_history = []
    new_edge_history = []
    recall_history = []

    count = ne(graph)
    push!(count_history, "iter_0"=>count)
    push!(diam_history, "iter_0"=>diameter_stats(graph))
    push!(new_edge_history, "iter_0"=>num_new_edges(graph))
    recall₀ = recall(indices, true_indices)
    push!(recall_history, "iter_0"=>recall₀)
    iter = 1
    while count > 0
        count, diams, nne, rec = do_iter(graph, data, metric, sample_rate, true_indices)
        push!(count_history, "iter_$iter"=>count)
        push!(diam_history, "iter_$iter"=>diams)
        push!(new_edge_history, "iter_$iter"=>nne)
        push!(recall_history, "iter_$iter"=>rec)
        iter += 1
    end
    return count_history, diam_history, new_edge_history, recall_history
end

mnist_x_arr = [col for col in
    eachcol(MNIST.convert2features(MNIST.traintensor(Float64, 1:10000)))]
k = 20
metric = SqEuclidean()

true_graph = NearestNeighborDescent.brute_knn(mnist_x_arr, metric, k)
true_indices = getindex.(true_graph, 1)

c_hist_100, d_hist_100, e_hist_100, r_hist_100 =
    join_with_history(mnist_x_arr, k, metric, 1, true_indices)

c_hist_25, d_hist_25, e_hist_25, r_hist_25 =
    join_with_history(mnist_x_arr, k, metric, .25, true_indices)


text:
using ForneyLab
using LinearAlgebra
import ForneyLab: SoftFactor, @ensureVariables, generateId, addNode!, associate!,
                  averageEnergy, differentialEntropy, Interface, Variable, slug, ProbabilityDistribution,
				  unsafeLogMean, unsafeMean, unsafeCov, unsafePrecision, unsafeMeanCov,
                  ultimatePartner, region, Region, isClamped, currentInferenceAlgorithm
export LatentAutoregressiveX, LARX

"""
Description:

    A Latent Autoregressive model with eXogenous input (LARX).

    Consider states x ∈ R, with y = [x_k, x_k-1, ... x_k-M+1] and
    z = [x_k-1, ... x_k-M]. The first element of y (x_k) is produced by an inner
    product between autoregression coefficients θ and M previous states (i.e. z),
    as well as a coefficient η times an input variable u. The remainder of y is
    produced by the first M-1 elements of z.

    The node function is a Gaussian with mean-precision parameterization:

    f(y, θ, z, η, u, γ) = 𝒩(y | A(θ)z + B(η)u, V(γ)),

    where A(θ) is a state transition matrix dependent on unknown coefficients θ.
    The matrix is split into a matrix

        S = |0 .. 0
              I   0|,

    which shifts the elements of the state vector downwards and drops the
    oldest state x_k-M, and a vector s = [1 .. 0]', which can be used to
    construct a matrix with the autoregression coefficients as top row;

        sθ' = | θ_1 .. θ_M
                 ⋮       ⋮
                 0  ..  0 |

    B(η)u is a scaled linear additive control and V(γ) is a covariance
    matrix based on process precision γ.

Interfaces:

    1. y (state vector at time k; i.e. [x_k, x_k-1, ... x_k-M+1])
    2. θ (autoregression coefficients; i.e. [θ_1, ... θ_M])
    3. z (state vector at time k-1; i.e. [x_k-1, ... x_k-M])
    4. η (control coefficient)
    5. u (control / input)
    6. γ (precision)

Construction:

    LatentAutoregressiveX(y, θ, z, η, u, γ, id=:some_id)
"""

mutable struct LatentAutoregressiveX <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function LatentAutoregressiveX(y, θ, z, η, u, γ; id=generateId(LatentAutoregressiveX))
        @ensureVariables(y, z, θ, η, u, γ)
        self = new(id, Array{Interface}(undef, 6), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:y] = self.interfaces[1] = associate!(Interface(self), y)
        self.i[:z] = self.interfaces[2] = associate!(Interface(self), z)
        self.i[:θ] = self.interfaces[3] = associate!(Interface(self), θ)
        self.i[:η] = self.interfaces[4] = associate!(Interface(self), η)
        self.i[:u] = self.interfaces[5] = associate!(Interface(self), u)
        self.i[:γ] = self.interfaces[6] = associate!(Interface(self), γ)
        return self
    end
end

slug(::Type{LatentAutoregressiveX}) = "LARX"

function averageEnergy(::Type{LatentAutoregressiveX},
                       marg_y::ProbabilityDistribution{Multivariate},
                       marg_z::ProbabilityDistribution{Multivariate},
                       marg_θ::ProbabilityDistribution{Multivariate},
                       marg_η::ProbabilityDistribution{Univariate},
                       marg_u::ProbabilityDistribution{Univariate},
                       marg_γ::ProbabilityDistribution{Univariate})

    # Expectations of marginal beliefs
    my, Vy = unsafeMeanCov(marg_y)
    mz, Vz = unsafeMeanCov(marg_z)
    mθ, Vθ = unsafeMeanCov(marg_θ)
    mη, vη = unsafeMeanCov(marg_η)
    mu, vu = unsafeMeanCov(marg_u)
    mγ = unsafeMean(marg_γ)

    error("Not implemented yet")

    # Compute
    Az = mθ*mz
    Eg2 = Eg*Eg' + Jx'*Vz*Jx + Jθ'*Vθ*Jθ

    # Expand square and pre-compute terms
    sq1 = my[1]^2 + Vy[1,1]
	sq2 = my[1]*(Eg + mη*mu)
	sq3 = Eg2 + 2*Eg*mη*mu + (mη^2 + vη)*(mu + vu)

	# Compute average energy
	AE = 1/2*log(2*π) -1/2*unsafeLogMean(marg_γ) +1/2*mγ*(sq1 -2*sq2 + sq3)

    # Correction
    AE += differentialEntropy(marg_y)
    marg_y1 = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=my[1], v=Vy[1,1])
    AE -= differentialEntropy(marg_y1)

    return AE
end

max_stars_repo_name: wmkouw/LARX

max_issues_repo_name: wmkouw/LARX

max_forks_repo_name: wmkouw/LARX


text:
module QuantumOperator

using PyCall

import ..@pyfunc

#qulacs = PyNULL()
quantum_operator = PyNULL()

const quantum_operator_functions = [
    :create_quantum_operator_from_openfermion_file,
    :create_quantum_operator_from_openfermion_text,
    :create_split_quantum_operator,
]

#=
for func in quantum_operator_functions
    @eval begin
        function $(func)(args...; kwargs...)
            qulacs.quantum_operator.$(func)(args...; kwargs...)
        end
        export $(func)
    end
end
=#

for func in quantum_operator_functions
    @eval begin
        @pyfunc quantum_operator $(func)
        export $(func)
    end
end

function __init__()
    #copy!(qulacs, pyimport("qulacs"))
    copy!(quantum_operator, pyimport("qulacs.quantum_operator"))
end

end # module QuantumOperator