text:
module TestConstraintNormInfinity

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_NormInfinity()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NormInfinity{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariables_NormInfinityCone",
            "test_basic_VectorAffineFunction_NormInfinityCone",
            "test_basic_VectorQuadraticFunction_NormInfinityCone",
        ],
    )
    return
end

function test_NormInfinity_1()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NormInfinity{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.5, 1.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[0.0, 1.0, 0.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [-1]],
        )
    MOI.Test.test_conic_NormInfinityCone_VectorOfVariables(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_NormInfinityCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    var_names = ["x", "y", "z"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
    zeros = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Zeros,
        }(),
    )
    @test length(zeros) == 2
    MOI.set(mock, MOI.ConstraintName(), zeros[1], "x_eq")
    MOI.set(mock, MOI.ConstraintName(), zeros[2], "y_eq")

    s = """
    variables: x, y, z
    nonneg: [x + -1.0y, x + -1.0z, x + y, x + z] in VecMathOptInterface.Nonnegatives(4)
    x_eq: [-1.0 + x] in VecMathOptInterface.Zeros(1)
    y_eq: [-0.5 + y] in VecMathOptInterface.Zeros(1)
    maxobjective: y + z
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names,
        ["nonneg", "x_eq", "y_eq"],
    )
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    norminf = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.NormInfinityCone,
        }(),
    )
    @test length(norminf) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), norminf[1], "norminf")
    zeros = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Zeros,
        }(),
    )
    @test length(zeros) == 2
    MOI.set(bridged_mock, MOI.ConstraintName(), zeros[1], "x_eq")
    MOI.set(bridged_mock, MOI.ConstraintName(), zeros[2], "y_eq")

    s = """
    variables: x, y, z
    norminf: [1.0x, y, z] in VecMathOptInterface.NormInfinityCone(3)
    x_eq: [-1.0 + x] in VecMathOptInterface.Zeros(1)
    y_eq: [-0.5 + y] in VecMathOptInterface.Zeros(1)
    maxobjective: y + z
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["norminf", "x_eq", "y_eq"],
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.NormInfinityCone,
            }(),
        ),
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [4.0, 1.0, -2.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            nonneg_value = Float64[3, 6, 5, 2]
        else
            nonneg_value = [0.25, 2.25, 1.25, 0.25]
        end
        @test MOI.get(mock, attr, nonneg[1]) ≈ nonneg_value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),),
    )
    return
end

function test_conic_NormInfinityCone_3()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NormInfinity{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [2, -1, -1, -1],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [vcat(fill(inv(3), 3), zeros(3)), fill(inv(3), 3)],
        )
    MOI.Test.test_conic_NormInfinityCone_3(bridged_mock, config)
    var_names = ["x", "y1", "y2", "y3"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    @test length(nonneg) == 2
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg1")
    MOI.set(mock, MOI.ConstraintName(), nonneg[2], "nonneg2")

    s = """
    variables: x, y1, y2, y3
    nonneg1: [x + -1.0y1 + -3.0, x + -1.0y2 + -3.0, x + -1.0y3 + -3.0, x + y1 + 1.0, x + y2 + 1.0, x + y3 + 1.0] in VecMathOptInterface.Nonnegatives(6)
    nonneg2: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in VecMathOptInterface.Nonnegatives(3)
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names,
        ["nonneg1", "nonneg2"],
    )
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    norminf = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.NormInfinityCone,
        }(),
    )
    @test length(norminf) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), norminf[1], "norminf")
    nonneg = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    @test length(nonneg) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), nonneg[1], "nonneg")

    s = """
    variables: x, y1, y2, y3
    norminf: [x + -1.0, y1 + 2.0, y2 + 2.0, y3 + 2.0] in VecMathOptInterface.NormInfinityCone(4)
    nonneg: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in VecMathOptInterface.Nonnegatives(3)
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["norminf", "nonneg"],
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.NormInfinityCone,
            }(),
        ),
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = Float64[5, 1, -2, 3]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            nonneg_value = Float64[4, 7, 2, 6, 3, 8]
        else
            nonneg_value = [-1, 11, -1, 5, -1, 17] / 6
        end
        @test MOI.get(mock, attr, nonneg[1]) ≈ nonneg_value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        4,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 1),),
    )
    return
end

function test_NormOne()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NormOne{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariables_NormOneCone",
            "test_basic_VectorAffineFunction_NormOneCone",
            "test_basic_VectorQuadraticFunction_NormOneCone",
        ],
    )
    return
end

function test_conic_NormOneCone_VectorOfVariables()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NormOne{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.5, 0.5, 0.5, 0.5],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[1.0, 1.0, 1.0, 0.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [0]],
        )
    MOI.Test.test_conic_NormOneCone_VectorOfVariables(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_NormOneCone_VectorAffineFunction(bridged_mock, config)
    var_names = ["x", "y", "z"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    u, v = MOI.get(mock, MOI.ListOfVariableIndices())[4:5]
    MOI.set(mock, MOI.VariableName(), u, "u")
    MOI.set(mock, MOI.VariableName(), v, "v")
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
    zeros = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Zeros,
        }(),
    )
    @test length(zeros) == 2
    MOI.set(mock, MOI.ConstraintName(), zeros[1], "x_eq")
    MOI.set(mock, MOI.ConstraintName(), zeros[2], "y_eq")

    s = """
    variables: x, y, z, u, v
    nonneg: [x + -1.0u + -1.0v, u + -1.0y, v + -1.0z, u + y, v + z] in VecMathOptInterface.Nonnegatives(5)
    x_eq: [-1.0 + x] in VecMathOptInterface.Zeros(1)
    y_eq: [-0.5 + y] in VecMathOptInterface.Zeros(1)
    maxobjective: y + z
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        [var_names; "u"; "v"],
        ["nonneg", "x_eq", "y_eq"],
    )
    normone = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.NormOneCone,
        }(),
    )
    @test length(normone) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), normone[1], "normone")
    zeros = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Zeros,
        }(),
    )
    @test length(zeros) == 2
    MOI.set(bridged_mock, MOI.ConstraintName(), zeros[1], "x_eq")
    MOI.set(bridged_mock, MOI.ConstraintName(), zeros[2], "y_eq")

    s = """
    variables: x, y, z
    normone: [1.0x, y, z] in VecMathOptInterface.NormOneCone(3)
    x_eq: [-1.0 + x] in VecMathOptInterface.Zeros(1)
    y_eq: [-0.5 + y] in VecMathOptInterface.Zeros(1)
    maxobjective: y + z
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["normone", "x_eq", "y_eq"],
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.NormOneCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [4.0, 1.0, -2.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), u) == 1
            @test MOI.get(mock, MOI.VariablePrimalStart(), v) == 2
            @test MOI.get(mock, attr, nonneg[1]) == Float64[1, 0, 4, 2, 0]
        else
            @test MOI.get(mock, attr, nonneg[1]) == Float64[4, 0, 2, 1, 0]
        end
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),),
    )
    return
end

function test_conic_NormOneCone()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NormOne{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [4, -1, -1, -1, 1, 1, 1],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [vcat(ones(4), zeros(3)), ones(3)],
        )
    MOI.Test.test_conic_NormOneCone(bridged_mock, config)
    var_names = ["x", "y1", "y2", "y3"]
    var_names_all = vcat(var_names, "z1", "z2", "z3")
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names_all,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    @test length(nonneg) == 2
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg1")
    MOI.set(mock, MOI.ConstraintName(), nonneg[2], "nonneg2")
    s = """
    variables: x, y1, y2, y3, z1, z2, z3
    nonneg1: [x + -1.0 + -1.0z1 + -1.0z2 + -1.0z3, z1 + -1.0y1 + -2.0, z2 + -1.0y2 + -2.0, z3 + -1.0y3 + -2.0, z1 + y1 + 2.0, z2 + y2 + 2.0, z3 + y3 + 2.0] in VecMathOptInterface.Nonnegatives(7)
    nonneg2: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in VecMathOptInterface.Nonnegatives(3)
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names_all,
        ["nonneg1", "nonneg2"],
    )
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    normone = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.NormOneCone,
        }(),
    )
    @test length(normone) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), normone[1], "normone")
    nonneg = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    @test length(nonneg) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), nonneg[1], "nonneg")
    s = """
    variables: x, y1, y2, y3
    normone: [x + -1.0, y1 + 2.0, y2 + 2.0, y3 + 2.0] in VecMathOptInterface.NormOneCone(4)
    nonneg: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in VecMathOptInterface.Nonnegatives(3)
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["normone", "nonneg"],
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.NormOneCone,
            }(),
        ),
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = (
            attr isa MOI.ConstraintPrimalStart ? vcat(3, ones(3)) :
            vcat(1, fill(-1, 3))
        )
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        nonneg_value = (
            attr isa MOI.ConstraintPrimalStart ?
            vcat(zeros(4), fill(2.0, 3)) : vcat(ones(4), zeros(3))
        )
        @test MOI.get(mock, attr, nonneg[1]) ≈ nonneg_value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        4,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 1),),
    )
    return
end

end  # module

TestConstraintNormInfinity.runtests()


text:
@testset "Generating base HMM model..." begin
        d = ChTreebank()
	m = KongYiji.HMM(d, KongYiji.dir("usrdict"), 1.)
        file = KongYiji.dir("hmm.jld2")
        save(file, "m", m)
        @time m2 = KongYiji.HMM()
        @assert m == m2

        x = KongYiji.rawsents(d)
        z = KongYiji.wordsents(d)
        y = m(x)
        tb = KongYiji.HmmScoreTable(z, y)
        println(tb)

end

#=
@testset "Cross validating HMM on CTB..." begin
        @time ctb = ChTreebank()
        nb = 10
        batches = KongYiji.foldbatch(ctb, nb)
        tbs = KongYiji.HmmScoreTable[]
        @showprogress 1 "Cross Validating HMM..." for (tr, te) in batches
                hmm = KongYiji.HMM(tr)
                hmm = KongYiji.HMM(;base=hmm)
                x = KongYiji.rawsents(te)
                z = KongYiji.wordsents(te)
                y = hmm(x)
                push!(tbs, KongYiji.HmmScoreTable(z, y))
                println(tbs[end])
        end
        hmmscoretable = sum(tbs)
        println(hmmscoretable)
end

@testset "Test HMM exported debug infos..." begin
        tk = KongYiji.HMM()
        @show postable(tk)
        @show hprtable(tk)
        @show h2htable(tk)
        @show h2vtable(tk, "NT")
        @show v2htable(tk, "中国")
end

@testset "Cross validating HMM ..." begin
        @time d = ZhuXian()
        nb = 10
        batches = KongYiji.foldbatch(d, nb)
        tbs = KongYiji.HmmScoreTable[]
        @showprogress 1 "Cross Validating HMM..." for (tr, te) in batches
                hmm = KongYiji.HMM(tr, .99, KongYiji.dir("usrdict"))
                x = KongYiji.rawsents(te)
                z = KongYiji.wordsents(te)
                y = hmm(x)
                push!(tbs, KongYiji.HmmScoreTable(z, y))
                println(tbs[end])
        end
        hmmscoretable = sum(tbs)
        println(hmmscoretable)
end

@testset "Generating HMM ..." begin
        @time d = ChTreebank()
        m = HMM(d, 1., KongYiji.dir("usrdict"))
        save(KongYiji.dir("hmm.jld2"), "m", m)
        m2 = HMM()
        @test m == m2
        x = KongYiji.rawsents(d)
        z = KongYiji.wordsents(d)
        y = m(x)
        tb = KongYiji.HmmScoreTable(z, y)
        println(tb)
end

@testset "Testing HMM on icwb..." begin
        m = HMM()
        testonicwb(m)
end

@testset "Testing HMM Tagger..." begin
	#x = split("我  扔  了  两颗  手榴弹  ，  他  一下子  出  溜  下去  。")
        x = "唐山东日新能源材料有限公司"
	m = HMM()
	println(m(x; withpos=true))
end

=#
@testset "Names..." begin
        m = HMM()
        is = joinpath("D:\\", "project-names.csv")
        os = joinpath("D:\\", "project-names-segmented.csv")
        xs = [x for x in eachline(is)]
        ys = m(xs; withpos=true)
        open(os, "w") do io
        	for y in ys
        		for w in y; print(io, w, ' '); end
        		println(io)
		end
	end
end

#=

=#

ext: jl


text:
using Wright
using Test

@testset "Wright function" begin
    @test isapprox(wright(1, 2, -1), 0.57672)
    @test isapprox(wright(1, 1/2, -1), -0.23479)
    @test isapprox(wright(1, 2, -9/16), 0.74392)
    @test isapprox(wright(1, 1/2, -9/16), 0.03991)
    @test isapprox(wright(1, 3/2, -121/100), -0.14475)
    @test isapprox(wright(-1/2, 1/2, -5/2), 0.11826)
    @test isapprox(wright(-1/2, 1/2, -1/2), 0.53001)
    @test isapprox(wright(-1/3, 2/3, -5/2), 0.10772)
    @test isapprox(wright(-1/3, 2/3, -1/2), 0.55633)
    @test isapprox(wright(1, 2, -2500), -0.002)
end


text:
using Base: @deprecate

# BEGIN TimeZones 0.9 deprecations

@deprecate build(version::AbstractString, regions) build(version) false

# END TimeZones 0.9 deprecations


text:
using GLPK

@testset "hydrovalleymodel" begin
    @testset "Test kwargs" begin
        params = create_param( stages                  = 3,
                            model_constructor_grid  = DCPPowerModel,
                            post_method             = PowerModels.post_opf,
                            optimizer                  = GLPK.Optimizer)
        # no hydro data
        @test_throws Exception hydrovalleymodel(Dict("powersystem" => Dict()), params)

    end
end


text:
# p32.jl - solve u_xx = exp(4x), u(-1)=0, u(1)=1 (compare p13.jl)

N = 16;
(D,x) = cheb(N);
D2 = D^2;
D2 = D2[2:N,2:N];                   # boundary conditions
f = exp(4*x[2:N]);
u = D2\f;                           # Poisson eq. solved here
u = [0;u;0] + (x+1)/2;
clf();
plot(x,u,".",markersize=10);
xx = -1:.01:1;
uu = polyval(polyfit(x,u),xx);
plot(xx,uu); grid(true);
exact = @. (exp(4*xx) - sinh(4)*xx - cosh(4))/16 + (xx+1)/2;
title("max err = $(signif(norm(uu-exact,Inf),4))",fontsize=12);


text:
export iauGc2gde
"""
Transform geocentric coordinates to geodetic for a reference
ellipsoid of specified form.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   a       double     equatorial radius (Notes 2,4)
   f       double     flattening (Note 3)
   xyz     double[3]  geocentric vector (Note 4)

Returned:
   elong   double     longitude (radians, east +ve)
   phi     double     latitude (geodetic, radians)
   height  double     height above ellipsoid (geodetic, Note 4)

Returned (function value):
            int        status:  0 = OK
                              -1 = illegal f
                              -2 = illegal a

Notes:

   1. This function is based on the GCONV2H Fortran subroutine by
   Toshio Fukushima (see reference).

   2. The equatorial radius, a, can be in any units, but meters is
   the conventional choice.

   3. The flattening, f, is (for the Earth) a value around 0.00335,
   i.e. around 1/298.

   4. The equatorial radius, a, and the geocentric vector, xyz,
   must be given in the same units, and determine the units of
   the returned height, height.

   5. If an error occurs (status < 0), elong, phi and height are
   unchanged.

   6. The inverse transformation is performed in the function
   iauGd2gce.

   7. The transformation for a standard ellipsoid (such as WGS84) can
   more conveniently be performed by calling iauGc2gd, which uses a
   numerical code to identify the required A and F values.

Reference:

   Fukushima, T., "Transformation from Cartesian to geodetic
   coordinates accelerated by Halley's method", J.Geodesy (2006)
   79: 689-693

This revision:  2014 November 7

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauGc2gde ( double a, double f, double xyz[3],
#                 double *elong, double *phi, double *height )

function iauGc2gde(a::Real, f::Real, xyz::AbstractVector{<:Real})
   ref_elong  = Ref{Float64}(0.0)
   ref_phi    = Ref{Float64}(0.0)
   ref_height = Ref{Float64}(0.0)

   status = ccall((:iauGc2gde, libsofa_c), Cint,
         (Cdouble, Cdouble, Ptr{Cdouble},
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
         convert(Float64, a), convert(Float64, f), xyz,
         ref_elong, ref_phi, ref_height)

   return status, ref_elong[], ref_phi[], ref_height[]
 end


text:
module Outer

baremodule Inner

import ..Outer: Outer
Core.include(Outer, "include1.jl")

end # module Inner

end # module Outer


text:
mutable struct RNumber{T<:Number,R} <: Number
    v::T
    record::R
end
recorded(::Type{E}, c::AbstractClock, n::Number) where {E<:AbstractEntry} =
    RNumber(n, Record{E}(c, n))

Base.parent(x::RNumber) = x.v
getrecord(x::RNumber) = x.record

mutable struct RReal{T<:Real,R} <: Real
    v::T
    record::R
end
recorded(::Type{E}, c::AbstractClock, n::Real) where {E<:AbstractEntry} =
    RReal(n, Record{E}(c, n))

Base.parent(x::RReal) = x.v
getrecord(x::RReal) = x.record

"""
    RecordedNumber{T}

A `Union` of recorded numbers with type `T`.

!!! info

    `RecordedNumber{S} <: T` will always return `false` where `T` is a subtype
    of `Number`, even if `S <: T`.
    There is a function `issubtype(RecordedNumber{S}, T)`  which would return `true`
    if `S <: T`. Besides, `isnum(T, n::RecordedNumber{S})` would return `true` if
    `S <: T`.

!!! note

    Store a value `v` to an `RecordedNumber` `x` by `x[] = v` or `x[1] = v` instead of
    `x = v`.
"""
const RecordedNumber{T} = Union{RNumber{T},RReal{T}}
const RN = RecordedNumber

"""
    getentries(x::RecordedNumber)

Get entries of a recorded number `x`.
"""
getentries(x::RecordedNumber) = parent(getrecord(x))

"""
    issubtype(S, T) -> Bool

Similar to `S <: T`, but for a `S` where `S <: RecordedNumber{P}`, return `true` if
`P <: T`.
"""
issubtype(::Type{S}, ::Type{T}) where {T,S} = S <: T
issubtype(::Type{<:RecordedNumber{S}}, ::Type{T}) where {T,S} = S <: T

"""
    isnum(T, x) -> Bool

Similar to `x is T`, but for `x` where `x::RecordedNumber{S}`, return `true` if `S <:T`.
"""
isnum(::Type{T}, t::Number) where {T} = t isa T
isnum(::Type{T}, t::RecordedNumber{S}) where {T,S} = S <: T

Base.convert(::Type{T}, n::RecordedNumber{T}) where {T<:Number} = state(n)
Base.convert(::Type{T}, n::RecordedNumber) where {T<:Number} = convert(T, state(n))
Base.convert(::Type{T}, n::RecordedNumber) where {P,T<:RecordedNumber{P}} =
    convert(P, state(n))

Base.promote_rule(::Type{<:RN{S}}, ::Type{<:RN{T}}) where {S<:Number,T<:Number} =
    promote_type(S, T)
Base.promote_rule(::Type{<:RN{S}}, ::Type{T}) where {S<:Number,T<:Number} =
    promote_type(S, T)

Base.show(io::IO, x::RecordedNumber) = show(io, state(x))
Base.show(io::IO, ::MIME"text/plain", x::RecordedNumber) = show(io, state(x))

_unpack(x::RecordedNumber) = _unpack(getentries(x))

@inline state(n::RecordedNumber{T}) where {T} = parent(n)::T

@inline Base.getindex(x::RecordedNumber) = state(x)
@inline function Base.getindex(x::RecordedNumber, i::Integer)
    @boundscheck isone(i) || throw(BoundsError(x, i))
    return parent(x)
end
@inline function Base.getindex(x::RecordedNumber, I::Integer...)
    @boundscheck all(isone, I) || throw(BoundsError(x, I))
    return parent(x)
end
@inline Base.getindex(x::RecordedNumber, ::CartesianIndex{0}) = state(x) # to fix ambiguity
Base.@propagate_inbounds Base.getindex(x::RecordedNumber, I::CartesianIndex) = x[I.I...]

@inline Base.setindex!(x::RecordedNumber, v) = (getrecord(x)[] = v; x.v = v; x)
@inline function Base.setindex!(x::RecordedNumber, v, i::Integer)
    @boundscheck isone(i) || throw(BoundsError(x, i))
    getrecord(x)[] = v
    x.v = v
    return x
end
@inline function Base.setindex!(x::RecordedNumber, v, I::Integer...)
    @boundscheck all(isone, I) || throw(BoundsError(x, I))
    getrecord(x)[] = v
    x.v = v
    return x
end
Base.@propagate_inbounds Base.setindex!(x::RecordedNumber, v, I::CartesianIndex) =
    setindex!(x, v, I.I...)

for op in (:+, :-, :conj, :real, :imag, :float)
    @eval @inline Base.$op(x::RecordedNumber) = $op(parent(x))
end

for op in (:+, :-, :*, :/, :\, :^, :(==))
    @eval @inline Base.$op(x::RecordedNumber{T}, y::RecordedNumber{T}) where {T} =
        $op(parent(x), parent(y))
end

Base.hash(x::RecordedNumber, h::UInt) = hash(state(x), h)


text:
@testset "_parse_usecols" begin
    dta = "$(@__DIR__)/../data/sample.dta"
    f = read_dta(dta)
    @test _parse_usecols(f, :dtime) == _parse_usecols(f, "dtime") == (4,)
    @test_throws ArgumentError _parse_usecols(f, :time)
    @test _parse_usecols(f, [:dtime]) == _parse_usecols(f, ["dtime"]) == [4]
    @test_throws ArgumentError _parse_usecols(f, [:time])
    @test _parse_usecols(f, 1) == (1,)
    @test_throws ArgumentError _parse_usecols(f, 8)
    @test _parse_usecols(f, [1, 2]) == 1:2
    @test_throws ArgumentError _parse_usecols(f, 1:8)
end

@testset "readstat dta" begin
    dta = "$(@__DIR__)/../data/sample.dta"
    d = readstat(dta)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,120)) == """
        5×7 ReadStatTable:
         Row │ mychar    mynum      mydate                dtime         mylabl           myord               mytime
             │ String  Float64       Date?            DateTime?  Labeled{Int8}  Labeled{Int8?}             DateTime
        ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────
           1 │      a      1.1  2018-05-06  2018-05-06T10:10:10           Male             low  1960-01-01T10:10:10
           2 │      b      1.2  1880-05-06  1880-05-06T10:10:10         Female          medium  1960-01-01T23:10:10
           3 │      c  -1000.3  1960-01-01  1960-01-01T00:00:00           Male            high  1960-01-01T00:00:00
           4 │      d     -1.4  1583-01-01  1583-01-01T00:00:00         Female             low  1960-01-01T16:10:10
           5 │      e   1000.3     missing              missing           Male         missing  2000-01-01T00:00:00"""

    @test length(varlabels(d)) == 7
    @test varlabels(d)[:myord] == "ordinal"
    @test length(varformats(d)) == 7
    @test varformats(d)[:mytime] == "%tcHH:MM:SS"
    @test val_label_keys(d)[:mylabl] == "mylabl"
    @test val_label_dict(d)["mylabl"] == d.mylabl.labels
    @test filelabel(d) == "A test file"
    ts = filetimestamp(d)
    @test minute(ts) == 36
    @test fileext(d) == ".dta"

    @test sprint(show, getmeta(d)) == "ReadStatMeta"
    ss = sprint(show, MIME("text/plain"), getmeta(d))
    @test first(ss, 40) == "ReadStatMeta:\n  variable labels:    Dict"
    @test last(ss, 27) == "\n  file extension:     .dta"

    df = DataFrame(d)
    @test all(n->isequal(df[!, n], getproperty(d, n)), columnnames(d))
    df = DataFrame(d, copycols=false)
    @test all(n->df[!, n] === getproperty(d, n), columnnames(d))

    d = readstat(dta, usecols=Int[])
    @test sprint(show, d) == "0×0 ReadStatTable"
    @test isempty(varlabels(d))
    @test length(val_label_dict(d)) == 2

    d = readstat(dta, usecols=1:3, convert_datetime=3)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,120)) == """
        5×3 ReadStatTable:
         Row │ mychar    mynum      mydate
             │ String  Float64       Date?
        ─────┼─────────────────────────────
           1 │      a      1.1  2018-05-06
           2 │      b      1.2  1880-05-06
           3 │      c  -1000.3  1960-01-01
           4 │      d     -1.4  1583-01-01
           5 │      e   1000.3     missing"""

    d = readstat(dta, usecols=[:dtime, :mylabl], convert_datetime=false, apply_value_labels=10:20)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,120)) == """
        5×2 ReadStatTable:
         Row │       dtime  mylabl
             │    Float64?    Int8
        ─────┼─────────────────────
           1 │  1.84122e12       1
           2 │  -2.5136e12       2
           3 │         0.0       1
           4 │ -1.18969e13       2
           5 │     missing       1"""

    d = readstat(dta, usecols=["dtime", "mylabl"], convert_datetime=:dtime, apply_value_labels=["mylabl"])
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,120)) == """
            5×2 ReadStatTable:
             Row │               dtime         mylabl
                 │           DateTime?  Labeled{Int8}
            ─────┼────────────────────────────────────
               1 │ 2018-05-06T10:10:10           Male
               2 │ 1880-05-06T10:10:10         Female
               3 │ 1960-01-01T00:00:00           Male
               4 │ 1583-01-01T00:00:00         Female
               5 │             missing           Male"""

    d = readstat(dta, usecols=:myord, missingvalue=0)
    @test d.myord.values[5] == 0
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,120)) == """
        5×1 ReadStatTable:
         Row │          myord
             │ Labeled{Int64}
        ─────┼────────────────
           1 │            low
           2 │         medium
           3 │           high
           4 │            low
           5 │              0"""

    @test_throws ArgumentError readstat("$(@__DIR__)/../data/README.md")
end

@testset "readstat sav" begin
    sav = "$(@__DIR__)/../data/sample.sav"
    d = readstat(sav)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,150)) == """
        5×7 ReadStatTable:
         Row │ mychar    mynum               mydate                dtime            mylabl             myord               mytime
             │ String  Float64            DateTime?            DateTime?  Labeled{Float64}  Labeled{Float64}            DateTime?
        ─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
           1 │      a      1.1  2018-05-06T00:00:00  2018-05-06T10:10:10              Male               low  1582-10-14T10:10:10
           2 │      b      1.2  1880-05-06T00:00:00  1880-05-06T10:10:10            Female            medium  1582-10-14T23:10:10
           3 │      c  -1000.3  1960-01-01T00:00:00  1960-01-01T00:00:00              Male              high  1582-10-14T00:00:00
           4 │      d     -1.4  1583-01-01T00:00:00  1583-01-01T00:00:00            Female               low  1582-10-14T16:10:10
           5 │      e   1000.3              missing              missing              Male               low              missing"""
end

@testset "readstat por" begin
    por = "$(@__DIR__)/../data/sample.por"
    d = readstat(por)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,150)) == """
        5×7 ReadStatTable:
         Row │ MYCHAR    MYNUM               MYDATE                DTIME            MYLABL             MYORD               MYTIME
             │ String  Float64            DateTime?            DateTime?  Labeled{Float64}  Labeled{Float64}            DateTime?
        ─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
           1 │      a      1.1  2018-05-06T00:00:00  2018-05-06T10:10:10              Male               low  1582-10-14T10:10:10
           2 │      b      1.2  1880-05-06T00:00:00  1880-05-06T10:10:10            Female            medium  1582-10-14T23:10:10
           3 │      c  -1000.3  1960-01-01T00:00:00  1960-01-01T00:00:00              Male              high  1582-10-14T00:00:00
           4 │      d     -1.4  1583-01-01T00:00:00  1583-01-01T00:00:00            Female               low  1582-10-14T16:10:10
           5 │      e   1000.3              missing              missing              Male               low              missing"""
end

@testset "readstat sas7bdat" begin
    sas7 = "$(@__DIR__)/../data/sample.sas7bdat"
    d = readstat(sas7)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,150)) == """
        5×7 ReadStatTable:
         Row │ mychar    mynum      mydate                dtime   mylabl    myord               mytime
             │ String  Float64       Date?            DateTime?  Float64  Float64            DateTime?
        ─────┼─────────────────────────────────────────────────────────────────────────────────────────
           1 │      a      1.1  2018-05-06  2018-05-06T10:10:10      1.0      1.0  1960-01-01T10:10:10
           2 │      b      1.2  1880-05-06  1880-05-06T10:10:10      2.0      2.0  1960-01-01T23:10:10
           3 │      c  -1000.3  1960-01-01  1960-01-01T00:00:00      1.0      3.0  1960-01-01T00:00:00
           4 │      d     -1.4  1583-01-01  1583-01-01T00:00:00      2.0      1.0  1960-01-01T16:10:10
           5 │      e   1000.3     missing              missing      1.0      1.0              missing"""
end

@testset "readstat xpt" begin
    xpt = "$(@__DIR__)/../data/sample.xpt"
    d = readstat(xpt)
    @test sprint(show, MIME("text/plain"), d, context=:displaysize=>(15,150)) == """
        5×7 ReadStatTable:
         Row │ MYCHAR    MYNUM      MYDATE                DTIME   MYLABL    MYORD               MYTIME
             │ String  Float64       Date?            DateTime?  Float64  Float64            DateTime?
        ─────┼─────────────────────────────────────────────────────────────────────────────────────────
           1 │      a      1.1  2018-05-06  2018-05-06T10:10:10      1.0      1.0  1960-01-01T10:10:10
           2 │      b      1.2  1880-05-06  1880-05-06T10:10:10      2.0      2.0  1960-01-01T23:10:10
           3 │      c  -1000.3  1960-01-01  1960-01-01T00:00:00      1.0      3.0  1960-01-01T00:00:00
           4 │      d     -1.4  1583-01-01  1583-01-01T00:00:00      2.0      1.0  1960-01-01T16:10:10
           5 │      e   1000.3     missing              missing      1.0      1.0              missing"""
end


text:
export perms, automaton_constructor, islight
using LightGraphs
"""
    perms(p)

Generates a vector of all sequences of length `p` made up of elements `-1` or `1`.
"""
function perms(p)
    vals = []
    for i in 0:p
        out = unique(permutations([ones(p-i); -ones(i)]))
        for j in 1:length(out)
            push!(vals, out[j])
        end
    end
    return vals
end

"""
    automaton_constructor(Σ::AbstractVector{<:AbstractMatrix})

Generates an automaton tracking admissible transitions between states. `Σ` is a vector of matrices corresponding to a TAR model, assuming a default order generated by [`perms`](@ref).

    automaton_constructor(Σ::AbstractVector{<:AbstractMatrix}, seqlist)

Generates an automaton tracking admissible transitions between states given a custom list of regime sequences.
"""
function automaton_constructor(Σ::AbstractVector{<:AbstractMatrix})
    G = LightAutomaton(length(Σ))
    p = Int(log2(length(Σ)))
    seqlist = perms(p)
    for i in 1:length(seqlist)
        for j in 1:length(seqlist)
            if seqlist[i][1:end-1] == seqlist[j][2:end]
                add_transition!(G, i, j, i)
            end
        end
    end
    return G
end

function automaton_constructor(Σ::AbstractVector{<:AbstractMatrix}, seqlist)
    G = LightAutomaton(length(Σ))
    for i in 1:length(seqlist)
        for j in 1:length(seqlist)
            if seqlist[i][1:end-1] == seqlist[j][2:end]
                add_transition!(G, i, j, i)
            end
        end
    end
    return G
end

function HybridSystems.event(A::AbstractAutomaton) end

function HybridSystems.event(A::LightAutomaton, q, r)
    if has_transition(A, q, r)
        if length(A.Σ[Edge(q, r)]) > 1
            @warn "Edge associated with more than one label! Reporting first label only."
        end
        return collect(values(A.Σ[Edge(q, r)]))[1]
    end
    nothing
end

function HybridSystems.event(A::OneStateAutomaton, q, r)
    @assert q == 1
    @assert r == 1
    collect(1:A.nt)
end

islight(A::AbstractAutomaton) = typeof(A) <: LightAutomaton


text:
module TEMPsolver_PDE

using ForwardDiff
using LinearAlgebra

include("FVMPDE_problem.jl")
export FVMPDEProblem

end


text:
type Arr <: Op; init; initialized; out; Arr(;init=nothing,o...)=new(init,false); end
Kenv.kdef(:arr,Arr)
infersize(a::Arr,ysize)=tuple(size(a.init))
ninputs(::Arr)=0
canoverwrite(::Arr)=false
back_reads_x(::Arr)=false
back_reads_y(::Arr)=false
back(::Arr,dy;o...)=nothing

function forw(a::Arr,y;o...)
    if !a.initialized
        a.out=copysync!(y, a.init)
        a.initialized=true
    end
    a.out===y || error("Constant modified")
    return y
end


text:
"""
    optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts,
        optimopts::Optim.Options, minimize = true)

Optimize the rotation angles of a particle collection for minimization or
maximization (depending on `minimize`) of the field intensity at `points`.
`optimopts` and `method` define the optimization method, convergence criteria,
and other optimization parameters.
`adjoint` dictates whether the gradient is calculated using the adjoint method
(faster) or directly.
Returns an object of type `Optim.MultivariateOptimizationResults`.
"""
function optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts,
                    optimopts::Optim.Options, method; minimize = true, adjoint = true)

    #stuff that is done once
    verify_min_distance(shapes, centers, ids, points) || error("Particles are too close.")
    mFMM, scatteringMatrices, scatteringLU, buf =
        prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopts)
    Ns = size(centers,1)
    H = optimizationHmatrix(points, centers, Ns, P, k0)
    α = u2α(k0, ui, centers, P)
    uinc_ = uinc(k0, points, ui)
    # Allocate buffers
    shared_var = OptimBuffer(Ns,P,size(points,1))
    initial_φs = copy(φs0)
    last_φs = similar(initial_φs)
    last_φs == initial_φs && (last_φs[1] += 1) #should never happen but has

    fopt = φs -> ifelse(minimize,1,-1)*optimize_φ_f(φs, shared_var, last_φs, α,
                                        H, points, P, uinc_, Ns, k0, centers,
                                        scatteringMatrices, ids, mFMM, fmmopts,
                                        buf)
    if adjoint
        gopt! = (grad_stor, φs) -> optimize_φ_adj_g!(grad_stor, φs, shared_var,
                                    last_φs, α, H, points, P, uinc_, Ns, k0,
                                    centers, scatteringMatrices, scatteringLU,
                                    ids, mFMM, fmmopts, buf, minimize)
    else
        gopt! = (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var,
                                    last_φs, α, H, points, P, uinc_, Ns, k0,
                                    centers, scatteringMatrices, scatteringLU,
                                    ids, mFMM, fmmopts, buf, minimize)
    end
    df = OnceDifferentiable(fopt, gopt!, initial_φs)
    optimize(df, initial_φs, method, optimopts)
end

function optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, ids, mFMM, opt, buf)
    if φs != last_φs
        copyto!(last_φs, φs)
        #do whatever common calculations and save to shared_var
        #construct rhs
        for ic = 1:Ns
            rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
            if φs[ic] == 0.0
                buf.rhs[rng] = scatteringMatrices[ids[ic]]*α[rng]
            else
                #rotate without matrix
                rotateMultipole!(view(buf.rhs,rng),view(α,rng),-φs[ic],P)
                buf.rhs[rng] = scatteringMatrices[ids[ic]]*buf.rhs[rng]
                rotateMultipole!(view(buf.rhs,rng),φs[ic],P)
            end
        end

        MVP = LinearMap{eltype(buf.rhs)}(
                (output_, x_) -> FMM_mainMVP_pre!(output_, x_,
                                    scatteringMatrices, φs, ids, P, mFMM,
                                    buf.pre_agg, buf.trans),
                Ns*(2*P+1), Ns*(2*P+1), ismutating = true)

        fill!(shared_var.β,0.0)
        shared_var.β,ch = gmres!(shared_var.β, MVP, buf.rhs,
                            restart = Ns*(2*P+1), tol = opt.tol, log = true,
							initially_zero = true) #no restart, preconditioning
        !ch.isconverged && error("FMM process did not converge")

        shared_var.f[:] = transpose(H)*shared_var.β
        shared_var.f .+= uinc_
    end
end

function optimize_φ_f(φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers,scatteringMatrices, ids, mFMM, opt, buf)
    optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns,
        k0, centers,scatteringMatrices, ids, mFMM, opt, buf)

    func = sum(abs2, shared_var.f)
end

function optimize_φ_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize)
    optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns,
        k0, centers,scatteringMatrices, ids, mFMM, opt, buf)

    MVP = LinearMap{eltype(buf.rhs)}(
            (output_, x_) -> FMM_mainMVP_pre!(output_, x_, scatteringMatrices,
                                            φs, ids, P, mFMM, buf.pre_agg,
                                            buf.trans),
            Ns*(2*P+1), Ns*(2*P+1), ismutating = true)
    #time for gradient
    shared_var.rhs_grad[:] = 0.0
    shared_var.∂β[:] = 0.0
    D = -1.0im*collect(-P:P)
    tempn = Array{Complex{Float64}}(undef, 2*P+1)
    for n = 1:Ns
        #compute n-th gradient
        rng = (n-1)*(2*P+1) .+ (1:2*P+1)
        rotateMultipole!(tempn, shared_var.β[rng], -φs[n], P)
        tempn[:] = scatteringLU[ids[n]]\tempn #LU decomp with pivoting
        tempn[:] .*= -D
        v = view(shared_var.rhs_grad, rng)
        mul!(v, scatteringMatrices[ids[n]], tempn)
        rotateMultipole!(v, φs[n], P)
        v[:] += D.*shared_var.β[rng]

        shared_var.∂β[:,n], ch = gmres!(view(shared_var.∂β,:,n), MVP,
                                    shared_var.rhs_grad, restart = Ns*(2*P+1),
                                    tol = 10*opt.tol, log = true,
									initially_zero = true)
        #@warn("using dbdn_tol = 10*opt.tol = $(10*opt.tol)")

        if ch.isconverged == false
            display("FMM process did not converge for partial derivative $n/$Ns. ")
            error("..")
        end
        #prepare for next one
        v[:] = 0.0im
    end

    grad_stor[:] = ifelse(minimize, 2, -2)*
                    real(transpose(shared_var.∂β)*(H*conj(shared_var.f)))
end

function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize)
    optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns,
        k0, centers,scatteringMatrices, ids, mFMM, opt, buf)

    MVP = LinearMap{eltype(buf.rhs)}(
            (output_, x_) -> FMM_mainMVP_transpose!(output_, x_,
                                scatteringMatrices, φs, ids, P, mFMM,
                                buf.pre_agg, buf.trans),
            Ns*(2*P+1), Ns*(2*P+1), ismutating = true)

    #build rhs of adjoint problem
    shared_var.rhs_grad[:] = -(H*conj(shared_var.f))
    #solve adjoint problem
    λadj, ch = gmres(MVP, shared_var.rhs_grad, restart = Ns*(2*P+1),
                    tol = opt.tol, log = true)
    # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1),
    #                 tol = opt.tol, log = true, initially_zero = true)
    if ch.isconverged == false
        display("FMM process did not converge for adjoint system.")
        error("..")
    end

    #note: this assumes that there are exactly 2P+1 non zero elements in ∂X. If
    #not, v must be (2P+1)Ns × 1.
    D = -1.0im*collect(-P:P)
    v = Array{Complex{Float64}}(undef, 2P+1) #TODO: minimize dynamic alloc
    for n = 1:Ns
        #compute n-th element of gradient
        rng = (n-1)*(2*P+1) .+ (1:2*P+1)
        rotateMultipole!(v, view(shared_var.β,rng), -φs[n], P)
        v[:] = scatteringLU[ids[n]]\v #LU decomp with pivoting
        v[:] .*= -D
        v[:] = scatteringMatrices[ids[n]]*v
        rotateMultipole!(v, φs[n], P)
        v[:] += D.*shared_var.β[rng]

        grad_stor[n] = ifelse(minimize, -2, 2)*real(transpose(view(λadj,rng))*v)
        #prepare for next one - #TODO: check why this is here
        v[:] .= 0.0im
    end
end

function optimizationHmatrix(points, centers, Ns, P, k0)
    points_moved = Array{Float64}(undef, 2)
    H = Array{Complex{Float64}}(undef, Ns*(2*P+1), size(points,1))
    for ic = 1:Ns, i = 1:size(points,1)
        points_moved[1] = points[i,1] - centers[ic,1]
        points_moved[2] = points[i,2] - centers[ic,2]
        r_angle = atan(points_moved[2], points_moved[1])
        kr = k0*hypot(points_moved[1], points_moved[2])

        ind = (ic-1)*(2*P+1) + P + 1
        H[ind,i] = besselh(0, kr)
        for l = 1:P
            b = besselh(l, kr)
            H[ind + l,i] = b*exp(1.0im*l*r_angle)
            H[ind - l,i] = b*(-1)^l*exp(-1.0im*l*r_angle)
        end
    end
    H
end

function prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopt)
    #setup FMM reusal
    Ns = size(centers,1)
    groups, boxSize = divideSpace(centers, fmmopt)
    P2, Q = FMMtruncation(fmmopt.acc, boxSize, k0)
    mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true)
    scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids)
    scatteringLU = [lu(scatteringMatrices[iid]) for iid = 1:length(shapes)]
    buf = FMMbuffer(Ns,P,Q,length(groups))
    return mFMM, scatteringMatrices, scatteringLU, buf
end


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
prodrank = prod(1:rank)

# Not possible to PROD a Char (and neither Int8 with OpenMPI)
for T in setdiff([Base.uniontypes(MPI.MPIDatatype)...], [Char, Int8, UInt8])
    A = ArrayType{T}(fill(T(rank+1), 4))
    B = similar(A)
    MPI.Exscan!(A, B, *, comm)
    if rank > 0
        @test B == ArrayType{T}(fill(T(prodrank), 4))
    end

    B = MPI.Exscan(A, *, comm)
    @test B isa ArrayType{T}
    if rank > 0
        @test B == ArrayType{T}(fill(T(prodrank), 4))
    end

    B = MPI.Exscan(T(rank+1), *, comm)
    if rank > 0
        @test B[1] == prodrank
    end
end

MPI.Finalize()
@test MPI.Finalized()


text:
function readoutput(calculation::Calculation{Wannier90}, file; kwargs...)
    return wan_read_output(file; kwargs...)
end

#THIS IS THE MOST HORRIBLE FUNCTION I HAVE EVER CREATED!!!
#extracts only atoms with projections
function extract_atoms(atoms_block, proj_block, cell::Mat3, spinors = false)
    if atoms_block.name == :atoms_cart
        cell = Mat3(Matrix(1.0Ang * I, 3, 3))
    end
    projections_ = proj_block.data
    atoms = atoms_block.data
    t_start = 1
    out_ats = Atom[]
    if projections_ !== nothing
        t_ats = Atom[]
        for (proj_at, projs) in projections_
            for (pos_at, pos) in atoms
                for ps in pos
                    if proj_at != pos_at
                        continue
                    end
                    for proj in projs
                        size = spinors ? 2 * length(proj) : length(proj)
                        push!(t_ats,
                              Atom(; name = pos_at, position_cart = cell * ps,
                                   position_cryst = ps,
                                   projections = [Projection(Structures.orbital(proj),
                                                             t_start, t_start + size - 1)]))
                        t_start += size
                    end
                end
            end
        end
        for (pos_at, pos) in atoms
            for ps in pos
                same_ats = Atom[]
                for at in t_ats
                    if at.position_cart == cell * ps
                        push!(same_ats, at)
                    end
                end
                if isempty(same_ats)
                    continue
                end
                if length(same_ats) > 1
                    for at in same_ats[2:end]
                        push!(same_ats[1].projections, at.projections[1])
                    end
                end
                push!(out_ats, same_ats[1])
            end
        end
    else
        for (pos_at, pos) in atoms
            for p in pos
                push!(out_ats,
                      Atom(pos_at; position_cart = cell * p, position_cryst = p,
                           projections = Projection[]))
            end
        end
    end
    return out_ats
end

function extract_structure(cell_block, atoms_block, projections_block,
                           spinors = false)
    if atoms_block == nothing || cell_block == nothing
        return nothing
    end
    if cell_block.option == :bohr
        cell = cell_block.data' .* 1bohr
    else
        cell = cell_block.data' .* 1Ang
    end

    atoms = extract_atoms(atoms_block, projections_block, cell, spinors)
    return Structure(Mat3(cell), atoms)
end

function wan_read_calculation(::Type{T}, f::IO) where {T}
    flags       = Dict{Symbol,Any}()
    data        = InputData[]
    atoms_block = nothing
    cell_block  = nothing
    proj_block  = nothing
    line        = readline(f)
    while !eof(f)
        @label start_label

        if occursin("!", line) ||
           line == "" ||
           occursin("end", lowercase(line)) ||
           occursin("#", line)
            line = readline(f)
            continue
        end

        if occursin("begin", lowercase(line))
            block_name = Symbol(split(lowercase(line))[end])

            if block_name == :projections
                proj_dict = Tuple{Symbol,Array{String,1}}[]
                line      = readline(f)
                while !occursin("end", lowercase(line))
                    if occursin("!", line) || line == ""
                        line = readline(f)
                        continue
                    end
                    if occursin("random", line)
                        proj_block = InputData(:projections, :random, nothing)
                        line = readline(f)
                        break
                    else
                        split_line  = strip_split(line, ':')
                        atom        = Symbol(split_line[1])
                        projections = [proj for proj in strip_split(split_line[2], ';')]
                        push!(proj_dict, (atom, projections))
                        line = readline(f)
                    end
                end
                proj_block = InputData(:projections, :none, proj_dict)
                @goto start_label

            elseif block_name == :kpoint_path
                line = readline(f)
                k_path_array = Array{Tuple{Symbol,Array{T,1}},1}()
                while !occursin("end", lowercase(line))
                    if occursin("!", line) || line == ""
                        line = readline(f)
                        continue
                    end
                    split_line = split(line)
                    push!(k_path_array,
                          (Symbol(split_line[1]), parse_string_array(T, split_line[2:4])))
                    push!(k_path_array,
                          (Symbol(split_line[5]), parse_string_array(T, split_line[6:8])))
                    line = readline(f)
                end
                push!(data, InputData(:kpoint_path, :none, k_path_array))
                @goto start_label

            elseif block_name == :unit_cell_cart
                line = readline(f)
                if length(split(line)) == 1
                    option = Symbol(strip(lowercase(line)))
                    line = readline(f)
                else
                    option = :ang
                end
                cell_param = Matrix{T}(undef, 3, 3)
                for i in 1:3
                    cell_param[i, :] = parse_line(T, line)
                    line = readline(f)
                end
                cell_block = InputData(:unit_cell_cart, option, Mat3(cell_param))
                # line = readline(f)
                @goto start_label

            elseif block_name == :atoms_frac || block_name == :atoms_cart
                line   = readline(f)
                atoms  = Dict{Symbol,Array{Point3{T},1}}()
                option = :ang
                while !occursin("end", lowercase(line))
                    if occursin("!", line) || line == ""
                        line = readline(f)
                        continue
                    end
                    if length(split(line)) == 1
                        option = Symbol(strip(line))
                        line = readline(f)
                        continue
                    end
                    split_line = strip_split(line)
                    atom       = Symbol(split_line[1])
                    position   = Point3(parse_string_array(T, split_line[2:4]))
                    if !haskey(atoms, atom)
                        atoms[atom] = [position]
                    else
                        push!(atoms[atom], position)
                    end
                    line = readline(f)
                end
                atoms_block = InputData(block_name, option, atoms)
                @goto start_label

            elseif block_name == :kpoints
                line     = readline(f)
                k_points = Array{NTuple{3,T},1}()
                while !occursin("end", lowercase(line))
                    if line == ""
                        line = readline(f)
                        continue
                    end
                    push!(k_points, (parse_line(T, line)...,))
                    line = readline(f)
                end
                push!(data, InputData(:kpoints, :none, k_points))
                @goto start_label
            end

        else
            if occursin("mp_grid", line)
                flags[:mp_grid] = parse_string_array(Int, split(line)[end-2:end])
            else
                flag, val = wan_parse_flag_line(line)
                flags[flag] = val
            end
        end
        line = readline(f)
    end
    return (flags=flags, data=data, atoms=atoms_block, cell=cell_block, projections=proj_block)
end

wan_read_calculation(f::IO) = wan_read_calculation(Float64, f)

"""
    wan_read_calculation(filename::String, T=Float64; run=true, exec=Exec(exec="wannier90.x"), structure_name="NoName")

Reads a `Calculation{Wannier90}` and the included `Structure` from a WANNIER90 calculation file.
"""
function wan_read_calculation(filename::String, T = Float64;
                              exec = Exec(; exec = "wannier90.x"), kwargs...)
    flags       = Dict{Symbol,Any}()
    data        = InputData[]
    atoms_block = nothing
    cell_block  = nothing
    proj_block  = nothing
    open(filename, "r") do f
        return flags, data, atoms_block, cell_block, proj_block = wan_read_calculation(T, f)
    end
    structure = extract_structure(cell_block, atoms_block, proj_block,
                                  get(flags, :spinors, false))
    dir, file = splitdir(filename)
    flags[:preprocess] = Calculations.hasflag(exec,
                                              :pp) ? true : false
    return Calculation{Wannier90}(; name = splitext(file)[1], flags = flags,
                                  data = data, exec = exec, kwargs...), structure
end

function wan_parse_array_value(eltyp, value_str)
    value_str = replace(value_str, "," => "")
    spl = split(value_str)
    arr = eltyp[]
    for i in 1:length(spl)
        if spl[i] == "t"
            spl[i] = "true"
        elseif spl[i] == "f"
            spl[i] = "false"
        end
        s = spl[i]
        if occursin("-", s)
            t = split(s, '-')
            append!(arr, collect(parse(eltyp, t[1]):parse(eltyp, t[2])))
        else
            append!(arr, parse(eltyp, s))
        end
    end
    return arr
end

function wan_parse_flag_line(line::String)
    split_line = strip_split(line, '=')
    flag       = Symbol(split_line[1])
    flagtyp    = Calculations.flagtype(Wannier90, flag)
    value      = strip(lowercase(split_line[2]), '.')
    if flagtyp != String
        value = replace(value, "d" => "e")
    end
    if flagtyp <: Vector
        parsed_val = wan_parse_array_value(eltype(flagtyp), value)
    elseif flagtyp == String
        parsed_val = value
    else
        parsed_val = parse.(eltype(flagtyp), split(value))
    end
    if length(parsed_val) == 1
        val = parsed_val[1]
    else
        val = parsed_val
    end
    return flag, val
end

function wan_write_projections(f::IO, atoms::Vector{Atom})
    write(f, "begin projections\n")
    uniats = unique(atoms)
    projs = map(x -> x.projections, uniats)
    if all(isempty.(projs))
        write(f, "random\n")
    else
        for (at, prjs) in zip(uniats, projs)
            if isempty(prjs)
                continue
            end
            write(f, Structures.projections_string(at))
            write(f, "\n")
        end
    end
    return write(f, "end projections\n")
end

"""
    save(calculation::Calculation{Wannier90}, structure, filename::String)

Writes the `Calculation{Wannier90}` and `structure` to a file, that can be interpreted by WANNIER90.
The atoms in the structure must have projections defined.
"""
function save(calculation::Calculation{Wannier90}, structure,
              filename::String)

    projs = vcat(map(structure.atoms) do x
                     ps = x.projections
                     @assert !isempty(ps) "Please first set projections for all atoms in the Structure."
                     return ps
                 end...)

    Structures.sanitize!(projs, Calculations.issoc(calculation))

    open(filename, "w") do f
        for (flag, value) in calculation.flags
            write_flag_line(f, flag, value)
        end
        write(f, "\n")

        if structure != nothing
            write(f, "begin unit_cell_cart\n")
            writedlm(f, ustrip.(structure.cell'))
            write(f, "end unit_cell_cart\n")
            write(f, "\n")
        end
        wan_write_projections(f, structure.atoms)

        write(f, "\n")
        write(f, "begin atoms_frac\n")
        for at in structure.atoms
            pos = round.(at.position_cryst, digits = 5)
            write(f, "$(at.name)  $(pos[1]) $(pos[2]) $(pos[3])\n")
        end
        write(f, "end atoms_frac\n")
        write(f, "\n")

        for block in calculation.data
            write(f, "begin $(block.name)\n")
            if block.name == :kpoint_path
                for i in 1:2:length(block.data)
                    letter1, k_points1 = block.data[i]
                    letter2, k_points2 = block.data[i+1]
                    write(f,
                          "$letter1 $(k_points1[1]) $(k_points1[2]) $(k_points1[3]) $letter2 $(k_points2[1]) $(k_points2[2]) $(k_points2[3])\n")
                end

            elseif block.name == :kpoints
                for k in block.data
                    write(f, "$(k[1]) $(k[2]) $(k[3])\n")
                end
            end
            write(f, "end $(block.name)\n\n")
        end
    end
end

function write_flag_line(f, flag, data, seperator = "=", i = "")
    flagstr = string(flag)
    if flagstr[end-1] == '_' && tryparse(Int, string(flagstr[end])) != nothing
        flagstr = flagstr[1:end-2] * "($(flagstr[end]))"
    end
    write(f, "  $flagstr$i $seperator ")

    if typeof(data) <: Array
        if length(data) % 3 == 0 && eltype(data) != Int
            write(f, "\n")
            for i in 1:3:length(data)
                write(f, "  $(data[i]) $(data[i + 1]) $(data[i + 2])\n")
            end
        else
            write(f, "  $(data[1])")
            for x in data[2:end]
                write(f, " $x")
            end
            write(f, "\n")
        end

    else #this should work for anything singular valued data such as bools, ''s and other types
        write(f, "$data\n")
    end
end

function wan_parse_disentanglement(out, line, f)
    DisTuple = NamedTuple{(:Iter, :Ω_i_1, :Ω_i, :δ, :Time),
                          Tuple{Int,Float64,Float64,Float64,Float64}}
    out[:disentanglement] = DisTuple[]
    for i in 1:4
        readline(f)
    end
    line = strip(readline(f))
    while !isempty(line)
        push!(out[:disentanglement], parse(DisTuple, split(line)[1:end-2]))
        line = strip(readline(f))
    end
end
function wan_parse_wannierise(out, line, f)
    WanTuple = NamedTuple{(:Iter, :δ_spread, :∇RMS, :Spread, :Time),
                          Tuple{Int,Float64,Float64,Float64,Float64}}
    WfcTuple = NamedTuple{(:index, :center, :Spread),Tuple{Int,Point3{Float64},Float64}}
    out[:wannierise] = WanTuple[]
    line = strip(readline(f))
    while !occursin("Final State", line)
        if occursin("CONV", line)
            push!(out[:wannierise], parse(WanTuple, split(line)[1:end-2]))
        end
        line = strip(readline(f))
    end
    out[:final_state] = WfcTuple[]
    line = strip(readline(f))
    while !occursin("Sum", line)
        line = replace(replace(replace(line, "(" => ""), ")" => ""), "," => "")
        push!(out[:final_state], parse(WfcTuple, split(line)[5:end]))
        line = strip(readline(f))
    end
end

const WAN_PARSE_FUNCS = ["Extraction of optimally-connected subspace" => wan_parse_disentanglement,
                         "Initial State" => wan_parse_wannierise]

"""
    wan_read_output(filename::AbstractString; parse_funcs::Vector{Pair{String,Function}}=Pair{String,Function}[])

Reads an outputfile for wannier.
Parsed info:
    :disentanglement,
    :wannierise,
    :final_state,
"""
function wan_read_output(filename::AbstractString;
                         parse_funcs::Vector{<:Pair{String}} = Pair{String}[])
    return parse_file(filename, WAN_PARSE_FUNCS; extra_parse_funcs = parse_funcs)
end

function outfiles(c::Calculation{Wannier90})
    files = [c.name]
    if get(c, :wannier_plot, false)
        push!(files, "UNK")
    end
    return unique(files)
end


text:
using Sched
using Sched: Priority
using Dates
using Test

# Time as Float64
_time = FloatTimeFunc

function print_time()
    println("From print_time $(time())")
end

function print_time_noparam()
    println("From print_time_noparam $(_time())")
end

function print_time_args(x)
    println("From print_time_args $(_time()) $x")
end

function print_time_kwargs(; a="default")
    println("From print_time_kwargs $(_time()) $a")
end

@testset "Priority Tests" begin

    @testset "Priority Equality" begin
        p1 = Priority(DateTime(2010, 1, 1), 0)
        p2 = Priority(DateTime(2010, 1, 1), 0)
        @test p1 == p2
    end;

    @testset "Priority Inequality" begin
        @testset "Same DateTime" begin
            p1 = Priority(DateTime(2010, 1, 1), 0)
            p2 = Priority(DateTime(2010, 1, 1), 1)
            @test p1 != p2
        end;

        @testset "Same priority level" begin
            p1 = Priority(DateTime(2010, 1, 1), 0)
            p2 = Priority(DateTime(2010, 1, 2), 0)
            @test p1 != p2
        end;

    end;

    @testset "Priority order" begin
        p1 = Priority(DateTime(2010, 1, 2), 0)
        p2 = Priority(DateTime(2010, 1, 1), 0)
        @test p2 > p1

        p1 = Priority(DateTime(2010, 1, 1), 1)
        p2 = Priority(DateTime(2010, 1, 1), 0)
        @test p2 > p1
    end;

    @testset "Sample" begin

        s = Scheduler()

        # Time as Float64
        #_time = FloatTimeFunc

        # Time as DateTime
        _time = UTCDateTimeFunc

        function print_some_times()
            t0 = _time()
            println(t0)
            enter(s, Dates.Second(10), 1, print_time_noparam)
            enter(s, Dates.Second(5), 2, print_time_args, ("positional, argument"))
            enter(s, Dates.Second(5), 1, print_time_kwargs; Dict(:a=>"keyword")...)
            run(s)
            t1 = _time()
            println(t1)
            t1 - t0
        end

        delta_t = print_some_times()
        @test delta_t >= Dates.Second(10)
        @test delta_t < Dates.Second(11)

    end;

    @testset "empty" begin
        s = Scheduler()
        @test isempty(s)
        enter(s, Dates.Second(10), 0, print_time)
        @test !isempty(s)
    end;

    @testset "cancel" begin
        s = Scheduler()
        event = enter(s, Dates.Second(10), 0, print_time)
        @test !isempty(s)
        cancel(s, event)
        @test isempty(s)
    end;

    @testset "queue" begin
        @testset "ms resolution" begin
            s = Scheduler()
            enter(s, Dates.Second(10), 1, print_time_noparam)
            enter(s, Dates.Second(5), 2, print_time_args, ("positional, argument"))
            enter(s, Dates.Second(5), 1, print_time_kwargs; Dict(:a=>"keyword")...)
            a = queue(s)
            @test length(a) == 3
            @test a[1].action == print_time_kwargs
            @test a[2].action == print_time_args
            @test a[3].action == print_time_noparam
        end;

        @testset "lower (ns) resolution" begin
            s = Scheduler(timefunc=FloatTimeFunc)
            enter(s, 10.0, 1, print_time_noparam)
            enter(s, 5.0, 2, print_time_args, ("positional, argument"))
            enter(s, 5.0, 1, print_time_kwargs; Dict(:a=>"keyword")...)
            a = queue(s)
            @test length(a) == 3
            @test a[1].action == print_time_args     # execution order is different!
            @test a[2].action == print_time_kwargs
            @test a[3].action == print_time_noparam
        end;
    end;
end;


text:
JSON.lower(hf::HavingSpec) = non_nothing_dict(hf)

"""
    EqualTo(aggregation::String, value::Real)

aggregation = value filter.
"""
struct EqualTo <: HavingSpec
    type::String
    aggregation::String
    value::Real
    EqualTo(aggregation, value) = new("equalTo", aggregation, value)
end

"""
    GreaterThan(aggregation::String, value::Real)

aggregation >= value filter.
"""
struct GreaterThan <: HavingSpec
    type::String
    aggregation::String
    value::Real
    GreaterThan(aggregation, value) = new("greaterThan", aggregation, value)
end

"""
    LessThan(aggregation::String, value::Real)

aggregation <= value filter.
"""
struct LessThan <: HavingSpec
    type::String
    aggregation::String
    value::Real
    LessThan(aggregation, value) = new("lessThan", aggregation, value)
end

"""
    DimSelector(dimension::String, value)

dimension = value filter.
"""
struct DimSelector <: HavingSpec
    type::String
    dimension::String
    value
    DimSelector(dimension, value) = new("dimSelector", dimension, value)
end

"""
    AndH(havingSpecs::Vector{<:HavingSpec})

Match all `havingSpecs`.
"""
struct AndH <: HavingSpec
    type::String
    havingSpecs::Vector{<:HavingSpec}
    AndH(havingSpecs) = new("and", havingSpecs)
end

"""
    OrH(havingSpecs::Vector{<:HavingSpec})

Match at leats one of the `havingSpecs`.
"""
struct OrH <: HavingSpec
    type::String
    havingSpecs::Vector{<:HavingSpec}
    OrH(havingSpecs) = new("or", havingSpecs)
end

"""
    NotH(havingSpec::HavingSpec)

Do not match the `havingSpec`.
"""
struct NotH <: HavingSpec
    type::String
    havingSpec::HavingSpec
    NotH(havingSpec) = new("not", havingSpec)
end


text:
using Test     #src
# # Continuous-time Controlled Invariant Set
#
#md # [![Binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/continuous_controlled.ipynb)
#md # [![nbviewer](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/continuous_controlled.ipynb)
#
# ## Introduction
#
# This reproduces reproduces the numerical result of the Section 4 of [LJ21].
#
# This example considers the continuous-time constrained linear control system:
# ```math
# \begin{aligned}
# \dot{x}_1(t) & = x_2(t)\\
# \dot{x}_2(t) & = u(t)
# \end{aligned}
# ```
# with state constraint $x \in [-1, 1]^2$ and input constraint $u \in [-1, 1]$.
#
# In order to compute controlled invariant sets for this system, we consider
# the projection onto the first two dimensions of controlled invariant sets of the
# following lifted system:
# ```math
# \begin{aligned}
# \dot{x}_1(t) & = x_2(t)\\
# \dot{x}_2(t) & = x_3(t)\\
# \dot{x}_3(t) & = u(t)
# \end{aligned}
# ```
# with state constraint $x \in [-1, 1]^3$.
#
# The matricial form of this system is given by $\dot{x}(t) = Ax(t) + Bu(t)$ where `A` and `B` are as defined below.
# As shown in Proposition 5 of [LJ21], a set is controlled invariant for this system if and only if it is invariant for the algebraic system
# ```math
# \begin{aligned}
# \dot{x}_1(t) & = x_2(t)\\
# \dot{x}_2(t) & = x_3(t)
# \end{aligned}
# ```
# The matricial form of this system is given by $E\dot{x}(t) = Cx(t)$ where
#
# [LJ21] B. Legat and R. M. Jungers.
# *Continuous-time controlled invariant sets, a geometric approach*.
# 7th IFAC Conference on Analysis and Design of Hybrid Systems ADHS 2021, **2021**.

A = [0.0 1.0 0.0
     0.0 0.0 1.0
     0.0 0.0 0.0]
B = reshape([0.0, 0.0, 1.0], 3, 1)
E = [1.0 0.0 0.0
     0.0 1.0 0.0]
C = A[1:2, :]

# The invariance of a set $S$ for this system is characterized by the following condition (see Proposition 7 of [LJ21]):
# ```math
# \forall x \in \partial S, \exists y \in T_S(x), Ey = Cx.
# ```
# The search for a set satisfying this condition can be formulated as the following set program; see [L20] for an intoduction to set programming.
#
# [L20] Legat, B. (2020). *Set programming : theory and computation*. Ph.D. thesis, UCLouvain.

using SetProg
function maximal_invariant(family, γ = nothing; dirs=dirs)
    model = Model(sdp_solver)
    @variable(model, S, family)
    @constraint(model, S ⊆ □_3)
    x = boundary_point(S, :x)
    @constraint(model, C * x in E * tangent_cone(S, x))
    S_2 = project(S, 1:2)
    if γ === nothing
        @variable(model, γ)
    end
    for point in dirs
        @constraint(model, γ * point in S_2)
    end
    @show γ
    @objective(model, Max, γ)
    @show JuMP.objective_function(model)
    JuMP.optimize!(model)
    @show solve_time(model)
    @show JuMP.termination_status(model)
    @show JuMP.objective_value(model)
    if JuMP.termination_status(model) == MOI.OPTIMAL
        return JuMP.value(S), JuMP.objective_value(model)
    else
        return
    end
end

import GLPK
lp_solver = optimizer_with_attributes(GLPK.Optimizer, MOI.Silent() => true, "presolve" => GLPK.GLP_ON)
import CSDP
sdp_solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
using Polyhedra
interval = HalfSpace([1.0], 1.0) ∩ HalfSpace([-1.0], 1.0)
lib = Polyhedra.DefaultLibrary{Float64}(lp_solver)
□_2 = polyhedron(interval * interval, lib)
□_3 = □_2 * interval
dirs = [[-1 + √3, -1 + √3], [-1, 1]]
all_dirs = [dirs; (-).(dirs)]
inner = polyhedron(vrep(all_dirs), lib)
outer = polar(inner)

# ## Ellipsoidal template
#
# We start with the ellipsoidal template. We consider two different objectives (see Section 4.2 of [L20]):
# * the volume of the set (which corresponds to $\log(\det(Q))$ or $\sqrt[n]{\det(Q)}$ in the objective function) and
# * the sum of the squares of the length of its semi-axes of the polar (which corresponds to the trace of $Q$ in the objective function).

sol_ell, γ_ell = maximal_invariant(Ellipsoid(symmetric=true))

using Plots
function hexcolor(rgb::UInt32)
    r = ((0xff0000 & rgb) >> 16) / 255
    g = ((0x00ff00 & rgb) >>  8) / 255
    b = ((0x0000ff & rgb)      ) / 255
    Plots.RGBA(r, g, b)
end # Values taken from http://www.toutes-les-couleurs.com/code-couleur-rvb.php
lichen = hexcolor(0x85c17e)
canard = hexcolor(0x048b9a)
aurore = hexcolor(0xffcb60)
frambo = hexcolor(0xc72c48)
cols = [canard, frambo]

x2 = range(0, stop=1, length=20)
x1 = 1 .- x2.^2 / 2
upper = [[[-1, 1]]; [[x1[i], x2[i]] for i in eachindex(x2)]]
mci = polyhedron(vrep([upper; (-).(upper)]), lib)
polar_mci = polar(mci)

SetProg.Sets.print_support_function(project(sol_ell, 1:2))

# We can plot the primal solution as follows:

function primal_plot(set, γ=nothing; npoints=256, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05), args...)
    plot(ratio=:equal, tickfont=Plots.font(12); xlim=xlim, ylim=ylim, args...)
    plot!(□_2, color=lichen)
    plot!(mci, color=aurore)
    plot!(set, color=canard, npoints=npoints)
    γ === nothing || plot!(γ * inner, color=frambo)
    plot!()
end
primal_plot(project(sol_ell, 1:2), γ_ell)

# and the dual plot as follows:

function polar_plot(set, γ; npoints=256, xlim=(-1.5, 1.5), ylim=(-1.5, 1.5), args...)
    plot(ratio=:equal, tickfont=Plots.font(12); xlim=xlim, ylim=ylim, args...)
    γ === nothing || plot!(inv(γ) * outer, color=frambo)
    plot!(polar(set), color=canard, npoints=npoints)
    plot!(polar_mci, color=aurore)
    plot!(polar(□_2), color=lichen)
end
polar_plot(project(sol_ell, 1:2), γ_ell)

# ## Polyset template

p4, γ4 = maximal_invariant(PolySet(symmetric=true, degree=4, convex=true), 0.91)
γ4

# Below is the primal plot:

primal_plot(project(p4, 1:2), γ4)

# and here is the polar plot:

polar_plot(project(p4, 1:2), γ4)

# ## Piecewise semi-ellipsoidal template

sol_piece_◇, γ_piece_◇ = maximal_invariant(Ellipsoid(symmetric=true, piecewise=polar(□_3)), dirs=all_dirs)
γ_piece_◇

# Below is the primal plot:

primal_plot(project(sol_piece_◇, 1:2), γ_piece_◇)

# and here is the polar plot:

polar_plot(project(sol_piece_◇, 1:2), γ_piece_◇)


text:
module Helper

using DocStringExtensions

# Manually reexport used names in DocStringExtensions
export DOCSTRING, TYPEDEF, FIELDS, METHODLIST, SIGNATURES, @template

export Address,
    emptymeta,
    wrap_vector,
    typeunion,
    push_to_dict,
    add_to_dict,
    rev_dict,
    rev_dict_safe,
    intern,
    rand_cartesian,
    dim_max,
    dim_min

# Convenience mapping
# TODO: Think about replacing this because it could cause confusion ...
const Address = CartesianIndex

macro SetupDocStringTemplates()
    return esc(
        quote
            # Set up default template.
            @template DEFAULT = """
                $(DOCSTRING)
                """

            # Template for type definitions.
            @template TYPES = """
                $(TYPEDEF)

                Fields
                ------
                $(FIELDS)

                Documentation
                -------------
                $(DOCSTRING)

                Method List
                -----------
                $(METHODLIST)
                """
        end,
    )
end

"Return an empty `Dict{String,Any}()` for `metadata` fields."
emptymeta() = Dict{String,Any}()

wrap_vector(v) = [v]
wrap_vector(v::Vector) = identity(v)

function typeunion(A::Array)
    types = unique(map(typeof, A))
    return Array{Union{types...}}(A)
end

"""
    push_to_dict(d, k, v)

Push value `v` to the vector found in dictionary `d` at `d[k]`. If `d[k]`
does not exist, create a new vector by `d[k] = [v]`.
"""
function push_to_dict(d::AbstractDict{K,V}, k, v) where {K,V}
    haskey(d, k) ? push!(d[k], v) : d[k] = V([v])
    return nothing
end

"""
    add_to_dict(d::Dict{K}, k::K, v = 1; b = 1) where K

Increment `d[k]` by `v`. If `d[k]` does not exist, initialize `d[k] = b`.
"""
add_to_dict(d, k, v = 1; b = 1) = haskey(d, k) ? (d[k] += v) : (d[k] = b)

"""
    rev_dict(d)

Reverse the keys and values of dictionary `d`. Behavior if multiple values
are equivalent is not defined.
"""
rev_dict(d) = Dict(v => k for (k, v) in d)

"""
    rev_dict_safe(d::Dict{K,V}) where {K,V}

Reverse the keys and values of dictionary `d`. Returns a dictionary of type
`Dict{V, Vector{K}}` to handle the case where the same value in `d` points to
multiple keys.
"""
function rev_dict_safe(d::Dict{K,V}) where {K,V}
    r = Dict{V,Vector{K}}()
    for (k, v) in d
        push_to_dict(r, v, k)
    end
    return r
end

"""
    intern(x)

Given a mutable collection `x`, make all equivalent values in `x` point to a
single instance in memory. If `x` is made up of many of the same arrays, this
can greatly decrease the amount of memory required to store `x`.
"""
function intern(x)
    d = Dict{eltype(x),eltype(x)}()
    for i in eachindex(x)
        haskey(d, x[i]) ? (x[i] = d[x[i]]) : (d[x[i]] = x[i])
    end
    return nothing
end

_dimop(x, op) = reduce(op, i.I for i in x)
dim_min(x) = _dimop(x, (a, b) -> min.(a, b))
dim_max(x) = _dimop(x, (a, b) -> max.(a, b))

# Documentation for the max and min functions generated above.
@doc """
    `dim_max(indices)` returns tuple of the minimum componentwise values
    from a collection of CartesianIndices.
    """ dim_max

@doc """
    `dim_min(indices)` returns tuple of the minimum componentwise values
    from a collection of CartesianIndices.
    """ dim_min

end # module Helper


text:
#=
Working from left-to-right if no digit is exceeded by the digit to its left it is called an increasing number; for example, 134468.

Similarly if no digit is exceeded by the digit to its right it is called a decreasing number; for example, 66420.

We shall call a positive integer that is neither increasing nor decreasing a "bouncy" number; for example, 155349.

Clearly there cannot be any bouncy numbers below one-hundred, but just over half of the numbers below one-thousand (525) are bouncy. In fact, the least number for which the proportion of bouncy numbers first reaches 50% is 538.

Surprisingly, bouncy numbers become more and more common and by the time we reach 21780 the proportion of bouncy numbers is equal to 90%.

Find the least number for which the proportion of bouncy numbers is exactly 99%.
=#
function isBouncy(x)
  if length(digits(x)) < 3
    return false
  end
  arr = reverse(digits(x))
  isIncreasing = issorted(arr)
  isDecreasing = issorted(arr, rev=true)
  return !isIncreasing && !isDecreasing
end

function calc()
  bouncycount, proportion, i = 0, 0, 1
  while proportion < .99
    bouncycount = isBouncy(i) ? bouncycount + 1 : bouncycount
    proportion = (bouncycount > 0) ? bouncycount / i : 0
    i = i + 1
  end
  i - 1
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
include("activation_functions.jl")

"""
Derivative of Sigmoid
"""
function sigmoid_backwards(δA, activated_cache)
    s = sigmoid(activated_cache).A
    s = [max(min(x, 1-(eps(1.0))), eps(1.0)) for x in s]
    δZ = [δA[i] * s[i] * (1 - s[i]) for i in 1:length(s)]

    #δZ = [log(abs( 1 / x - 1 )) for x in δA]
    return δZ
end # function

"""
Derivative of ReLU
"""
function relu_backwards(δA, activated_cache)
    return [δA[i] * (activated_cache[i] > 0) for i in 1:length(δA)]
end # function

"""
Derivative of Tanh
"""
function tanh_backwards(δA, activated_cache)
    return δA[i] .* (1 .- (tanh(activated_cache[i])).^2)
end # function

"""
Derivative of Softmax
"""
function softmax_backwards(δA, activated_cache)
    return δA * ( - exp.(activated_cache) .* sum(exp.(activated_cache)) )
end # function

"""
Derivative of Swish
"""
function swish_backwards(δA, activated_cache)
    # β = 1
    return δA .* (swish(activated_cache).A .+ sigmoid(activated_cache).A .* (1 .- swish(activated_cache).A))
end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

"""
3-ary Cartesian product
"""

tripschema = @acset LabeledGraph begin
    V = 2
    E = 3
    vlabel = [:s, :s3]
    elabel = [:p1, :p2, :p3]
    src = [2,2,2]
    tgt = [1,1,1]
end

tripd = @acset LabeledGraph begin
    V = 4
    vlabel = [:s, :s, :s,]
end

trips = FLS(:trips, tripschema, [
    Cone(lconed, :s2, [1=>:p1,2=>:p2, 3=>:p3])], [])


text:
# ## A simple example
# Let's have another look at the example from the [introduction](README.md).
#
# First we create some dummy data of observations and targets:
x = rand(128, 4096)  # 4096 samples of size (128,)
y = rand(1, 4096)    # 4096 samples of size (1,)

# We then create a `DataLoader` with batch size 16:
dataloader = DataLoader((x, y), 16)

# Et voilá, a training loop:
for (xs, ys) in dataloader  # size(xs) == (128, 16)
    ## optimization step
end

# !!! note "Threading"
#
#     If you want to benefit from parallel data loading, start Julia with multiple threads
#     by setting the `-t n` flag on startup.
#
#
# `DataLoaders` supports many different data containers by building off the interface of
# [`MLDataPattern.jl`](https://mldatautilsjl.readthedocs.io/en/latest/data/pattern.html).
# Above, we pass in a tuple of datasets, hence the batch is also tuple.
#
# ---
#
# Let's [next](/docs/datacontainers.md) look at a realistic use case
# and show how to support custom data containers.


text:
#------USEARCH WRAPPERS---------

"""
    usearch_filter(inFastqPath, outFastXPath;
                   errorRate=0.01, minLength=100,
                   labelPrefix="seq", errorOut=true)

Julia wrapper for usearch error rate filter. Input file may be .fastq or .fasta.
"""
function usearch_filter(inFastqPath, outFastXPath;
                        errorRate=0.01, minLength=100,
                        labelPrefix="seq", errorOut=true)
    if last(outFastXPath) == 'q'
        outputTypeString = "-fastqout"
    else
        if last(outFastXPath)=='a'
            outputTypeString = "-fastaout"
            else error("Output filepath must end in 'q' (.fastq) or 'a' (.fasta)")
        end
    end
    errorOutString = ""
    if errorOut == true
        errorOutString = "-fastq_eeout"
    end
    run(`$(PATHS.usearch) -fastq_filter $inFastqPath -fastq_maxee_rate $errorRate -fastq_qmax 99 -fastq_minlen $minLength $outputTypeString $outFastXPath -relabel $labelPrefix $errorOutString`)
end

"""
    usearch_ref(inFastXPath, outFastaPath, referenceString;
                     divergence=0.9)

Julia wrapper for usearch.
"""
function usearch_ref(inFastXPath, outFastaPath, referenceString;
                     divergence=0.9)
    tempRefPath = "/tmp/usearch9tempRefFile"
    tempOutPath = "/tmp/usearch9tempOutFile"
    write_fasta(tempRefPath,[DNASequence(referenceString)])
    run(`$(PATHS.usearch) -usearch_global $inFastXPath -db $tempRefPath -id $divergence -fastapairs $tempOutPath -strand both`)
    names, seqs = read_fasta_with_names(tempOutPath)
    names2nd = [names[i * 2 - 1] for i in 1:Int64(round(length(names) / 2))]
    seqs2nd = [degap(seqs[i * 2 - 1]) for i in 1:Int64(round(length(names) / 2))]
    write_fasta(outFastaPath, seqs2nd, names=names2nd)
    return names2nd, seqs2nd
end

"""
    usearch_ref_database(inFastXPath, outFastaPath, referencePath;
                         divergence=0.9)

Julia wrapper for usearch.
"""
function usearch_ref_database(inFastXPath, outFastaPath, referencePath;
                              divergence=0.9)
    tempOutPath = "/tmp/usearch9tempOutFile"
    run(`$(PATHS.usearch) -usearch_global $inFastXPath -db $referencePath -id $divergence -fastapairs $tempOutPath -strand both -top_hit_only`)
    names, seqs = read_fasta_with_names(tempOutPath)
    names2nd = [names[i*2-1] for i in 1:Int64(round(length(names)/2))]
    seqs2nd = [degap(seqs[i*2-1]) for i in 1:Int64(round(length(names)/2))]
    write_fasta(outFastaPath, seqs2nd, names=names2nd)
    return names2nd, seqs2nd
end

function substring_coords(query,ref)
    start_coord = searchindex(ref,query)
    if start_coord == 0
        return false
    end
    return (start_coord,start_coord+length(query)-1)
end

function trim_fastq_entry(trimmed_seq,fastq_read,phred_arr)
    trimmed_seq = uppercase(trimmed_seq)
    fastq_read = uppercase(fastq_read)
    #first look for seq in current orientation
    revcomp = false
    try_coords = substring_coords(trimmed_seq,fastq_read)
    #then try the reverse complement
    if try_coords==false
        revcomp = true
        fastq_read = reverse_complement(fastq_read)
        try_coords = substring_coords(trimmed_seq,fastq_read)
        if try_coords==false
            return false
        else
            phred_arr = reverse(phred_arr)
        end
    end
    return fastq_read[try_coords[1]:try_coords[2]],phred_arr[try_coords[1]:try_coords[2]]
end

function array2index_dict(arr)
    dic = Dict{String,Int64}()
    for i in 1:length(arr)
        dic[arr[i]]=i
    end
    return dic
end

function usearch_trim_fastq_with_phreds(inFastqPath, referencePath; divergence=0.9, outFastqPath="")
    #READ IN THE FASTQ FILE
    fastqseqs, fastqphreds, fastqnames = read_fastq(inFastqPath)
    fastq_name_dic = array2index_dict(fastqnames)

    #USEARCH TO TRIM, BUT ONLY STORES .FASTA OUTPUT
    tempOutPath = "/tmp/usearch9tempOutFile"
    run(`$(PATHS.usearch) -usearch_global $inFastqPath -db $referencePath -id $divergence -fastapairs $tempOutPath -strand both -top_hit_only`)
    names, seqs = read_fasta_with_names(tempOutPath)
    ref_names = [names[i*2] for i in 1:Int64(round(length(names)/2))]

    names2nd = [names[i*2-1] for i in 1:Int64(round(length(names)/2))]
    seqs2nd = [degap(seqs[i*2-1]) for i in 1:Int64(round(length(names)/2))]

    #trim_fastq_entry returns read,phred
    trimmed_fastq = [(fastqnames[fastq_name_dic[names2nd[i]]],trim_fastq_entry(seqs2nd[i],fastqseqs[fastq_name_dic[names2nd[i]]],fastqphreds[fastq_name_dic[names2nd[i]]])) for i in 1:length(names2nd)]

    if outFastqPath != ""
        write_fastq(outFastqPath, [i[2][1] for i in trimmed_fastq], [i[2][2] for i in trimmed_fastq],names=[i[1] for i in trimmed_fastq])
    end

    return ref_names,trimmed_fastq
end

#-------MAFFT WRAPPER--------
"""
    mafft(inpath, outpath; path="", flags::Vector{String}=String[], kwargs...)

Julia wrapper for mafft.
"""
function mafft(inpath, outpath; path="", flags::Vector{String}=String[], kwargs...)
    mafft = length(path) == 0 ? PATHS.mafft : path
    flagstrings = []
    for flag in flags
        push!(flagstrings, "--$flag")
    end
    args = []
    for (k, v) in kwargs
        push!(args, "--$k")
        push!(args, "$v")
    end
    run(`$mafft $flagstrings $args --quiet --adjustdirection --progress /tmp/mafft.progress --out $outpath $inpath`)
end

"""
    mafft_consensus{T<:BioSequence}(seqs::Vector{T}; kwargs...)

Julia wrapper for mafft.
"""
function mafft_consensus{T<:BioSequence}(seqs::Vector{T}; kwargs...)
    mktempdir() do mydir
        seqfile = string(mydir, "/sequences.fasta")
        mafftout = string(mydir, "/mafft.fasta")

        write_fasta(string(mydir, "/sequences.fasta"), seqs)
        mafft(seqfile, mafftout; kwargs...)
        aligned = Array{String}(read_fasta(mafftout))
        return alignment_consensus(aligned)
    end
end