text:
using Pkg
Pkg.activate("/home/Morgan/FrankWolfe.jl")

using Test
using FrankWolfe
# include("/home/Morgan/FrankWolfe.jl/src/FrankWolfe.jl")
using LinearAlgebra
import SparseArrays

import FrankWolfe: compute_extreme_point, LpNormLMO, KSparseLMO
# include("/home/Morgan/FrankWolfe.jl/src/norm_oracles.jl")

import GLPK
import MathOptInterface
const MOI = MathOptInterface

import Clp
using Random
import Hypatia

@testset "Simplex LMOs" begin
    n = 6
    direction = zeros(6)
    rhs = 10 * rand()
    lmo_prob = FrankWolfe.ProbabilitySimplexOracle(rhs)
    lmo_unit = FrankWolfe.UnitSimplexOracle(rhs)
    @testset "Choosing improving direction" for idx in 1:n
        direction .= 0
        direction[idx] = -1
        res_point_prob = FrankWolfe.compute_extreme_point(lmo_prob, direction)
        res_point_unit = FrankWolfe.compute_extreme_point(lmo_unit, direction)
        for j in eachindex(res_point_prob)
            if j == idx
                @test res_point_prob[j] == res_point_unit[j] == rhs
            else
                @test res_point_prob[j] == res_point_unit[j] == 0
            end
        end
        # computing dual solutions and testing complementarity
        dual, redC = FrankWolfe.compute_dual_solution(lmo_prob, direction, res_point_prob)
        @test sum((redC .* res_point_prob)) + (dual[1] * (rhs - sum(res_point_prob))) == 0
        dual, redC = FrankWolfe.compute_dual_solution(lmo_unit, direction, res_point_unit)
        @test sum((redC .* res_point_unit)) + (dual[1] * (rhs - sum(res_point_unit))) == 0
    end
    @testset "Choosing least-degrading direction" for idx in 1:n
        # all directions worsening, must pick idx
        direction .= 2
        direction[idx] = 1
        res_point_prob = FrankWolfe.compute_extreme_point(lmo_prob, direction)
        res_point_unit = FrankWolfe.compute_extreme_point(lmo_unit, direction)
        for j in eachindex(res_point_unit)
            @test res_point_unit[j] == 0
            if j == idx
                @test res_point_prob[j] == rhs
            else
                @test res_point_prob[j] == 0
            end
        end
        # computing dual solutions and testing complementarity
        dual, redC = FrankWolfe.compute_dual_solution(lmo_prob, direction, res_point_prob)
        @test sum((redC .* res_point_prob)) + (dual[1] * (rhs - sum(res_point_prob))) == 0
        dual, redC = FrankWolfe.compute_dual_solution(lmo_unit, direction, res_point_unit)
        @test sum((redC .* res_point_unit)) + (dual[1] * (rhs - sum(res_point_unit))) == 0
    end
end

@testset "Lp-norm epigraph LMO" begin
    for n in (1, 2, 5, 10)
        τ = 5 + 3 * rand()
        # tests that the "special" p behaves like the "any" p, i.e. 2.0 and 2
        @testset "$p-norm" for p in (1, 1.0, 1.5, 2, 2.0, Inf, Inf32)
            lmo = LpNormLMO{Float64,p}(τ)
            for _ in 1:100
                c = 5 * randn(n)
                v = FrankWolfe.compute_extreme_point(lmo, c)
                @test norm(v, p) ≈ τ
            end
            c = zeros(n)
            v = FrankWolfe.compute_extreme_point(lmo, c)
            @test !any(isnan, v)
        end
        @testset "K-Norm ball $K" for K in 1:n
            lmo_ball = FrankWolfe.KNormBallLMO(K, τ)
            for _ in 1:20
                c = 5 * randn(n)
                v = FrankWolfe.compute_extreme_point(lmo_ball, c)
                v1 = FrankWolfe.compute_extreme_point(FrankWolfe.LpNormLMO{1}(τ), c)
                v_inf = FrankWolfe.compute_extreme_point(FrankWolfe.LpNormLMO{Inf}(τ / K), c)
                # K-norm is convex hull of union of the two norm epigraphs
                # => cannot do better than the best of them
                @test dot(v, c) ≈ min(dot(v1, c), dot(v_inf, c))
                # test according to original norm definition
                # norm constraint must be tight
                K_sum = 0.0
                for vi in sort!(abs.(v), rev=true)[1:K]
                    K_sum += vi
                end
                @test K_sum ≈ τ
            end
        end
    end
    # testing issue on zero direction
    for n in (1, 5)
        lmo = FrankWolfe.LpNormLMO{Float64,2}(1.0)
        x0 = FrankWolfe.compute_extreme_point(lmo, zeros(n))
        @test all(!isnan, x0)
    end
end

@testset "K-sparse polytope LMO" begin
    @testset "$n-dimension" for n in (1, 2, 10)
        τ = 5 + 3 * rand()
        for K in 1:n
            lmo = KSparseLMO(K, τ)
            x = 10 * randn(n) # dense vector
            v = compute_extreme_point(lmo, x)
            # K-sparsity
            @test count(!iszero, v) == K
            @test sum(abs.(v)) ≈ K * τ
            xsort = sort!(10 * rand(n))
            v = compute_extreme_point(lmo, xsort)
            @test all(iszero, v[1:n-K])
            @test all(abs(vi) ≈ abs(τ * sign(vi)) for vi in v[K:end])
            reverse!(xsort)
            v = compute_extreme_point(lmo, xsort)
            @test all(iszero, v[K+1:end])
            @test all(abs(vi) ≈ abs(τ * sign(vi)) for vi in v[1:K])
        end
    end
    # type stability of K-sparse polytope LMO
    lmo = KSparseLMO(3, 2.0)
    x = 10 * randn(10) # dense vector
    @inferred compute_extreme_point(lmo, x)
end

@testset "Caching on simplex LMOs" begin
    n = 6
    direction = zeros(6)
    rhs = 10 * rand()
    lmo_unit = FrankWolfe.UnitSimplexOracle(rhs)
    lmo_never_cached = FrankWolfe.SingleLastCachedLMO(lmo_unit)
    lmo_cached = FrankWolfe.SingleLastCachedLMO(lmo_unit)
    lmo_multicached = FrankWolfe.MultiCacheLMO{3}(lmo_unit)
    lmo_veccached = FrankWolfe.VectorCacheLMO(lmo_unit)
    @testset "Forcing no cache remains nothing" for idx in 1:n
        direction .= 0
        direction[idx] = -1
        res_point_unit = FrankWolfe.compute_extreme_point(lmo_unit, direction)
        res_point_cached = FrankWolfe.compute_extreme_point(lmo_cached, direction, threshold=0)
        res_point_cached_multi =
            FrankWolfe.compute_extreme_point(lmo_multicached, direction, threshold=-1000)
        res_point_cached_vec =
            FrankWolfe.compute_extreme_point(lmo_veccached, direction, threshold=-1000)
        res_point_never_cached =
            FrankWolfe.compute_extreme_point(lmo_cached, direction, store_cache=false)
        @test res_point_never_cached == res_point_unit
        @test lmo_never_cached.last_vertex === nothing
        @test length(lmo_never_cached) == 0
        empty!(lmo_never_cached)
        @test lmo_cached.last_vertex !== nothing
        @test length(lmo_cached) == 1
        @test count(!isnothing, lmo_multicached.vertices) == min(3, idx) == length(lmo_multicached)
        @test length(lmo_veccached.vertices) == idx == length(lmo_veccached)
        # we set the cache at least at the first iteration
        if idx == 1
            @test lmo_cached.last_vertex == res_point_unit
        end
        # whatever the iteration, last vertex is always the one returned
        @test lmo_cached.last_vertex == res_point_cached
    end
    empty!(lmo_multicached)
    @test length(lmo_multicached) == 0
    empty!(lmo_veccached)
    @test length(lmo_veccached) == 0
end

function _is_doubly_stochastic(m)
    for col in eachcol(m)
        @test sum(col) == 1
    end
    for row in eachrow(m)
        @test sum(row) == 1
    end
end

@testset "Birkhoff polytope" begin
    Random.seed!(42)
    lmo = FrankWolfe.BirkhoffPolytopeLMO()
    for n in (1, 2, 10)
        cost = rand(n, n)
        res = FrankWolfe.compute_extreme_point(lmo, cost)
        _is_doubly_stochastic(res)
        res2 = FrankWolfe.compute_extreme_point(lmo, vec(cost))
        @test res ≈ res2
    end
    cost_mat = [
        2 3 3
        3 2 3
        3 3 2
    ]
    res = FrankWolfe.compute_extreme_point(lmo, cost_mat)
    @test res == I
    @test sum(cost_mat .* res) == 6
    cost_mat = [
        3 2 3
        2 3 3
        3 3 2
    ]
    res = FrankWolfe.compute_extreme_point(lmo, cost_mat)
    @test sum(cost_mat .* res) == 6
    @test res == Bool[
        0 1 0
        1 0 0
        0 0 1
    ]
    cost_mat = [
        3 2 3
        3 3 2
        2 3 3
    ]
    res = FrankWolfe.compute_extreme_point(lmo, cost_mat)
    @test dot(cost_mat, res) == 6
    @test res == Bool[
        0 1 0
        0 0 1
        1 0 0
    ]
end

@testset "Matrix completion and nuclear norm" begin
    nfeat = 50
    nobs = 100
    r = 5
    Xreal = Matrix{Float64}(undef, nobs, nfeat)
    X_gen_cols = randn(nfeat, r)
    X_gen_rows = randn(r, nobs)
    svals = 100 * rand(r)
    for i in 1:nobs
        for j in 1:nfeat
            Xreal[i, j] = sum(X_gen_cols[j, k] * X_gen_rows[k, i] * svals[k] for k in 1:r)
        end
    end
    @test rank(Xreal) == r
    missing_entries = unique!([(rand(1:nobs), rand(1:nfeat)) for _ in 1:1000])
    f(X) =
        0.5 *
        sum((X[i, j] - Xreal[i, j])^2 for i in 1:nobs, j in 1:nfeat if (i, j) ∉ missing_entries)
    function grad!(storage, X)
        storage .= 0
        for i in 1:nobs
            for j in 1:nfeat
                if (i, j) ∉ missing_entries
                    storage[i, j] = X[i, j] - Xreal[i, j]
                end
            end
        end
        return nothing
    end
    # TODO value of radius?
    lmo = FrankWolfe.NuclearNormLMO(sum(svdvals(Xreal)))
    x0 = FrankWolfe.compute_extreme_point(lmo, zero(Xreal))
    gradient = similar(x0)
    grad!(gradient, x0)
    v0 = FrankWolfe.compute_extreme_point(lmo, gradient)
    @test dot(v0 - x0, gradient) < 0
    xfin, vmin, _ = FrankWolfe.frank_wolfe(
        f,
        grad!,
        lmo,
        x0;
        epsilon=1e-6,
        max_iteration=400,
        print_iter=100,
        trajectory=false,
        verbose=false,
        linesearch_tol=1e-7,
        line_search=FrankWolfe.Backtracking(),
        emphasis=FrankWolfe.memory,
    )
    @test 1 - (f(x0) - f(xfin)) / f(x0) < 1e-3
    svals_fin = svdvals(xfin)
    @test sum(svals_fin[r+1:end]) / sum(svals_fin) ≤ 2e-2
    xfin, vmin, _ = FrankWolfe.lazified_conditional_gradient(
        f,
        grad!,
        lmo,
        x0;
        epsilon=1e-6,
        max_iteration=400,
        print_iter=100,
        trajectory=false,
        verbose=false,
        linesearch_tol=1e-7,
        line_search=FrankWolfe.Backtracking(),
        emphasis=FrankWolfe.memory,
    )
end

@testset "Spectral norms" begin
    Random.seed!(42)
    o = Hypatia.Optimizer()
    MOI.set(o, MOI.Silent(), true)
    optimizer = MOI.Bridges.full_bridge_optimizer(
            MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            o,
        ),
        Float64,
    )
    radius = 5.0
    @testset "Spectraplex $n" for n in (2, 10)
        lmo = FrankWolfe.SpectraplexLMO(radius, n)
        direction = Matrix{Float64}(undef, n, n)
        lmo_moi = FrankWolfe.convert_mathopt(lmo, optimizer; side_dimension=n)
        for _ in 1:10
            Random.randn!(direction)
            v = @inferred FrankWolfe.compute_extreme_point(lmo, direction)
            vsym = @inferred FrankWolfe.compute_extreme_point(lmo, direction + direction')
            vsym2 = @inferred FrankWolfe.compute_extreme_point(lmo, Symmetric(direction + direction'))
            @test v ≈ vsym atol=1e-6
            @test v ≈ vsym2 atol=1e-6
            @testset "Vertex properties" begin
                eigen_v = eigen(Matrix(v))
                @test eigmax(Matrix(v)) ≈ radius
                @test norm(eigen_v.values[1:end-1]) ≈ 0 atol=1e-7
                # u can be sqrt(r) * vec or -sqrt(r) * vec
                case_pos = ≈(norm(eigen_v.vectors[:,n] * sqrt(eigen_v.values[n]) - v.u), 0, atol=1e-9)
                case_neg = ≈(norm(eigen_v.vectors[:,n] * sqrt(eigen_v.values[n]) + v.u), 0, atol=1e-9)
                @test case_pos || case_neg
            end
            @testset "Comparison with SDP solution" begin
                v_moi = FrankWolfe.compute_extreme_point(lmo_moi, direction)
                @test norm(v - v_moi) <= 1e-6
            end
        end
    end
    @testset "Unit spectrahedron $n" for n in (2, 10)
        lmo = FrankWolfe.UnitSpectrahedronLMO(radius, n)
        direction = Matrix{Float64}(undef, n, n)
        lmo_moi = FrankWolfe.convert_mathopt(lmo, optimizer; side_dimension=n)
        direction_sym = similar(direction)
        for _ in 1:10
            Random.randn!(direction)
            @. direction_sym = direction + direction'
            v = @inferred FrankWolfe.compute_extreme_point(lmo, direction)
            vsym = @inferred FrankWolfe.compute_extreme_point(lmo, direction_sym)
            @test v ≈ vsym atol=1e-6
            @testset "Vertex properties" begin
                emin = eigmin(direction_sym)
                if emin ≥ 0
                    @test norm(Matrix(v)) ≈ 0
                else
                    eigen_v = eigen(Matrix(v))
                    @test eigmax(Matrix(v)) ≈ radius
                    @test norm(eigen_v.values[1:end-1]) ≈ 0 atol=1e-7
                    # u can be sqrt(r) * vec or -sqrt(r) * vec
                    case_pos = ≈(norm(eigen_v.vectors[:,n] * sqrt(eigen_v.values[n]) - v.u), 0, atol=1e-9)
                    case_neg = ≈(norm(eigen_v.vectors[:,n] * sqrt(eigen_v.values[n]) + v.u), 0, atol=1e-9)
                    @test case_pos || case_neg
                    # make direction PSD
                    direction_sym += 1.1 * abs(emin) * I
                    @assert isposdef(direction_sym)
                    v2 = FrankWolfe.compute_extreme_point(lmo, direction_sym)
                    @test norm(Matrix(v2)) ≈ 0
                end
            end
            @testset "Comparison with SDP solution" begin
                v_moi = FrankWolfe.compute_extreme_point(lmo_moi, direction)
                @test norm(v - v_moi) <= 1e-6
                # forcing PSD direction to test 0 matrix case
                @. direction_sym = direction + direction'
                direction_sym += 1.1 * abs(eigmin(direction_sym)) * I
                @assert isposdef(direction_sym)
                v_moi2 = FrankWolfe.compute_extreme_point(lmo_moi, direction_sym)
                v_lmo2 = FrankWolfe.compute_extreme_point(lmo_moi, direction_sym)
                @test norm(v_moi2 - v_lmo2) <= 1e-6
                @test norm(v_moi2) <= 1e-6
            end
        end
    end
end

@testset "MOI oracle consistency" begin
    Random.seed!(42)
    o = GLPK.Optimizer()
    MOI.set(o, MOI.Silent(), true)
    @testset "MOI oracle consistent with unit simplex" for n in (1, 2, 10)
        MOI.empty!(o)
        x = MOI.add_variables(o, n)
        for xi in x
            MOI.add_constraint(o, xi, MOI.Interval(0.0, 1.0))
        end
        MOI.add_constraint(
            o,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
            MOI.LessThan(1.0),
        )
        lmo = FrankWolfe.MathOptLMO(o)
        lmo_ref = FrankWolfe.UnitSimplexOracle(1.0)
        lmo_moi_ref = FrankWolfe.convert_mathopt(lmo_ref, GLPK.Optimizer(), dimension=n)
        direction = Vector{Float64}(undef, n)
        for _ in 1:5
            Random.randn!(direction)
            vref = compute_extreme_point(lmo_ref, direction)
            v = compute_extreme_point(lmo, direction)
            v_moi = compute_extreme_point(lmo_moi_ref, direction)
            @test vref ≈ v
            @test vref ≈ v_moi
        end
    end
    @testset "MOI consistent probability simplex" for n in (1, 2, 10)
        MOI.empty!(o)
        x = MOI.add_variables(o, n)
        for xi in x
            MOI.add_constraint(o, xi, MOI.Interval(0.0, 1.0))
        end
        MOI.add_constraint(
            o,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
            MOI.EqualTo(1.0),
        )
        lmo = FrankWolfe.MathOptLMO(o)
        lmo_ref = FrankWolfe.ProbabilitySimplexOracle(1.0)
        lmo_moi_ref = FrankWolfe.convert_mathopt(lmo_ref, GLPK.Optimizer(), dimension=n)
        direction = Vector{Float64}(undef, n)
        for _ in 1:5
            Random.randn!(direction)
            vref = compute_extreme_point(lmo_ref, direction)
            v = compute_extreme_point(lmo, direction)
            v_moi = compute_extreme_point(lmo_moi_ref, direction)
            @test vref ≈ v
            @test vref ≈ v_moi
        end
    end
    @testset "Direction with coefficients" begin
        n = 5
        MOI.empty!(o)
        x = MOI.add_variables(o, n)
        for xi in x
            MOI.add_constraint(o, xi, MOI.Interval(0.0, 1.0))
        end
        MOI.add_constraint(
            o,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
            MOI.EqualTo(1.0),
        )
        lmo = FrankWolfe.MathOptLMO(o)
        direction = [MOI.ScalarAffineTerm(-2.0i, x[i]) for i in 2:3]
        v = compute_extreme_point(lmo, direction)
        @test v ≈ [0, 1]
    end
    @testset "Non-settable optimizer with cache" begin
        n = 5
        o = MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            Clp.Optimizer(),
        )
        MOI.set(o, MOI.Silent(), true)
        x = MOI.add_variables(o, 5)
        for xi in x
            MOI.add_constraint(o, xi, MOI.Interval(0.0, 1.0))
        end
        MOI.add_constraint(
            o,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
            MOI.EqualTo(1.0),
        )
        lmo = FrankWolfe.MathOptLMO(o)
        lmo_ref = FrankWolfe.ProbabilitySimplexOracle(1.0)
        direction = Vector{Float64}(undef, n)
        for _ in 1:5
            Random.randn!(direction)
            vref = compute_extreme_point(lmo_ref, direction)
            v = compute_extreme_point(lmo, direction)
            @test vref ≈ v
        end
    end
    @testset "Nuclear norm" for n in (5, 10)
        o = Hypatia.Optimizer()
        MOI.set(o, MOI.Silent(), true)
        optimizer = MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            o,
        )
        MOI.set(optimizer, MOI.Silent(), true)
        nrows = 3n
        ncols = n
        direction = Matrix{Float64}(undef, nrows, ncols)
        τ = 10.0
        lmo = FrankWolfe.NuclearNormLMO(τ)
        lmo_moi =
            FrankWolfe.convert_mathopt(lmo, optimizer, row_dimension=nrows, col_dimension=ncols)
        for _ in 1:10
            randn!(direction)
            v_r = FrankWolfe.compute_extreme_point(lmo, direction)
            flattened = collect(vec(direction))
            push!(flattened, 0)
            v_moi = FrankWolfe.compute_extreme_point(lmo_moi, flattened)
            if v_moi === nothing
                # ignore non-terminating MOI solver results
                continue
            end
            v_moi_mat = reshape(v_moi[1:end-1], nrows, ncols)
            @test v_r ≈ v_moi_mat rtol = 1e-2
        end
    end
end

@testset "MOI oracle on Birkhoff polytope" begin
    o = GLPK.Optimizer()
    o_ref = GLPK.Optimizer()
    for n in (1, 2, 10)
        MOI.empty!(o)
        (x, _) = MOI.add_constrained_variables(o, fill(MOI.Interval(0.0, 1.0), n * n))
        xmat = reshape(x, n, n)
        for idx in 1:n
            # column constraint
            MOI.add_constraint(
                o,
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), xmat[:, idx]), 0.0),
                MOI.EqualTo(1.0),
            )
            # row constraint
            MOI.add_constraint(
                o,
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), xmat[idx, :]), 0.0),
                MOI.EqualTo(1.0),
            )
        end
        direction_vec = Vector{Float64}(undef, n * n)
        lmo_bkf = FrankWolfe.BirkhoffPolytopeLMO()
        lmo_moi = FrankWolfe.MathOptLMO(o)
        lmo_moi_ref = FrankWolfe.convert_mathopt(lmo_bkf, o_ref, dimension=n)
        for _ in 1:10
            randn!(direction_vec)
            direction_mat = reshape(direction_vec, n, n)
            v_moi = FrankWolfe.compute_extreme_point(lmo_moi, direction_vec)
            v_moi_mat = reshape(v_moi, n, n)
            v_bfk = FrankWolfe.compute_extreme_point(lmo_bkf, direction_mat)
            @test all(isapprox.(v_moi_mat, v_bfk, atol=1e-4))
            v_moi_mat2 = FrankWolfe.compute_extreme_point(lmo_moi, direction_mat)
            @test all(isapprox.(v_moi_mat2, v_moi_mat))
        end
    end
end

@testset "MOI oracle and KSparseLMO" begin
    o_base = GLPK.Optimizer()
    cached = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        o_base,
    )
    o = MOI.Bridges.full_bridge_optimizer(cached, Float64)
    o_ref = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            GLPK.Optimizer(),
        ),
        Float64,
    )
    for n in (1, 2, 5, 10)
        for K in 1:3:n
            τ = 10 * rand()
            MOI.empty!(o)
            x = MOI.add_variables(o, n)
            tinf = MOI.add_variable(o)
            MOI.add_constraint(o, MOI.VectorOfVariables([tinf; x]), MOI.NormInfinityCone(n + 1))
            MOI.add_constraint(o, tinf, MOI.LessThan(τ))
            t1 = MOI.add_variable(o)
            MOI.add_constraint(o, MOI.VectorOfVariables([t1; x]), MOI.NormOneCone(n + 1))
            MOI.add_constraint(o, t1, MOI.LessThan(τ * K))
            direction = Vector{Float64}(undef, n)
            lmo_moi = FrankWolfe.MathOptLMO(o)
            lmo_ksp = FrankWolfe.KSparseLMO(K, τ)
            lmo_moi_convert = FrankWolfe.convert_mathopt(lmo_ksp, o_ref, dimension=n)
            for _ in 1:20
                randn!(direction)
                v_moi =
                    FrankWolfe.compute_extreme_point(lmo_moi, MOI.ScalarAffineTerm.(direction, x))
                v_ksp = FrankWolfe.compute_extreme_point(lmo_ksp, direction)
                v_moi_conv = FrankWolfe.compute_extreme_point(
                    lmo_moi_convert,
                    MOI.ScalarAffineTerm.(direction, x),
                )
                for idx in eachindex(v_moi)
                    @test isapprox(v_moi[idx], v_ksp[idx], atol=1e-4)
                    @test isapprox(v_moi_conv[idx], v_ksp[idx], atol=1e-4)
                end
            end
            # verifying absence of a bug
            if n == 5
                direction .= (
                    -0.07020498519126772,
                    0.4298929981513661,
                    -0.8678437699266819,
                    -0.08899938054920563,
                    1.160622285477465,
                )
                v_moi =
                    FrankWolfe.compute_extreme_point(lmo_moi, MOI.ScalarAffineTerm.(direction, x))
                v_ksp = FrankWolfe.compute_extreme_point(lmo_ksp, direction)
                for idx in eachindex(v_moi)
                    @test isapprox(v_moi[idx], v_ksp[idx], atol=1e-4)
                end
            end
        end
    end
end

@testset "Product LMO" begin
    lmo = FrankWolfe.ProductLMO(FrankWolfe.LpNormLMO{Inf}(3.0), FrankWolfe.LpNormLMO{1}(2.0))
    dinf = randn(10)
    d1 = randn(5)
    vtup = FrankWolfe.compute_extreme_point(lmo, (dinf, d1))
    @test length(vtup) == 2
    (vinf, v1) = vtup
    @test sum(abs, vinf) ≈ 10 * 3.0
    @test sum(!iszero, v1) == 1

    vvec = FrankWolfe.compute_extreme_point(lmo, [dinf; d1]; direction_indices=(1:10, 11:15))
    @test vvec ≈ [vinf; v1]
end

@testset "Scaled L-1 norm polytopes" begin
    lmo = FrankWolfe.ScaledBoundL1NormBall(-ones(10), ones(10))
    # equivalent to LMO
    lmo_ref = FrankWolfe.LpNormLMO{1}(1)
    # all coordinates shifted up
    lmo_shifted = FrankWolfe.ScaledBoundL1NormBall(zeros(10), 2 * ones(10))
    lmo_scaled = FrankWolfe.ScaledBoundL1NormBall(-2 * ones(10), 2 * ones(10))
    for _ in 1:100
        d = randn(10)
        v = FrankWolfe.compute_extreme_point(lmo, d)
        vref = FrankWolfe.compute_extreme_point(lmo_ref, d)
        @test v ≈ vref
        vshift = FrankWolfe.compute_extreme_point(lmo_shifted, d)
        @test v .+ 1 ≈ vshift
        v2 = FrankWolfe.compute_extreme_point(lmo_scaled, d)
        @test v2 ≈ 2v
    end
    d = zeros(10)
    v = FrankWolfe.compute_extreme_point(lmo, d)
    vref = FrankWolfe.compute_extreme_point(lmo_ref, d)
    @test v ≈ vref
    @test norm(v) == 1
    # non-uniform scaling
    # validates bugfix
    lmo_nonunif = FrankWolfe.ScaledBoundL1NormBall([-1.0,-1.0], [3.0,1.0])
    direction = [-0.8272727272727383, -0.977272727272718]
    v = FrankWolfe.compute_extreme_point(lmo_nonunif, direction)
    @test v ≈ [3, 0]
    v = FrankWolfe.compute_extreme_point(lmo_nonunif, -direction)
    @test v ≈ [-1, 0]
end

@testset "Scaled L-inf norm polytopes" begin
    # tests ScaledBoundLInfNormBall for the standard hypercube, a shifted one, and a scaled one
    lmo = FrankWolfe.ScaledBoundLInfNormBall(-ones(10), ones(10))
    lmo_ref = FrankWolfe.LpNormLMO{Inf}(1)
    lmo_shifted = FrankWolfe.ScaledBoundLInfNormBall(zeros(10), 2 * ones(10))
    lmo_scaled = FrankWolfe.ScaledBoundLInfNormBall(-2 * ones(10), 2 * ones(10))
    bounds = collect(1.0:10)
    # tests another ScaledBoundLInfNormBall with unequal bounds against a MOI optimizer
    lmo_scaled_unequally = FrankWolfe.ScaledBoundLInfNormBall(-bounds, bounds)
    o = GLPK.Optimizer()
    MOI.set(o, MOI.Silent(), true)
    x = MOI.add_variables(o, 10)
    MOI.add_constraint.(o, x, MOI.GreaterThan.(-bounds))
    MOI.add_constraint.(o, x, MOI.LessThan.(bounds))
    scaled_unequally_opt = FrankWolfe.MathOptLMO(o)
    for _ in 1:100
        d = randn(10)
        v = FrankWolfe.compute_extreme_point(lmo, d)
        vref = FrankWolfe.compute_extreme_point(lmo_ref, d)
        @test v ≈ vref
        vshift = FrankWolfe.compute_extreme_point(lmo_shifted, d)
        @test v .+ 1 ≈ vshift
        v2 = FrankWolfe.compute_extreme_point(lmo_scaled, d)
        @test v2 ≈ 2v
        v3 = FrankWolfe.compute_extreme_point(lmo_scaled_unequally, d)
        v3_test = compute_extreme_point(scaled_unequally_opt, d)
        @test v3 ≈ v3_test
    end
    d = zeros(10)
    v = FrankWolfe.compute_extreme_point(lmo, d)
    vref = FrankWolfe.compute_extreme_point(lmo_ref, d)
    @test v ≈ vref
    @test norm(v, Inf) == 1
end


text:
# abstract type AbstractOrderedVector{T} <: AbstractVector{T} end

# struct OrderedVector{T} <: AbstractOrderedVector{T}
#     data::Vector{T}
# end

# as(::Type{OrderedVector}, transformation::TV.AbstractTransform, dim::Int) =
#     Ordered(transformation, dim)

using MeasureTheory

export Ordered
struct Ordered{T <: TV.AbstractTransform} <: TV.VectorTransform
    transformation::T
    dim::Int
end

TV.dimension(t::Ordered) = t.dim

addlogjac(ℓ, Δℓ) = ℓ + Δℓ
addlogjac(::TV.NoLogJac, Δℓ) = TV.NoLogJac()

using MappedArrays

bounds(t::TV.ShiftedExp{true}) = (t.shift, TV.∞)
bounds(t::TV.ShiftedExp{false}) = (-TV.∞, t.shift)
bounds(t::TV.ScaledShiftedLogistic) = (t.shift, t.scale + t.shift)
bounds(::TV.Identity) = (-TV.∞, TV.∞)

const OrderedΔx = -8.0

# See https://mc-stan.org/docs/2_27/reference-manual/ordered-vector.html
function TV.transform_with(flag::TV.LogJacFlag, t::Ordered, x, index::T) where {T}
    transformation, len = (t.transformation, t.dim)
    @assert dimension(transformation) == 1
    y = similar(x, len)

    (lo,hi) = bounds(transformation)

    x = mappedarray(xj -> xj + OrderedΔx, x)

    @inbounds (y[1], ℓ, _) = TV.transform_with(flag, as(Real, lo, hi), x, index)
    index += 1

    @inbounds for i in 2:len
        (y[i], Δℓ, _) =  TV.transform_with(flag, as(Real, y[i-1], hi), x, index)
        ℓ = addlogjac(ℓ, Δℓ)
        index += 1
    end

    return (y, ℓ, index)
end

TV.inverse_eltype(t::Ordered, y::AbstractVector) = TV.extended_eltype(y)

Ordered(n::Int) = Ordered(asℝ, n)

function TV.inverse_at!(x::AbstractVector, index, t::Ordered, y::AbstractVector)
    (lo,hi) = bounds(t.transformation)

    @inbounds x[index] = inverse(as(Real, lo, hi), y[1]) - OrderedΔx
    index += 1

    @inbounds for i in 2:length(y)
        x[index] = inverse(as(Real, y[i-1], hi), y[i]) - OrderedΔx
        index += 1
    end
    return index
end

export Sorted
struct Sorted{M} <: AbstractMeasure
    μ::M
    n::Int
end

logdensity(s::Sorted, x) = logdensity(s.μ ^ s.n, x)

TV.as(s::Sorted) = Ordered(as(s.μ), s.n)

function Random.rand!(rng::AbstractRNG, d::Sorted, x::AbstractArray)
    rand!(rng, d.μ ^ d.n, x)
    sort!(x)
    return x
end

function Base.rand(rng::AbstractRNG, T::Type, d::Sorted)
    # TODO: Use `sampletype` for this
    elT = typeof(rand(rng, T, d.μ))
    x = Vector{elT}(undef, d.n)
    rand!(rng, d, x)
end

# logdensity(d, rand(d))

# TV.transform_with(TV.LogJac(), Ordered(asℝ, 4), zeros(4), 1)
# TV.transform_with(TV.LogJac(), Ordered(asℝ, 4), randn(4), 1)

# d = Pushforward(Ordered(10), Normal()^10, false)
# logdensity(Lebesgue(ℝ)^10, rand(d))

ext: jl


text:
const sources = [
    "im2col.cpp",
    "pooling.cpp"]

if Sys.iswindows()
    builddir = dirname(Base.source_path())
    println("Build directory is $(builddir)")

    import WinRPM

    println("Installing gcc-c++.")
    WinRPM.install("gcc-c++"; yes=true)
    WinRPM.install("gcc"; yes=true)
    WinRPM.install("headers"; yes=true)

    gpp = Pkg.dir("WinRPM","deps","usr","x86_64-w64-mingw32","sys-root","mingw","bin","g++")
    RPMbindir = Pkg.dir("WinRPM","deps","usr","x86_64-w64-mingw32","sys-root","mingw","bin")
    incdir = Pkg.dir("WinRPM","deps","usr","x86_64-w64-mingw32","sys-root","mingw","include")

    push!(Base.Libdl.DL_LOAD_PATH,RPMbindir)
    ENV["PATH"] = ENV["PATH"] * ";" * RPMbindir

    run(`$gpp --version`)
    cmd = `$gpp -Wall -O3 -fPIC -shared -I $incdir -o libmerlin.dll $sources`
elseif Sys.isapple()
    cmd = `g++ -Wall -O3 -fPIC -shared -o libmerlin.so $sources`
elseif Sys.islinux()
    cmd = `g++ -Wall -O3 -fPIC -shared -o libmerlin.so $sources`
else
    throw("Unknown OS.")
end

println(cmd)
run(cmd)
println()
println("Build success.")


text:
using Test
using MARY_fRG.Triangulated
using MARY_fRG.Fermi
using MARY_fRG.Fermi.Surface
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Kagome

#@testset "正方格子的费米面" begin
#    qs = common_square_lattice(0.)
#    ltris, ladjs = split_square(qs.brillouin, 100)
#    sf = const_energy_triangle(ltris, 0.2, qs.dispersion[1])
#    plt = draw_lines(qs.brillouin.edges)
#    draw_polygon!(plt, sf; pcidx=ones(Int64, length(sf)))
#    #draw_polygon!(plt, ltris)
#    sf2 = const_energy_triangle(ltris, -0.2, qs.dispersion[1])
#    draw_polygon!(plt, sf2; pcidx=2ones(Int64, length(sf2)))
#    @savepng plt "surfacet1"
#end
#
#
#@testset "三角格子的费米面" begin
#    ts = common_triangle_lattice(0.)
#    ltris, ladjs = split_hexagon(ts.brillouin, 100)
#    sf = const_energy_triangle(ltris, 0., ts.dispersion[1])
#    plt = draw_lines(ts.brillouin.edges)
#    draw_polygon!(plt, sf; pcidx=ones(Int64, length(sf)))
#    sf2 = const_energy_triangle(ltris, 2.0, ts.dispersion[1])
#    draw_polygon!(plt, sf2; pcidx=2ones(Int64, length(sf2)))
#    @savepng plt "surfacet2"
#end

@testset "kagome格子的费米面" begin
    ks = upperband_kagome_lattice(0.)
    ltris, ladjs = split_hexagon(ks.brillouin, 100)
    sf = const_energy_triangle(ltris, 0., ks.dispersion[1])
    plt = draw_lines(ks.brillouin.edges)
    draw_polygon!(plt, sf; pcidx=ones(Int64, length(sf)))
    sf2 = const_energy_triangle(ltris, 2.0, ks.dispersion[1])
    draw_polygon!(plt, sf2; pcidx=2ones(Int64, length(sf2)))
    @savepng plt "surfacet3"
end


text:
"""
    ImplicitBernsteinForm{T, PT<:AbstractPolynomialLike{T},
                          D, N, TC, NC, SC} <: AbstractBernsteinForm{T}

Type defining an implicit Bernstein form of a polynomial over a hyperrectangular
domain.

### Fields

- `pol`         -- polynomial in one or several variables
- `dom`         -- hyperrectangular set of type `IntervalBox` representing the
                   domain of the expansion
- `numvars`     -- number of variables in the polynomial
- `coeffs`      -- vector of arrays that holds the Bernstein coefficients implicitly
- `numcoeffs`   -- integer holding the number of Bernstein coefficients that are
                   implicitly defined
- `degs`        -- vector of degrees of the Bernstein polynomial for each variable

### Examples
"""
struct ImplicitBernsteinForm{T, PT<:AbstractPolynomialLike{T},
                             D, N, TC, NC, SC, VD<:AbstractVector{Int}} <: AbstractBernsteinForm{T}
    pol::PT
    dom::IntervalBox{D, N}
    numvars::Int
    coeffs::VectorOfArray{TC, NC, SC}
    numcoeffs::Int
    degs::VD
end

# univariate case is not special cased
function ImplicitBernsteinForm(pol::PT, dom::Interval, args...) where {T, PT<:AbstractPolynomialLike}
    ImplicitBernsteinForm(pol, IntervalBox(dom), args...)
end

#=
# constructor from
function ImplicitBernsteinForm(pol::PT, dom::IntervalBox, args...) where {T, PT<:AbstractPolynomialLike}
    numvars = nvariables(pol)

end
=#

function getcoeff(bf::ImplicitBernsteinForm, ids::NTuple{D, Int}) where {D}
    return prod(bf.coeffs[k, i] for (i, k) in enumerate(ids))
end

function assemble(array::Vector{VN}) where {N, VN<:AbstractVector{N}}
    n = length(array)
    if n == 1
        return hcat(first(array))
    elseif n == 2
        return array[1] * array[2]'
    else
        error("not implemented yet")
    end
    # preallocate grid ...
end

#=
# mutating versions
function _assemble!(..)
    ...
end
=#

max_stars_repo_name: JuliaReach/BernsteinExpansions

max_issues_repo_name: JuliaReach/BernsteinExpansions

max_forks_repo_name: JuliaReach/BernsteinExpansions


text:
module GCloud

using JSON

import Base: show, addprocs

export external_ip, instances, create_instance, delete_instance, instance

function input(prompt::String="")
    print(prompt)
    chomp(readline())
end

## ------------- ##
#- Machine Types -#
## ------------- ##

## NOTE: This data was taken from https://cloud.google.com/compute/pricing
#        on 2014-10-14 00:50:09

immutable MachineType
    configuration::String
    virtual_cores::Int
    memory::Real
    GCEU::Union(Real, String)
    local_disk::Real
    price_per_hour::Real
end

function show(io::IO, mt::MachineType)
    kind = split(mt.configuration, '-')[2]
    msg = "$kind machine ($(mt.virtual_cores) cores, $(mt.memory) GB memory, "
    msg *= "\$$(mt.price_per_hour)/hr.)"
    print(io, msg)
end

const machine_types = Dict{String, MachineType}()

# standard
machine_types["n1-standard-1"] = MachineType("n1-standard-1", 1, 3.75, 2.75, 0,
                                             0.063)
machine_types["n1-standard-2"] = MachineType("n1-standard-2", 2, 7.50, 5.50, 0,
                                             0.126)
machine_types["n1-standard-8"] = MachineType("n1-standard-8", 8, 30, 22, 0,
                                             0.504)
machine_types["n1-standard-4"] = MachineType("n1-standard-4", 4, 15, 11, 0,
                                             0.252)
machine_types["n1-standard-16"] = MachineType("n1-standard-16", 16, 60, 44, 0,
                                              1.008)

# bursting
machine_types["f1-micro"] = MachineType("f1-micro", 1, 0.60,
                                        "Shared CPU, not guaranteed", 0, 0.012)
machine_types["g1-small"] = MachineType("g1-small", 1, 1.70, 1.38, 0, 0.032)

# high-memory
machine_types["n1-highmem-2"] = MachineType("n1-highmem-2", 2, 13, 5.50, 0,
                                            0.148)
machine_types["n1-highmem-4"] = MachineType("n1-highmem-4", 4, 26, 11, 0,
                                            0.296)
machine_types["n1-highmem-8"] = MachineType("n1-highmem-8", 8, 52, 22, 0,
                                            0.592)
machine_types["n1-highmem-16"] = MachineType("n1-highmem-16", 16, 104, 44, 0,
                                             1.184)

# high-cpu
machine_types["n1-highcpu-2"] = MachineType("n1-highcpu-2", 2, 1.80,  5.50, 0,
                                            0.080)
machine_types["n1-highcpu-4"] = MachineType("n1-highcpu-4", 4, 3.60, 11, 0,
                                            0.160)
machine_types["n1-highcpu-8"] = MachineType("n1-highcpu-8", 8, 7.20, 22, 0,
                                            0.320)
machine_types["n1-highcpu-16"] = MachineType("n1-highcpu-16", 16, 14.40, 44, 0,
                                             0.640)

## --------- ##
#- Instances -#
## --------- ##

immutable Instance
    name::String
    zone::String
    kind::MachineType
    external_ip::String
    status::String
end

function show(io::IO, i::Instance)
    kind = split(i.kind.configuration, '-')[2]
    msg = "GCloud compute Instance\n"
    msg *= "  * name: $(i.name)\n"
    msg *= "  * zone: $(i.zone)\n"
    msg *= "  * machine type: $(i.kind)\n"
    msg *= "  * external ip address: $(i.external_ip)\n"
    msg *= "  * status: $(i.status)"
    print(io, msg)
end

function instances()
    data = JSON.parse(readchomp(`gcloud compute instances list --format json`))

    out = Array(Instance, length(data))
    for (i, machine) in enumerate(data)
        out[i] = Instance(machine["name"],
                          machine["zone"],
                          machine_types[machine["machineType"]],
                          external_ip(machine),
                          machine["status"])
    end
    return out
end

function instance(name::String)
    all_instances = instances()

    for i in all_instances
        if i.name == name
            return i
        end
    end

    error("Instance with name $name not found")
end

function create_instance()

    name = input("Instance name: ")
    # for now let google handle the rest with their own menus (will just have
    # to pick zone)
    run(`gcloud compute instances create $name`)
end

function create_instance(name::String,
                         m_type::MachineType,
                         zone::String;
                         options...)
    # Parse options

    opts1 = isempty(options) ? "" :
                    join(map(x->"--$(x[1]) $(x[2])", options), " ")
    opts2 = isempty(unary_options) ? "" : join(unary_options, " ")
    opts = join(Any[opts1, opts2], " ")

    # construct command
    cmd = `gcloud compute instances create $name --zone $zone
           --machine-type $(m_type.configuration)
           --network "default"
           --maintenance-policy "MIGRATE"
           --scopes 'https://www.googleapis.com/auth/devstorage.read_only'
           $opts
           `

    return cmd
    # run command
    # run(cmd)

    # return instance(name)
end

function create_instance(name::String,
                         m_type::String,
                         zone::String,
                         project::String,
                         image::String)
    # construct command
    cmd = `gcloud compute instances create $name --zone $zone
           --machine-type $m_type
           --network 'default'
           --maintenance-policy 'MIGRATE'
           --scopes 'https://www.googleapis.com/auth/devstorage.read_only'
           --image $image
           --project $project
           `
    # run command
    run(cmd)

    return instance(name)
end

function create_instance{T <: String}(name::String,
                                      m_type::String,
                                      zone::String,
                                      options::Dict{T, T}=Dict{T, T}(),
                                      unary_options::Array{T}=T[])
    m = machine_types[m_type]
    return create_instance(name, m, zone, options, unary_options)
end

delete_instance(name::String) = delete(filter(x->x.name==name, instances())[1])

function delete(instance::Instance)
    run(`gcloud compute instances delete $(instance.name)
         --zone $(instance.zone)`)
end

# accepts machine dict from parsed JSON
external_ip(d::Dict) = d["networkInterfaces"][1]["accessConfigs"][1]["natIP"]

external_ip(i::Instance) = i.external_ip

function external_ip(name::String)
    data = JSON.parse(readchomp(`gcloud compute instances list --format json`))

    # filter data to find data for instance with given name
    data = filter(x->x["name"] == name, data)

    if isempty(data)
        throw(ArgumentError("no instance with name: $name"))
    end

    data[1]["networkInterfaces"][1]["accessConfigs"][1]["natIP"]
end

function addprocs(i::Instance, dir::String,
                  username::String=readchomp(`whoami`); n::Int=1)
    ip = external_ip(i)
    # for some reason I have to do it in a loop or ssh borks with high n
    for i=1:n
        addprocs(["$(username)@$(ip)"], tunnel=true, dir=dir)
    end
    nothing
end

## ------- ##
#- Storage -#
## ------- ##

# TODO: work with storage

end # module


text:
module BenchGroupBy

using BenchmarkTools
using SplitApplyCombine: groupreduce
using Transducers

suite = BenchmarkGroup()

n = 10_000
table = map(rand(0:9, n), rand(n)) do key, value
    (key = key, value = value)
end

suite["sum"] = BenchmarkGroup()
suite["sum"]["sac"] = @benchmarkable(groupreduce(x -> x.key, x -> x.value, +, $table))
suite["sum"]["xf-without-init"] = @benchmarkable(foldl(
    right,
    GroupBy(x -> x.key, Map(x -> last(x).value), +),  # no init
    $table,
))
suite["sum"]["xf-with-init"] = @benchmarkable(foldl(
    right,
    GroupBy(x -> x.key, Map(x -> last(x).value), +, 0.0),  # with init
    $table,
))

end  # module
BenchGroupBy.suite


text:
# solve the KS function and store data for plotting

include("kuramoto-sivashinsky.jl")
using kuramoto_sivashinsky

# set the parameters
tp = Float64            # floating point type to use
num_nodes = 127         # number of (interior) nodes to use
num_steps = 5000        # number of time steps to use
Lx = convert(tp, 128)   # domain size
order = 16              # order of accuracy
Time = convert(tp, 500) # time period for simulation

# set initial condition and set-up the KSData structure
dx = Lx/num_nodes
x = Array{tp}(0.0:dx:Lx-dx)
u = exp.(-((x-64).^2)./512)
ks = kuramoto_sivashinsky.buildKSDataPeriodic(order, num_nodes, tp)

# storage for solution over all time steps; for plotting
dt = Time/num_steps
num_steps = round(Int, Time/dt, RoundUp)
sol = zeros(num_nodes+1, num_steps+1)
sol[1:end-1,1] = u
sol[end,1] = u[1]

# solve using midpoint; don't do this in practice!!!
t = convert(tp, 0.0)
Jac = zeros(ks.linear_op)
A = zeros(ks.linear_op)
u_mid = zeros(u)
u_old = zeros(u)
r = zeros(u)
newt_tol = 1e-10 # We will need to ensure this is sufficiently small
max_newt_iter = 50
for n = 1:num_steps
    println("Iteration ",n,": time ",t)

    # Newton iterations
    u_old = u
    for k = 1:max_newt_iter
        u_mid = 0.5*(u + u_old)
        kuramoto_sivashinsky.getKSFuncPeriodic!(ks, t, u_mid, r)
        r = u - u_old - dt*r
        println("\tNewton iter ",k,": res norm = ",norm(r))
        if norm(r) < newt_tol break end
        kuramoto_sivashinsky.getKSJacPeriodic!(ks, t, u_mid, Jac)
        Jac *= 0.5*dt
        Jac -= speye(num_nodes)
        du = Jac\r
        u += du
    end

    t += dt
    sol[1:end-1,n+1] = u
    sol[end,n+1] = u[1]
end

# functionals; here I use only trapezoidal for the temporal integration
u_avg = zero(tp)
u2_avg = zero(tp)
for n = 1:num_steps+1
    wt = dx*dt
    if n == 1 || n == num_steps+1 wt *= 0.5 end
    for val in sol[1:end-1,n]
        u_avg += val*wt
        u2_avg += val*val*wt
    end
end
u_avg /= (Time*Lx)
u2_avg /= (Time*Lx)
println("u averaged = ", u_avg)
println("u^2 averaged = ", u2_avg)

# output to file for plotting
file = open("ks_solution.dat", "w")
for n = 1:num_steps+1
    for val in sol[:,n]
        print(file, val, " ")
    end
    println(file)
end
close(file)

max_stars_repo_name: jehicken/high-order-and-chaos

max_issues_repo_name: jehicken/high-order-and-chaos

max_forks_repo_name: jehicken/high-order-and-chaos


text:
n = NamedArray(rand(2,4))
setnames!(n, ["one", "two"], 1)
setnames!(n, ["a", "b", "c", "d"], 2)

v = NamedArray(rand(1:100, 6), (["a", "b", "c", "d", "e", "f"],), (:index,))

m = NamedArray(rand(2,3,4,3,2,3))


text:
struct MatSpread
    mat::Vector{Vector{Float64}}
    spread::Float64
end

struct MorrisSensitivity
    means
    variances
    elementary_effects
end

function generate_design_matrix(p_range,p_steps;k = 10)
    ps = [range(p_range[i][1], stop=p_range[i][2], length=p_steps[i]) for i in 1:length(p_range)]
    indices = [rand(1:i) for i in p_steps]
    all_idxs = Vector{typeof(indices)}(undef,k)

    for i in 1:k
        flag = 0
        j = rand(1:length(p_range))
        indices[j] += (rand() < 0.5 ? -1 : 1)
        if indices[j] > p_steps[j]
            indices[j] -= 2
        elseif indices[j] < 1.0
            indices[j] += 2
        end
        all_idxs[i] = copy(indices)
    end

    B = Array{Array{Float64}}(undef,k)
    for j in 1:k
        cur_p = [ps[u][(all_idxs[j][u])] for u in 1:length(p_range)]
        B[j] = cur_p
    end
    B
end

function calculate_spread(matrix)
    spread = 0.0
    for i in 2:length(matrix)
        spread += sqrt(sum(abs2.(matrix[i] - matrix[i-1])))
    end
    spread
end

function sample_matrices(p_range,p_steps;len_trajectory=10,num_trajectory=10,total_num_trajectory=5*num_trajectory)
    matrix_array = []
    if total_num_trajectory<num_trajectory
        error("total_num_trajectory should be greater than num_trajectory preferably atleast 3-4 times higher")
    end
    for i in 1:total_num_trajectory
        mat = generate_design_matrix(p_range,p_steps;k=len_trajectory)
        spread = calculate_spread(mat)
        push!(matrix_array,MatSpread(mat,spread))
    end
    sort!(matrix_array,by = x -> x.spread,rev=true)
    matrices = [i.mat for i in matrix_array[1:num_trajectory]]
    matrices
end

function morris_sensitivity(prob::DiffEqBase.DEProblem,alg,t,p_range,p_steps;kwargs...)
    f = function (p)
      prob1 = remake(prob;p=p)
      Array(solve(prob1,alg;saveat=t))
    end
    morris_sensitivity(f,p_range,p_steps;kwargs...)
end

function morris_sensitivity(f,p_range,p_steps;relative_scale=false,kwargs...)
    design_matrices = sample_matrices(p_range,p_steps;kwargs...)
    y1 = f(design_matrices[1][1])
    effects = [[]]
    for i in design_matrices
        y1 = f(i[1])
        for j in 1:length(i)-1
            y2 = y1
            del = i[j+1] - i[j]
            change_index = 0
            for k in 1:length(del)
                if abs(del[k]) > 0
                    change_index = k
                    break
                end
            end
            del = sum(del)
            y1 = f(i[j+1])
            if relative_scale == false
                elem_effect = @. abs((y1-y2)/(del))
            else
                if del > 0
                    elem_effect = @. abs((y1-y2)/(y2*del))
                else
                    elem_effect = @. abs((y1-y2)/(y1*del))
                end
            end
            if length(effects) >= change_index && change_index > 0
                push!(effects[change_index],elem_effect)
            elseif change_index > 0
                push!(effects,[elem_effect])
            end
        end
    end
    means = []
    variances = []
    for k in effects
        push!(means,mean(convert(Array, VectorOfArray(k)),dims=2))
        push!(variances,var(convert(Array, VectorOfArray(k)),dims=2))
    end
    MorrisSensitivity(means,variances,effects)
end


text:
function EventShifter(arg0::EventDetector, arg1::jboolean, arg2::jdouble, arg3::jdouble)
    return EventShifter((EventDetector, jboolean, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function g(obj::EventShifter, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_decreasing_time_shift(obj::EventShifter)
    return jcall(obj, "getDecreasingTimeShift", jdouble, ())
end

function get_increasing_time_shift(obj::EventShifter)
    return jcall(obj, "getIncreasingTimeShift", jdouble, ())
end

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


text:
module fock

export FockBasis, number, destroy, create, displace, fockstate, coherentstate, coherentstate!

import Base: ==

using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
using SparseArrays

"""
    FockBasis(N)

Basis for a Fock space where `N` specifies a cutoff, i.e. what the highest
included fock state is. Note that the dimension of this basis then is N+1.
"""
mutable struct FockBasis <: Basis
    shape::Vector{Int}
    N::Int
    function FockBasis(N::Int)
        if N < 0
            throw(DimensionMismatch())
        end
        new([N+1], N)
    end
end

==(b1::FockBasis, b2::FockBasis) = b1.N==b2.N

"""
    number(b::FockBasis)

Number operator for the specified Fock space.
"""
function number(b::FockBasis)
    diag = complex.(0.:b.N)
    data = spdiagm(0 => diag)
    SparseOperator(b, data)
end

"""
    destroy(b::FockBasis)

Annihilation operator for the specified Fock space.
"""
function destroy(b::FockBasis)
    diag = complex.(sqrt.(1.:b.N))
    data = spdiagm(1 => diag)
    SparseOperator(b, data)
end

"""
    create(b::FockBasis)

Creation operator for the specified Fock space.
"""
function create(b::FockBasis)
    diag = complex.(sqrt.(1.:b.N))
    data = spdiagm(-1 => diag)
    SparseOperator(b, data)
end

"""
    displace(b::FockBasis, alpha)

Displacement operator ``D(α)`` for the specified Fock space.
"""
displace(b::FockBasis, alpha::Number) = exp(dense(alpha*create(b) - conj(alpha)*destroy(b)))

"""
    fockstate(b::FockBasis, n)

Fock state ``|n⟩`` for the specified Fock space.
"""
function fockstate(b::FockBasis, n::Int)
    @assert n <= b.N
    basisstate(b, n+1)
end

"""
    coherentstate(b::FockBasis, alpha)

Coherent state ``|α⟩`` for the specified Fock space.
"""
function coherentstate(b::FockBasis, alpha::Number)
    result = Ket(b, Vector{ComplexF64}(undef, length(b)))
    coherentstate!(result, b, alpha)
    return result
end

"""
    coherentstate!(ket::Ket, b::FockBasis, alpha)

Inplace creation of coherent state ``|α⟩`` for the specified Fock space.
"""
function coherentstate!(ket::Ket, b::FockBasis, alpha::Number)
    alpha = complex(alpha)
    data = ket.data
    data[1] = exp(-abs2(alpha)/2)
    @inbounds for n=1:b.N
        data[n+1] = data[n]*alpha/sqrt(n)
    end
end

end # module


text:
using DataFrames, Missings;

"""
   $(SIGNATURES)

Retrieve molecular identifier from other databases, `targetDb`, for single or mulitple query IDs, `queryId`,
and moreover information on Ensembl gene, transcript and peptide IDs, such as ID and genomic loation.

### Arguments
- `conn::Neo4j.Connection` : a valid connection to a Neo4j graph DB instance.
- `cypher::String` : Cypher `MATCH` query returning tabular data.
- `params::Pair` : parameters which are passed on to the cypher query.
- `elTypes::Vector{Type}` : column types can be provided manually as a Vector{Type}
- `nRowsElTypeCheck::Int` : Number of rows which are used to determine column datatypes (defaults to 1000)

### Examples
```julia-repl
julia> cypherQuery(
         Neo4j.Connection("localhost"),
         "MATCH (p :Person {name: {name}}) RETURN p.name AS Name, p.age AS Age;",
         "name" => "John Doe")

```
"""
function cypherQuery(
      conn::Connection,
      cypher::AbstractString,
      params::Pair...;
      elTypes::Vector{DataType} = Vector{DataType}(),
      nRowsElTypeCheck::Int = 1000)::DataFrames.DataFrame

   url = connurl(conn, "transaction/commit")
   headers = connheaders(conn)
   body = Dict("statements" => [Statement(cypher, Dict(params))])

   resp = HTTP.post(url; headers=headers, body=JSON.json(body))

   if resp.status != 200
     error("Failed to commit transaction ($(resp.status)): $(txn)\n$(resp)")
   end
   respdata = JSON.parse(String(resp.body))

   if !isempty(respdata["errors"])
      error(join(map(i -> (i * ": " * respdata["errors"][1][i]), keys(respdata["errors"][1])), "\n"));
   end
   # parse results into data sink
   # Result(respdata["results"], respdata["errors"])
   if !isempty(respdata["results"][1]["data"])
      return parseResults(respdata["results"][1], elTypes = elTypes, nRowsElTypeCheck = nRowsElTypeCheck);
   else
      return DataFrames.DataFrame();
   end
end

# Currently only supports DataFrames.DataFrame objects
#  -> Future: Allow different data sink types, such as tables from JuliaDB
function parseResults(res::Dict{String, Any}; elTypes::Vector{DataType} = Vector(), nRowsElTypeCheck::Int = 100)::DataFrames.DataFrame
   # Get elementary types from a column where there is no NA value (nothing)
   if isempty(elTypes)
      elTypes = getElTypes(res["data"], nRowsElTypeCheck);
   end
   colNames = Symbol.(collect(res["columns"])) # collect(Symbol, res["columns"]);
   nRows = length(res["data"]);

   x = DataFrames.DataFrame(elTypes, colNames, nRows);

   for (rowIdx, rowVal) in enumerate(res["data"])
      for (colIdx,colVal) in enumerate(rowVal["row"])
         if colVal != nothing
            x[rowIdx,colIdx] = colVal;
         end
      end
   end

   return x;
end

function getElTypes(x::Vector{Any}, nRowsElTypeCheck::Int = 0)::Vector{Type}
   nRecords = length(x);
   elTypes::Vector{Type} = Type[Union{Nothing, Missings.Missing} for i in 1:length(x[1]["row"])];
   nMaxRows = nRecords;
   # elTypes = Type[Union{Nothing, Missings.Missing} for i in 1:length(x[1]["row"])];
   nMaxRows = (nRowsElTypeCheck != 0 && nRowsElTypeCheck <= nMaxRows) ? nRowsElTypeCheck : nRecords;
   checkIdx = trues(length(x[1]["row"]));
   for i in 1:nMaxRows
      # check each column individually
      for el in findall(checkIdx)
         if !(x[i]["row"][el] == nothing)
            if !(typeof(x[i]["row"][el]) === Array{Any,1})
               elTypes[el] = i > 1 ?
                     Union{typeof(x[i]["row"][el]), Missings.Missing} :
                     typeof(x[i]["row"][el]);
            else
               elTypes[el] = i > 1 ?
                     Union{Vector{typeof(x[i]["row"][el][1])}, Missings.Missing} :
                     Vector{typeof(x[i]["row"][el][1])};
            end
            checkIdx[el] = false;
         end
      end
      if isempty(findall(checkIdx))
         break;
      end
   end

   return elTypes;
end


text:
using Dash, DashBootstrapComponents

app = dash(external_stylesheets = [dbc_themes.GRID]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# Moreau envelope of box indicator

f = IndBox(-1, 1)
g = MoreauEnvelope(f, 1e-2)

predicates_test(g)

@test ProximalOperators.is_smooth(g) == true
@test ProximalOperators.is_quadratic(g) == false
@test ProximalOperators.is_set(g) == false

x = [1.0, 2.0, 3.0, 4.0, 5.0]

grad_g_x, g_x = gradient(g, x)

# Moreau envelope of L2 norm = (circulant) Huber loss (for some parameters)

rho = 1e-1
mu = 1e0
f = NormL2(mu)
g = MoreauEnvelope(f, rho)
h = HuberLoss(rho, mu/rho)

predicates_test(g)

@test ProximalOperators.is_smooth(g) == true
@test ProximalOperators.is_quadratic(g) == false
@test ProximalOperators.is_set(g) == false

x = [1.0, 2.0, 3.0, 4.0, 5.0]

@test abs(g(x) - h(x)) <= 1e-12

grad_g_x, g_x = gradient(g, x)
grad_h_x, h_x = gradient(h, x)

@test abs(g_x - g(x)) <= 1e-12
@test abs(h_x - h(x)) <= 1e-12
@test norm(grad_g_x - grad_h_x, Inf) <= 1e-12


text:
function trinary_to_decimal(str::AbstractString)
    typeof(match(r"^[0-2]+$", str)) == Void && return 0
    mapreduce(i->(i[2]-'0')*3^i[1], +, zip(0:length(str), reverse(str)))
end


text: incy = incx


text:
const NOPIVOT = VERSION >= v"1.7" ? NoPivot : Val{false}

function Base.getproperty(LU::LU{<:Any,<:HMatrix},s::Symbol)
    H = getfield(LU,:factors) # the underlying hierarchical matrix
    if s == :L
        return UnitLowerTriangular(H)
    elseif s == :U
        return UpperTriangular(H)
    else
        return getfield(LU,s)
    end
end

"""
    lu!(M::HMatrix,comp)

Hierarhical LU facotrization of `M`, using `comp` to generate the compressed
blocks during the multiplication routines.
"""
function lu!(M::HMatrix,compressor;threads=false)
    # perform the lu decomposition of M in place
    @timeit_debug "lu factorization" begin
        _lu!(M,compressor,threads)
    end
    # wrap the result in the LU structure
    return LU(M,LinearAlgebra.BlasInt[],LinearAlgebra.BlasInt(0))
end

"""
    lu!(M::HMatrix;atol=0,rank=typemax(Int),rtol=atol>0 ||
    rank<typemax(Int) ? 0 : sqrt(eps(Float64)))

Hierarhical LU facotrization of `M`, using the `PartialACA(;atol,rtol;rank)` compressor.
"""
function lu!(M::HMatrix;atol=0,rank=typemax(Int),rtol=atol>0 || rank<typemax(Int) ? 0 : sqrt(eps(Float64)),kwargs...)
    compressor = PartialACA(atol,rank,rtol)
    lu!(M,compressor)
end

"""
    lu(M::HMatrix,args...;kwargs...)

Hierarchical LU factorization. See [`lu!`](@ref) for the available options.
"""
lu(M::HMatrix,args...;kwargs...) = lu!(deepcopy(M),args...;kwargs...)

function _lu!(M::HMatrix,compressor,threads)
    if isleaf(M)
        d = data(M)
        @assert d isa Matrix
        @timeit_debug "dense lu factorization" begin
            lu!(d,NOPIVOT())
        end
    else
        @assert !hasdata(M)
        chdM = children(M)
        m,n = size(chdM)
        for i=1:m
            _lu!(chdM[i,i],compressor,threads)
            for j=i+1:n
                @sync begin
                    @timeit_debug "ldiv! solution" begin
                        if threads
                            Threads.@spawn ldiv!(UnitLowerTriangular(chdM[i,i]),chdM[i,j],compressor)
                        else
                            ldiv!(UnitLowerTriangular(chdM[i,i]),chdM[i,j],compressor)
                        end
                    end
                    @timeit_debug "rdiv! solution" begin
                        if threads
                            Threads.@spawn rdiv!(chdM[j,i],UpperTriangular(chdM[i,i]),compressor)
                        else
                            rdiv!(chdM[j,i],UpperTriangular(chdM[i,i]),compressor)
                        end
                    end
                end
            end
            for j in i+1:m
                for k in i+1:n
                    @timeit_debug "hmul!" begin
                        hmul!(chdM[j,k],chdM[j,i],chdM[i,k],-1,1,compressor)
                    end
                end
            end
        end
    end
    return M
end

function ldiv!(A::LU{<:Any,<:HMatrix},y::AbstractVector;global_index=true)
    p         = A.factors # underlying data
    ctree     = coltree(p)
    rtree     = rowtree(p)
    # permute input
    global_index && permute!(y,loc2glob(ctree))
    L,U = A.L, A.U
    # solve LUx = y through:
    # (a) L(z) = y
    # (b) Ux   = z
    ldiv!(L,y)
    ldiv!(U,y)
    global_index && invpermute!(y,loc2glob(rtree))
    return y
end

function ldiv!(L::HUnitLowerTriangular, y::AbstractVector)
    H = parent(L)
    if isleaf(H)
        d = data(H)
        ldiv!(UnitLowerTriangular(d), y) # B <-- L\B
    else
        @assert !hasdata(H) "only leaves are allowed to have data when using `ldiv`!"
        shift = pivot(H) .- 1
        chdH  = children(H)
        m, n   = size(chdH)
        @assert m === n
        for i = 1:m
            irows  = colrange(chdH[i,i]) .- shift[2]
            bi     = view(y, irows)
            for j = 1:(i - 1)# j<i
                jrows  = colrange(chdH[i,j]) .- shift[2]
                bj     = view(y, jrows)
                _mul131!(bi, chdH[i,j], bj, -1)
            end
            # recursion stage
            ldiv!(UnitLowerTriangular(chdH[i,i]), bi)
        end
    end
    return y
end

function ldiv!(U::HUpperTriangular, y::AbstractVector)
    H = parent(U)
    if isleaf(H)
        d = data(H)
        ldiv!(UpperTriangular(d), y) # B <-- L\B
    else
        @assert !hasdata(H) "only leaves are allowed to have data when using `ldiv`!"
        shift = pivot(H) .- 1
        chdH  = children(H)
        m, n   = size(chdH)
        @assert m === n
        for i = m:-1:1
            irows  = colrange(chdH[i,i]) .- shift[2]
            bi     = view(y, irows)
            for j = i+1:n # j>i
                jrows  = colrange(chdH[i,j]) .- shift[2]
                bj     = view(y, jrows)
                _mul131!(bi, chdH[i,j], bj, -1)
            end
            # recursion stage
            ldiv!(UpperTriangular(chdH[i,i]), bi)
        end
    end
    return y
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
module ComposeAtoms

using Compose, Colors, JuLIP

import Base.display
import Compose: compose, context

export compose, compose_layers, display_layers

plotschemes = Dict(
   :default =>
   Dict( :atcol    => "tomato",
         :qmcol    => "cyan",
         :bufcol   => "purple",
         :bondcol  => "darkblue",
         :bdrycol => "darkgrey",
         :atradius => 0.2,
         :bdryradius => 0.1,
         :atborder => 0.3,
         :bondwidth => 0.3,
         :axisbuffer => 0.05,
         :atcol2 => "aquamarine2"
         ),

   :presentation =>
      Dict( :atcol    => RGB(78/255,115/255,174/255),    # blue
            :qmcol    => RGB(194/255, 79/255, 84/255),   # red
            :bufcol   => "purple",   # RGB(129/255,116/255,176/255),   # purple
            :bondcol  => RGB(0.35, 0.35, 0.35),     # grey
            :bdrycol => RGB(0.7, 0.7, 0.7),     # grey
            :atradius => 0.2,
            :bdryradius => 0.1,
            :atborder => 0.3,
            :bondwidth => 0.3,
            :axisbuffer => 0.05
            ),

   :P2 =>
      Dict( :atcol    => RGB(100/255,130/255,220/255),    # blue
            :qmcol    => RGB(220/255, 120/255, 120/255),   # red
            :bufcol   => RGB(170/255,150/255,210/255),   # RGB(129/255,116/255,176/255),   # purple
            :bondcol  => RGB(0.35, 0.35, 0.35),     # grey
            :bdrycol => RGB(0.7, 0.7, 0.7),     # grey
            :atradius => 0.2,
            :bdryradius => 0.1,
            :atborder => 0.3,
            :bondwidth => 0.3,
            :axisbuffer => 0.05
            ),

   :layers =>
      Dict( :atcol    => "tomato",
            :qmcol    => "cyan",
            :bufcol   => "purple",
            :bondcol  => RGB(0.35, 0.35, 0.35),
            :bdrycol => "darkgrey",
            :atradius => 0.2,
            :atborder => 0.3,
            :bondwidth => 0.3,
            :axisbuffer => 0.05,
            :atcolors => ["tomato", "aquamarine3", "purple"],
            :atradii => [0.4, 0.35, 0.3],
            :bdryradius => 0.1
   ),
)

plotscheme(s) = plotschemes[s]

function analyse_kwargs!(scheme; kwargs...)
   for (key, val) in kwargs
      scheme[key] = val
   end
end

coordinates(X::Matrix) = X[1,:][:], X[2,:][:]
coordinates(X::Vector{JVecF}) = coordinates(mat(X))

"return a reasonable choice for the plotting axis"
function autoaxis(X, scheme)
   x, y = coordinates(X)
    xLim = [extrema(x)...]
    width = xLim[2]-xLim[1]
    yLim = [extrema(y)...]
    height = yLim[2]-yLim[1]
    buffer_factor = scheme[:axisbuffer]
    buffer = min( 1.0, 0.05 * width, 0.05 * height )
    xLim[1] -= buffer; xLim[2] += buffer
    yLim[1] -= buffer; yLim[2] += buffer
    return [xLim[1], xLim[2], yLim[1], yLim[2]]
end

"turn axis into UnitBox for creating contexts"
_ub_(axis) = UnitBox(axis[1], axis[3], axis[2]-axis[1], axis[4]-axis[3])

context(axis::Vector) = context(units=_ub_(axis))

relative_height(ax, width) = (ax[4]-ax[3])/(ax[2]-ax[1]) * width
relative_width(ax, height) = (ax[2]-ax[1])/(ax[4]-ax[3]) * height

function compose_atoms(X, axis, r, atcol, atborder, bordercol)
   # return an empty context if we have an empty set of atoms
   if length(X) == 0
      return (context(axis),)
   end
   x, y = coordinates(X)
   r = ones(length(x)) * r
   return ( context(axis), circle(x, y, r),
              stroke(bordercol), linewidth(atborder), fill(atcol) )
end

compose_bg(axis, bg) = (
   bg ? (compose(context(axis), rectangle()), fill("white")) :
   (context(axis),) )

function unique_bonds(at::AbstractAtoms, rcut; X = nothing)
   bc = pbc(at)
   set_pbc!(at, false)
   if X != nothing
      X0 = positions(at)
      set_positions!(at, X)
   end
   nlist = neighbourlist(at, rcut)
   ij = sort([nlist.i nlist.j], 2)
   ij = unique(ij, 1)
   if X != nothing
      set_positions!(at, X0)
   end
   set_pbc!(at, bc)
   return (ij[:,1], ij[:, 2])
end

# edges2path(X, B) = insert_nans( [X[:, B[1]]; X[:, B[2]]] )
# insert_nans(P) = reshape( [P; NaN * ones(2, size(P, 2))], 2, 3 * size(P, 2) )

function compose_bonds(X, B, scheme, axis)
   i, j = B[1], B[2]
   x, y = coordinates(X)
   P = NTuple{2, Float64}[]
   for n = 1:length(i)
      append!(P, [ (x[i[n]], y[i[n]]), (x[j[n]], y[j[n]]), (NaN, NaN)])
   end
   # P = mat2points(reshape(P, 2, length(P) ÷ 2))
   # return the required tuple to compose
   return ( context(axis), line(P), stroke(scheme[:bondcol]), linewidth(scheme[:bondwidth]) )
end

function compose(at::AbstractAtoms;
                  X=nothing, scheme = plotscheme[:default],
                  axis = autoaxis(X, scheme),
                  rcut = rnn(at) * 1.2, bg = false,
                  Iqm = Int[], Ibdry = Int[], Ibuf = Int[],
                  kwargs... )

   analyse_kwargs!(scheme; kwargs...)
   i, j = unique_bonds(at, rcut; X=X)
   if X == nothing
      X = positions(at)
   end

   Idom = collect(1:length(at))
   Imm = setdiff(Idom, union(Iqm, Ibdry, Ibuf))

   ctx = compose( context(axis),
      compose_atoms( X[Imm], axis, scheme[:atradius], scheme[:atcol], scheme[:atborder], scheme[:bondcol] ),
      compose_atoms( X[Iqm], axis, scheme[:atradius], scheme[:qmcol], scheme[:atborder], scheme[:bondcol] ),
      compose_atoms( X[Ibuf], axis, scheme[:atradius], scheme[:atcol], scheme[:atborder], scheme[:qmcol] ),
      compose_atoms( X[Ibdry], axis, scheme[:bdryradius], scheme[:bdrycol], scheme[:atborder], scheme[:bondcol] ),
      compose_bonds( X, (i, j), scheme, axis )
                  # compose_bg(axis, bg)
               )
   return ctx
end

function compose_layers( at::AbstractAtoms, zlayers;
                        X=nothing, scheme = plotscheme[:layers],
                        axis = autoaxis(X, scheme),
                        rcut = rnn(at) * 1.2, bg = false,
                        Iqm = Int[], Ibdry = Int[], Ibuf = Int[],   # ignore all except Ibrdy for now
                        kwargs... )

   analyse_kwargs!(scheme; kwargs...)
   i, j = unique_bonds(at, rcut; X=X)
   if X == nothing
      X = positions(at)
   end

   # add atoms to one of the two layers
   z = mat(X)[3,:]
   II = [ Int[] for j = 1:length(zlayers) ]
   for j = 1:length(z)
      _, i = findmin(abs(z[j] - zlayers))
      push!(II[i], j)
   end
   for Ij in II
      Ij = setdiff(Ij, Ibdry)
   end

   cat = [ compose_atoms( X[II[j]], axis, scheme[:atradii][j], scheme[:atcolors][j],
                          scheme[:atborder], scheme[:bondcol] )  for j = 1:length(II) ]

   ctx = compose( context(axis), cat...,
      compose_atoms( X[Ibdry], axis, scheme[:bdryradius], scheme[:bdrycol], scheme[:atborder], scheme[:bondcol] ),
      compose_bonds( X, (i, j), scheme, axis )
               )
   return ctx
end

function _display(ctx, filename, display, plotwidth, plotheight, Img)
   if (filename == nothing) || (display == true)
      draw( SVG(plotwidth, plotheight), ctx )
   end
   if filename != nothing
      draw( Img(filename, plotwidth, plotheight), ctx )
   end
end

"""
`display(at::AbstractAtoms; kwargs...)`

## Keyword arguments
* X : positions (default, reference positions)
* axis : default is an autoaxis
* plotwidth: default is 12cm
* Img : image type, default is SVG
* rnn : bond-length, default given by `rnn(at)` (cf JuLIP)
"""
display(at::AbstractAtoms;
               X=positions(at), scheme = plotscheme(:default),
               axis = autoaxis(X, scheme), display = true,
               plotwidth = 12cm, plotheight= relative_height(axis, plotwidth),
               Img=SVG, filename = nothing, kwargs... ) =
   _display( compose(at; X=X, axis=axis, scheme=scheme, kwargs...),
             filename, display, plotwidth, plotheight, Img )

display_layers(at::AbstractAtoms, zlayers;
                  X=positions(at), scheme = plotscheme(:layers),
                  axis = autoaxis(X, scheme), display = true,
                  plotwidth = 12cm, plotheight= relative_height(axis, plotwidth),
                  Img=SVG, filename = nothing, kwargs... ) =
   _display( compose_layers(at, zlayers; X=X, axis=axis, scheme=scheme, kwargs...),
            filename, display, plotwidth, plotheight, Img )

end


text:
using Random
using MultiAgentPathFinding
using MultiAgentAllocationTransit
using BenchmarkTools
using Logging

global_logger(SimpleLogger(stderr, Logging.Warn))

## Hard-code a bunch of depots and many more sites
const DEPOT1 = LatLonCoords((lat = 37.762892, lon = -122.472193))
const DEPOT2 = LatLonCoords((lat = 37.751751, lon = -122.410654))
const DEPOT3 = LatLonCoords((lat = 37.718779, lon = -122.462401))

const SITE1 = LatLonCoords((lat = 37.789290, lon = -122.426797)) # Top right
const SITE2 = LatLonCoords((lat = 37.739611, lon = -122.492203)) # Way left
const SITE3 = LatLonCoords((lat = 37.780593, lon = -122.434555)) # Under SITE 1
const SITE4 = LatLonCoords((lat = 37.739011, lon = -122.430059)) # Near the bottom
const SITE5 = LatLonCoords((lat = 37.749018, lon = -122.462966)) # 2nd most left
const SITE1DUMMY = LatLonCoords((lat = 37.788042, lon = -122.429579))

const TRANSIT_CAP_RANGE = (1, 1)
const ECBS_WEIGHT = 1.05

stop_coords_file = "data/sfmta/stop_to_coords.json"
trips_file = "data/sfmta/trips.json"
drone_params_file = "data/drone_params.toml"

depots = [DEPOT1, DEPOT2, DEPOT3]
sites = [SITE1, SITE2, SITE3, SITE4, SITE5, SITE1DUMMY]

agent_tasks = [AgentTask((origin = 1, site = 1, dest = 2)),
               AgentTask((origin = 2, site = 6, dest = 1)),
               AgentTask((origin = 1, site = 3, dest = 2)),
               AgentTask((origin = 2, site = 4, dest = 3)),
               AgentTask((origin = 3, site = 5, dest = 1))
               ]

rng = MersenneTwister(1234)

# Load drone params and setup ENV!!
drone_params = parse_drone_params(drone_params_file)

# Generate OTG and load TG and drone params
otg = OffTransitGraph(depots = depots, sites = sites)
tg = load_transit_graph_latlong(stop_coords_file, trips_file, TRANSIT_CAP_RANGE, rng)

tg, stop_idx_to_trips = transit_graph_preprocessing(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean, drone_params)

aug_trips_fws_dists = augmented_trip_meta_graph_fws_dists(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                                                      length(depots), length(sites),
                                                      vcat(depots, sites),
                                                      drone_params)

state_graph, depot_sites_to_vtx, trip_to_vtx_range = setup_state_graph(tg, otg)

agent_states = [AgentState(task=agt_task) for agt_task in agent_tasks]

env = MAPFTransitEnv(off_transit_graph = otg, transit_graph = tg, state_graph = state_graph,
                     agent_states = agent_states, depot_sites_to_vtx = depot_sites_to_vtx, trip_to_vtx_range = trip_to_vtx_range,
                     stop_idx_to_trips = stop_idx_to_trips, aug_trips_fws_dists = aug_trips_fws_dists,
                     plan_ref_times = zeros(length(agent_tasks)),
                     drone_params = drone_params, dist_fn = MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                     curr_site_points = [0 for _ = 1:length(agent_tasks)])

# Define ECBS Solver; Run search
initial_states = Vector{MAPFTransitVertexState}(undef, length(agent_tasks))
for i = 1:length(agent_tasks)
    initial_states[i] = env.state_graph.vertices[env.depot_sites_to_vtx[string("d-", agent_tasks[i].origin)]]
end

solver = ECBSSolver{MAPFTransitVertexState,MAPFTransitAction,Float64,Makespan,MAPFTransitConflict,MAPFTransitConstraints,MAPFTransitEnv}(env = env, weight = ECBS_WEIGHT)

solution = search!(solver, initial_states)
#
# # NEED TO RE_INIT SOLVER!!!!
agent_states = [AgentState(task=agt_task) for agt_task in agent_tasks]

env = MAPFTransitEnv(off_transit_graph = otg, transit_graph = tg, state_graph = state_graph,
                     agent_states = agent_states, depot_sites_to_vtx = depot_sites_to_vtx, trip_to_vtx_range = trip_to_vtx_range,
                     stop_idx_to_trips = stop_idx_to_trips, aug_trips_fws_dists = aug_trips_fws_dists,
                     plan_ref_times = zeros(length(agent_tasks)), valid_transit_options = zeros(length(agent_tasks)),
                     drone_params = drone_params, dist_fn = MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                     curr_site_points = [0 for _ = 1:length(agent_tasks)])

solver = ECBSSolver{MAPFTransitVertexState,MAPFTransitAction,Float64,Makespan,MAPFTransitConflict,MAPFTransitConstraints,MAPFTransitEnv}(env = env, weight = ECBS_WEIGHT)
@time solution = search!(solver, initial_states)


text:
using MonteCarlo, Test
import MPI

#=

t1 = MPI.Initialized()
m = HubbardModelAttractive(2, 2)
scheduler = SimpleScheduler(LocalSweep(), MPIReplicaPull())
mc = DQMC(m, beta = 1.0, scheduler = scheduler)
t2 = MPI.Initialized()

io = open("worker$(MPI.Comm_rank(MPI.COMM_WORLD)).log", "w")
old_stdout = stdout
old_stderr = stderr
redirect_stdout(io)
redirect_stderr(io)

append!(MonteCarlo.connected_ids, filter(i -> i != MPI.Comm_rank(MPI.COMM_WORLD), 0:3))
run!(mc)

=#

MPI.Initialized() || MPI.Init()

io = open("worker$(MPI.Comm_rank(MPI.COMM_WORLD)).log", "w")
old_stdout = stdout
old_stderr = stderr
redirect_stdout(io)
redirect_stderr(io)

# maps 0 1 2 3 ...
#       X   X
# to   0 1 2 3 ...
partner1(rank = MPI.Comm_rank(MPI.COMM_WORLD)) = rank - isodd(rank) + iseven(rank)
# maps N 0 1 2 ...
#       X   X
# to   N 0 1 2 ...
function partner2(rank = MPI.Comm_rank(MPI.COMM_WORLD), max = MPI.Comm_size(MPI.COMM_WORLD) - 1)
    mod1(rank + isodd(rank) - iseven(rank) + 1, max + 1) - 1
end

m = HubbardModelAttractive(2, 2)
scheduler = SimpleScheduler(
    LocalSweep(), MPIReplicaExchange(partner1()), MPIReplicaExchange(partner2())
)
println(MPI.Comm_rank(MPI.COMM_WORLD), " < - > ", partner1(), ", ", partner2())
show(scheduler)
mc = DQMC(m, beta = 1.0, scheduler = scheduler)

run!(mc)

println("Simulations finished, good.")
# @testset "MPI Updates" begin
#     # Check if init is working as expected
#     @test t1
#     @test t2
# end

redirect_stdout(old_stdout)
redirect_stderr(old_stderr)

MPI.Barrier(MPI.COMM_WORLD)

exit()


text:
import Term.style: apply_style

@testset "\e[31mMACROS" begin
    @test (@green "string") == apply_style("[green]string[/green]")
    @test (@blue "string") == apply_style("[blue]string[/blue]")
    @test (@red "string") == apply_style("[red]string[/red]")

    @test (@bold "string") == apply_style("[bold]string[/bold]")
    @test (@italic "string") == apply_style("[italic]string[/italic]")
    @test (@underline "string") == apply_style("[underline]string[/underline]")

    @test (@style "string" bold underline) ==
          apply_style("[bold underline]string[/bold underline]")
    @test (@style "string" bold greeen) == apply_style("[bold greeen]string[/bold greeen]")
    @test (@style "string" red on_blue underline) ==
          apply_style("[red on_blue underline]string[/red on_blue underline]")
    @test (@style "string" bold italic) == apply_style("[bold italic]string[/bold italic]")
end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

using Test
using MGVI

using Distributions
using LinearAlgebra
using Random
using ForwardDiff

if :ModelPolyfit ∉ names(Main)
    include("test_models/model_polyfit.jl")
    import .ModelPolyfit
end

Test.@testset "test_jacobians_consistency" begin
    let
        A = rand(4,6)
        f = let A=A; x -> (A*(x.^2).^2) end
        x = rand(6); l = rand(4); r = rand(6)
        J_ref = ForwardDiff.jacobian(f, x)

        J1 = @inferred FullJacobianFunc(f)(x)
        J2 = @inferred FwdRevADJacobianFunc(f)(x)
        J3 = @inferred FwdDerJacobianFunc(f)(x)

        for J in (J1, J2, J3)
            @test @inferred(Array(J)) ≈ J_ref
            @test @inferred(J' * l) ≈ J_ref' * l
            @test @inferred(J * r) ≈ J_ref * r
        end
    end

    epsilon = 1E-5
    Random.seed!(145)

    _flat_model = MGVI.flat_params ∘ ModelPolyfit.model
    true_params = ModelPolyfit.true_params

    full_jac = FullJacobianFunc(_flat_model)(true_params)
    fwdder_jac = FwdDerJacobianFunc(_flat_model)(true_params)
    fwdrevad_jac = FwdRevADJacobianFunc(_flat_model)(true_params)

    for i in 1:min(size(full_jac)...)

        vec = rand(size(full_jac, 2))

        Test.@test norm(fwdder_jac*vec - full_jac*vec) < epsilon

        Test.@test norm(fwdrevad_jac*vec - full_jac*vec) < epsilon

    end

end


text:
"""
    rich_club(g,k)

Return the non-normalised [rich-club coefficient](https://en.wikipedia.org/wiki/Rich-club_coefficient) of graph `g`,
with degree cut-off `k`.

```jldoctest
julia> using LightGraphs
julia> g = star_graph(5)
julia> rich_club(g,1)
0.4
```
"""
function rich_club(g::AbstractGraph{T},k::Int) where T
    E = zero(T)
    for e in edges(g)
        if (degree(g,src(e)) >= k) && (degree(g,dst(e)) >= k )
            E +=1
        end
    end
    N = count(degree(g) .>= k)
    return 2*E / (N*(N-1))
end


text:
function testconnectnodes()
  n1 = Node(1,2,3)
  n2 = Node(4,5,6)
  n3 = Node(7,8,9)
  sp = SegmentParameters(w=10,h=2,sigma=3)
  segmentarray = connectnodes([n1,n2,n3],sp)
  @test length(segmentarray) == 2
  g6 = Group([n1,n2,n3,segmentarray...])
  verified =
  """
  N_1 x=1.000000000e+00 y=2.000000000e+00 z=3.000000000e+00
  N_2 x=4.000000000e+00 y=5.000000000e+00 z=6.000000000e+00
  N_3 x=7.000000000e+00 y=8.000000000e+00 z=9.000000000e+00
  E_4 N_1 N_2
  +  w=1.000000000e+01 h=2.000000000e+00
  +  sigma=3.000000000e+00
  E_5 N_2 N_3
  +  w=1.000000000e+01 h=2.000000000e+00
  +  sigma=3.000000000e+00
  """
  testelement(g6,verified)
end
testconnectnodes()

function testviagroup()
  @test_throws(ArgumentError,viagroup(height=20, h=3))
  @test_throws(ArgumentError,viagroup(radius=10, h=3))
  @test_throws(ArgumentError,viagroup(radius=10, height=20))
  @test_throws(ArgumentError,viagroup(radius=10, height=20, h=3,n=1))
  via = viagroup(radius=10, height=20, h=3)
  verified =
  """
  N_1 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00
  N_2 x=7.071067812e+00 y=7.071067812e+00 z=0.000000000e+00
  N_3 x=6.123233996e-16 y=1.000000000e+01 z=0.000000000e+00
  N_4 x=-7.071067812e+00 y=7.071067812e+00 z=0.000000000e+00
  N_5 x=-1.000000000e+01 y=1.224646799e-15 z=0.000000000e+00
  N_6 x=-7.071067812e+00 y=-7.071067812e+00 z=0.000000000e+00
  N_7 x=-1.836970199e-15 y=-1.000000000e+01 z=0.000000000e+00
  N_8 x=7.071067812e+00 y=-7.071067812e+00 z=0.000000000e+00
  N_9 x=1.000000000e+01 y=0.000000000e+00 z=-2.000000000e+01
  N_10 x=7.071067812e+00 y=7.071067812e+00 z=-2.000000000e+01
  N_11 x=6.123233996e-16 y=1.000000000e+01 z=-2.000000000e+01
  N_12 x=-7.071067812e+00 y=7.071067812e+00 z=-2.000000000e+01
  N_13 x=-1.000000000e+01 y=1.224646799e-15 z=-2.000000000e+01
  N_14 x=-7.071067812e+00 y=-7.071067812e+00 z=-2.000000000e+01
  N_15 x=-1.836970199e-15 y=-1.000000000e+01 z=-2.000000000e+01
  N_16 x=7.071067812e+00 y=-7.071067812e+00 z=-2.000000000e+01
  N_17 x=0.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  N_18 x=0.000000000e+00 y=0.000000000e+00 z=-2.000000000e+01
  .equiv N_17 N_1 N_2 N_3 N_4 N_5 N_6 N_7 N_8
  .equiv N_18 N_9 N_10 N_11 N_12 N_13 N_14 N_15 N_16
  E_19 N_1 N_9
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=6.123233996e-17 wy=1.000000000e+00 wz=0.000000000e+00
  E_20 N_2 N_10
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=-7.071067812e-01 wy=7.071067812e-01 wz=0.000000000e+00
  E_21 N_3 N_11
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=-1.000000000e+00 wy=1.224646799e-16 wz=0.000000000e+00
  E_22 N_4 N_12
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=-7.071067812e-01 wy=-7.071067812e-01 wz=0.000000000e+00
  E_23 N_5 N_13
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=-1.836970199e-16 wy=-1.000000000e+00 wz=0.000000000e+00
  E_24 N_6 N_14
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=7.071067812e-01 wy=-7.071067812e-01 wz=0.000000000e+00
  E_25 N_7 N_15
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=1.000000000e+00 wy=-2.449293598e-16 wz=0.000000000e+00
  E_26 N_8 N_16
  +  w=8.284271247e+00 h=3.000000000e+00 nhinc=5 nwinc=1
  +  wx=7.071067812e-01 wy=7.071067812e-01 wz=0.000000000e+00
  """
  testelement(via, verified)
  @test length(via[:alltop]) == 8
  @test length(via[:allbot]) == 8
  @test typeof(via[:top]) == Node
  @test typeof(via[:top]) == Node
  xyz(via[:top]) == [0.0,0.0,0.0]
end
testviagroup()

function testhelixnodes()
  radius = 10.0
  pitch = 2.0
  radians = 4π
  radiansperpoint = π/4
  hp1 = helixnodes(radius, pitch, radians, radiansperpoint)
  verified =
  """
  N_1 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00
  N_2 x=7.071067812e+00 y=7.071067812e+00 z=2.500000000e-01
  N_3 x=6.123233996e-16 y=1.000000000e+01 z=5.000000000e-01
  N_4 x=-7.071067812e+00 y=7.071067812e+00 z=7.500000000e-01
  N_5 x=-1.000000000e+01 y=1.224646799e-15 z=1.000000000e+00
  N_6 x=-7.071067812e+00 y=-7.071067812e+00 z=1.250000000e+00
  N_7 x=-1.836970199e-15 y=-1.000000000e+01 z=1.500000000e+00
  N_8 x=7.071067812e+00 y=-7.071067812e+00 z=1.750000000e+00
  N_9 x=1.000000000e+01 y=-2.449293598e-15 z=2.000000000e+00
  N_10 x=7.071067812e+00 y=7.071067812e+00 z=2.250000000e+00
  N_11 x=3.061616998e-15 y=1.000000000e+01 z=2.500000000e+00
  N_12 x=-7.071067812e+00 y=7.071067812e+00 z=2.750000000e+00
  N_13 x=-1.000000000e+01 y=3.673940397e-15 z=3.000000000e+00
  N_14 x=-7.071067812e+00 y=-7.071067812e+00 z=3.250000000e+00
  N_15 x=-4.286263797e-15 y=-1.000000000e+01 z=3.500000000e+00
  N_16 x=7.071067812e+00 y=-7.071067812e+00 z=3.750000000e+00
  N_17 x=1.000000000e+01 y=-4.898587197e-15 z=4.000000000e+00
  """
  testelement(Group(hp1), verified)
  hp2 = helixnodes(radius, pitch, radians+0.1, radiansperpoint)
  verified =
  """
  N_1 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00
  N_2 x=7.071067812e+00 y=7.071067812e+00 z=2.500000000e-01
  N_3 x=6.123233996e-16 y=1.000000000e+01 z=5.000000000e-01
  N_4 x=-7.071067812e+00 y=7.071067812e+00 z=7.500000000e-01
  N_5 x=-1.000000000e+01 y=1.224646799e-15 z=1.000000000e+00
  N_6 x=-7.071067812e+00 y=-7.071067812e+00 z=1.250000000e+00
  N_7 x=-1.836970199e-15 y=-1.000000000e+01 z=1.500000000e+00
  N_8 x=7.071067812e+00 y=-7.071067812e+00 z=1.750000000e+00
  N_9 x=1.000000000e+01 y=-2.449293598e-15 z=2.000000000e+00
  N_10 x=7.071067812e+00 y=7.071067812e+00 z=2.250000000e+00
  N_11 x=3.061616998e-15 y=1.000000000e+01 z=2.500000000e+00
  N_12 x=-7.071067812e+00 y=7.071067812e+00 z=2.750000000e+00
  N_13 x=-1.000000000e+01 y=3.673940397e-15 z=3.000000000e+00
  N_14 x=-7.071067812e+00 y=-7.071067812e+00 z=3.250000000e+00
  N_15 x=-4.286263797e-15 y=-1.000000000e+01 z=3.500000000e+00
  N_16 x=7.071067812e+00 y=-7.071067812e+00 z=3.750000000e+00
  N_17 x=1.000000000e+01 y=-4.898587197e-15 z=4.000000000e+00
  N_18 x=9.950041653e+00 y=9.983341665e-01 z=4.031830989e+00
  """
  testelement(Group(hp2), verified)
end
testhelixnodes()

function testcoilcraft1010vsgroup()
  inductor = coilcraft1010vsgroup("1010VS-111ME")
  @test length(elements(inductor)) == 85
  @test length(keys(terms(inductor))) == 2
  verified =
  """
  .units mm
  N_1 x=3.550000000e+00 y=3.550000000e+00 z=3.250000000e-01
  N_2 x=3.550000000e+00 y=-3.550000000e+00 z=3.250000000e-01
  N_3 x=3.550000000e+00 y=-3.550000000e+00 z=2.025000000e+00
  N_4 x=3.550000000e+00 y=0.000000000e+00 z=2.025000000e+00
  N_5 x=2.510229073e+00 y=2.510229073e+00 z=2.106250000e+00
  N_6 x=2.173748068e-16 y=3.550000000e+00 z=2.187500000e+00
  N_7 x=-2.510229073e+00 y=2.510229073e+00 z=2.268750000e+00
  N_8 x=-3.550000000e+00 y=4.347496137e-16 z=2.350000000e+00
  N_9 x=-2.510229073e+00 y=-2.510229073e+00 z=2.431250000e+00
  N_10 x=-6.521244205e-16 y=-3.550000000e+00 z=2.512500000e+00
  N_11 x=2.510229073e+00 y=-2.510229073e+00 z=2.593750000e+00
  N_12 x=3.550000000e+00 y=-8.694992274e-16 z=2.675000000e+00
  N_13 x=2.510229073e+00 y=2.510229073e+00 z=2.756250000e+00
  N_14 x=1.086874034e-15 y=3.550000000e+00 z=2.837500000e+00
  N_15 x=-2.510229073e+00 y=2.510229073e+00 z=2.918750000e+00
  N_16 x=-3.550000000e+00 y=1.304248841e-15 z=3.000000000e+00
  N_17 x=-2.510229073e+00 y=-2.510229073e+00 z=3.081250000e+00
  N_18 x=-1.521623648e-15 y=-3.550000000e+00 z=3.162500000e+00
  N_19 x=2.510229073e+00 y=-2.510229073e+00 z=3.243750000e+00
  N_20 x=3.550000000e+00 y=-1.738998455e-15 z=3.325000000e+00
  N_21 x=2.510229073e+00 y=2.510229073e+00 z=3.406250000e+00
  N_22 x=1.956373262e-15 y=3.550000000e+00 z=3.487500000e+00
  N_23 x=-2.510229073e+00 y=2.510229073e+00 z=3.568750000e+00
  N_24 x=-3.550000000e+00 y=2.173748068e-15 z=3.650000000e+00
  N_25 x=-2.510229073e+00 y=-2.510229073e+00 z=3.731250000e+00
  N_26 x=-8.697189655e-15 y=-3.550000000e+00 z=3.812500000e+00
  N_27 x=2.510229073e+00 y=-2.510229073e+00 z=3.893750000e+00
  N_28 x=3.550000000e+00 y=-2.608497682e-15 z=3.975000000e+00
  N_29 x=2.510229073e+00 y=2.510229073e+00 z=4.056250000e+00
  N_30 x=-3.480194291e-15 y=3.550000000e+00 z=4.137500000e+00
  N_31 x=-2.510229073e+00 y=2.510229073e+00 z=4.218750000e+00
  N_32 x=-3.550000000e+00 y=3.043247296e-15 z=4.300000000e+00
  N_33 x=-2.510229073e+00 y=-2.510229073e+00 z=4.381250000e+00
  N_34 x=-9.566688883e-15 y=-3.550000000e+00 z=4.462500000e+00
  N_35 x=2.510229073e+00 y=-2.510229073e+00 z=4.543750000e+00
  N_36 x=3.550000000e+00 y=-3.477996910e-15 z=4.625000000e+00
  N_37 x=2.510229073e+00 y=2.510229073e+00 z=4.706250000e+00
  N_38 x=-2.610695063e-15 y=3.550000000e+00 z=4.787500000e+00
  N_39 x=-2.510229073e+00 y=2.510229073e+00 z=4.868750000e+00
  N_40 x=-3.550000000e+00 y=3.912746523e-15 z=4.950000000e+00
  N_41 x=-3.550000000e+00 y=-3.550000000e+00 z=4.950000000e+00
  N_42 x=-3.550000000e+00 y=-3.550000000e+00 z=3.250000000e-01
  N_43 x=-3.550000000e+00 y=3.550000000e+00 z=3.250000000e-01
  E_44 N_1 N_2
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_45 N_2 N_3
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_46 N_3 N_4
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_47 N_4 N_5
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_48 N_5 N_6
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_49 N_6 N_7
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_50 N_7 N_8
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_51 N_8 N_9
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_52 N_9 N_10
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_53 N_10 N_11
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_54 N_11 N_12
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_55 N_12 N_13
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_56 N_13 N_14
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_57 N_14 N_15
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_58 N_15 N_16
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_59 N_16 N_17
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_60 N_17 N_18
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_61 N_18 N_19
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_62 N_19 N_20
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_63 N_20 N_21
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_64 N_21 N_22
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_65 N_22 N_23
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_66 N_23 N_24
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_67 N_24 N_25
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_68 N_25 N_26
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_69 N_26 N_27
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_70 N_27 N_28
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_71 N_28 N_29
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_72 N_29 N_30
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_73 N_30 N_31
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_74 N_31 N_32
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_75 N_32 N_33
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_76 N_33 N_34
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_77 N_34 N_35
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_78 N_35 N_36
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_79 N_36 N_37
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_80 N_37 N_38
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_81 N_38 N_39
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_82 N_39 N_40
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_83 N_40 N_41
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_84 N_41 N_42
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  E_85 N_42 N_43
  +  w=1.700000000e+00 h=6.000000000e-01
  +  sigma=5.850000000e+04
  .units mm
  """
  testelement(inductor, verified)
end
testcoilcraft1010vsgroup()


text:
struct Dataset
    A::Matrix{Int}          # Adjacency matrix of an observed graph
    n::Int                  # Number of nodes in the graph
    m::Int                  # Number of edges in the graph
    n_communities::Int      # Number of communities to be recovered
    k::Vector{Int}          # Degrees of each node

    # Constructor function for struct Dataset
    function Dataset(A::Matrix{Int}, q::Int)
        if any(A .< 0)
            throw(DomainError("Adjacency matrix A must have only non-negative integer elements."))
        end
        if q <= 0
            throw(DomainError("Number of communities/clusters q must be a positive integer value."))
        end
        if ~all([A[i,i] for i=1:size(A,1)] .% 2 .== 0)
            throw(ArgumentError("Self-edges should count as 2. Please check the elements in the diagonal of the adjacency matrix."))
        end
        n = size(A,1)
        m = Int(sum(A)/2)
        k = vec(sum(A, dims=1))
        return Dataset(A, n, m, q, k)
    end

    function Dataset(dataset_path::String)
        return loadDataset(dataset_path)
    end

    function Dataset(A::Matrix{Int}, n::Int, m::Int, q::Int, k::Vector{Int})
        if any(A .< 0)
            throw(DomainError("Adjacency matrix A must have only non-negative integer elements."))
        end
        if (size(A,1) != n)
            throw(ArgumentError("The size of the adjacency matrix A must match the specified number of nodes n=$n"))
        end
        if (sum(A)/2 != m)
            throw(ArgumentError("Number of edges in adjacency matrix A does not match the specified value of m=$m."))
        end
        if (A != A')
            throw(ArgumentError("Adjacency matrix must be symmetric."))
        end
        if q <= 0
            throw(DomainError("Number of communities/clusters q must be a positive integer value."))
        end
        # TODO: check if matrix A follows the convention that self edges count as 2
        return new(A, n, m, q, k)
    end
end

function loadDataset(dataset_path::String)::Dataset
    lines = readlines(dataset_path)
    n,m,q = map(x->parse(Int,x), split(lines[1],","))
    A = Matrix(zeros(Int,n,n))
    for line in lines[2:end]
        i,j = map(x->parse(Int,x), split(line,","))
        A[i,j] += 1
    end
    A = A + A'
    k = vec(sum(A, dims=1)) # get node degree of each node
    dataset = Dataset(A, n, m, q, k)
    return dataset
end

function saveDataset(dataset::Dataset, path::String)
    firstline = string(join([dataset.n, dataset.m, dataset.n_communities], ","),"\n")
    open(path, "w") do f
        # Write fist line to file
        write(f, firstline)
        # Write edges as edge list
        for i=1:dataset.n, j=i:dataset.n
            n_edges = dataset.A[i,j]
            if i == j
                n_edges = Int(n_edges / 2)
            end
            for w in 1:n_edges
                write(f, "$i,$j\n")
            end
        end
    end
end

function loadKarate()
    dataset_path = "../instances/zachary.in"
    return Dataset(dataset_path)
end

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM