text:
ROS.init("testNodeHandle")
nh = ROS.NodeHandle()

@test ROS.param(nh, "float_param", 16.7) == 16.7
@test ROS.param(nh, "int_param", Int32(16)) == 16
@test ROS.param(nh, "int_param", 16) == 16
@test ROS.param(nh, "string_param", "s") == "s"
@test ROS.param(nh, "bool_param", true)
@test !ROS.hasParam(nh, "float_param")
@test !ROS.hasParam(nh, "int_param")
@test !ROS.hasParam(nh, "string_param")
@test !ROS.hasParam(nh, "bool_param")
ROS.setParam(nh, "float_param", 16.7)
ROS.setParam(nh, "int_param", Int32(16))
ROS.setParam(nh, "int_param", 16)
ROS.setParam(nh, "string_param", "s")
ROS.setParam(nh, "bool_param", true)
@test ROS.hasParam(nh, "float_param")
@test ROS.hasParam(nh, "int_param")
@test ROS.hasParam(nh, "string_param")
@test ROS.hasParam(nh, "bool_param")

ROS.deleteParam(nh, "float_param")
ROS.deleteParam(nh, "int_param")
ROS.deleteParam(nh, "string_param")
ROS.deleteParam(nh, "bool_param")

@test !ROS.hasParam(nh, "float_param")
@test !ROS.hasParam(nh, "int_param")
@test !ROS.hasParam(nh, "string_param")
@test !ROS.hasParam(nh, "bool_param")

ROS.shutdown(nh)

println("All $(basename(@__FILE__)) tests passed.")

ext: jl


text:
module HiddenMarkovModel

include("../src/Jaynes.jl")
using .Jaynes
using Distributions

function kernel(prev_latent::Float64)
    z = rand(:z, Normal(prev_latent, 1.0))
    if z > 2.0
        x = rand(:x, Normal(z, 3.0))
    else
        x = rand(:x, Normal(z, 10.0))
    end
    return z
end

# Specialized Markov call site informs tracer of dependency information.
kernelize = n -> markov(:k, kernel, n, 5.0)

simulation = () -> begin
    ps = initialize_filter(selection(), 1000, kernelize, (1, ))
    for i in 2:100
        sel = selection([(:k, i, :x) => 5.0])

        # Complexity of filter step is constant as a size of the trace.
        @time filter_step!(sel, ps, NoChange(), (i,))
    end
    return ps
end

ps = simulation()

end # module


text:
# AtmosTurbulenceStats
#
# Computes average kinetic energy and dissipation.

using ..Atmos
using ..Mesh.Topologies
using ..Mesh.Grids

struct TurbulenceStatsParams <: DiagnosticsGroupParams
    nor::Float64
    iter::Float64
end

"""
    setup_atmos_turbulence_stats(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String,
        nor::Float64,
        iter::Float64;
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create the "AtmosTurbulenceStats" `DiagnosticsGroup` which contains the
following diagnostic variables:

- E_k: volumetrically-averaged dimensionless kinetic energy
- dE: volumetrically-averaged kinetic energy dissipation

These are scalar variables computed on the DG grid (`interpol` may _not_
be specified), output on the (unlimited) `time` dimension at the specified
`interval`.
"""
function setup_atmos_turbulence_stats(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String,
    nor::Float64,
    iter::Float64;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    @assert isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosTurbulenceStats",
        Diagnostics.atmos_turbulence_stats_init,
        Diagnostics.atmos_turbulence_stats_fini,
        Diagnostics.atmos_turbulence_stats_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        TurbulenceStatsParams(nor, iter),
    )
end

# store average kinetic energy and dissipation
Base.@kwdef mutable struct AtmosTurbulenceStats{FT}
    E_k::FT = 0.0
end
const AtmosTurbulenceStatsCollected = AtmosTurbulenceStats()

function atmos_turbulence_stats_init(dgngrp, currtime)
    FT = eltype(Settings.Q)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    if mpirank == 0
        # outputs are scalars
        dims = OrderedDict()

        # set up the variables we're going to be writing
        vars = OrderedDict(
            "E_k" => ((), FT, Variables["E_k"].attrib),
            "dE" => ((), FT, Variables["dE"].attrib),
        )

        # create the output file
        dprefix = @sprintf(
            "%s_%s_%s",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        init_data(dgngrp.writer, dfilename, dims, vars)
    end

    return nothing
end

function average_kinetic_energy_and_dissipation(
    Q,
    vgeo,
    E_0,
    nor,
    iter,
    mpicomm,
)
    u₀ = Q.ρu ./ Q.ρ
    u_0 = u₀[:, 1, :] ./ nor
    v_0 = u₀[:, 2, :] ./ nor
    w_0 = u₀[:, 3, :] ./ nor

    M = vgeo[:, Grids._M, 1:size(u_0, 2)]
    SM = sum(M)

    E_k = 0.5 * sum((u_0 .^ 2 .+ v_0 .^ 2 .+ w_0 .^ 2) .* M) / SM
    E_k = MPI.Reduce(E_k, +, 0, mpicomm)

    mpirank = MPI.Comm_rank(mpicomm)
    nranks = MPI.Comm_size(mpicomm)
    if mpirank == 0
        E_k /= nranks
        dE = -(E_k - E_0) / iter
        return (E_k, dE)
    end

    return (0.0, 0.0)
end

function atmos_turbulence_stats_collect(dgngrp, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    nranks = MPI.Comm_size(mpicomm)

    E_k, dE = average_kinetic_energy_and_dissipation(
        Q,
        dg.grid.vgeo,
        AtmosTurbulenceStatsCollected.E_k,
        dgngrp.params.nor,
        dgngrp.params.iter,
        mpicomm,
    )

    if mpirank == 0
        AtmosTurbulenceStatsCollected.E_k = E_k

        varvals = OrderedDict("E_k" => E_k, "dE" => dE)

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

end

function atmos_turbulence_stats_fini(dgngrp, currtime) end

max_forks_repo_name: jm-c/CLIMA


text:
immutable AbsoluteError <: Loss
end

@inline function loss_value(loss::AbsoluteError, p::Real, y::Real)
    return abs(p - y)
end

Base.show(io::IO, loss::AbsoluteError) = print(io, "Absolute Error")


text:
using Manifolds: decorator_transparent_dispatch
using Base: decode_overlong

include("../utils.jl")
include("group_utils.jl")

@testset "General group tests" begin
    @test length(methods(has_biinvariant_metric)) == 1
    @test length(methods(has_invariant_metric)) == 1
    @test length(methods(has_biinvariant_metric)) == 1
    @testset "Not implemented operation" begin
        G = GroupManifold(NotImplementedManifold(), NotImplementedOperation())
        @test repr(G) ==
              "GroupManifold(NotImplementedManifold(), NotImplementedOperation())"
        x = [1.0, 2.0]
        v = [2.0, 3.0]
        eg = Identity(G, [0.0, 0.0])
        @test repr(eg) === "Identity($(G), $([0.0, 0.0]))"
        @test number_eltype(eg) == Bool
        @test is_point(G, eg) # identity transparent
        p = similar(x)
        copyto!(p, eg)
        @test p == eg.p
        @test isapprox(G, eg, p)
        @test isapprox(G, p, eg)
        @test isapprox(G, eg, eg)
        @test length(methods(is_group_decorator)) == 1

        @test Manifolds.is_group_decorator(G)
        @test Manifolds.decorator_group_dispatch(G) === Val{true}()
        @test Manifolds.default_decorator_dispatch(G) === Val{false}()
        @test !Manifolds.is_group_decorator(NotImplementedManifold())
        @test Manifolds.decorator_group_dispatch(NotImplementedManifold()) === Val{false}()

        @test Manifolds.decorator_transparent_dispatch(compose, G, x, x, x) ===
              Val{:intransparent}()
        @test Manifolds.decorator_transparent_dispatch(compose!, G, x, x, x) ===
              Val{:intransparent}()
        @test Manifolds.decorator_transparent_dispatch(group_exp, G, x, x) ===
              Val{:intransparent}()
        @test Manifolds.decorator_transparent_dispatch(group_log, G, x, x) ===
              Val{:intransparent}()
        @test Manifolds.decorator_transparent_dispatch(
            translate_diff!,
            G,
            x,
            x,
            x,
            x,
            x,
        ) === Val{:intransparent}()
        @test base_group(G) === G
        z = similar(x)
        copyto!(z, eg)
        @test z == eg.p
        @test NotImplementedOperation(NotImplementedManifold()) === G
        @test (NotImplementedOperation())(NotImplementedManifold()) === G

        @test_throws ErrorException allocate_result(
            G,
            get_vector,
            Identity(SpecialOrthogonal(3), x),
            v,
        )
        @test_throws ErrorException allocate_result(
            G,
            get_coordinates,
            Identity(SpecialOrthogonal(3), x),
            v,
        )
        @test_throws ErrorException allocate_result(
            ValidationManifold(NotImplementedManifold()),
            get_coordinates,
            Identity(SpecialOrthogonal(3), x),
            v,
        )
        @test_throws ErrorException base_group(
            MetricManifold(Euclidean(3), EuclideanMetric()),
        )
        @test_throws ErrorException hat(Rotations(3), eg, [1, 2, 3])
        @test_throws ErrorException hat(
            GroupManifold(Rotations(3), NotImplementedOperation()),
            eg,
            [1, 2, 3],
        )
        @test_throws ErrorException vee(Rotations(3), eg, [1, 2, 3])
        @test_throws ErrorException vee(
            GroupManifold(Rotations(3), NotImplementedOperation()),
            eg,
            [1, 2, 3],
        )
        @test_throws ErrorException Identity(Euclidean(3), [0, 0, 0])

        @test_throws ErrorException inv!(G, x, x)
        @test_throws ErrorException inv!(G, x, eg)
        @test_throws ErrorException inv(G, x)

        @test copyto!(x, eg) === x
        @test isapprox(G, x, eg)
        @test_throws ErrorException identity!(G, x, x)
        @test_throws ErrorException identity(G, x)

        @test_throws ErrorException compose(G, x, x)
        @test_throws ErrorException compose(G, x, eg)
        @test_throws ErrorException compose!(G, x, eg, x)
        @test_throws ErrorException compose!(G, x, x, eg)
        @test_throws ErrorException compose!(G, x, x, x)
        @test_throws ErrorException compose!(G, x, eg, eg)

        @test_throws ErrorException translate(G, x, x)
        @test_throws ErrorException translate(G, x, x, LeftAction())
        @test_throws ErrorException translate(G, x, x, RightAction())
        @test_throws ErrorException translate!(G, x, x, x)
        @test_throws ErrorException translate!(G, x, x, x, LeftAction())
        @test_throws ErrorException translate!(G, x, x, x, RightAction())

        @test_throws ErrorException inverse_translate(G, x, x)
        @test_throws ErrorException inverse_translate(G, x, x, LeftAction())
        @test_throws ErrorException inverse_translate(G, x, x, RightAction())
        @test_throws ErrorException inverse_translate!(G, x, x, x)
        @test_throws ErrorException inverse_translate!(G, x, x, x, LeftAction())
        @test_throws ErrorException inverse_translate!(G, x, x, x, RightAction())

        @test_throws ErrorException translate_diff(G, x, x, v)
        @test_throws ErrorException translate_diff(G, x, x, v, LeftAction())
        @test_throws ErrorException translate_diff(G, x, x, v, RightAction())
        @test_throws ErrorException translate_diff!(G, v, x, x, v)
        @test_throws ErrorException translate_diff!(G, v, x, x, v, LeftAction())
        @test_throws ErrorException translate_diff!(G, v, x, x, v, RightAction())

        @test_throws ErrorException inverse_translate_diff(G, x, x, v)
        @test_throws ErrorException inverse_translate_diff(G, x, x, v, LeftAction())
        @test_throws ErrorException inverse_translate_diff(G, x, x, v, RightAction())
        @test_throws ErrorException inverse_translate_diff!(G, v, x, x, v)
        @test_throws ErrorException inverse_translate_diff!(G, v, x, x, v, LeftAction())
        @test_throws ErrorException inverse_translate_diff!(G, v, x, x, v, RightAction())

        @test_throws ErrorException group_exp(G, v)
        @test_throws ErrorException group_exp!(G, x, v)
        @test_throws ErrorException group_log(G, x)
        @test_throws ErrorException group_log!(G, v, x)

        for f in [translate, translate!]
            @test Manifolds.decorator_transparent_dispatch(f, G) === Val{:intransparent}()
        end
        for f in [inverse_translate_diff!, inverse_translate_diff]
            @test Manifolds.decorator_transparent_dispatch(f, G) === Val{:transparent}()
        end
        for f in [group_exp!, group_exp, group_log, group_log!]
            @test Manifolds.decorator_transparent_dispatch(f, G, x, x) ===
                  Val{:intransparent}()
        end
        for f in [get_vector, get_coordinates]
            @test Manifolds.decorator_transparent_dispatch(f, G) === Val{:parent}()
        end
        @test Manifolds.decorator_transparent_dispatch(identity!, G, x, x) ===
              Val{:intransparent}()
        @test Manifolds.decorator_transparent_dispatch(isapprox, G, eg, x) ===
              Val{:transparent}()
        @test Manifolds.decorator_transparent_dispatch(isapprox, G, x, eg) ===
              Val{:transparent}()
        @test Manifolds.decorator_transparent_dispatch(isapprox, G, eg, eg) ===
              Val{:transparent}()
    end

    @testset "Action direction" begin
        @test switch_direction(LeftAction()) == RightAction()
        @test switch_direction(RightAction()) == LeftAction()

        @test Manifolds._action_order(1, 2, LeftAction()) === (1, 2)
        @test Manifolds._action_order(1, 2, RightAction()) === (2, 1)
    end

    @testset "Addition operation" begin
        G = GroupManifold(NotImplementedManifold(), Manifolds.AdditionOperation())
        test_group(G, [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], [], [[1.0, 2.0]])

        @test_throws DomainError is_point(
            G,
            Identity(
                GroupManifold(NotImplementedManifold(), NotImplementedOperation()),
                [0.0, 0.0],
            ),
            true,
        )

        x = [1.0, 2.0]
        v = [3.0, 4.0]
        ge = Identity(G, [0.0, 0.0])
        @test zero(ge) === ge
        @test number_eltype(ge) == Bool
        @test copyto!(ge, ge) === ge
        y = allocate(x)
        copyto!(y, ge)
        @test y ≈ zero(x)
        @test ge - x == -x
        @test x - ge === x
        @test ge - ge === ge
        @test ge + x ≈ x
        @test x + ge ≈ x
        @test ge + ge === ge
        @test ge + Identity(G, 1) === ge
        @test -(ge) === ge
        @test +(ge) === ge
        @test ge - Identity(G, 1) === ge
        @test ge * 1 === ge
        @test 1 * ge === ge
        @test ge * ge === ge
        @test ge.p ≈ zero(x)
        @test zero(ge) == ge
        @test inv(G, x) ≈ -x
        @test inv(G, ge) === ge
        @test identity(G, x) ≈ zero(x)
        @test identity(G, ge) === ge
        y = allocate(x)
        identity!(G, y, x)
        @test y ≈ zero(x)
        @test compose(G, x, x) ≈ x + x
        @test compose(G, x, ge) ≈ x
        @test compose(G, ge, x) ≈ x
        @test compose(G, ge, ge) == ge
        compose!(G, y, x, x)
        @test y ≈ x + x
        compose!(G, y, x, ge)
        @test y ≈ x
        compose!(G, y, ge, x)
        @test y ≈ x
        @test group_exp(G, v) === v
        @test group_log(G, x) === x

        y = identity(G, x)
        @test isapprox(y, ge; atol=1e-10)
        @test isapprox(ge, y; atol=1e-10)
        @test isapprox(ge, ge)
    end

    @testset "Multiplication operation" begin
        G = GroupManifold(NotImplementedManifold(), Manifolds.MultiplicationOperation())
        test_group(
            G,
            [[2.0 1.0; 3.0 4.0], [3.0 2.0; 4.0 5.0], [4.0 3.0; 5.0 6.0]],
            [],
            [[1.0 2.0; 3.0 4.0]];
            test_group_exp_log=true,
        )

        x = [2.0 1.0; 2.0 3.0]
        ge = Identity(G, [1.0 0.0; 0.0 1.0])
        @test number_eltype(ge) == Bool
        @test copyto!(ge, ge) === ge
        y = allocate(x)
        copyto!(y, ge)
        @test y ≈ one(x)
        @test one(ge) === ge
        @test transpose(ge) === ge
        @test det(ge) == 1
        @test ge * x ≈ x
        @test x * ge ≈ x
        @test ge * ge === ge
        @test inv(ge) === ge
        @test *(ge) === ge

        @test x / ge ≈ x
        @test ge \ x ≈ x
        @test ge / ge === ge
        @test ge \ ge === ge
        @test ge / x ≈ inv(x)
        @test x \ ge ≈ inv(x)
        y = allocate(x)
        @test LinearAlgebra.mul!(y, x, ge) === y
        @test y ≈ x
        y = allocate(x)
        @test LinearAlgebra.mul!(y, ge, x) === y
        @test y ≈ x
        y = allocate(x)
        @test LinearAlgebra.mul!(y, ge, ge) === y
        @test y ≈ one(y)

        @test ge.p ≈ one(x)
        @test inv(G, x) ≈ inv(x)
        @test inv(G, ge) === ge
        @test identity(G, x) ≈ one(x)
        @test identity(G, ge) === ge
        y = allocate(x)
        identity!(G, y, x)
        @test y ≈ one(x)
        z = allocate(x)
        copyto!(G, z, x)
        z2 = allocate(x)
        copyto!(G.manifold, z2, x)
        @test z == z2
        X = zeros(2, 2)
        Y = allocate(X)
        copyto!(G, Y, x, X)
        Y2 = allocate(X)
        copyto!(G.manifold, Y2, x, X)
        @test Y == Y2

        @test_throws ErrorException identity!(G, [0.0], ge)
        @test compose(G, x, x) ≈ x * x
        @test compose(G, x, ge) ≈ x
        @test compose(G, ge, x) ≈ x
        @test compose(G, ge, ge) == ge
        compose!(G, y, x, x)
        @test y ≈ x * x
        compose!(G, y, x, ge)
        @test y ≈ x
        compose!(G, y, ge, x)
        @test y ≈ x
        X = [1.0 2.0; 3.0 4.0]
        @test group_exp!(G, y, X) === y
        @test_throws ErrorException group_exp!(G, y, :a)
        @test y ≈ exp(X)
        Y = allocate(X)
        @test group_log!(G, Y, y) === Y
        @test Y ≈ log(y)

        @testset "identity optimization" begin
            x2 = copy(x)
            identity!(G, x2, x)
            x3 = copy(x)
            invoke(
                identity!,
                Tuple{
                    AbstractGroupManifold{ℝ,Manifolds.MultiplicationOperation},
                    Any,
                    AbstractMatrix,
                },
                G,
                x3,
                x,
            )
            @test isapprox(G, x2, x3)
        end
    end

    @testset "Identity on Group Manifolds" begin
        G = TranslationGroup(3)
        e = Identity(G, zeros(3))
        @test get_vector(G, e, ones(3), DefaultOrthogonalBasis()) == ones(3)
        @test e - e == e
        @test ones(3) + e == ones(3)
    end

    @testset "Transparency tests" begin
        G = DefaultTransparencyGroup(Euclidean(3), AdditionOperation())
        p = ones(3)
        q = 2 * p
        X = zeros(3)
        Y = similar(X)
        for f in
            [vector_transport_along!, vector_transport_direction!, vector_transport_to!]
            @test ManifoldsBase.decorator_transparent_dispatch(
                f,
                G,
                Y,
                p,
                X,
                q,
                ParallelTransport(),
            ) == Val(:intransparent)
        end
    end
end

struct NotImplementedAction <: AbstractGroupAction{LeftAction} end

@testset "General group action tests" begin
    @testset "Not implemented operations" begin
        A = NotImplementedAction()
        x = [1.0, 2.0]
        a = [1.0, 2.0]
        v = [1.0, 2.0]

        @test_throws ErrorException base_group(A)
        @test_throws ErrorException g_manifold(A)
        @test_throws ErrorException apply(A, a, x)
        @test_throws ErrorException apply!(A, x, a, x)
        @test_throws ErrorException inverse_apply(A, a, x)
        @test_throws ErrorException inverse_apply!(A, x, a, x)
        @test_throws ErrorException apply_diff(A, a, x, v)
        @test_throws ErrorException apply_diff!(A, v, x, a, v)
        @test_throws ErrorException inverse_apply_diff(A, a, x, v)
        @test_throws ErrorException inverse_apply_diff!(A, v, x, a, v)
        @test_throws ErrorException compose(A, a, a)
        @test_throws ErrorException compose!(A, a, a, a)
        @test_throws ErrorException optimal_alignment(A, x, x)
        @test_throws ErrorException optimal_alignment!(A, a, x, x)
    end
end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, LinearAlgebra
using Random, RandomNumbers
seed=100; Random.seed!(seed)
@testset "Two noise processes for different m" begin
  for m=1:4
    rand_prototype = zeros(m)
    rand_prototype2 = similar(rand_prototype,Int(m*(m-1)/2))
    rand_prototype2 .= false
    Random.seed!(seed)
    W = WienerProcess!(0.0,rand_prototype,rand_prototype2,
                     save_everystep=true,
                     rng = Xorshifts.Xoroshiro128Plus(seed))
    prob = NoiseProblem(W,(0.0,1.0))
    sol = solve(prob;dt=0.2)

    Random.seed!(seed)
    Woop = WienerProcess(0.0,rand_prototype,rand_prototype2,
                     save_everystep=true,
                     rng = Xorshifts.Xoroshiro128Plus(seed))
    proboop = NoiseProblem(Woop,(0.0,1.0))
    soloop = solve(proboop;dt=0.2)

    @test norm(soloop.W - sol.W) == 0.0
    @test norm(soloop.Z - sol.Z) == 0.0

  end
end


text:
using Plots
using WriteVTK
gr()

include( "SemsFDTD.jl" )

testsim = SemsFDTD.new_sim( Float32, 1000, 1000, 25 )

SemsFDTD.timestep( testsim )

@time begin
    for i = 1:50
        SemsFDTD.timestep( testsim )
        print( ' ', i )
    end
end

heatmap( testsim.coeff_e_xe[:,:,2], aspect_ratio = 1 )
# vtk_write_array( "test_file_out", rand( 100, 100, 100 )

max_stars_repo_name: asandwhich/SemsFDTD

max_issues_repo_name: asandwhich/SemsFDTD

max_forks_repo_name: asandwhich/SemsFDTD


text:
"""
"""
struct RestrictedTriangulation{Dc,Dp,G,A} <: Triangulation{Dc,Dp}
  parent_trian::G
  cell_to_parent_cell::A
  @doc """
  """
  function RestrictedTriangulation(
    parent_trian::Triangulation,
    cell_to_parent_cell::AbstractVector{<:Integer})

    Dc = num_cell_dims(parent_trian)
    Dp = num_point_dims(parent_trian)
    G = typeof(parent_trian)
    A = typeof(cell_to_parent_cell)
    new{Dc,Dp,G,A}(parent_trian,cell_to_parent_cell)
  end
end

# Constructors

function RestrictedTriangulation(
  parent_trian::Triangulation,
  parent_cell_to_mask::AbstractArray{Bool})

  cell_to_parent_cell = findall(collect1d(parent_cell_to_mask))
  RestrictedTriangulation(parent_trian,cell_to_parent_cell)
end

function RestrictedTriangulation(
  parent_trian::Triangulation,
  parent_cell_to_mask::AbstractVector{Bool})

  cell_to_parent_cell = findall(parent_cell_to_mask)
  RestrictedTriangulation(parent_trian,cell_to_parent_cell)
end

function Triangulation(
  parent_trian::Triangulation,
  cell_to_parent_cell::AbstractVector{<:Integer})

  RestrictedTriangulation(parent_trian,cell_to_parent_cell)
end

function Triangulation(
  parent_trian::Triangulation,
  parent_cell_to_mask::AbstractArray{Bool})

  RestrictedTriangulation(parent_trian,parent_cell_to_mask)
end

function Triangulation(
  parent_model::DiscreteModel,
  cell_to_parent_cell::AbstractVector{<:Integer})

  parent_trian = Triangulation(parent_model)
  RestrictedTriangulation(parent_trian,cell_to_parent_cell)
end

function Triangulation(
  parent_model::DiscreteModel,
  parent_cell_to_mask::AbstractArray{Bool})

  parent_trian = Triangulation(parent_model)
  RestrictedTriangulation(parent_trian,parent_cell_to_mask)
end

function Triangulation(parent_model::DiscreteModel, labels::FaceLabeling; tags)
  parent_trian = Triangulation(parent_model)
  parent_cell_to_mask = get_face_mask(labels,tags,num_cell_dims(parent_model))
  RestrictedTriangulation(parent_trian,parent_cell_to_mask)
end

# Triangulation API

TriangulationStyle(::Type{<:RestrictedTriangulation}) = SubTriangulation()

get_background_triangulation(trian::RestrictedTriangulation) = get_background_triangulation(trian.parent_trian)

get_node_coordinates(trian::RestrictedTriangulation) = get_node_coordinates(trian.parent_trian)

get_reffes(trian::RestrictedTriangulation) = get_reffes(trian.parent_trian)

function get_cell_coordinates(trian::RestrictedTriangulation)
  parent_cell_data = get_cell_coordinates(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end

function get_cell_type(trian::RestrictedTriangulation)
  parent_cell_data = get_cell_type(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end

function get_cell_node_ids(trian::RestrictedTriangulation)
  parent_cell_data = get_cell_node_ids(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end

function get_cell_map(trian::RestrictedTriangulation)
  parent_cell_data = get_cell_map(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end

function get_cell_to_bgcell(trian::RestrictedTriangulation)
  parent_cell_data = get_cell_to_bgcell(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end

function get_facet_normal(trian::RestrictedTriangulation)
  parent_cell_data = get_facet_normal(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end

function get_cell_ref_map(trian::RestrictedTriangulation)
  parent_cell_data = get_cell_ref_map(trian.parent_trian)
  lazy_map(Reindex(parent_cell_data),trian.cell_to_parent_cell)
end


text:
using StatisticalTests
using Test, RDatasets, Statistics, StatsModels, Distributions

test_show(x) = show(IOBuffer(), x)

include("test_t.jl")
include("test_f.jl")
include("test_ks.jl")
include("test_wilcox.jl")


text:
using NBody
using ArgParse

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--position_only", "-p"
      help = "Position only"
      action = :store_true
  end

  parse_args(s)
end

function read_nb()
  str = readline(STDIN)
  NBody.read_nbody_json(str)
end

function main(args)
  nb_a = read_nb()
  nb_b = read_nb()

  nb_c = NBody.diff(nb_a, nb_b)

  distance = if args["position_only"]
    NBody.abs_pos(nb_c)
  else
    NBody.abs(nb_c)
  end

  @printf("Phase Space Distance: %.16e\n", distance)
end

main(parse_commandline())


text:
using Nettle

labels = ["\"\" (empty string)", "\"a\"", "\"abc\"",
        "\"message digest\"", "\"a...z\"",
        "\"abcdbcde...nopq\"", "\"A...Za...z0...9\"",
        "8 times \"1234567890\"", "1 million times \"a\""]
texts = ["", "a", "abc", "message digest", "abcdefghijklmnopqrstuvwxyz",
        "abcdbcdecdefdefgefghfghighijhijkijkljklmklmnlmnomnopnopq",
        "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789",
        "1234567890" ^ 8, "a" ^ 1_000_000]
expects = ["9c1185a5c5e9fc54612808977ee8f548b2258d31",
        "0bdc9d2d256b3ee9daae347be6f4dc835a467ffe",
        "8eb208f7e05d987a9b044a8e98c6b087f15a0bfc",
        "5d0689ef49d2fae572b881b123a85ffa21595f36",
        "f71c27109c692c1b56bbdceb5b9d2865b3708dbc",
        "12a053384a9c0c88e405a06c27dcf49ada62eb2b",
        "b0e20b6e3116640286ed3a87a5713079b21f5189",
        "9b752e45573d4b39f4dbd3323cab82bf63326bfb",
        "52783243c1697bdbe16d37f97f68f08325dc1528"]

for (lab, text, expect) in zip(labels, texts, expects)
    digest = hexdigest("ripemd160", text)
    println("# $lab\n -> digest: $digest\n -> expect: $expect")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Hs::Vector{MPO}
psi0::MPS
Hs = permute.(Hs, Ref((linkind, siteinds, linkind)))
PH = ExtendedProjMPOSum(Hs)
orthogonalize!(psi0, 1)
position!(PH, psi0, 1)
dmrg(PH, psi, sweeps; kwargs...)  # Keeping track of PHK allows one to restart without extra cost
"""
module ExtendedProjMPOs

using ITensors

import ITensors: product, position!, noiseterm, lproj, rproj

"""
Edge tensor with added on MPO tensor; srw 4/30/21
o--o--o-
|  |  |  |
o--o--o--o-
|  |  |  |
o--o--o-
"""
mutable struct ExtendedProjMPO
  projmpo::ProjMPO
  EL::ITensor
  ER::ITensor
  ExtendedProjMPO(H::MPO) = new(ProjMPO(H), ITensor(), ITensor())
end

nsite(p::ExtendedProjMPO) = ITensors.nsite(p.projmpo)
Base.length(p::ExtendedProjMPO) = length(p.projmpo)
(p::ExtendedProjMPO)(v::ITensor) = product(p.projmpo, v)
Base.eltype(p::ExtendedProjMPO) = ITensors.eltype(p.projmpo)
Base.size(p::ExtendedProjMPO) = ITensors.size(p.projmpo)

lproj(p::ExtendedProjMPO) = ITensors.lproj(p.projmpo)
rproj(p::ExtendedProjMPO) = ITensors.rproj(p.projmpo)

function position!(p::ExtendedProjMPO, psi::MPS, pos::Int)
  ITensors.position!(p.projmpo, psi, pos)
  L, R = lproj(p), rproj(p)
  isnothing(L) && (L = 1.0)
  isnothing(R) && (R = 1.0)
  r = p.projmpo.rpos
  p.EL = L * p.projmpo.H[pos]
  p.ER = R * p.projmpo.H[r - 1]
  return p
end

function product(p::ExtendedProjMPO, v::ITensor)::ITensor
  Hv = p.EL * v * p.ER
  if order(Hv) != order(v)
    error(
      string(
        "The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, ",
        "this is probably due to an index mismatch.\nCommon reasons for this error: \n",
        "(1) You are trying to multiply the ProjMPO with the $(nsite(p))-site wave-function at the wrong position.\n",
        "(2) `orthognalize!` was called, changing the MPS without updating the ProjMPO.\n\n",
        "P*v inds: $(inds(Hv)) \n\n",
        "v inds: $(inds(v))",
      ),
    )
  end
  return noprime(Hv)
end

function noiseterm(p::ExtendedProjMPO, phi::ITensor, ortho::String) # ITensors.noiseterm(p.projmpo,phi,ortho)
  if nsite(p) != 2
    error("noise term only defined for 2-site ProjMPO")
  end
  LR = (ortho == "left" ? p.EL : p.ER)
  U, S, V = svd(phi, commoninds(LR, phi)...; cutoff=1e-4)
  nt = LR * U * S
  return nt * dag(noprime(nt))
end

mutable struct ExtendedProjMPOSum
  pm::Vector{ExtendedProjMPO}
end

function ExtendedProjMPOSum(mpos::Vector{MPO})
  return ExtendedProjMPOSum([ExtendedProjMPO(M) for M in mpos])
end

ExtendedProjMPOSum(Ms::MPO...) = ExtendedProjMPOSum([Ms...])

nsite(P::ExtendedProjMPOSum) = nsite(P.pm[1])

Base.length(P::ExtendedProjMPOSum) = length(P.pm[1])

function product(P::ExtendedProjMPOSum, v::ITensor)::ITensor
  Pvs = fill(ITensor(), Threads.nthreads())
  Threads.@threads for Ppmn in P.pm
    Pvs[Threads.threadid()] += product(Ppmn, v)
  end
  return sum(Pvs)
end

function Base.eltype(P::ExtendedProjMPOSum)
  elT = eltype(P.pm[1])
  for n in 2:length(P.pm)
    elT = promote_type(elT, eltype(P.pm[n]))
  end
  return elT
end

(P::ExtendedProjMPOSum)(v::ITensor) = product(P, v)

Base.size(P::ExtendedProjMPOSum) = size(P.pm[1])

function position!(P::ExtendedProjMPOSum, psi::MPS, pos::Int)
  Threads.@threads for M in P.pm
    position!(M, psi, pos)
  end
  return P
end

function noiseterm(P::ExtendedProjMPOSum, phi::ITensor, dir::String)
  nts = fill(ITensor(), Threads.nthreads())
  Threads.@threads for pp in P.pm
    nts[Threads.threadid()] += noiseterm(pp, phi, dir)
  end
  return sum(nts)
end

export ExtendedProjMPO, ExtendedProjMPOSum

end


text:
"""
A forge is an online platform for Git repositories.
The most common example is [GitHub](https://github.com).

`Forge` subtypes can access their respective web APIs.
"""
abstract type Forge end

"""
    Endpoint(
        method::Symbol,
        url::$AStr;
        headers::Vector{<:Pair}=HTTP.Header[],
        query::Dict=Dict(),
        allow_404::Bool=false,
    ) -> Endpoint

Contains information on how to call an endpoint.

## Arguments
- `method::Symbol`: HTTP request method to use.
- `url::$AStr`: Endpoint URL, relative to the base URL.

## Keywords
- `headers::Vector{<:Pair}=HTTP.Header[]`: Request headers to add.
- `query::Dict=Dict()`: Query string parameters to add.
- `allow_404::Bool=false`: Sends responses  with 404 statuses to the postprocessor.
"""
struct Endpoint
    method::Symbol
    url::String
    headers::Vector{<:Pair}
    query::Dict
    allow_404::Bool

    function Endpoint(
        method::Symbol, url::AStr;
        headers::Vector{<:Pair}=HTTP.Header[],
        query::Dict=Dict(),
        allow_404::Bool=false,
    )
        return new(method, url, headers, query, allow_404)
    end
end

"""
    base_url(::Forge) -> String

Returns the base URL of all API endpoints.
"""
base_url(::Forge) = ""

"""
    request_headers(::Forge, ::Function) -> Vector{Pair}

Returns the headers that should be added to each request.
"""
request_headers(::Forge, ::Function) = HEADERS

"""
    request_query(::Forge, ::Function) -> Dict

Returns the query string parameters that should be added to each request.
"""
request_query(::Forge, ::Function) = Dict()

"""
    request_kwargs(::Forge, ::Function) -> Dict{Symbol}

Returns the extra keyword arguments that should be passed to `HTTP.request`.
"""
request_kwargs(::Forge, ::Function) = Dict()

"""
    postprocessor(::Forge, ::Function) -> PostProcessor

Returns the [`PostProcessor`](@ref) to be used.
"""
postprocessor(::Forge, ::Function) = DoNothing()

"""
    into(::Forge, ::Function) -> Type

Returns the type that the [`PostProcessor`](@ref) should create from the response.

"""
into(::Forge, ::Function) = Nothing

"""
    endpoint(::Forge, ::Function, args...) -> Endpoint

Returns an [`Endpoint`](@ref) for a given function.
Trailing arguments are usually important for routing.
For example, [`get_user`](@ref) can take some ID parameter which becomes part of the URL.
"""
endpoint(f::T, ::Function, args...) where T <: Forge =
    error("$T has not implimented this function")

"""
    has_rate_limits(::Forge, ::Function) -> Bool

Returns whether or not the forge server uses rate limiting.
"""
has_rate_limits(::Forge, ::Function) = false

"""
    rate_limit_check(::Forge, ::Function) -> Bool

Returns whether or not there is an active rate limit.
If one is found, [`on_rate_limit`](@ref) is called to determine how to react.
"""
rate_limit_check(::Forge, ::Function) = false

"""
    on_rate_limit(::Forge, ::Function) -> OnRateLimit

Returns an [`OnRateLimit`](@ref) that determines how to react to an exceeded rate limit.
"""
on_rate_limit(::Forge, ::Function) = ORL_RETURN

"""
    rate_limit_wait(::Forge, ::Function)

Wait for a rate limit to expire.
"""
rate_limit_wait(::Forge, ::Function) = nothing

"""
    rate_limit_period(::Forge, ::Function) -> Period

Compute the amount of time until a rate limit expires.
"""
rate_limit_period(::Forge, ::Function) = nothing

"""
    rate_limit_update!(::Forge, ::Function, ::HTTP.Response)

Update the rate limiter with a new response.
"""
rate_limit_update!(::Forge, ::Function, ::HTTP.Response) = nothing


text:
module JSON

using HTTP
using JSON2

import HTTP: handle

export JSONHandler

"""
    JSONHandler{T,F} <: HTTP.RequestHandler

HTTP.RequestHandler specialization holding a function of type `F`. The wrapped
function receives the payload of type `T`, which is parsed from the HTTP request
body in JSON format. The function's return value will be converted to JSON and
written to the HTTP response body
"""
struct JSONHandler{T,F} <: HTTP.RequestHandler
    func::F # func(reqPayload::T)
end

"""
    JSONHandler{T}(func::F) where {T, F}

Outer constructor for JSONHandler, inferring `T` and `F` from parameters.
"""
JSONHandler(::Type{T}, func::F) where {T, F} = JSONHandler{T, F}(func)

"""
    HTTP.handle(h::JSONHandler{T}, ::HTTP.Request) where {T}

Overloaded method of `HTTP.handle` for `JSONHandler`. Parses the request body
into an instance of `T` and passes it to `h`. The return value of `h` is
converted to JSON and written to the response body.

Registering a `JSONHandler` with `HTTP.@register` will make the `HTTP` module
call this method for a matching request.
"""
function handle(h::JSONHandler{T}, req::HTTP.Request) where {T}
    res = req.response

    local reqPayload
    try
        reqPayload = JSON2.read(IOBuffer(HTTP.payload(req)), T)
    catch e
        # JSON parsing failed, return "bad request"
        res.status = 400
        return res
    end

    resPayload = h.func(reqPayload)
    HTTP.setheader(res, "Content-Type" => "application/json")
    JSON2.write(IOBuffer(res.body, write = true), resPayload)
    return res
end

end


text:
using Test
using ExamplePkg

@testset "square" begin
    @test ExamplePkg.square(4.0) == 16.0
    @test ExamplePkg.square(0.0) == 0.0
    @test ExamplePkg.square(-4.0) == 16.0
    @test ExamplePkg.square(1e-300) == 0.0
end;


text:
using .Zygote

Zygote.@adjoint function (ev::Eval)(args...)
    Z = ev.fwd(args...)
    Z, Δ -> begin
        isnothing(ev.rev) && error("no gradient definition here!")
        tuple(nothing, ev.rev(Δ, Z, args...)...)
    end
end

Tullio.promote_storage(T::Type, ::Type{<:Zygote.Fill}) = T
Tullio.promote_storage(::Type{<:Zygote.Fill}, T::Type) = T


text:
"""
    flux_variability_analysis(
        model::MetabolicModel,
        reactions::Vector{Int},
        optimizer;
        modifications = [],
        workers = [myid()],
        bounds = z -> (z,z),
        ret = objective_value,
    )::Matrix{Float64}

Flux variability analysis solves a pair of optimization problems in `model` for
each flux listed in `reactions`:
```
 min,max xᵢ
s.t. S x = b
    xₗ ≤ x ≤ xᵤ
     cᵀx ≥ bounds(Z₀)[1]
     cᵀx ≤ bounds(Z₀)[2]
```
where Z₀:= cᵀx₀ is the objective value of an optimal solution of the associated
FBA problem (see [`flux_balance_analysis`](@ref)). See "Gudmundsson, S., Thiele,
I. Computationally efficient flux variability analysis. BMC Bioinformatics 11,
489 (2010). https://doi.org/10.1186/1471-2105-11-489" for more information.

The `bounds` is a user-supplied function that specifies the objective bounds
for the variability optimizations, by default it restricts the flux objective
value to the precise optimum reached in FBA. It can return `-Inf` and `Inf` in
first and second pair to remove the limit. Use [`gamma_bounds`](@ref) and
[`objective_bounds`](@ref) for simple bounds.

`optimizer` must be set to a `JuMP`-compatible optimizer. The computation of
the individual optimization problems is transparently distributed to `workers`
(see `Distributed.workers()`).

`ret` is a function used to extract results from optimized JuMP models of the
individual reactions. By default, it calls and returns the value of
`JuMP.objective_value`. More information can be extracted e.g. by setting it to
a function that returns a more elaborate data structure; such as `m ->
(JuMP.objective_value(m), JuMP.value.(m[:x]))`.

Returns a matrix of extracted `ret` values for minima and maxima, of total size
(`length(reactions)`,2). The optimizer result status is checked with
[`is_solved`](@ref); `nothing` is returned if the optimization failed for any
reason.
"""
function flux_variability_analysis(
    model::MetabolicModel,
    reactions::Vector{Int},
    optimizer;
    modifications = [],
    workers = [myid()],
    bounds = z -> (z, z),
    ret = objective_value,
)
    if any(reactions .< 1) || any(reactions .> n_reactions(model))
        throw(DomainError(reactions, "Index exceeds number of reactions."))
    end

    Z = bounds(
        objective_value(
            flux_balance_analysis(model, optimizer; modifications = modifications),
        ),
    )

    # store a JuMP optimization model at all workers
    save_model = :(
        begin
            model = $model
            opt_model = $COBREXA.make_optimization_model(model, $optimizer)
            for mod in $modifications
                mod(model, opt_model)
            end
            $COBREXA._FVA_add_constraint(opt_model, $(objective(model)), opt_model[:x], $Z)
            opt_model
        end
    )
    map(fetch, save_at.(workers, :cobrexa_parfva_model, Ref(save_model)))
    save_model = nothing # this has some volume, free it again

    # schedule FVA parts parallely using pmap
    fluxes = dpmap(
        rid -> :($COBREXA._FVA_optimize_reaction(cobrexa_parfva_model, $rid, $ret)),
        CachingPool(workers),
        [-reactions reactions],
    )

    # free the data on workers
    map(fetch, remove_from.(workers, :cobrexa_parfva_model))

    return fluxes
end

"""
    flux_variability_analysis(
        model::MetabolicModel,
        optimizer;
        kwargs...
    )

A simpler version of [`flux_variability_analysis`](@ref) that maximizes and minimizes all reactions in the model. Arguments are forwarded.
"""
function flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...)
    n = n_reactions(model)
    return flux_variability_analysis(model, collect(1:n), optimizer; kwargs...)
end

"""
    flux_variability_analysis_dict(
        model::MetabolicModel,
        optimizer;
        kwargs...
    )

A variant of [`flux_variability_analysis`](@ref) that returns the individual
maximized and minimized fluxes of all reactions as two dictionaries (of
dictionaries). All keyword arguments except `ret` are passed through.
"""
function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
    vs = flux_variability_analysis(
        model,
        optimizer;
        kwargs...,
        ret = m -> JuMP.value.(m[:x]),
    )
    rxns = reactions(model)

    return (
        Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 1]])),
        Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 2]])),
    )
end

"""
    _FVA_add_constraint(model, c, x, Z)

Internal helper function for adding constraints to a model. Exists mainly
because for avoiding namespace problems on remote workers.
"""
function _FVA_add_constraint(model, c, x, Z)
    Z[1] > -Inf && @constraint(model, c' * x >= Z[1])
    Z[2] < Inf && @constraint(model, c' * x <= Z[2])
end

"""
    _FVA_get_opt(model, rid)

Internal helper for creating the optimized model on a remote worker, for
avoiding namespace problems.
"""
function _FVA_optimize_reaction(model, rid, ret)
    sense = rid > 0 ? MAX_SENSE : MIN_SENSE
    var = all_variables(model)[abs(rid)]

    @objective(model, sense, var)
    optimize!(model)

    if is_solved(model)
        return ret(model)
    else
        return nothing
    end
end


text:
"""
a direct radixsort on strings
    10m - 3-4 seconds slightly slower than radixsort! which only takes 2.5
    100m - 75 seconds
"""
function direct_radixsort_string!(strvec::AbstractVector{S}, sizeofstr = sizeof.(strvec), sim_strvec = similar(strvec)) where {S <: String}
    l = length(strvec)

    # Init
    iters = maximum(sizeofstr)
    bin = zeros(UInt32, 256, iters)

    # Histogram for each element, radix
    @time for i = 1:l
         @inbounds sz = sizeofstr[i]
        for j = 1:sz
            @inbounds idx = Int(codeunit(strvec[i], j))+1
            @inbounds bin[idx,j] += 1
        end
        for j = sz+1:iters
            @inbounds bin[1,j] += 1
        end
    end

    # Sort!
    swaps = 0
    for j = iters:-1:1
        # Unroll first data iteration, check for degenerate case
        @inbounds idx = (sizeofstr[l] >= j) ? Int(codeunit(strvec[l], j))+1 : 1

        # are all values the same at this radix?
        if bin[idx,j] == l;  continue;  end

        cbin = cumsum(bin[:,j])
        ci = cbin[idx]
        sim_strvec[ci] = strvec[l]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in l-1:-1:1
            @inbounds idx = (sizeofstr[i] >= j) ? Int(codeunit(strvec[i], j))+1 : 1

            ci = cbin[idx]
            # println(ci)
            sim_strvec[ci] = strvec[i]
            # println("hello")
            cbin[idx] -= 1
            # println("hello2")
        end
        sim_strvec, strvec = strvec, sim_strvec
        # println("hello3")

        # try
        sizeofstr = sizeof.(strvec)
        swaps += 1
    end

    if isodd(swaps)
        sim_strvec,strvec = strvec,sim_strvec
        for i = 1:l
            @inbounds strvec[i] = sim_strvec[i]
        end
    end
    strvec
end

# const M=1000; const K=100
# srand(1)
# @time svec1 = rand([string(rand(Char.(1:255), rand(1:8))...) for k in 1:M÷K], M)
# @time radixsort_string!(svec1)
# issorted(svec1)
#
# # test equal length
# const M=10_000_000; const K=100
# srand(1)
# @time svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
# # @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
# @time radixsort_string!(svec1) # 3.9 seconds for
# issorted(svec1)
#
# const M=100_000_000; const K=100
# srand(1)
# @time svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
# # @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
# @time radixsort_string!(svec1) # 7 seconds for
# issorted(svec1)


text:
#= ######################################################
model structs
###################################################### =#

mutable struct ElementRaw{T1<:Integer, T2<:AbstractFloat}
    node1::T1
    node2::T1
    com::Point2D{T2}
    nvec::Point2D{T2}
    area::T2
end

struct Element{T1<:Integer, T2<:AbstractFloat}
    node1::T1
    node2::T1
    com::Point2D{T2}
    nvec::Point2D{T2}
    area::T2
end

struct Part{T1<:Integer, T2<:AbstractFloat}
    nodes::Vector{Point2D{T2}}
    elem::Vector{ElementRaw{T1,T2}}
    nnodes::T1
    nelem::T1
end

struct ElementAssign{T<:Integer}
    first::T
    last::T
end

mutable struct ModelRaw{T1<:Integer, T2<:AbstractFloat}
    nodes::Vector{Point2D{T2}}
    elem::Vector{ElementRaw{T1,T2}}
    elem2par::Vector{ElementAssign{T1}}
    nnodes::T1
    nelem::T1
    npar::T1
end

struct Model{T1<:Integer, T2<:AbstractFloat}
    nodes::Vector{Point2D{T2}}
    elem::Vector{Element{T1,T2}}
    elem2par::Vector{ElementAssign{T1}}
    nnodes::T1
    nelem::T1
    npar::T1
end

#= ######################################################
Functions for working with structs
###################################################### =#

function create_empty_model()
    return ModelRaw(Vector{Point2D{Float64}}(), Vector{ElementRaw{Int64,Float64}}(), Vector{ElementAssign{Int64}}(), 0, 0, 0)
end

function norm(v::Point2D)
    l = sqrt(v.x^2 + v.y^2)
    return l
end

function normit(v::Point2D)::Point2D
    vn = v / norm(v)
    return vn
end

function add_offset!(p::Part, offset::Integer)
    for i = 1:p.nelem
        p.elem[i].node1 += offset
        p.elem[i].node2 += offset
    end
end

function get_com(p1::Point2D{T}, p2::Point2D{T})::Point2D{T} where T<:AbstractFloat
    # calculate center of line (center of mass)
    c = p1 + ((p2 - p1) / 2)
    return c
end

function get_nvec(p1::Point2D{T}, p2::Point2D{T}, dir::String)::Point2D{T} where T<:AbstractFloat
    # calculate normal vector
    n = Point2D((-1) * (p2.y - p1.y), (p2.x - p1.x))
    (dir == "neg") && (n = n * (-1.0))
    n_norm = normit(n)
    return n_norm
end

function get_length(p1::Point2D, p2::Point2D)
    # calculate length of line element
    # is used as area here
    l = sqrt((p1.x-p2.x)^2 + (p1.y-p2.y)^2)
    return l
end

function add!(m::ModelRaw, p::Part)
    append!(m.nodes, p.nodes)
    add_offset!(p, m.nnodes)
    append!(m.elem, p.elem)
    push!(m.elem2par, ElementAssign(m.nelem+1, m.nelem+p.nelem))
    m.nnodes += p.nnodes
    m.nelem += p.nelem
    m.npar += 1
end

function get_nodes_min_and_max(m)
    # get minimum and maximum value of nodes
    # as Vector of Point2D
    xmin = minimum(m.nodes[i].x for i = 1:m.nnodes)
    xmax = maximum(m.nodes[i].x for i = 1:m.nnodes)
    ymin = minimum(m.nodes[i].y for i = 1:m.nnodes)
    ymax = maximum(m.nodes[i].y for i = 1:m.nnodes)
    return Point2D(xmin, ymin), Point2D(xmax, ymax)
end

function offset_model!(m::ModelRaw)
    # offset all nodes to only positiv coords
    nmin, nmax = get_nodes_min_and_max(m)
    for i = 1:m.nnodes
        m.nodes[i] = m.nodes[i] - nmin
    end
    for i = 1:m.nelem
        m.elem[i].com = m.elem[i].com - nmin
    end
end

function make_model_immutable(m::ModelRaw)::Model
    nodes = m.nodes
    elem = Vector{Element{Int64,Float64}}(undef,m.nelem)
    for i = 1:m.nelem
        i1 = m.elem[i].node1
        i2 = m.elem[i].node2
        c = m.elem[i].com
        nv = m.elem[i].nvec
        a = m.elem[i].area
        elem[i] = Element(i1, i2, c, nv, a)
    end
    elem2par = m.elem2par
    nnodes = m.nnodes
    nelem = m.nelem
    npar = m.npar
    return Model(nodes, elem, elem2par, nnodes, nelem, npar)
end

#= ######################################################
Define discretized geometries
###################################################### =#

function edge(p1::Point2D, p2::Point2D; seed::T = 10, dir::String = "pos") where T<:Integer
    """
    creating edge between two points
    # println("edge between ", p1, " and ",p2)
    """
    dx = (p2.x - p1.x) / seed
    dy = (p2.y - p1.y) / seed
    nnodes = seed+1
    nodes = Vector{Point2D{Float64}}(undef,nnodes)
    for i = 1:nnodes
        nodes[i] = Point2D(p1.x + (i-1) * dx, p1.y + (i-1) * dy)
    end
    # nodes = [Point2D(p1.x + (i-1) * dx, p1.y + (i-1) * dy) for i = 1:nnodes]
    nelements = seed
    elements = Vector{ElementRaw{Int64,Float64}}(undef,nelements)
    for i = 1:nelements
        i1 = i
        i2 = i + 1
        c = get_com(nodes[i1], nodes[i2])
        nv = get_nvec(nodes[i1], nodes[i2], dir)
        a = get_length(nodes[i1], nodes[i2])
        elements[i] = ElementRaw(i1, i2, c, nv, a)
    end
    return Part(nodes, elements, nnodes, nelements)
end

function rectangle(x::T1, y::T1, c::Point2D; seedx::T2 = 10, seedy::T2 = 10, dir::String = "pos") where {T1<:Real, T2<:Integer}
    # creating rectangle with center
    p1 = c + Point2D(-0.5*x, -0.5*y)
    p2 = c + Point2D(-0.5*x, 0.5*y)
    p3 = c + Point2D(0.5*x, 0.5*y)
    p4 = c + Point2D(0.5*x, -0.5*y)
    dx = (p3.x - p2.x) / seedx
    dy = (p2.y - p1.y) / seedy
    nnodes = 2*seedx + 2*seedy
    nodes = Vector{Point2D{Float64}}(undef,nnodes)
    for i = 1:seedy
        nodes[i] = Point2D(p1.x, p1.y + (i-1) * dy)
    end
    offset = seedy
    for i = 1:seedx
        nodes[offset + i] = Point2D(p2.x + (i-1) * dx, p2.y)
    end
    offset = seedy + seedx
    for i = 1:seedy
        nodes[offset + i] = Point2D(p3.x, p3.y - (i-1) * dy)
    end
    offset = seedy + seedx + seedy
    for i = 1:seedx
        nodes[offset + i] = Point2D(p4.x - (i-1) * dx, p4.y)
    end
    nelements = nnodes
    elements = Vector{ElementRaw{Int64,Float64}}(undef,nelements)
    for i = 1:nelements-1
        i1 = i
        i2 = i + 1
        c = get_com(nodes[i1], nodes[i2])
        nv = get_nvec(nodes[i1], nodes[i2], dir)
        a = get_length(nodes[i1], nodes[i2])
        elements[i] = ElementRaw(i1, i2, c, nv, a)
    end
    i1 = nnodes
    i2 = 1
    c = get_com(nodes[i1], nodes[i2])
    nv = get_nvec(nodes[i1], nodes[i2], dir)
    a = get_length(nodes[i1], nodes[i2])
    elements[end] = ElementRaw(i1, i2, c, nv, a)
    return Part(nodes, elements, nnodes, nelements)
end

function circle(d::T1, c::Point2D; seed::T2 = 12, dir::String = "pos") where {T1<:Real, T2<:Integer}
    # creating circle based on diameter and center
    r = d/2
    phi = 2 * pi / seed
    nnodes = seed
    nodes = Vector{Point2D{Float64}}(undef,nnodes)
    for i = 1:nnodes
        nodes[i] = Point2D(c.x + r * cos(i*phi), c.y + r * sin(i*phi))
    end
    nelements = seed
    elements = Vector{ElementRaw{Int64,Float64}}(undef,nelements)
    for i = 1:nelements-1
        i1 = i
        i2 = i + 1
        c = get_com(nodes[i1], nodes[i2])
        nv = get_nvec(nodes[i1], nodes[i2], dir)
        a = get_length(nodes[i1], nodes[i2])
        elements[i] = ElementRaw(i1, i2, c, nv, a)
    end
    i1 = nnodes
    i2 = 1
    c = get_com(nodes[i1], nodes[i2])
    nv = get_nvec(nodes[i1], nodes[i2], dir)
    a = get_length(nodes[i1], nodes[i2])
    elements[end] = ElementRaw(i1, i2, c, nv, a)
    return Part(nodes, elements, nnodes, nelements)
end

function circle_open(d::T1, c::Point2D; seed::T2 = 12, leftout::T2 = 2, open::String = "right", dir::String = "pos") where {T1<:Real, T2 <: Integer}
    # creating circle based on diameter and center which is open at one side
    r = d/2
    phi = 2 * pi / seed
    # use half of leftout as start
    phio = 2 * pi / seed * (leftout/2) # phi offset
    if open == "left"
        phio = phio + 1 * pi
    end
    nnodes = seed-leftout+1
    nodes = Vector{Point2D{Float64}}(undef,nnodes)
    for i = 1:nnodes
        nodes[i] = Point2D(c.x + r * cos((i-1) * phi + phio), c.y + r * sin((i-1) * phi + phio))
    end
    nelements = seed - leftout
    elements = Vector{ElementRaw{Int64,Float64}}(undef,nelements)
    for i = 1:nelements
        i1 = i
        i2 = i + 1
        c = get_com(nodes[i1], nodes[i2])
        nv = get_nvec(nodes[i1], nodes[i2], dir)
        a = get_length(nodes[i1], nodes[i2])
        elements[i] = ElementRaw(i1, i2, c, nv, a)
    end
    return Part(nodes, elements, nnodes, nelements)
end

function cosinus(a::T1, b::T1, p1::Point2D; seed::T2 = 30, dir::String = "pos") where {T1<:Real, T2<:Integer}
    # creating cosinus wave
    nnodes = seed+1
    nodes = Vector{Point2D{Float64}}(undef,nnodes)
    for i = 1:nnodes
        nodes[i] = Point2D(p1.x + b*(i-1)*2*pi/(seed), p1.y + a*cos((i-1)*2*pi/(seed)))
    end
    nelements = seed
    elements = Vector{ElementRaw{Int64,Float64}}(undef,nelements)
    for i = 1:nelements
        i1 = i
        i2 = i + 1
        c = get_com(nodes[i1], nodes[i2])
        nv = get_nvec(nodes[i1], nodes[i2], dir)
        a = get_length(nodes[i1], nodes[i2])
        elements[i] = ElementRaw(i1, i2, c, nv, a)
    end
    return Part(nodes, elements, nnodes, nelements)
end

#= ######################################################
Other functions
###################################################### =#

function element_analysis(m::Model; printit = true)
    # analysis of elements
    area = [m.elem[i].area for i = m.elem2par[1].first:m.elem2par[end].last]
    e_length_mean = sum(area) / m.nelem
    # e_length_mean = round(e_length_mean, digits=2)
    e_length_min = minimum(area)
    e_length_max = maximum(area)
    if printit
        println("Element length:")
        println("    Mean: ", e_length_mean)
        println("    Min: ", e_length_min)
        println("    Max: ", e_length_max)
    else
        return t_max, t_min, t_mean_occ, n_occ/n_all
    end
end


text:
function get_prime_factors(channel::Channel, i::Integer)
    max_f = i÷2
    for f = 2:max_f
        n = 0
        while i % f == 0
            n += 1
            i ÷= f
        end
        if n > 0
            put!(channel,(f,n))
        end
    end
    if i > 1
        put!(channel, (i,1))
    end
end

function get_primes_below(channel::Channel, n::Integer)
    max_n::Int = n-1
    sieve_prime::Array{Bool} = trues(max_n)
    for i = 2:max_n
        if sieve_prime[i-1]
            put!(channel,i)
            for j = i ^ 2:i:max_n
                sieve_prime[j-1] = false
            end
        end
    end
end

function get_up_to_nth_prime(channel::Channel, n::Integer)
    max_n::Int = n <= 6 ? 13 : Int(floor(n*(log(n)+log(log(n)))))
    have_primes = 0
    for x in Channel(c -> get_primes_below(c,max_n+1))
        put!(channel,x)
        have_primes += 1
        if have_primes == n
            return
        end
    end
end

function get_nth_prime(n::Integer)
    res = 0
    for x in Channel(c -> get_up_to_nth_prime(c,n))
        res = x
    end
    return res
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
### Congestion management with renewable energy sources and storages
function CongestionManagementPHS80PtG(hours,
							  nodes::Nodes,
							  lines::Lines,
							  powerplants::PowerPlants,
							  renewables::Renewables,
							  storages::Storages,
							  price,
							  P_opt,
							  P_R_opt,
							  P_S_opt,
							  D_S_opt,
							  L_S_opt,
							  η_E,
							  η_M,
						      )

	h,b = calc_h_b(lines, nodes)

	T = hours
	G = powerplants.unit
	R = renewables.unit
	S = storages.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	# Efficiencies
	# η_E = 0.76		# Electrolysis
	# η_M = 0.83		# Methanation

	# Subset of power plants: GFPP (Fuel = "NaturalGas")
	G_NG = getKeyVector(powerplants.fuel, "NaturalGas")
	G_other = setdiff(G, G_NG)
	R = renewables.unit
	R_SolarPV = getKeyVector(renewables.fuel, :SolarPV)
	R_WindOffshore = getKeyVector(renewables.fuel, :WindOffshore)
	R_WindOnshore = getKeyVector(renewables.fuel, :WindOnshore)

	R_PtG = [R_SolarPV; R_WindOffshore; R_WindOnshore]
	R_other = setdiff(R, R_PtG)

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	L = lines.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, ΔP_R_up[T, R] == 0)
	@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_flow[T, L])
	@variable(CM_mod, P_load_lost[T, J] >= 0)
	@variable(CM_mod, P_gen_lost[T, J] >= 0)

	# Storage voriables
	@variable(CM_mod, P_S_up[T, S] >= 0)	# Power generation from storage unit
	# @variable(ED_mod, D_S[T, S] >= 0)	# Electricity demand from storage unit
	@variable(CM_mod, L_S[T, S] >= 0)

	# PtG extensions
	@variable(CM_mod, P_syn[T, G_NG] >= 0)	# Power generation from syngas
	@variable(CM_mod, D_PtG[T, R_PtG] >= 0)	# Electricity demand from PtG units
 	@variable(CM_mod, L_syn[T] >=0)

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	for s in S
		fix(L_S[T[1], s], 0.2*storages.storage[s]/Sbase; force = true)
	end

	# Fix SNG parameters in first period
	fix(L_syn[T[1]], 0; force = true)

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G_other, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g][t]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G_other, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g][t]/Sbase);

	# Upper generation limit for NG subset
	# @constraint(CM_mod, GenerationLimitUpNG[e = G_NG, t = T],
	#     (P_opt[t, e]/Sbase + ΔP_up[t, e] - ΔP_dn[t, e]) +
	# 	P_syn[t, e] <= powerplants.pmax[e][t]/Sbase);

	@constraint(CM_mod, GenerationLimitUpNG[e = G_NG, t = T],
		P_opt[t, e]/Sbase + ΔP_up[t, e] + P_syn[t, e] <= powerplants.pmax[e][t]/Sbase);

	# Lower generation limit for NG subset
	# @constraint(CM_mod, GenerationLimitDownNG[e = G_NG, t = T],
	# 	(P_opt[t, e]/Sbase + ΔP_up[t, e] - ΔP_dn[t, e]) +
	# 	P_syn[t, e] >= powerplants.pmin[e][t]/Sbase);

	@constraint(CM_mod, GenerationLimitDownNG[e = G_NG, t = T],
		P_opt[t, e]/Sbase - ΔP_dn[t, e]  >= powerplants.pmin[e][t]/Sbase);

	# Upper generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitUp[r = R_other, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
			renewables.infeed[r][t]/Sbase);

	# Lower generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitDown[r = R_other, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
			0);

	# # Ramp-up limit
	# @constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
	# 	(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
	# 		powerplants.rup[g]/Sbase)
	#
	# # Ramp-down limit
	# @constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
	# 	-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
	# 		powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) +
		sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r])
			for r in getKeyVector(renewables.node, j)) +
		sum((P_S_opt[t, s]/Sbase + P_S_up[t, s]) for s in getKeyVector(storages.node, j)) -
		sum((D_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) +
		sum(P_syn[t, e] for e in intersect(G_NG, getKeyVector(powerplants.node, j))) -
		sum(D_PtG[t, r] for r in intersect(R_PtG, getKeyVector(renewables.node, j))) -
		P_gen_lost[t, j] +
		sum(b[j,jj] * Θ[t, jj] for jj in J) ==
		nodes.load[j][t]/Sbase - nodes.exchange[j][t]/Sbase - P_load_lost[t, j]);

		@constraint(CM_mod, LinePowerFlowMax[l = L, t = T],
			P_flow[t, l] <= lines.pmax[l]/Sbase);

		@constraint(CM_mod, LinePowerFlowMin[l = L, t = T],
			P_flow[t, l] >= -lines.pmax[l]/Sbase);

		@constraint(CM_mod, LineFlow[l = L, t = T],
			P_flow[t, l] == sum(h[l,j] * Θ[t,j] for j in J));

			### PHS constraints
			@constraint(CM_mod, StoragePowerOutput[s = S, t = T],
		 		P_S_opt[t, s]/Sbase + P_S_up[t, s] <= storages.power[s]/Sbase);

			# Maximum storage power input
			# @constraint(CM_mod, StoragePowerInput[s = S, t = T],
			# 	D_S[t, s] <= storages.power[s]);

			# Generation CM_mod determined by remaining storage level
			@constraint(CM_mod, StorageLevelGen[s = S, t = T],
				P_S_up[t, s] <= L_S[t, s]);

			# Storage
			@constraint(CM_mod, Storage[s = S, t = T; t > T[1]],
				L_S[t-1, s] - P_S_up[t-1, s] == L_S[t, s]);

		### PtG limits
		# @constraint(CM_mod, PowerToGasLimitUp[r = R_PtG, t = T],
		# 	D_PtG[t, r] + (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
		# 	renewables.infeed[r][t]/Sbase);
		#
		# @constraint(CM_mod, PowerToGasLimitDown[r = R_PtG, t = T],
		# 	D_PtG[t, r] + (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
		# 	0);

		@constraint(CM_mod, PowerToGasLimitUp[r = R_PtG, t = T],
			P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r] - D_PtG[t, r] <=
			renewables.infeed[r][t]/Sbase);

		@constraint(CM_mod, PowerToGasLimitDown[r = R_PtG, t = T],
			P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r] - D_PtG[t, r] >=
			0);

		@constraint(CM_mod, PtGStorageGen[e = G_NG, t = T],
			sum((1/powerplants.efficiency[e]) * P_syn[t, e] for e in G_NG) <= L_syn[t]);

		# Virtual synthetic methane storage
		@constraint(CM_mod, PtGStorage[t = T; t > T[1]],
			L_syn[t-1] - sum((1/powerplants.efficiency[e]) * P_syn[t-1, e] for e in G_NG) + η_E * η_M * sum(D_PtG[t-1, r] for r in R_PtG) == L_syn[t]);

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T) +
		sum(0 * (ΔP_R_up[t, r]) * Sbase +
			(price[t] - 0) * ΔP_R_dn[t, r] *
				Sbase for r in R, t in T) +
		sum((price[t]/storages.efficiency[s]) * P_S_up[t, s] * Sbase for s in S, t in T) +
		1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
		1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T) +
		sum(price[t] * (1/(η_E*η_M)) * D_PtG[t, r] * Sbase for r in R_PtG, t in T) +
		sum(powerplants.varcost[e][1]* P_syn[t, e] * Sbase for e in G_NG, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
	ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	price = JuMP.dual.(PowerInjectionBal)
	P_gen_lost_opt = JuMP.value.(P_gen_lost) * Sbase
	P_load_lost_opt = JuMP.value.(P_load_lost) * Sbase
	P_S_up_opt = JuMP.value.(P_S_up) * Sbase
	P_syn_opt = JuMP.value.(P_syn)*Sbase
	D_PtG_opt = JuMP.value.(D_PtG)*Sbase
	L_syn_opt = JuMP.value.(L_syn)*Sbase

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
		   ΔP_R_up_opt,
		   ΔP_R_dn_opt,
	       Θ_opt,
	       P_flow_opt,
		   price,
		   P_gen_lost_opt,
		   P_load_lost_opt,
		   P_S_up_opt,
		   P_syn_opt,
   		   D_PtG_opt,
   	       L_syn_opt
	       )
end

### Congestion management with renewable energy sources and storages
function CongestionManagementPHS80PtG(hours,
							  nodes::Nodes,
							  lines::Lines,
							  powerplants::PowerPlants,
							  renewables::Renewables,
							  storages::Storages,
							  price,
							  P_opt,
							  P_R_opt,
							  P_S_opt,
							  D_S_opt,
							  L_S_opt,
							  P_syn_opt_prev,
							  D_PtG_opt_prev,
							  L_syn_opt_prev,
							  η_E,
							  η_M,
						      )

	h,b = calc_h_b(lines, nodes)

	T = hours
	G = powerplants.unit
	R = renewables.unit
	S = storages.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	# Efficiencies
	# η_E = 0.76		# Electrolysis
	# η_M = 0.83		# Methanation

	# Subset of power plants: GFPP (Fuel = "NaturalGas")
	G_NG = getKeyVector(powerplants.fuel, "NaturalGas")
	G_other = setdiff(G, G_NG)
	R = renewables.unit
	R_SolarPV = getKeyVector(renewables.fuel, :SolarPV)
	R_WindOffshore = getKeyVector(renewables.fuel, :WindOffshore)
	R_WindOnshore = getKeyVector(renewables.fuel, :WindOnshore)

	R_PtG = [R_SolarPV; R_WindOffshore; R_WindOnshore]
	R_other = setdiff(R, R_PtG)

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	L = lines.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, ΔP_R_up[T, R] == 0)
	@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_flow[T, L])
	@variable(CM_mod, P_load_lost[T, J] >= 0)
	@variable(CM_mod, P_gen_lost[T, J] >= 0)

	# Storage voriables
	@variable(CM_mod, P_S_up[T, S] >= 0)	# Power generation from storage unit
	# @variable(ED_mod, D_S[T, S] >= 0)	# Electricity demand from storage unit
	@variable(CM_mod, L_S[T, S] >= 0)

	# PtG extensions
	@variable(CM_mod, P_syn[T, G_NG] >= 0)	# Power generation from syngas
	@variable(CM_mod, D_PtG[T, R_PtG] >= 0)	# Electricity demand from PtG units
 	@variable(CM_mod, L_syn[T] >=0)

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	for s in S
		fix(L_S[T[1], s], 0.2*storages.storage[s]/Sbase; force = true)
	end

	# Fix PtG demand
	for r in R_PtG
		fix(D_PtG[T[1], r], D_PtG_opt_prev[(T[1]), r]/Sbase; force = true)
	end

	# Fix SNG generation through GfG power plants
	for e in G_NG
		fix(P_syn[T[1], e], P_syn_opt_prev[(T[1]), e]/Sbase; force = true)
	end

	# Fix SNG storage level
	fix(L_syn[T[1]], L_syn_opt_prev[T[1]]/Sbase; force = true)

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G_other, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g][t]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G_other, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g][t]/Sbase);

	# Upper generation limit for NG subset
	# @constraint(CM_mod, GenerationLimitUpNG[e = G_NG, t = T],
	#     (P_opt[t, e]/Sbase + ΔP_up[t, e] - ΔP_dn[t, e]) +
	# 	P_syn[t, e] <= powerplants.pmax[e][t]/Sbase);

	@constraint(CM_mod, GenerationLimitUpNG[e = G_NG, t = T],
		P_opt[t, e]/Sbase + ΔP_up[t, e] + P_syn[t, e] <= powerplants.pmax[e][t]/Sbase);

	# Lower generation limit for NG subset
	# @constraint(CM_mod, GenerationLimitDownNG[e = G_NG, t = T],
	# 	(P_opt[t, e]/Sbase + ΔP_up[t, e] - ΔP_dn[t, e]) +
	# 	P_syn[t, e] >= powerplants.pmin[e][t]/Sbase);

	@constraint(CM_mod, GenerationLimitDownNG[e = G_NG, t = T],
		P_opt[t, e]/Sbase - ΔP_dn[t, e]  >= powerplants.pmin[e][t]/Sbase);

	# Upper generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitUp[r = R_other, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
			renewables.infeed[r][t]/Sbase);

	# Lower generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitDown[r = R_other, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
			0);

	# # Ramp-up limit
	# @constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
	# 	(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
	# 		powerplants.rup[g]/Sbase)
	#
	# # Ramp-down limit
	# @constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
	# 	-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
	# 		powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) +
		sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r])
			for r in getKeyVector(renewables.node, j)) +
		sum((P_S_opt[t, s]/Sbase + P_S_up[t, s]) for s in getKeyVector(storages.node, j)) -
		sum((D_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) +
		sum(P_syn[t, e] for e in intersect(G_NG, getKeyVector(powerplants.node, j))) -
		sum(D_PtG[t, r] for r in intersect(R_PtG, getKeyVector(renewables.node, j))) -
		P_gen_lost[t, j] +
		sum(b[j,jj] * Θ[t, jj] for jj in J) ==
		nodes.load[j][t]/Sbase - nodes.exchange[j][t]/Sbase - P_load_lost[t, j]);

		@constraint(CM_mod, LinePowerFlowMax[l = L, t = T],
			P_flow[t, l] <= lines.pmax[l]/Sbase);

		@constraint(CM_mod, LinePowerFlowMin[l = L, t = T],
			P_flow[t, l] >= -lines.pmax[l]/Sbase);

		@constraint(CM_mod, LineFlow[l = L, t = T],
			P_flow[t, l] == sum(h[l,j] * Θ[t,j] for j in J));

			### PHS constraints
			@constraint(CM_mod, StoragePowerOutput[s = S, t = T],
		 		P_S_opt[t, s]/Sbase + P_S_up[t, s] <= storages.power[s]/Sbase);

			# Maximum storage power input
			# @constraint(CM_mod, StoragePowerInput[s = S, t = T],
			# 	D_S[t, s] <= storages.power[s]);

			# Generation CM_mod determined by remaining storage level
			@constraint(CM_mod, StorageLevelGen[s = S, t = T],
				P_S_up[t, s] <= L_S[t, s]);

			# Storage
			@constraint(CM_mod, Storage[s = S, t = T; t > T[1]],
				L_S[t-1, s] - P_S_up[t-1, s] == L_S[t, s]);

		### PtG limits
		# @constraint(CM_mod, PowerToGasLimitUp[r = R_PtG, t = T],
		# 	D_PtG[t, r] + (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
		# 	renewables.infeed[r][t]/Sbase);
		#
		# @constraint(CM_mod, PowerToGasLimitDown[r = R_PtG, t = T],
		# 	D_PtG[t, r] + (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
		# 	0);

		@constraint(CM_mod, PowerToGasLimitUp[r = R_PtG, t = T],
			P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r] - D_PtG[t, r] <=
			renewables.infeed[r][t]/Sbase);

		@constraint(CM_mod, PowerToGasLimitDown[r = R_PtG, t = T],
			P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r] - D_PtG[t, r] >=
			0);

		@constraint(CM_mod, PtGStorageGen[e = G_NG, t = T],
			sum((1/powerplants.efficiency[e]) * P_syn[t, e] for e in G_NG) <= L_syn[t]);

		# Virtual synthetic methane storage
		@constraint(CM_mod, PtGStorage[t = T; t > T[1]],
			L_syn[t-1] - sum((1/powerplants.efficiency[e]) * P_syn[t-1, e] for e in G_NG) + η_E * η_M * sum(D_PtG[t-1, r] for r in R_PtG) == L_syn[t]);

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T) +
		sum(0 * (ΔP_R_up[t, r]) * Sbase +
			(price[t] - 0) * ΔP_R_dn[t, r] *
				Sbase for r in R, t in T) +
		sum((price[t]/storages.efficiency[s]) * P_S_up[t, s] * Sbase for s in S, t in T) +
		1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
		1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T) +
		sum(price[t] * (1/(η_E*η_M)) * D_PtG[t, r] * Sbase for r in R_PtG, t in T) +
		sum(powerplants.varcost[e][1]* P_syn[t, e] * Sbase for e in G_NG, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
	ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	price = JuMP.dual.(PowerInjectionBal)
	P_gen_lost_opt = JuMP.value.(P_gen_lost) * Sbase
	P_load_lost_opt = JuMP.value.(P_load_lost) * Sbase
	P_S_up_opt = JuMP.value.(P_S_up) * Sbase
	P_syn_opt = JuMP.value.(P_syn)*Sbase
	D_PtG_opt = JuMP.value.(D_PtG)*Sbase
	L_syn_opt = JuMP.value.(L_syn)*Sbase

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
		   ΔP_R_up_opt,
		   ΔP_R_dn_opt,
	       Θ_opt,
	       P_flow_opt,
		   price,
		   P_gen_lost_opt,
		   P_load_lost_opt,
		   P_S_up_opt,
		   P_syn_opt,
   		   D_PtG_opt,
   	       L_syn_opt
	       )
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
@testset "FileIO" begin
    img = testimage("cameraman")
    ref = jpeg_decode(jpeg_encode(img))
    tmpfile = File{format"JPEG"}(joinpath(tmpdir, "tmp.jpg"))

    JpegTurbo.fileio_save(tmpfile, img)
    data = JpegTurbo.fileio_load(tmpfile)
    @test data == ref

    open(tmpfile, "w") do s
        JpegTurbo.fileio_save(s, img)
    end
    data = open(tmpfile) do s
        JpegTurbo.fileio_load(s)
    end
    @test data == ref
end


text:
# Arbitrary GM-PHD example
let
    Random.seed!(0) # reset seed

    x0 = GaussianMixture([1.0],[randn(2)],[Matrix(1.0I,2,2)])
    Z = [randn(1)]
    k(z) = 0
    gamma = GaussianMixture([1.0],[randn(2)],[Matrix(1.0I,2,2)])
    sp_dyn = [Dynamics(Matrix(1.0I,2,2),Matrix(1.0I,2,2),[0.0, 0.0])]
    sp = Spawn(GaussianMixture([1.0],[randn(2)],[Matrix(1.0I,2,2)]), sp_dyn)
    dyn = Dynamics(Matrix(1.0I,2,2),Matrix(1.0I,2,2))
    mea = Measurement(randn(1,2),ones(1,1))

    # Should construct identically
    phd1 = PHDFilter(gamma, sp, dyn, mea, 0.5, 0.5, k)
    phd2 = PHDFilter(gamma, sp, [dyn], mea, 0.5, 0.5, k) # should be same
    @test phd1.γ.N == 1
    @test phd2.γ.N == 1
    @test phd1.dyn == phd2.dyn

    # Predict/Measure and test update
    xp = predict(phd1, x0)
    s = measure(phd1, xp, Z)
    @test s.N == 6
    array_isapprox(s.w, [0.5, 0.5, 0.25, 0.425012, 0.383326, 0.191663], atol=0.002)
    @test length(s.Σ) == s.N
    @test length(s.μ) == s.N

    # Prune and test pruning
    T = 10^-5 # Truncation threshold
    U = 0.1 # Merging threshold
    J_max = 10; # Max number of Gaussian terms
    p = prune(s, T, U, J_max)
    @test p.N == 2
    array_isapprox(p.w, [0.925012, 1.32499], atol=0.002)
    @test length(p.μ) == p.N

    # Test single step updating
    sp = update(phd1, x0, Z, T, U, J_max)
    @test sp.w == p.w
    @test sp.μ == p.μ

    # Test state extraction at different thresholds
    states1 = multiple_target_state_extraction(sp, 0.8)
    @test length(states1) == 2
    array_isapprox(states1[1], [-0.132261, 0.598847], atol=0.002)

    states2 = multiple_target_state_extraction(sp, 2.0)
    @test length(states2) == 0
end

# Surveillance example
let
    Random.seed!(0)

    # Dynamics
    Δ = 1 # Sampling Period
    σ_p = 5
    F = [ Matrix{Float64}(I,2,2) Δ*Matrix{Float64}(I,2,2) ;
        zeros(Float64,2,2) Matrix{Float64}(I,2,2) ]
    Q  = σ_p^2 * [ (Δ^4)/4*Matrix{Float64}(I,2,2) (Δ^3)/4*Matrix{Float64}(I,2,2) ;
                (Δ^3)/4*Matrix{Float64}(I,2,2) (Δ^2)*Matrix{Float64}(I,2,2) ]
    Dyns = Dynamics(F,Q);

    # Measurement
    σ_m = 10
    C = [ Matrix{Float64}(I,2,2) zeros(Float64,2,2) ]
    R = σ_m^2 * Matrix{Float64}(I,2,2)
    Meas = Measurement(C,R);

    # Birth
    w1 = 0.1
    w2 = 0.1
    m_γ1 = [ 250 , 250 , 0 , 0]
    m_γ2 = [ -250 , -250 , 0 , 0 ]
    P_γ = Matrix(Diagonal([100 , 100 , 25 , 25]));

    # Spawn
    Q_β = Matrix(Diagonal([100.0 , 100.0 , 400.0 , 400.0]))
    β = GaussianMixture([w1, w2] , [m_γ1, m_γ2] , [P_γ, P_γ])
    spawn = Spawn(β , [Dyns,Dyns]);

    # Clutter
    function κ(z)
        V = 4*10^6 #Surveillance volume
        λ_c = 12.5*10^-6 #Average clutter returns per unit volume
        return λ_c*V*(1/(2000^2))
    end;

    # Misc
    Ps = 0.99 # Probability of survival
    Pd = 0.98 # Probability of detection
    T = 10^-5 # Truncation threshold
    U = 4 # Merging threshold
    J_max = 10; # Max number of Gaussian terms

    # Instantiation
    w_target1 = 1;
    w_target2 = 1;
    mu_target1 = [ -300.0 , -300.0 , 0.0 , 0.0 ]
    mu_target2 = [ 300.0 , 300.0 , 0.0 , 0.0 ]
    P_agent1 = Matrix(Diagonal([100.0 , 100.0 , 25.0 , 25.0]))
    P_agent2 = P_agent1
    γ = GaussianMixture([w_target1, w_target2] , [mu_target1, mu_target2] ,
                        [P_agent1, P_agent2]);

    # Simulation
    xsim = [γ.μ]
    Tf = 20
    for t = 1 : Δ : Tf
        MVD = MvNormal(Dyns.Q)
        x_new = [  Dyns.A*xsim[t][i] + rand(MVD,1)[:] for i=1:length(xsim[t]) ]
        if t == 13
            MVDspawn = MvNormal(Q_β)
            x1 = xsim[t][1]
            x_spawn = Dyns.A*x1 + rand(MVDspawn,1)[:]
            push!(x_new,x_spawn)
        end
        if t == 8
            x_birth = m_γ2
            push!(x_new,x_birth)
        end
        push!(xsim,x_new)
    end

    # Filter
    phd = PHDFilter(γ,spawn,Dyns,Meas,Ps,Pd,κ)
    x =  [γ]
    for t = 1:Δ:Tf
        MVD = MvNormal(Meas.R)
        z = [ Meas.C*xsim[t][i] + rand(MVD,1)[:] for i=1:length(xsim[t])]
        x_new_pruned = update(phd,x[t],z,T,U,J_max)
        push!(x, x_new_pruned)
    end

    sf = x[end]
    @test sf.N == 6
    array_isapprox(sf.w, [0.999298, 0.983239, 0.95002,
                        1.03055, 0.0204876, 0.000549536], atol=0.002)
    @test length(sf.Σ) == 6

    esf = multiple_target_state_extraction(sf,0.5)
    @test length(esf) == 4
    array_isapprox(esf[1], [-282.13, -724.875, -1.9644, -46.7133], atol=0.002)
    array_isapprox(esf[4], [-313.933, -301.799, -3.09298, -1.9735], atol=0.002)
end


text:
using DataStructures
using IterableTables: getiterator
using DataFrames

abstract type AbstractBranchTree{RT, N, B, LI, ND} <:
    AbstractTree{OneTree, RT, String, N, B}
end

import Phylo.API: _hasinbound
function _hasinbound(tree::AbstractBranchTree{<: Rooted}, name::String)
    return tree.nodes[name].inbound !== nothing
end

import Phylo.API: _getinbound
function _getinbound(tree::AbstractBranchTree{<: Rooted}, name::String)
    return tree.nodes[name].inbound
end

import Phylo.API: _addinbound!
function _addinbound!(tree::AbstractBranchTree{<: Rooted, N, B},
                      name::String,
                      inbound::B) where {N, B <: Branch}
    _hasinbound(tree, name) &&
        error("Node already has an inbound connection")
    tree.nodes[name].inbound = inbound
end

import Phylo.API: _removeinbound!
function _removeinbound!(tree::AbstractBranchTree{<: Rooted, N, B},
                         name::String,
                         inbound::B) where {N, B <: Branch}
    _hasinbound(tree, name) || error("Node has no inbound connection")
    node = tree.nodes[name]
    node.inbound === inbound ||
        error("Node has no inbound connection from branch $inbound")
    node.inbound = nothing
end

import Phylo.API: _leafinfotype
_leafinfotype(::Type{T}) where {RT, N, B, LI,
                                T <: AbstractBranchTree{RT, N, B, LI}} = LI

import Phylo.API: _nodedatatype
_nodedatatype(::Type{<: AbstractBranchTree{RT, N, B, LI, ND}}) where {RT, N, B, LI, ND} = ND

import Phylo.API: _branchdatatype
_branchdatatype(::Type{<: AbstractBranchTree}) = Nothing

import Phylo.API: _getnodes
_getnodes(bt::AbstractBranchTree) = keys(bt.nodes)

import Phylo.API: _getnodenames
_getnodenames(bt::AbstractBranchTree) = keys(bt.nodes)

import Phylo.API: _getbranches
_getbranches(bt::AbstractBranchTree) = values(bt.branches)

import Phylo.API: _getbranchnames
_getbranchnames(bt::AbstractBranchTree) = keys(bt.branches)

import Phylo.API: _getleafinfo
_getleafinfo(bt::AbstractBranchTree) = bt.leafinfos
function _getleafinfo(bt::AbstractBranchTree, leafname)
    return [info for info in getiterator(bt.leafinfos) if info[1] == leafname]
end

import Phylo.API: _setleafinfo!
function _setleafinfo!(bt::AbstractBranchTree, info)
    bt.leafinfos = info
end

import Phylo.API: _resetleaves!
function _resetleaves!(bt::AbstractBranchTree)
    # bt.leafinfos = empty!(bt.leafinfos)
    return true
end

import Phylo.API: _getnodedata
function _getnodedata(bt::AbstractBranchTree, name::String)
    return bt.nodedata[name]
end

import Phylo.API: _setnodedata!
function _setnodedata!(bt::AbstractBranchTree{RT, N, B, LI, ND},
                       name::String, value::ND) where {RT, N, B, LI, ND}
    bt.nodedata[name] = value
end

import Phylo.API: _hasnode
_hasnode(tree::AbstractBranchTree, name::String) = haskey(tree.nodes, name)

import Phylo.API: _getnode
@traitfn _getnode(tree::T,
                  nodename::String) where {T <: AbstractBranchTree;
                                           PreferNodeObjects{T}} =
                                               tree.nodes[name]

import Phylo.API: _createnode!
function _createnode!(tree::AbstractBranchTree{RT, N, B, LI, ND},
                      name::Union{String, Missing},
                      data::ND = ND()) where {RT, N, B, LI, ND}
    nodename = ismissing(name) ? _newnodelabel(tree) : name
    _hasnode(tree, nodename) && error("Node $name already present in tree")
    tree.nodes[nodename] = N(nodename)
    _setnodedata!(tree, nodename, data)
    return nodename
end

import Phylo.API: _deletenode!
function _deletenode!(tree::AbstractBranchTree, node::AbstractNode)
    name = getnodename(tree, node)
    if _hasinbound(tree, node)
        deletebranch!(tree, _getinbound(tree, node))
    end
    for b in _getoutbounds(tree, node)
        deletebranch!(tree, b)
    end
    delete!(tree.nodes, name)
    delete!(tree.nodedata, name)
    return true
end
function _deletenode!(tree::AbstractBranchTree, name::String)
    node = getnode(tree, name)
    if _hasinbound(tree, node)
        deletebranch!(tree, _getinbound(tree, node))
    end
    for b in _getoutbounds(tree, node)
        deletebranch!(tree, b)
    end
    delete!(tree.nodes, name)
    delete!(tree.nodedata, name)
    return true
end

import Phylo.API: _hasbranch
_hasbranch(tree::AbstractBranchTree, name::Int) = haskey(tree.branches, name)

import Phylo.API: _getbranch
_getbranch(tree::AbstractBranchTree, name::Int) = tree.branches[name]

import Phylo.API: _createbranch!
function _createbranch!(tree::AbstractBranchTree{RT}, source, destination,
                        len::Union{Float64, Missing}, name::Union{Int, Missing},
                        data::Nothing = nothing) where RT
    # Add the new branch
    id = ismissing(name) ? _newbranchlabel(tree) : name
    branch = Branch{RT}(id, _getnodename(tree, source),
                        _getnodename(tree, destination), len)
    tree.branches[id] = branch

    # Update the associated source and destination nodes
    _addoutbound!(tree, _getnode(tree, source), branch)
    _addinbound!(tree, _getnode(tree, destination), branch)

    # Return updated tree
    return branch
end

import Phylo.API: _deletebranch!
function _deletebranch!(tree::AbstractBranchTree, name::Int)
    # Find the branch
    branch = _getbranch(tree, name)
    # Remove branch reference from its source node
    _removeoutbound!(tree, _getnode(tree, _src(tree, branch)), branch)
    # Remove branch reference from its destination node
    _removeinbound!(tree, _getnode(tree, _dst(tree, branch)), branch)
    # Remove branch itself
    delete!(tree.branches, name)
    # Return the branch name
    return true
end
function _deletebranch!(tree::AbstractBranchTree, branch::Branch)
    # Find the branch
    name = _getbranchname(tree, branch)
    # Remove branch reference from its source node
    _removeoutbound!(tree, _getnode(tree, _src(tree, branch)), branch)
    # Remove branch reference from its destination node
    _removeinbound!(tree, _getnode(tree, _dst(tree, branch)), branch)
    # Remove branch itself
    delete!(tree.branches, name)
    # Return the branch name
    return true
end

import Phylo.API: _validate!
function _validate!(tree::TREE) where {RT, TREE <: AbstractBranchTree{RT}}
    if _leafinfotype(TREE) != Nothing && length(getiterator(tree.leafinfos)) > 0
        if Set(info[1] for info in getiterator(tree.leafinfos)) !=
            Set(_getleafnames(tree, anyorder))
            @warn "LeafInfo names do not match actual leaves of tree"
            return false
        end
    end

    if Set(keys(tree.nodedata)) != Set(_getnodenames(tree))
        @warn "Node names do not match node records of tree"
        return false
    end

    nr = nroots(tree)
    if RT == OneRoot
        if nr != 1
            @warn "Wrong number of roots for $RT tree ($nr)"
            return false
        end
    elseif RT == ManyRoots
        if nr < 1
            @warn "Wrong number of roots for $RT tree ($nr)"
            return false
        end
    end

    return true
end

import Phylo.API: _hasrootheight
function _hasrootheight(tree::AbstractBranchTree)
    return !ismissing(tree.rootheight)
end

import Phylo.API: _getrootheight
function _getrootheight(tree::AbstractBranchTree)
    return tree.rootheight
end

import Phylo.API: _setrootheight!
function _setrootheight!(tree::AbstractBranchTree, height::Float64)
    tree.rootheight = height
    return height
end

import Phylo.API: _clearrootheight!
function _clearrootheight!(tree::AbstractBranchTree)
    tree.rootheight = missing
end

"""
    BinaryTree

Binary phylogenetic tree object with known leaves and per node data
"""
mutable struct BinaryTree{RT <: Rooted, LI, ND} <:
    AbstractBranchTree{RT, BinaryNode{RT, String, Branch{RT, String}},
                       Branch{RT, String}, LI, ND}
    nodes::OrderedDict{String, BinaryNode{RT, String, Branch{RT, String}}}
    branches::Dict{Int, Branch{RT, String}}
    leafinfos::LI
    nodedata::OrderedDict{String, ND}
    rootheight::Union{Float64, Missing}
end

function BinaryTree(lt::BinaryTree{RT, LI, ND};
                    copyinfo=false) where {RT, LI, ND}
    validate!(lt) || error("Tree to copy is not valid")
    leafnames = getleafnames(lt)
    # Leaf records are conserved across trees, as could be invariant?
    leafinfos = copyinfo ? deepcopy(lt.leafinfos) : lt.leafinfos
    nodes = OrderedDict(leaf => BinaryNode{RT, String, Branch{RT, String}}(leaf)
                        for leaf in leafnames)
    branches = Dict{Int, Branch{RT, String}}()
    nodedata = OrderedDict(leaf => ND() for leaf in leafnames)
    return BinaryTree{RT, LI, ND}(nodes, branches, leafinfos, nodedata,
                                  lt.rootheight)
end

function BinaryTree{RT, LI, ND}(leafnames::Vector{String},
                                treetype::Type{BinaryTree{RT, LI, ND}} =
                                BinaryTree{RT, LI, ND};
                                rootheight::Union{Float64, Missing} =
                                missing) where {RT, LI, ND}
    nodes = OrderedDict(leaf => BinaryNode{RT, String, Branch{RT, String}}(leaf)
                        for leaf in leafnames)
    leafinfos = LI()
    nodedata = OrderedDict(leaf => ND() for leaf in leafnames)
    return BinaryTree{RT, LI, ND}(nodes, Dict{Int, Branch{RT, String}}(),
                                  leafinfos, nodedata, rootheight)
end

function BinaryTree{RT, LI, ND}(numleaves::Int = 0,
                                treetype::Type{BinaryTree{RT, LI, ND}} =
                                BinaryTree{RT, LI, ND};
                                rootheight::Union{Float64, Missing} =
                                missing) where {RT, LI, ND}
    leafnames = ["Leaf $num" for num in Base.OneTo(numleaves)]
    return BinaryTree{RT, LI, ND}(leafnames, treetype; rootheight = rootheight)
end

function BinaryTree{RT, LI, ND}(leafinfos::LI,
                                treetype::Type{BinaryTree{RT, LI, ND}} =
                                BinaryTree{RT, LI, ND};
                                rootheight::Union{Float64, Missing} =
                                missing) where {RT, LI, ND}
    leafnames = unique(info[1] for info in getiterator(leafinfos))
    nodes = OrderedDict(leaf => BinaryNode{RT, String, Branch{RT, String}}(leaf)
                        for leaf in leafnames)
    branches = Dict{Int, Branch{RT, String}}()
    nodedata = OrderedDict(leaf => ND() for leaf in leafnames)
    return BinaryTree{RT, LI, ND}(nodes, branches, leafinfos,
                                  nodedata, rootheight)
end

BinaryTree{RT}(leafinfos::LI;
               rootheight::Union{Float64, Missing} =
               missing) where {RT <: Rootedness, LI} =
               BinaryTree{RT, LI, Dict{String, Any}}(leafinfos;
                                                     rootheight = rootheight)

BinaryTree(leafinfos::LI;
           rootheight::Union{Float64, Missing} = missing) where LI =
    BinaryTree{OneRoot}(leafinfos; rootheight = rootheight)

"""
    NamedBinaryTree

Binary phylogenetic tree object with known leaves
"""
const NamedBinaryTree = BinaryTree{ManyRoots, DataFrame, Dict{String, Any}}

"""
    PolytomousTree

Phylogenetic tree object with polytomous branching, and known leaves
and per node data
"""
mutable struct PolytomousTree{RT <: Rooted, LI, ND} <:
    AbstractBranchTree{RT, Node{RT, String, Branch{RT, String}},
                       Branch{RT, String}, LI, ND}
    nodes::OrderedDict{String, Node{RT, String, Branch{RT, String}}}
    branches::Dict{Int, Branch{RT, String}}
    leafinfos::LI
    nodedata::OrderedDict{String, ND}
    rootheight::Union{Float64, Missing}
end

function PolytomousTree(lt::PolytomousTree{RT, LI, ND};
                        copyinfo=false) where {RT, LI, ND}
    validate!(lt) || error("Tree to copy is not valid")
    leafnames = getleafnames(lt)
    # Leaf records may be conserved across trees, as could be invariant?
    leafinfos = copyinfo ? deepcopy(lt.leafinfos) : lt.leafinfos
    nodes = OrderedDict(leaf => Node{RT, String, Branch{RT, String}}(leaf)
                        for leaf in leafnames)
    branches = Dict{Int, Branch{RT, String}}()
    nodedata = OrderedDict(leaf => ND() for leaf in leafnames)
    return PolytomousTree{RT, LI, ND}(nodes, branches, leafinfos,
                                      nodedata, lt.rootheight)
end

function PolytomousTree{RT, LI, ND}(leafnames::Vector{String},
                                    treetype::Type{PolytomousTree{RT, LI, ND}} =
                                    PolytomousTree{RT, LI, ND};
                                    rootheight::Union{Float64, Missing} =
                                    missing) where {RT, LI, ND}
    nodes = OrderedDict(leaf => Node{RT, String, Branch{RT, String}}(leaf)
                        for leaf in leafnames)
    leafinfos = LI()
    nodedata = OrderedDict(leaf => ND() for leaf in leafnames)
    return PolytomousTree{RT, LI, ND}(nodes, Dict{Int, Branch{RT, String}}(),
                                      leafinfos, nodedata, rootheight)
end

function PolytomousTree{RT, LI, ND}(numleaves::Int = 0,
                                    treetype::Type{PolytomousTree{RT, LI, ND}} =
                                    PolytomousTree{RT, LI, ND};
                                    rootheight::Union{Float64, Missing} =
                                    missing) where {RT, LI, ND}
    leafnames = ["Leaf $num"  for num in Base.OneTo(numleaves)]
    return PolytomousTree{RT, LI, ND}(leafnames, treetype;
                                      rootheight = rootheight)
end

function PolytomousTree{RT, LI, ND}(leafinfos::LI,
                                    treetype::Type{PolytomousTree{RT, LI, ND}} =
                                    PolytomousTree{RT, LI, ND};
                                    rootheight::Union{Float64, Missing} =
                                    missing) where {RT, LI, ND}
    leafnames = unique(info[1] for info in getiterator(leafinfos))
    nodes = OrderedDict(leaf => Node{RT, String, Branch{RT, String}}(leaf)
                        for leaf in leafnames)
    branches = Dict{Int, Branch{RT, String}}()
    nodedata = OrderedDict(leaf => ND() for leaf in leafnames)
    return PolytomousTree{RT, LI, ND}(nodes, branches, leafinfos,
                                      nodedata, rootheight)
end

PolytomousTree{RT}(leafinfos::LI;
                   rootheight::Union{Float64, Missing} = missing) where {RT, LI} =
    PolytomousTree{RT, LI, Dict{String, Any}}(leafinfos; rootheight = rootheight)

PolytomousTree(leafinfos::LI;
               rootheight::Union{Float64, Missing} = missing) where LI =
    PolytomousTree{OneRoot}(leafinfos; rootheight = rootheight)

"""
    NamedPolytomousTree

Polytomous phylogenetic tree object with known leaves
"""
const NamedTree = NamedPolytomousTree =
    PolytomousTree{ManyRoots, DataFrame, Dict{String, Any}}

import Phylo.API: _outdegree
function _outdegree(tree::BinaryTree{<: Rooted}, name::String)
    node = tree.nodes[name]
    return (node.outbounds[1] ≡ nothing ? 0 : 1) +
        (node.outbounds[2] ≡ nothing ? 0 : 1)
end
function _outdegree(tree::PolytomousTree{<: Rooted}, name::String)
    node = tree.nodes[name]
    return length(node.outbounds)
end

import Phylo.API: _hasoutboundspace
_hasoutboundspace(tree::BinaryTree{<: Rooted}, node::String) =
    _outdegree(tree, node) < 2

import Phylo.API: _getoutbounds
function _getoutbounds(tree::BinaryTree{RT}, name::String) where RT <: Rooted
    node = tree.nodes[name]
    return (node.outbounds[1] ≡ nothing ?
            (node.outbounds[2] ≡ nothing ? Branch{RT, String}[] :
             [node.outbounds[2]]) :
            (node.outbounds[2] ≡ nothing ? [node.outbounds[1]] :
             [node.outbounds[1], node.outbounds[2]]))
end
function _getoutbounds(tree::PolytomousTree{<: Rooted}, name::String)
    return tree.nodes[name].outbounds
end

import Phylo.API: _addoutbound!
function _addoutbound!(tree::BinaryTree{RT}, name::String, branch::B) where
    {RT <: Rooted, N, B <: Branch{RT, String}}
    node = tree.nodes[name]
    node.outbounds[1] ≡ nothing ?
        node.outbounds = (branch, node.outbounds[2]) :
        (node.outbounds[2] ≡ nothing ?
         node.outbounds = (node.outbounds[1], branch) :
         error("BinaryNode already has two outbound connections"))
end

function _addoutbound!(tree::PolytomousTree{<: Rooted}, name::String,
                       branch::Branch)
    node = tree.nodes[name]
    push!(node.outbounds, branch)
end

import Phylo.API: _removeoutbound!
function _removeoutbound!(tree::BinaryTree{<: Rooted}, name::String,
                          branch::Branch)
    node = tree.nodes[name]
    node.outbounds[1] ≡ branch ?
        node.outbounds = (node.outbounds[2], nothing) :
        (node.outbounds[2] ≡ branch ?
         node.outbounds = (node.outbounds[1], nothing) :
         error("BinaryNode does not have outbound connection to branch " *
               "$branch"))
end

_removeoutbound!(tree::PolytomousTree{<: Rooted}, name::String,
                 branch::Branch) =
                     filter!(x -> x ≢ branch, tree.nodes[name].outbounds)