text:
using AstroTime
using Test
using ERFA

AstroTime.update()

const speed_of_light = 299792458.0 # m/s
const astronomical_unit = 149597870700.0 # m

@timescale GMT ep tai_offset(UTC, ep)

@timescale SCET ep distance begin
    tai_offset(UTC, ep) + distance / speed_of_light
end

@testset "AstroTime" begin
    include("periods.jl")
    include("astrodates.jl")
    include("offsets.jl")
    include("epochs.jl")
    @testset "Custom Time Scales" begin
        utc = UTCEpoch(2000, 1, 1, 0, 0, 0.1)
        gmt = GMTEpoch(2000, 1, 1, 0, 0, 0.1)
        @test utc == gmt

        @test string(GMT) == "GMT"
        @test typeof(GMT) == GMTType
        @test string(gmt) == "2000-01-01T00:00:00.100 GMT"

        @test string(SCET) == "SCET"
        @test typeof(SCET) == SCETType
        scet = SCETEpoch(2000, 1, 1, 0, 8, 19.10478383615643, astronomical_unit)
        @test TAIEpoch(scet) ≈ TAIEpoch(utc)
        @test SCETEpoch(j2000(scet), 0.0days, astronomical_unit) ≈ scet
    end
end


text:
import Dates
import Distributions
import Random
import TimeZones

const _kwargs_docstring_for_sequence_and_random_date_shift = """
- `round_to::Dates.Period`: Resolution to which all intervals should be rounded.
- `day::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Day` shift amount will be sampled.
- `hour::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Hour` shift amount will be sampled.
- `minute::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Minute` shift amount will be sampled.
- `second::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Second` shift amount will be sampled.
- `millisecond::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Millisecond` shift amount will be sampled.
- `microsecond::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Microsecond` shift amount will be sampled.
- `nanosecond::::Union{Distributions.Sampleable, Nothing}`: Probability distribution from which the `Nanosecond` shift amount will be sampled.
"""

"""
    sequence_and_random_date_shift(rng::Random.AbstractRNG,
                                   input_dt_list::Vector{Dates.DateTime};
                                   kwargs...)

## Arguments
- `rng::Random.AbstractRNG`: Random number generator that will be used for sampling.
- `input_dt_list::Vector{Dates.DateTime}`: A vector of `DateTime`s.

## Keyword Arguments
- `time_zone::Dates.TimeZone`: Time zone for the input dates.
$(_kwargs_docstring_for_sequence_and_random_date_shift)

## Example

```jldoctest; setup = :(import Random; Random.seed!(123))
julia> using Dates

julia> using DateShifting

julia> using Distributions

julia> using Random

julia> using TimeZones

julia> dates = [
           DateTime("2000-01-01T00:00:00"),
           DateTime("2000-02-01T00:00:00"),
           DateTime("2000-01-05T00:00:00"),
           DateTime("2000-01-02T04:05:06"),
           DateTime("2000-01-02T01:02:03"),
       ]
5-element Array{DateTime,1}:
 2000-01-01T00:00:00
 2000-02-01T00:00:00
 2000-01-05T00:00:00
 2000-01-02T04:05:06
 2000-01-02T01:02:03

julia> sequence, shifted_dates = sequence_and_random_date_shift(
           Random.GLOBAL_RNG,
           dates;
           round_to = Day(1),
           time_zone = TimeZone("America/New_York"),
           day = DiscreteUniform(-31, 31),
           hour = DiscreteUniform(-24, 24),
           minute = DiscreteUniform(-60, 60),
           second = DiscreteUniform(-60, 60),
           millisecond = DiscreteUniform(-1000, 1000),
           microsecond = DiscreteUniform(-1000, 1000),
           nanosecond = DiscreteUniform(-1000, 1000),
       )
([1, 5, 4, 3, 2], TimeZones.ZonedDateTime[ZonedDateTime(1999, 12, 5, tz"America/New_York"), ZonedDateTime(2000, 1, 5, tz"America/New_York"), ZonedDateTime(1999, 12, 9, tz"America/New_York"), ZonedDateTime(1999, 12, 6, tz"America/New_York"), ZonedDateTime(1999, 12, 6, tz"America/New_York")])

julia> sequence
5-element Array{Int64,1}:
 1
 5
 4
 3
 2

julia> shifted_dates
5-element Array{ZonedDateTime,1}:
 1999-12-05T00:00:00-05:00
 2000-01-05T00:00:00-05:00
 1999-12-09T00:00:00-05:00
 1999-12-06T00:00:00-05:00
 1999-12-06T00:00:00-05:00
```
"""
function sequence_and_random_date_shift(rng::Random.AbstractRNG,
                                        input_dt_list::Vector{Dates.DateTime};
                                        time_zone::Dates.TimeZone,
                                        kwargs...)
    input_zdt_list = [TimeZones.ZonedDateTime(x, time_zone) for x in input_dt_list]
    return sequence_and_random_date_shift(rng, input_zdt_list; time_zone = time_zone, kwargs...)
end

"""
    sequence_and_random_date_shift(rng::Random.AbstractRNG,
                                   input_zdt_list::Vector{TimeZones.ZonedDateTime};
                                   kwargs...)

## Arguments
- `rng::Random.AbstractRNG`: Random number generator that will be used for sampling.
- `input_zdt_list::Vector{TimeZones.ZonedDateTime}`: A vector of `ZonedDateTime`s.

## Keyword Arguments
- `time_zone::Dates.TimeZone`: Time zone to which all dates should be converted.
$(_kwargs_docstring_for_sequence_and_random_date_shift)

## Example

```jldoctest; setup = :(import Random; Random.seed!(123))
julia> using Dates

julia> using DateShifting

julia> using Distributions

julia> using Random

julia> using TimeZones

julia> dates = [
           ZonedDateTime(DateTime("2000-01-01T00:00:00"), tz"America/New_York"),
           ZonedDateTime(DateTime("2000-02-01T00:00:00"), tz"America/New_York"),
           ZonedDateTime(DateTime("2000-01-05T00:00:00"), tz"America/New_York"),
           ZonedDateTime(DateTime("2000-01-02T03:05:06"), tz"America/Chicago"),
           ZonedDateTime(DateTime("2000-01-02T01:02:03"), tz"America/New_York"),
       ]
5-element Array{ZonedDateTime,1}:
 2000-01-01T00:00:00-05:00
 2000-02-01T00:00:00-05:00
 2000-01-05T00:00:00-05:00
 2000-01-02T03:05:06-06:00
 2000-01-02T01:02:03-05:00

julia> sequence, shifted_dates = sequence_and_random_date_shift(
           Random.GLOBAL_RNG,
           dates;
           round_to = Day(1),
           time_zone = TimeZone("America/New_York"),
           day = DiscreteUniform(-31, 31),
           hour = DiscreteUniform(-24, 24),
           minute = DiscreteUniform(-60, 60),
           second = DiscreteUniform(-60, 60),
           millisecond = DiscreteUniform(-1000, 1000),
           microsecond = DiscreteUniform(-1000, 1000),
           nanosecond = DiscreteUniform(-1000, 1000),
       )
([1, 5, 4, 3, 2], TimeZones.ZonedDateTime[ZonedDateTime(1999, 12, 5, tz"America/New_York"), ZonedDateTime(2000, 1, 5, tz"America/New_York"), ZonedDateTime(1999, 12, 9, tz"America/New_York"), ZonedDateTime(1999, 12, 6, tz"America/New_York"), ZonedDateTime(1999, 12, 6, tz"America/New_York")])

julia> sequence
5-element Array{Int64,1}:
 1
 5
 4
 3
 2

julia> shifted_dates
5-element Array{ZonedDateTime,1}:
 1999-12-05T00:00:00-05:00
 2000-01-05T00:00:00-05:00
 1999-12-09T00:00:00-05:00
 1999-12-06T00:00:00-05:00
 1999-12-06T00:00:00-05:00
```
"""
function sequence_and_random_date_shift(rng::Random.AbstractRNG,
                                        input_zdt_list::Vector{TimeZones.ZonedDateTime};
                                        round_to::Dates.Period,
                                        time_zone::Dates.TimeZone,
                                        day::Union{Distributions.Sampleable, Nothing},
                                        hour::Union{Distributions.Sampleable, Nothing},
                                        minute::Union{Distributions.Sampleable, Nothing},
                                        second::Union{Distributions.Sampleable, Nothing},
                                        millisecond::Union{Distributions.Sampleable, Nothing},
                                        microsecond::Union{Distributions.Sampleable, Nothing},
                                        nanosecond::Union{Distributions.Sampleable, Nothing})
    sequence = sortperm(input_zdt_list)

    input_zdt_list_sametimezone = [TimeZones.astimezone(x, time_zone) for x in input_zdt_list]
    always_assert(sortperm(input_zdt_list_sametimezone) == sequence)

    day_shift_amount = _draw_sample(rng, Dates.Day, day)
    hour_shift_amount = _draw_sample(rng, Dates.Hour, hour)
    minute_shift_amount = _draw_sample(rng, Dates.Minute, minute)
    second_shift_amount = _draw_sample(rng, Dates.Second, second)
    millisecond_shift_amount = _draw_sample(rng, Dates.Millisecond, millisecond)
    microsecond_shift_amount = _draw_sample(rng, Dates.Microsecond, microsecond)
    nanosecond_shift_amount = _draw_sample(rng, Dates.Nanosecond, nanosecond)

    total_shift_amount = day_shift_amount + hour_shift_amount +
                                            minute_shift_amount +
                                            second_shift_amount +
                                            millisecond_shift_amount +
                                            microsecond_shift_amount +
                                            nanosecond_shift_amount

    shifted_dates_nonrounded = [x + total_shift_amount for x in input_zdt_list_sametimezone]
    shifted_dates = [round(x, round_to) for x in shifted_dates_nonrounded]

    always_assert(all(TimeZones.timezone.(shifted_dates_nonrounded) .== time_zone))
    always_assert(all(TimeZones.timezone.(shifted_dates) .== time_zone))

    return sequence, shifted_dates
end

ext: jl


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
ENV["MOCHA_USE_NATIVE_EXT"] = "true"

using Mocha

backend = CPUBackend()
init(backend)
img_width, img_height, img_channels = (224, 224, 3)
batch_size = 64  # could be larger if you want to classify a bunch of images at a time

layers = [
  MemoryDataLayer(name="data", tops=[:data], batch_size=batch_size,
      data = Array[zeros(img_width, img_height, img_channels, batch_size)])
  ConvolutionLayer(name="conv1", tops=[:conv1], bottoms=[:data],
      kernel=(3,3), stride=(1,1), n_filter=64, neuron=Neurons.ReLU())
  PoolingLayer(name="pool1", tops=[:pool1], bottoms=[:conv1],
      kernel=(2,2), stride=(2,2), pooling=Pooling.Max())

  ConvolutionLayer(name="conv2", tops=[:conv2], bottoms=[:pool1],
      kernel=(3,3), pad=(1,1), n_filter=128, neuron=Neurons.ReLU())
  PoolingLayer(name="pool2", tops=[:pool2], bottoms=[:conv2],
      kernel=(2,2), stride=(2,2), pooling=Pooling.Max())

  ConvolutionLayer(name="conv3_1", tops=[:conv3_1], bottoms=[:pool2],
      kernel=(3,3), pad=(1,1), n_filter=256, neuron=Neurons.ReLU())
  ConvolutionLayer(name="conv3_2", tops=[:conv3_2], bottoms=[:conv3_1],
      kernel=(3,3), pad=(1,1), n_filter=256, neuron=Neurons.ReLU())
  PoolingLayer(name="pool3", tops=[:pool3], bottoms=[:conv3_2],
      kernel=(2,2), stride=(2,2), pooling=Pooling.Max())

  ConvolutionLayer(name="conv4_1", tops=[:conv4_1], bottoms=[:pool3],
      kernel=(3,3), pad=(1,1), n_filter=512, neuron=Neurons.ReLU())
  ConvolutionLayer(name="conv4_2", tops=[:conv4_2], bottoms=[:conv4_1],
      kernel=(3,3), pad=(1,1), n_filter=512, neuron=Neurons.ReLU())
  PoolingLayer(name="pool4", tops=[:pool4], bottoms=[:conv4_2],
      kernel=(2,2), stride=(2,2), pooling=Pooling.Max())

  ConvolutionLayer(name="conv5_1", tops=[:conv5_1], bottoms=[:pool4],
      kernel=(3,3), pad=(1,1), n_filter=512, neuron=Neurons.ReLU())
  ConvolutionLayer(name="conv5_2", tops=[:conv5_2], bottoms=[:conv5_1],
      kernel=(3,3), pad=(1,1), n_filter=512, neuron=Neurons.ReLU())
  PoolingLayer(name="pool5", tops=[:pool5], bottoms=[:conv5_2],
      kernel=(2,2), stride=(2,2), pooling=Pooling.Max())

  InnerProductLayer(name="fc6", tops=[:fc6], bottoms=[:pool5],
      output_dim=4096)
  InnerProductLayer(name="fc7", tops=[:fc7], bottoms=[:fc6],
      output_dim=4096)
  InnerProductLayer(name="fc8", tops=[:fc8], bottoms=[:fc7],
      output_dim=1000)
]

net = Net("imagenet", backend, layers)

init(net)

for i = 1:3
    forward(net)
    backward(net)
end

num_trials = 10

forward_time = 0.0
backward_time = 0.0
for i = 1:num_trials
    tic()
    forward(net)
    forward_time += toq()
    tic()
    backward(net)
    backward_time += toq()
end
println("Avg forward time for $num_trials runs: $(forward_time / num_trials * 1000.0)ms")
println("Avg backward time for $num_trials runs: $(backward_time / num_trials * 1000.0)ms")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ExternalMetricStatus indicates the current value of a global metric not associated with any Kubernetes object.

    IoK8sApiAutoscalingV2beta1ExternalMetricStatus(;
        currentAverageValue=nothing,
        currentValue=nothing,
        metricName=nothing,
        metricSelector=nothing,
    )

    - currentAverageValue::IoK8sApimachineryPkgApiResourceQuantity : currentAverageValue is the current value of metric averaged over autoscaled pods.
    - currentValue::IoK8sApimachineryPkgApiResourceQuantity : currentValue is the current value of the metric (as a quantity)
    - metricName::String : metricName is the name of a metric used for autoscaling in metric system.
    - metricSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : metricSelector is used to identify a specific time series within a given metric.
"""
mutable struct IoK8sApiAutoscalingV2beta1ExternalMetricStatus <: SwaggerModel
    currentAverageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: currentAverageValue
    currentValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: currentValue
    metricName::Any # spec type: Union{ Nothing, String } # spec name: metricName
    metricSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: metricSelector

    function IoK8sApiAutoscalingV2beta1ExternalMetricStatus(;currentAverageValue=nothing, currentValue=nothing, metricName=nothing, metricSelector=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricStatus, Symbol("currentAverageValue"), currentAverageValue)
        setfield!(o, Symbol("currentAverageValue"), currentAverageValue)
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricStatus, Symbol("currentValue"), currentValue)
        setfield!(o, Symbol("currentValue"), currentValue)
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricStatus, Symbol("metricName"), metricName)
        setfield!(o, Symbol("metricName"), metricName)
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricStatus, Symbol("metricSelector"), metricSelector)
        setfield!(o, Symbol("metricSelector"), metricSelector)
        o
    end
end # type IoK8sApiAutoscalingV2beta1ExternalMetricStatus

const _property_map_IoK8sApiAutoscalingV2beta1ExternalMetricStatus = Dict{Symbol,Symbol}(Symbol("currentAverageValue")=>Symbol("currentAverageValue"), Symbol("currentValue")=>Symbol("currentValue"), Symbol("metricName")=>Symbol("metricName"), Symbol("metricSelector")=>Symbol("metricSelector"))
const _property_types_IoK8sApiAutoscalingV2beta1ExternalMetricStatus = Dict{Symbol,String}(Symbol("currentAverageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("currentValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("metricName")=>"String", Symbol("metricSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1ExternalMetricStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1ExternalMetricStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1ExternalMetricStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1ExternalMetricStatus)
    (getproperty(o, Symbol("currentValue")) === nothing) && (return false)
    (getproperty(o, Symbol("metricName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricStatus }, name::Symbol, val)
end


text:
using Underscores
using Utils

const MASK_PREFIX = "mask = "

function run_decoder(f)
    mem::Dict{UInt16,UInt64} = Dict()
    set_mask::UInt64 = 0
    unset_mask = typemax(UInt64)

    for line in eachline(f)
        if startswith(line, MASK_PREFIX)
            set_mask = 0
            unset_mask = typemax(typeof(unset_mask))
            for bit in length(MASK_PREFIX)+1:length(line)
                set_mask <<= 1
                unset_mask <<= 1
                if line[bit] != '0'
                    unset_mask |= 1
                end
                if line[bit] == '1'
                    set_mask |= 1
                end
            end
        else
            m = match(r"^mem\[(\d+)\] = (\d+)$", line)
            addr = parse(keytype(mem), m[1])
            value = parse(valtype(mem), m[2])
            mem[addr] = (value & unset_mask) | set_mask
        end
    end

    sum(values(mem))
end

function set_all(mem, addr, mask, value, lookat)
    if lookat == 0
        mem[addr] = value
    else
        while lookat != 0 && (mask & lookat) == 0
            lookat >>= 1
        end

        set_all(mem, addr | lookat, mask, value, lookat >> 1)
        set_all(mem, addr & ~lookat, mask, value, lookat >> 1)
    end
end

function run_decoder2(f)
    mem::Dict{UInt64,UInt64} = Dict()
    set_mask::UInt64 = 0
    floating_mask::UInt64 = 0

    for (i, line) in enumerate(eachline(f))
        if startswith(line, MASK_PREFIX)
            set_mask = 0
            floating_mask = 0
            for bit in length(MASK_PREFIX)+1:length(line)
                set_mask <<= 1
                floating_mask <<= 1
                if line[bit] == '1'
                    set_mask |= 1
                elseif line[bit] == 'X'
                    floating_mask |= 1
                end
            end
        else
            m = match(r"^mem\[(\d+)\] = (\d+)$", line)
            addr = parse(keytype(mem), m[1]) | set_mask
            value = parse(valtype(mem), m[2])
            set_all(mem, addr, floating_mask, value, 2^35)
        end
    end

    sum(values(mem))
end

function run_decoder()
    open(joinpath(@__DIR__, "14_docking-data.txt"), "r") do f
        run_decoder(f)
    end
end

function run_decoder2()
    open(joinpath(@__DIR__, "14_docking-data.txt"), "r") do f
        run_decoder2(f)
    end
end

ex1() = run_decoder()
ex2() = run_decoder2()

println("Final memory sum: ", ex1())
println("Final floating memory sum: ", ex2())

using Test

@testset "Docking Data" begin
    @testset "example 1" begin
        input = """
        mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X
        mem[8] = 11
        mem[7] = 101
        mem[8] = 0
        """

        example = run_decoder(IOBuffer(input))

        @test example == 165
    end

    @testset "example 2" begin
        input = """
        mask = 000000000000000000000000000000X1001X
        mem[42] = 100
        mask = 00000000000000000000000000000000X0XX
        mem[26] = 1
        """

        example = run_decoder2(IOBuffer(input))

        @test example == 208
    end

    @testset "results" begin
        @test ex1() == 18630548206046
        @test ex2() == 4254673508445
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using MuZero
import Zygote

include("game.jl")

learner = Learner(Connect4)

# Play a quick game to make sure
# the basic game dynamics are working
#learner |> play

learner |> tests |> svg
learner |> test_results |> svg

selfplay!(learner; limit=10)

learner |> randgame |> svg

# still a fence post problem here
s = learner |> randtrainingsample
learner |> randtrainingsample |> svg

learner |> randtrainingbatch |> show

using Base.Threads: @spawn
for i = 1:num_actors(learner)
    @spawn selfplay!(learner)
end

length(learner.replay_buffer)

@time fit!(learner, steps=1)


text:
using AbstractPlotting, CairoMakie
using MakieRecipes
using MakieRecipes.RecipesBase
using MarketData, TimeSeries

# ## The simplest example model

const i = Ref(0)

macro test_and_save(arg)
    return quote
        i[] += 1
        @test_nowarn save("test_$(i[]).png", $arg)
    end
end

struct T end

RecipesBase.@recipe function plot(::T, n = 1; customcolor = :green)
    markershape --> :auto        # if markershape is unset, make it :auto
    markercolor :=  customcolor       # force markercolor to be customcolor
    xrotation   --> 45           # if xrotation is unset, make it 45
    zrotation   --> 90           # if zrotation is unset, make it 90
    rand(10,n)                   # return the arguments (input data) for the next recipe
end

@test_and_save recipeplot(T(); seriestype = :path)

# ## Testing out series decomposition

sc = Scene()
@test_and_save recipeplot!(sc, rand(10, 2); seriestype = :scatter)
@test_and_save recipeplot!(sc, 1:10, rand(10, 1); seriestype = :path)

# ## Distributions

# ###

# ## Differential Equations

using OrdinaryDiffEq, StochasticDiffEq, DiffEqNoiseProcess

# ### A simple exponential growth model

f(u,p,t) = 1.01.*u
u0 = [1/2, 1]
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
@test_and_save recipeplot(sol)

# ### Matrix DiffEq

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
@test_and_save recipeplot(sol)

# ### Stochastic DiffEq

f(du,u,p,t) = (du .= u)
g(du,u,p,t) = (du .= u)
u0 = rand(4,2)

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRIW1())
@test_and_save recipeplot(sol)

# ## Phylogenetic tree
using Phylo
assetpath = joinpath(dirname(pathof(MakieRecipes)), "..", "docs", "src", "assets")
hummer = open(t -> parsenewick(t, NamedPolytomousTree), joinpath(assetpath, "hummingbirds.tree"))
evolve(tree) = Phylo.map_depthfirst((val, node) -> val + randn(), 0., tree, Float64)
trait = evolve(hummer)

scp = recipeplot!(
    Scene(scale_plot = false, show_axis = false),
    hummer;
    treetype = :fan,
    line_z = trait,
    linewidth = 5,
    showtips = false,
    cgrad = :RdYlBu,
    seriestype = :path
)

# ## GraphRecipes
using GraphRecipes

# ### Julia AST with GraphRecipes

code = quote
    function mysum(list)
        out = 0
        for value in list
            out += value
        end
        out
    end
end

@test_and_save recipeplot(code; fontsize = 12, shorten = 0.01, axis_buffer = 0.15, nodeshape = :rect)

# ### Type tree with GraphRecipes

@test_and_save recipeplot(AbstractFloat; method = :tree, fontsize = 10)

# Timeseries with market data
@test_and_save recipeplot(MarketData.ohlc; seriestype = :path)


text:
@testset "Cfunction" begin
	lib = Clibrary()

	if !Sys.iswindows()
		f1 = Cfunction{Clong, Tuple{Ptr{Clong}}}(lib, :time)
		@test eltype(f1) <: Cfunction{Clong, Tuple{Ptr{Clong}}}
		@test typeof(f1(C_NULL)) === Clong
		@test f1(C_NULL) isa Clong
		@test_throws MethodError f1(0)

		before = f1(C_NULL)
		sleep(3)
		after = f1(C_NULL)
		@test before < after
	end

	f2 = Cfunction{Clong, Tuple{}}(lib, :jl_gc_total_bytes)
	@test eltype(f2) <: Cfunction{Clong, Tuple{}}
	@test typeof(f2()) === Clong
	@test f2() isa Clong
	@test_throws MethodError f2("no arguments, please!")

	if !Sys.iswindows()
		f3 = Cfunction{Cint, Tuple{Ptr{Cchar}, Cstring, Vararg}}(lib, :sprintf)
		@test eltype(f3) <: Cfunction{Cint, Tuple{Ptr{Cchar}, Cstring, Vararg}}
		str = zeros(Cchar, 100)
		@test typeof(f3(str, "")) === Cint
		@test f3(str, "%s %ld\n", "testing printf", 1234) == 20
		@test unsafe_string(pointer(str)) == "testing printf 1234\n"
		@test f3(str, "%s i%c %ld great test of CBinding.jl v%3.1lf%c\n", "This", 's', 1, 0.1, '!') == length("This is 1 great test of CBinding.jl v0.1!\n")
		@test unsafe_string(pointer(str)) == "This is 1 great test of CBinding.jl v0.1!\n"
		@test_throws MethodError f3(1234)
		@test_throws MethodError f3(1234, "still wrong")
	end

	(Cadd, add) = Cfunction{Cint, Tuple{Cint, Cint}}() do val1, val2
		return val1+val2
	end
	@test eltype(Cadd) <: Cfunction{Cint, Tuple{Cint, Cint}}
	@test eltype(Cadd) === Cfunction{Cint, Tuple{Cint, Cint}, CBinding.default_convention(Tuple{Cint, Cint})}
	@test typeof(add.f) <: Function
	@test typeof(Cadd(Cint(10), Cint(3))) === typeof(add.f(Cint(10), Cint(3)))
	@test Cadd(Cint(10), Cint(3)) == add.f(Cint(10), Cint(3))

	if !Sys.iswindows()
		@cstruct time_t {
			val::Clong
		}
		@cstruct tm {
			sec::Cint
			min::Cint
			hour::Cint
			mday::Cint
			mon::Cint
			year::Cint
			wday::Cint
			yday::Cint
			isdst::Cint
		}
		time = Cfunction{time_t, Tuple{Ptr{time_t}}}(lib, :time)
		localtime = Cfunction{Ptr{tm}, Tuple{Ptr{time_t}}}(lib, :localtime)

		t = Ref(time_t())
		@test t[].val == 0
		time(t)
		@test t[].val != 0

		p = localtime(t)
		@test p isa Ptr
		@test p != C_NULL

		x = unsafe_load(p)
		@test 2000 < 1900+x.year < 2100  # if this package is still around in 2100...

		mktime = Cfunction{time_t, Tuple{Ptr{tm}}}(lib, :mktime)
		u = mktime(p)
		@test u.val == t[].val
	end
end


text:
using PackageCompiler

build_shared_lib("example.jl", init_shared=true, verbose=true)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Phylo.API
import LightGraphs: src, dst, indegree, outdegree, degree
using SimpleTraits

# AbstractTree/Node/Branch type methods
"""
    treetype(::Type{AbstractTree})

Returns tree number (OneTree, ManyTrees) from a tree type.
"""
treetype(::Type{T}) where {TT <: TreeType, RT, NL, N, B,
                           T <: AbstractTree{TT, RT, NL, N, B}} = TT

"""
    roottype(::Type{AbstractTree})
    roottype(::Type{AbstractNode})
    roottype(::Type{AbstractBranch})

Returns root type from a tree type.
"""
roottype(::Type{T}) where {TT, RT <: Rootedness, NL, N, B,
                           T <: AbstractTree{TT, RT, NL, N, B}} = RT
roottype(::Type{N}) where {RT <: Rootedness, NL,
                           N <: AbstractNode{RT, NL}} = RT
roottype(::Type{B}) where {RT <: Rootedness, NL,
                           B <: AbstractBranch{RT, NL}} = RT

"""
    nodenametype(::Type{AbstractTree})
    nodenametype(::Type{AbstractNode})
    nodenametype(::Type{AbstractBranch})

Returns type of node names from a tree type.
"""
nodenametype(::Type{T}) where {TT, RT, NL, N, B,
                               T <: AbstractTree{TT, RT, NL, N, B}} = NL
nodenametype(::Type{N}) where {RT <: Rootedness, NL,
                               N <: AbstractNode{RT, NL}} = NL
nodenametype(::Type{B}) where {RT <: Rootedness, NL,
                               B <: AbstractBranch{RT, NL}} = NL

"""
    nodetype(::Type{AbstractTree})

Returns type of nodes from a tree type.
"""
nodetype(::Type{T}) where {TT, RT, NL, N <: AbstractNode, B,
                           T <: AbstractTree{TT, RT, NL, N, B}} = N

"""
    branchnametype(::AbstractTree)

Returns type of branch names from a branch type.
"""
branchnametype(::Type{<: AbstractTree}) = Int
branchnametype(::Type{<: AbstractNode}) = Int
branchnametype(::Type{<: AbstractBranch}) = Int

"""
    branchtype(::Type{AbstractTree})

Returns type of branches from a tree type.
"""
branchtype(::Type{T}) where {TT, RT, NL, N, B <: AbstractBranch,
                             T <: AbstractTree{TT, RT, NL, N, B}} = B

"""
    treenametype(::Type{AbstractTree})

Returns the name type for a tree type.
"""
treenametype(treetype::Type{<: AbstractTree}) = _treenametype(treetype)

"""
    leafinfotype(::Type{<: AbstractTree})

retrieve the leaf info type of a tree.
"""
leafinfotype(t::Type{<: AbstractTree}) = _leafinfotype(t)

"""
    nodedatatype(::Type{<: AbstractTree})

retrieve the node info type of a tree.
"""
nodedatatype(t::Type{<: AbstractTree}) = _nodedatatype(t)

"""
    branchdatatype(::Type{<: AbstractTree})

retrieve the branch info type of a tree.
"""
branchdatatype(t::Type{<: AbstractTree}) = _branchdatatype(t)

"""
    branchdims(::Type{<: AbstractTree})

retrieve the dimensions of the branch lengths for the tree.
"""
branchdims(t::Type{<: AbstractTree}) = _branchdims(t)

# AbstractTree methods
"""
    ntrees(tree::AbstractTree)

Returns the number of trees in a tree object, 1 for a OneTree tree type, and
the count of trees for a ManyTrees type.
"""
ntrees(tree::AbstractTree) = _ntrees(tree)

"""
    gettrees(tree::AbstractTree)

Returns an iterable of trees.
"""
gettrees(tree::AbstractTree) = _gettrees(tree)

"""
    gettree(tree::AbstractTree, label)

Returns a single OneTree object tree corresponding to the label given.
"""
function gettree end
gettree(tree::AbstractTree{OneTree}) = tree
function gettree(trees::AbstractTree{ManyTrees})
    @assert _ntrees(trees) == 1 "Must be only one tree, found $(_ntrees(trees))"
    return first(_gettrees(tree))
end
gettree(tree::AbstractTree, label) = _gettree(tree, label)

"""
    gettreenames(tree::AbstractTree)

Returns the names of the trees.
"""
function gettreenames end
gettreenames(tree::AbstractTree{OneTree}) = [_gettreename(tree)]
gettreenames(trees::AbstractTree{ManyTrees}) = _gettreenames(trees)

"""
    gettreeinfo(tree::AbstractTree)
    gettreeinfo(tree::AbstractTree, treename)

Returns the info data associated with the tree(s).
"""
function gettreeinfo end
gettreeinfo(tree::AbstractTree) = _gettreeinfo(tree)
gettreeinfo(tree::AbstractTree, treename) = _gettreeinfo(tree, treename)

"""
    gettreename(tree::AbstractTree)

Returns the name of the single tree.
"""
gettreename(tree::AbstractTree{OneTree}) = _gettreename(tree)
function gettreename(trees::AbstractTree{ManyTrees})
    @assert _ntrees(trees) == 1 "Must be only one tree. found $(_ntrees(tree))"
    return first(_gettreenames(trees))
end

"""
    nleaves(::AbstractTree)

Returns the number of leaves (tips) in a tree.
"""
nleaves(tree::AbstractTree) = _nleaves(tree)

"""
    nroots(::AbstractTree)

Returns the number of roots in a tree. For OneTree types, Unrooted trees will
return 0, OneRoot trees should return 1, and manyroots tree (ones with multiple
subtrees) will return the number of subtrees. ManyTrees types will return a
Dict of counts of the number of roots for each tree in the set.
"""
function nroots end
nroots(::AbstractTree{OneTree, Unrooted}) = 0
nroots(tree::AbstractTree{OneTree, <: Rooted}) = _nroots(tree)
nroots(trees::AbstractTree{ManyTrees}, name) = nroots(gettree(trees, name))
nroots(trees::AbstractTree{ManyTrees}) =
    Dict(name => nroots(gettree(trees, name)) for name in _gettreenames(trees))

"""
    getroots(::AbstractTree)
    getroots(::AbstractTree, id)

Returns a vector containing the root(s) of a single (OneTree) tree
or a set of (ManyTrees) trees.
"""
function getroots end
getroots(tree::AbstractTree{OneTree}) = _getroots(tree)
getroots(tree::AbstractTree{OneTree, Unrooted, NL, N}) where {NL, N} = N[]
getroots(trees::AbstractTree{ManyTrees}, name) = getroots(gettree(trees, name))
getroots(trees::AbstractTree{ManyTrees}) =
    Dict(name => getroots(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    getroot(::AbstractTree)

Returns the root of a single tree (must be only one tree for a ManyTrees tree).
"""
function getroot end
getroot(tree::AbstractTree{OneTree, <: Rooted}) = _getroot(tree)
getroot(trees::AbstractTree{ManyTrees}, name) = getroot(gettree(trees, name))
getroot(trees::AbstractTree{ManyTrees}) =
    Dict(name => getroot(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    getnodes(::AbstractTree[, ::TraversalOrder])

Returns the vector of nodes of a single tree, or a Dict of vectors of nodes
for multiple trees.
"""
function getnodes end
getnodes(tree::AbstractTree{OneTree}, order::TraversalOrder = preorder) =
    _getnodes(tree, order)
getnodes(trees::AbstractTree{ManyTrees}, name,
         order::TraversalOrder = preorder) =
             getnodes(gettree(trees, name), order)
getnodes(trees::AbstractTree{ManyTrees}, order::TraversalOrder = preorder) =
    Dict(name => getnodes(gettree(trees, name), order)
         for name in _gettreenames(trees))

"""
    nnodes(::AbstractTree)

Returns the number of nodes of a single tree, or a Dict of numbers of nodes
for multiple trees.
"""
function nnodes end
nnodes(tree::AbstractTree{OneTree}) = _nnodes(tree)
nnodes(trees::AbstractTree{ManyTrees}, name) = nnodes(gettree(trees, name))
nnodes(trees::AbstractTree{ManyTrees}) =
    Dict(name => nnodes(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    nbranches(::AbstractTree)

Returns the number of branches of a single tree, or a Dict of numbers of
branches for multiple trees.
"""
function nbranches end
nbranches(tree::AbstractTree{OneTree}) = _nbranches(tree)
nbranches(trees::AbstractTree{ManyTrees}, name) =
    nbranches(gettree(trees, name))
nbranches(trees::AbstractTree{ManyTrees}) =
    Dict(name => nbranches(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    getnodenames(::AbstractTree[, ::TraversalOrder])

Return a vector of node names of a single tree (identified by id for a
ManyTrees tree), or a Dict of vectors of node names for multiple trees.
"""
function getnodenames end
getnodenames(tree::AbstractTree{OneTree}, order::TraversalOrder = preorder) =
    _getnodenames(tree, order)
getnodenames(trees::AbstractTree{ManyTrees}, name,
             order::TraversalOrder = preorder) =
    getnodenames(gettree(trees, name), order)
getnodenames(trees::AbstractTree{ManyTrees}, order::TraversalOrder = preorder) =
    Dict(name => getnodenames(gettree(trees, name), order)
         for name in _gettreenames(trees))

"""
    getbranches(::AbstractTree)

Returns the vector of branches of a single tree, or a Dict of vectors of
branches for multiple trees.
"""
function getbranches end
getbranches(tree::AbstractTree{OneTree}) = _getbranches(tree)
getbranches(trees::AbstractTree{ManyTrees}, name) =
    getbranches(gettree(trees, name))
getbranches(trees::AbstractTree{ManyTrees}) =
    Dict(name => getbranches(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    getbranchnames(tree::AbstractTree)

Return a vector of branch names of a single tree, or a Dict of vectors of
branch names for multiple trees.
"""
function getbranchnames end
getbranchnames(tree::AbstractTree{OneTree}) = _getbranchnames(tree)
getbranchnames(trees::AbstractTree{ManyTrees}, name) =
    getbranchnames(gettree(trees, name))
getbranchnames(trees::AbstractTree{ManyTrees}) =
    Dict(name => getbranchnames(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    createbranch!(tree::AbstractTree, src, dst[, len::Number];
                  data)

Add a branch from `src` to `dst` on `tree` with optional length and
data. `source` and `destination` can be either nodes or nodenames.
"""
function createbranch! end
@traitfn function createbranch!(tree::T, src::N1, dst::N2, len = missing;
                                name = missing, data = missing) where
    {T <: AbstractTree{OneTree, <: Rooted}, N1, N2;
     !MatchNodeTypes{T, N1, N2}}
    hasnode(tree, src) ||
        error("Tree does not have an available source node - $src")
    hasnode(tree, dst) ||
        error("Tree does not have a destination node - $dst")
    sn = getnode(tree, src)
    dn = getnode(tree, dst)
    _hasoutboundspace(tree, sn) ||
        error(_getnodename(tree, sn) *
              " already has maximum number of outbound connections " *
              "($(_outdegree(tree, sn)))")
    _hasinboundspace(tree, dn) ||
        error("Tree does not have an available destination node called " *
              _getnodename(tree, dn))
    sn ≢ dn || error("Branch must connect different nodes")

    return ismissing(data) ?
        _createbranch!(tree, sn, dn, len, name) :
        _createbranch!(tree, sn, dn, len, name, data)
end
@traitfn function createbranch!(tree::T, src::N1, dst::N2, len = missing;
                                name = missing, data = missing) where
    {T <: AbstractTree{OneTree, <: Rooted}, N1, N2; MatchNodeTypes{T, N1, N2}}
    _hasnode(tree, src) ||
        error("Tree does not have an available source node - $src")
    _hasnode(tree, dst) ||
        error("Tree does not have a destination node - $dst")
    _hasoutboundspace(tree, src) ||
        error(_getnodename(tree, src) *
              " already has maximum number of outbound connections " *
              "($(_outdegree(tree, src)))")
    _hasinboundspace(tree, dst) ||
        error("Tree does not have an available destination node called " *
              _getnodename(tree, dst))
    src ≢ dst || error("Branch must connect different nodes")

    return ismissing(data) ?
        _createbranch!(tree, src, dst, len, name) :
        _createbranch!(tree, src, dst, len, name, data)
end
@traitfn function createbranch!(tree::T, src::N1, dst::N2, len = missing;
                                name = missing, data = missing) where
    {T <: AbstractTree{OneTree, Unrooted}, N1, N2;
     !MatchNodeTypes{T, N1, N2}}
    hasnode(tree, src) ||
        error("Tree does not have an available source node - $src")
    hasnode(tree, dst) ||
        error("Tree does not have a destination node - $dst")
    sn = getnode(tree, src)
    dn = getnode(tree, dst)
    _hasspace(tree, sn) ||
        error(_getnodename(tree, sn) *
              " already has maximum number of outbound connections " *
              "($(_outdegree(tree, sn)))")
    _hasspace(tree, dn) ||
        error("Tree does not have an available destination node called " *
              _getnodename(tree, dn))
    sn ≢ dn || error("Branch must connect different nodes")

    return ismissing(data) ?
        _createbranch!(tree, sn, dn, len, name) :
        _createbranch!(tree, sn, dn, len, name, data)
end
@traitfn function createbranch!(tree::T,
                                src::N1, dst::N2, len = missing;
                                name = missing, data = missing) where
    {T <: AbstractTree{OneTree, Unrooted}, N1, N2; MatchNodeTypes{T, N1, N2}}
    _hasnode(tree, src) ||
        error("Tree does not have an available source node - $src")
    _hasnode(tree, dst) ||
        error("Tree does not have a destination node - $dst")
    _hasspace(tree, src) ||
        error(_getnodename(tree, src) *
              " already has maximum number of outbound connections " *
              "($(_outdegree(tree, src)))")
    _hasspace(tree, dst) ||
        error("Tree does not have an available destination node called " *
              _getnodename(tree, dst))
    src ≢ dst || error("Branch must connect different nodes")

    return ismissing(data) ?
        _createbranch!(tree, src, dst, len, name) :
        _createbranch!(tree, src, dst, len, name, data)
end

"""
    deletebranch!(tree::AbstractTree, branch)
    deletebranch!(tree::AbstractTree, src, dst)

Delete the branch `branch` from `tree`, or branch connecting `src` node to
`dst` node.
"""
function deletebranch! end
@traitfn function deletebranch!(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) ||
        error("Tree does not have a branch $branch")
    return _deletebranch!(tree, getbranch(tree, branch))
end
@traitfn function deletebranch!(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) ||
        error("Tree does not have a branch called " *
              "$(_getbranchname(tree, branch))")
    return _deletebranch!(tree, branch)
end
@traitfn function deletebranch!(tree::T, src::N1, dst::N2) where
    {T <: AbstractTree{OneTree}, N1, N2; !MatchNodeTypes{T, N1, N2}}
    hasbranch(tree, src, dst) ||
        error("Tree does not have a branch between " *
              "$(_getnodename(tree, src)) and $(_getnodename(tree, dst))")
    return _deletebranch!(tree, getbranch(tree, src, dst))
end
@traitfn function deletebranch!(tree::T, src::N1, dst::N2) where
    {T <: AbstractTree{OneTree}, N1, N2; MatchNodeTypes{T, N1, N2}}
    hasbranch(tree, src, dst) ||
        error("Tree does not have a branch between " *
              "$(_getnodename(tree, src)) and $(_getnodename(tree, dst))")
    return _deletebranch!(tree, _getbranch(tree, src, dst))
end

"""
    createnode!(tree::AbstractTree[, nodename]; data)

Create a node on a tree with optional node info.
"""
function createnode!(tree::AbstractTree{OneTree, RT, NL, N, B},
                     nodename = missing; data = missing) where {RT, NL, N, B}
    return ismissing(data) ? _createnode!(tree, nodename) :
        _createnode!(tree, nodename, data)
end

"""
    createnodes!(tree::AbstractTree, count::Integer)
    createnodes!(tree::AbstractTree, nodenames)
    createnodes!(tree::AbstractTree, nodedict)

Add a number of nodes, a vector with given names, or a Dict with node names
and associated node info to a tree.
"""
function createnodes! end
createnodes!(tree::AbstractTree{OneTree}, count::Int) =
    [createnode!(tree) for count in Base.OneTo(count)]
function createnodes!(tree::AbstractTree{OneTree, RT, NL, N, B},
                      nodenames::AbstractVector{NL}) where {RT, NL, N, B}
    here = [name for name in nodenames if _hasnode(tree, name)]
    isempty(here) || error("Nodes $here already present in tree")
    return [_createnode!(tree, name) for name in nodenames]
end
@traitfn function createnodes!(tree::T, nodedata::NDS) where
    {RT, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B},
     ND, NDS <: Dict{NL, ND}; HoldsNodeData{T, ND}}
    here = [name for name in keys(nodedata) if _hasnode(tree, name)]
    isempty(here) || error("Nodes $here already present in tree")
    return [_createnode!(tree, info.first, info.second) for info in nodedata]
end

"""
    deletenode!(tree::AbstractTree, node)

Delete a node (or a name) from a tree
"""
function deletenode! end
@traitfn deletenode!(tree::T, node::NL) where
{NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    deletenode!(tree, getnode(tree, node))
@traitfn function deletenode!(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Tree does not have this node")
    return _deletenode!(tree, node)
end

"""
    hasnode(tree::AbstractTree, node)

Returns whether a tree has a given node (or node name) or not.
"""
hasnode(tree::AbstractTree{OneTree}, node) = _hasnode(tree, node)

"""
    getnode(tree::AbstractTree, nodename)

Returns a node from a tree.
"""
function getnode(tree::AbstractTree{OneTree}, node)
    hasnode(tree, node) || error("Node $node does not exist")
    return _getnode(tree, node)
end

"""
    getnodename(::AbstractTree, node)

Returns the node name associated with a node from a tree. For some
node types, it will be able to extract the node name without reference to
the tree.
"""
function getnodename end
@traitfn function getnodename(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _getnodename(tree, node)
end
@traitfn function getnodename(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getnodename(tree, node)
end

"""
    hasbranch(tree::AbstractTree, branch)
    hasbranch(tree::AbstractTree, source, dest)

Does `tree` have a branch `branch` or a branch from `source` to `dest`?
"""
function hasbranch end
@traitfn hasbranch(tree::T, branch::B) where
{T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}} =
    _hasbranch(tree, branch)
@traitfn hasbranch(tree::T, branch::B) where
{T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}} =
    _hasbranch(tree, _getbranch(tree, branch))
@traitfn hasbranch(tree::T, src::N1, dst::N2) where
{T <: AbstractTree{OneTree}, N1, N2; !MatchNodeTypes{T, N1, N2}} =
    hasnode(tree, src) && hasnode(tree, dst) &&
    _hasbranch(tree, _getnode(tree, src), _getnode(tree, dst))
@traitfn hasbranch(tree::T, src::N1, dst::N2) where
{T <: AbstractTree{OneTree}, N1, N2; MatchNodeTypes{T, N1, N2}} =
    _hasnode(tree, src) && _hasnode(tree, dst) && _hasbranch(tree, src, dst)

"""
    getbranch(tree::AbstractTree, branch)
    getbranch(tree::AbstractTree, source, dest)

Returns a branch from a tree by name or by source and destination node.
"""
function getbranch end
@traitfn function getbranch(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getbranch(tree, branch)
end
@traitfn function getbranch(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return branch
end
@traitfn getbranch(tree::T, src::N1, dst::N2) where
{T <: AbstractTree{OneTree}, N1, N2; !MatchNodeTypes{T, N1, N2}} =
    getbranch(tree, getnode(tree, src), getnode(tree, dst))
@traitfn getbranch(tree::T, src::N1, dst::N2) where
{T <: AbstractTree{OneTree}, N1, N2; MatchNodeTypes{T, N1, N2}} =
    hasnode(tree, src) && hasnode(tree, dst) && _getbranch(tree, src, dst)

"""
    getbranchname(::AbstractTree, branch)
    getbranchname(branch)

Returns the branch name associated with a branch from a tree. For some
branch types, it will be able to extract the branch name without reference to
the tree.
"""
function getbranchname end
@traitfn function getbranchname(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getbranchname(tree, branch)
end
@traitfn function getbranchname(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getbranchname(tree, branch)
end
@traitfn function getbranchname(tree::T, branch::Int) where
    {T <: AbstractTree{OneTree}; MatchBranchType{T, Int}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return branch
end

"""
    hasrootheight(tree::AbstractTree)

Does the tree have an explicit root height
"""
function hasrootheight end
hasrootheight(tree::AbstractTree{OneTree, <: Rooted}) = _hasrootheight(tree)
hasrootheight(tree::AbstractTree{OneTree, Unrooted}) = false
hasrootheight(trees::AbstractTree{ManyTrees}) =
    all(hasrootheight(tree) for tree in gettrees(trees))

"""
    getrootheight(tree::AbstractTree)

Get the tree's root height.
"""
getrootheight(tree::AbstractTree{OneTree, <: Rooted}) = _getrootheight(tree)
getrootheight(trees::AbstractTree{ManyTrees}, name) =
    getrootheight(gettree(trees, name))
getrootheight(trees::AbstractTree{ManyTrees}) =
    Dict(name => getrootheight(gettree(trees, name))
         for name in _gettreenames(trees))

"""
    setrootheight!(tree::AbstractTree, height)

Set the tree's root height.
"""
function setrootheight! end
function setrootheight!(tree::AbstractTree{OneTree, <: Rooted}, height)
    return _setrootheight!(tree, height)
end
function setrootheight!(trees::AbstractTree{ManyTrees},
                        heights::AbstractVector)
    for (tree, height) in zip(gettrees(tree), heights)
        setrootheight!(tree, height)
    end
end

# AbstractNode methods
"""
    isleaf(tree::AbstractTree, node)

Is the node (referenced by name or node object) a leaf of the tree?
"""
function isleaf end
@traitfn function isleaf(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _isleaf(tree, _getnode(tree, node))
end
@traitfn function isleaf(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _isleaf(tree, node)
end

"""
    isroot(tree::AbstractTree, node)

Is the node (referenced by name or node object) a root of the tree?
"""
function isroot end
@traitfn function isroot(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _isroot(tree, _getnode(tree, node))
end
@traitfn function isroot(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _isroot(tree, node)
end

"""
    isinternal(tree::AbstractTree, node)

Is the node (referenced by name or node object) internal to the tree (neither
root nor leaf)?
"""
function isinternal end
@traitfn function isinternal(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _isinternal(tree, _getnode(tree, node))
end
@traitfn function isinternal(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _isinternal(tree, node)
end

"""
    isunattached(tree::AbstractTree, node)

Is the node (referenced by name or node object) unattached (i.e. not connected
to other nodes)?
"""
function isunattached end
@traitfn function isunattached(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _isunattached(tree, _getnode(tree, node))
end
@traitfn function isunattached(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _isunattached(tree, node)
end

"""
    indegree(tree::AbstractTree, node)

Return in degree of node.
"""
function indegree end
@traitfn function indegree(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _indegree(tree, _getnode(tree, node))
end
@traitfn function indegree(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _indegree(tree, node)
end

"""
    outdegree(tree::AbstractTree, node)

Return out degree of node.
"""
function outdegree end
@traitfn function outdegree(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _outdegree(tree, _getnode(tree, node))
end
@traitfn function outdegree(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _outdegree(tree, node)
end

"""
    degree(tree::AbstractTree, node)

Return the degree of a node including all connections.
"""
function degree end
@traitfn function degree(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _degree(tree, _getnode(tree, node))
end
@traitfn function degree(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _degree(tree, node)
end

"""
    hasoutboundspace(tree::AbstractTree, node)

Does the node have space for an[other] outbound connection?
"""
function hasoutboundspace end
@traitfn function hasoutboundspace(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _hasoutboundspace(tree, _getnode(tree, node))
end
@traitfn function hasoutboundspace(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _hasoutboundspace(tree, node)
end

"""
    hasinbound(tree::AbstractTree, node)

Does the node have an inbound connection?
"""
function hasinbound end
@traitfn function hasinbound(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _hasinbound(tree, _getnode(tree, node))
end
@traitfn function hasinbound(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _hasinbound(tree, node)
end

"""
    hasinboundspace(tree::AbstractTree, node)

Does the node have space for an inbound connection?
"""
function hasinboundspace end
@traitfn function hasinboundspace(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _hasinboundspace(tree, _getnode(tree, node))
end
@traitfn function hasinboundspace(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _hasinboundspace(tree, node)
end

"""
    getinbound(tree::AbstractTree, node)

return the inbound branch to this node (returns name for node name, branch for
node).
"""
function getinbound end
@traitfn function getinbound(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _getinbound(tree, _getnode(tree, node))
end
@traitfn function getinbound(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getinbound(tree, node)
end

"""
    getparent(tree::AbstractTree, node)

Return [the name of] the parent node for this node [name]. Second method may
not be implemented for some node types.
"""
function getparent end
@traitfn function getparent(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _getnodename(tree, _getparent(tree, _getnode(tree, node)))
end
@traitfn function getparent(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getparent(tree, node)
end

"""
    getoutbounds(tree::AbstractTree, nodename)

Return the names of the outbound branches from this node.
"""
function getoutbounds end
@traitfn function getoutbounds(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _getoutbounds(tree, _getnode(tree, node))
end
@traitfn function getoutbounds(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getoutbounds(tree, node)
end

"""
    getchildren(tree::AbstractTree, node)

Return the [name(s) of] the child node(s) for this node [name].
"""
function getchildren end
@traitfn function getchildren(tree::T, node::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _getnodename.(tree, _getchildren(tree, _getnode(tree, node)))
end
@traitfn function getchildren(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getchildren(tree, node)
end

"""
    getconnections(tree::AbstractTree, nodee)

Returns all of the branches connected to a node.
"""
function getconnections end
@traitfn function getconnections(tree::T, node::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _getconnections(tree, _getnode(tree, node))
end
@traitfn function getconnections(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getconnections(tree, node)
end

"""
    getsiblings(tree::AbstractTree, node)

Returns all of the siblings of a node. Must be implemented for any unrooted
AbstractNode subtype, can be inferred from _getparent and _getchildren for
a rooted node.
"""
function getsiblings end
@traitfn getsiblings(tree::T, node::NL) where
{NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    _getsiblings(tree, getnode(tree, node))
@traitfn function getsiblings(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getsiblings(tree, node)
end

"""
    hasheight(tree::AbstractTree, node)

Does the node have a height defined?
"""
function hasheight end
@traitfn hasheight(tree::T, node::NL) where
{NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}} =
    hasheight(tree,  getnode(tree, node))
@traitfn function hasheight(tree::AbstractTree{OneTree, <: Rooted},
                            node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    hasnode(tree, node) || error("Node $node does not exist")
    return _hasheight(tree, node) ||
        (_hasrootheight(tree) &&
         mapreduce(b -> haslength(tree, b), &, branchhistory(tree, node);
                   init = _hasrootheight(tree)))
end

"""
    getheight(tree::AbstractTree, node)

Return the height of the node.
"""
function getheight end
@traitfn getheight(tree::T, node::NL) where
{NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}} =
    getheight(tree,  getnode(tree, node))
@traitfn function getheight(tree::T, node::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _hasheight(tree, node) ? _getheight(tree, node) :
        mapreduce(b -> getlength(tree, b), +, branchhistory(tree, node);
                  init = hasrootheight(tree) ?
                  getrootheight(tree) :
                  0.0 * upreferred(branchdims(T)))
end

"""
    setheight!(tree::AbstractTree, nodename, height)

Set the height of the node.
"""
function setheight! end
@traitfn setheight!(tree::T, node::NL, height::Number) where
{NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}} =
    _setheight!(tree, getnode(tree, node), height)
@traitfn function setheight!(tree::T, node::N, height::Number) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _setheight!(tree, node, height)
end

# AbstractBranch methods
"""
    src(tree::AbstractTree, branch)

Return the source node for this branch.
"""
function src end
@traitfn src(tree::T, branch::B) where
{T <: AbstractTree{OneTree, <: Rooted}, B; !MatchBranchType{T, B}} =
    _src(tree, getbranch(tree, branch))
@traitfn function src(tree::T, branch::B) where
    {T <: AbstractTree{OneTree, <: Rooted}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _src(tree, branch)
end

"""
    dst(tree::AbstractTree, branch)

Return the destination node for this branch.
"""
function dst end
@traitfn function dst(tree::T, branch::B) where
    {T <: AbstractTree{OneTree, <: Rooted}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _dst(tree, getbranch(tree, branch))
end
@traitfn function dst(tree::T, branch::B) where
    {T <: AbstractTree{OneTree, <: Rooted}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _dst(tree, branch)
end

"""
    conns(tree::AbstractTree, branch)

Return the nodes connected to `branch`.
"""
function conns end
@traitfn function conns(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _conns(tree, getbranch(tree, branch))
end
@traitfn function conns(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _conns(tree, branch)
end

"""
    conn(tree::AbstractTree, branch, exclude)

Return the other node connected to `branch` that is not `exclude`.
"""
function conn end
@traitfn function conn(tree::T, branch::B, exclude::N) where
    {T <: AbstractTree{OneTree}, B, N; !MatchBranchNodeType{T, B, N}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    hasnode(tree, exclude) || error("Node $exclude does not exist")
    return _conn(tree, _getbranch(tree, branch), _getnode(tree, exclude))
end
@traitfn function conn(tree::T, branch::B, exclude::N) where
    {T <: AbstractTree{OneTree}, B, N; MatchBranchNodeType{T, B, N}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    _hasnode(tree, exclude) || error("Node $exclude does not exist")
    return _conn(tree, branch, exclude)
end

"""
    getlength(tree::AbstractTree, branch)

Return the length of this branch.
"""
function getlength end
@traitfn function getlength(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B, N; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getlength(tree, _getbranch(tree, branch))
end
@traitfn function getlength(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B, N; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getlength(tree, branch)
end

"""
    getleafnames(::AbstractTree[, ::TraversalOrder])

Retrieve the leaf names from the tree (in some specific order).
"""
getleafnames(tree::AbstractTree, order::TraversalOrder = preorder) =
    _getleafnames(tree, order)

"""
    getleaves(::AbstractTree[, ::TraversalOrder])

Retrieve the leaves from the tree.
"""
getleaves(tree::AbstractTree{OneTree}, order::TraversalOrder = preorder) =
    _getleaves(tree, order)
getleaves(trees::AbstractTree{ManyTrees}, name,
          order::TraversalOrder = preorder) =
    _getleaves(gettree(trees, name), order)
getleaves(trees::AbstractTree{ManyTrees}, order::TraversalOrder = preorder) =
    Dict(name => _getleaves(gettree(trees, name), order)
         for name in _gettreenames(trees))

"""
    getleafinfo(::AbstractTree[, label])

retrieve the leaf info for a leaf of the tree.
"""
function getleafinfo end
getleafinfo(tree::AbstractTree) = _getleafinfo(tree)
@traitfn getleafinfo(tree::T, leaf::NL) where
{RT, NL, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    _getleafinfo(tree, getnode(tree, leaf))
@traitfn function getleafinfo(tree::T, leaf::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, leaf) || error("Node $node does not exist")
    return _getleafinfo(tree, leaf)
end

"""
    setleafinfo!(::AbstractTree, table)

Set the leaf info for the leaves of the tree.
"""
function setleafinfo!(tree::AbstractTree, table)
    _setleafinfo!(tree, table)
    validate!(tree) || error("Leaf info not consistent with tree")
end

"""
    getnodedata(::AbstractTree, node)

retrieve the node data for a node of the tree.
"""
function getnodedata end
@traitfn getnodedata(tree::T, node::NL) where
{NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    _getnodedata(tree, getnode(tree, node))
@traitfn function getnodedata(tree::T, node::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getnodedata(tree, node)
end
@traitfn getnodedata(tree::T, node::NL, label) where
{NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    _getnodedata(tree, getnode(tree, node), label)
@traitfn function getnodedata(tree::T, node::N, label) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getnodedata(tree, node, label)
end

"""
    setnodedata!(::AbstractTree, node, label, value)
    setnodedata!(::AbstractTree, node, data)

Set the node data for a node of the tree.
"""
function setnodedata! end
@traitfn setnodedata!(tree::T, node::NL, label, value) where
{NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    _setnodedata!(tree, getnode(tree, node), label, value)
@traitfn function setnodedata!(tree::T, node::N, label, value) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _setnodedata!(tree, node, label, value)
end
@traitfn setnodedata!(tree::T, node::NL, data) where
{NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}} =
    _setnodedata!(tree, getnode(tree, node), data)
@traitfn function setnodedata!(tree::T, node::N, data) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, node) || error("Node $node does not exist")
    return _setnodedata!(tree, node, data)
end

"""
    getbranchdata(::AbstractTree, label)

retrieve the branch data for a leaf of the tree.
"""
function getbranchdata end
@traitfn function getbranchdata(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getbranchdata(tree, _getbranch(tree, branch))
end
@traitfn function getbranchdata(tree::T, branch::B) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getbranchdata(tree, branch)
end

"""
    setbranchdata!(::AbstractTree, branch, label, value)
    setbranchdata!(::AbstractTree, branch, data)

Set the branch data for a branch of the tree.
"""
function setbranchdata! end
@traitfn function setbranchdata!(tree::T, branch::B, label, value) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _setbranchdata!(tree, _getbranch(tree, branch), label, value)
end
@traitfn function setbranchdata!(tree::T, branch::B, label, value) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _setbranchdata!(tree, branch, label, value)
end
@traitfn function setbranchdata!(tree::T, branch::B, data) where
    {T <: AbstractTree{OneTree}, B; !MatchBranchType{T, B}}
    hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _setbranchdata!(tree, _getbranch(tree, branch), data)
end
@traitfn function setbranchdata!(tree::T, branch::B, data) where
    {T <: AbstractTree{OneTree}, B; MatchBranchType{T, B}}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _setbranchdata!(tree, branch, data)
end

"""
    validate!(tree::AbstractTree)

Validate the tree by making sure that it is connected up correctly.
"""
function validate!(tree::T) where
    {TT, RT, NL, N, B, T <: AbstractTree{TT, RT, NL, N, B}}
    _resetleaves!(tree)
    nodes = _getnodes(tree)
    nodenames = _getnodenames(tree)
    branches = _getbranches(tree)
    branchnames = _getbranchnames(tree)
    if !isempty(nodes) || !isempty(branches)
        # We need to validate the connections
        if Set(_getinbound(tree, node) for node in nodes
               if _hasinbound(tree, node)) != Set(branches)
            @warn "Inbound branches must exactly match Branch labels"
            return false
        end

        if Set(mapreduce(node -> _getoutbounds(tree, node), append!,
                         nodes; init = B[])) != Set(branches)
            @warn "Node outbound branches must exactly match Branch labels"
            return false
        end

        if !(Set(_getnodename(tree, _src(tree, branch))
                 for branch in branches) ⊆ Set(nodenames))
            @warn "Branch sources must be nodes"
            return false
        end

        if !(Set(_getnodename(tree, _dst(tree, branch))
                for branch in branches) ⊆ Set(nodenames))
            @warn "Branch destinations must be nodes"
            return false
        end
    end

    return _validate!(tree)
end

"""
    traversal(::AbstractTree, ::TraversalOrder)
    traversal(::AbstractTree, ::TraversalOrder, init)

Return an iterable object for a tree containing nodes in given order -
preorder, inorder, postorder or breadthfirst - optionally starting from init.
"""
function traversal end
traversal(tree::AbstractTree{OneTree}, order::TraversalOrder = preorder) =
    _traversal(tree, order)
@traitfn function traversal(tree::T, order::TraversalOrder, init::NL) where
    {NL, RT, T <: AbstractTree{OneTree, RT, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, init) || error("Node $init does not exist")
    return _getnodename.(tree, _traversal(tree, order, [getnode(tree, init)]))
end
@traitfn function traversal(tree::T, order::TraversalOrder, init::N) where
    {T <: AbstractTree{OneTree}, N; MatchNodeType{T, N}}
    _hasnode(tree, init) || error("Node $init does not exist")
    return _traversal(tree, order, [init])
end

"""
    getancestors(tree::AbstractTree, node)

Return the name of all of the nodes that are ancestral to this node.
"""
function getancestors end
@traitfn function getancestors(tree::T, init::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, init) || error("Node $init does not exist")
    return getnodename.(tree, _treehistory(tree, getnode(tree, init))[2][2:end])
end
@traitfn function getancestors(tree::T, init::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    hasnode(tree, init) || error("Node $init does not exist")
    return _treehistory(tree, init)[2][2:end]
end

"""
    getdescendants(tree::AbstractTree, node)

Return the names of all of the nodes that descend from this node.
"""
function getdescendants end
@traitfn function getdescendants(tree::T, init::NL) where
    {NL, T <: AbstractTree{OneTree, <: Rooted, NL}; !MatchNodeType{T, NL}}
    hasnode(tree, init) || error("Node $init does not exist")
    return getnodename.(tree, _treefuture(tree, getnode(tree, init))[2][2:end])
end
@traitfn function getdescendants(tree::T, init::N) where
    {T <: AbstractTree{OneTree, <: Rooted}, N; MatchNodeType{T, N}}
    hasnode(tree, init) || error("Node $init does not exist")
    return _treefuture(tree, init)[2][2:end]
end


text:
module LearnVanishingIdeal
#TODO How to make sure that if z in
import HomotopyContinuation: @polyvar
import LinearAlgebra: norm, rank

export approximateVanishingIdeal,
       approximateVanishingIdeal_maxDegree,
       leastSquaresListOfEquations,
       leastSquaresListOfEquations_quick,
	   affineVeronese,
	   projVeronese,
	   addNoise,
	   calculateMeanDistanceToVariety

include("auxiliaryFunctions.jl")
using .auxiliaryFunctions

"""
  @input The function approximateVanishingIdeal takes a list of points,
  a list of Degrees of the generators of the vanishing ideal and a
  boolean quick as input that indicates, which algorithm to use.

  @return a coefficient vector and the error (Sampson Distance or Mean Distance)
  of the best-approximating polynomial
"""
function approximateVanishingIdeal(points, listOfDegrees; quick=false, affine=true)
	degreeList = findEqListOfDegrees(listOfDegrees)
	return( leastSquaresListOfEquations( points, degreeList, affine; quick=quick))
end

"""
  @input a maximum occuring degree n among the generators of the vanishing ideal,
  a list of points, the amount of generators (codimension of the variety, as complete
  intersection is assumed) and a boolean quick.

  @return a coefficient vector and the error (Sampson Distance or Mean Distance)
  of the best-approximating polynomial
"""
function approximateVanishingIdeal_maxDegree(n, points, numEq, quick=false, affine=true)
	listOfDegrees = [n for i in 1:numEq]
	return(approximateVanishingIdeal(points, listOfDegrees, quick, affine))
end

"""
  This method is a more accurate variant of the method
  leastSquaresListOfEquations_quick.
  It uses the Mean Distance (HomotopyContinuation) for determining the
  best-approximating polynomial. Also, it uses a threshold of tau=2.o and
  400 gradient descent steps.
"""
function leastSquaresListOfEquations(data, listOfDegrees, affine; TOL = 1e-8, quick=false)
	time1 = time()
	if quick == true
		maxiter, epochs, threshold = 200, 2, 1.5
	else
		maxiter, epochs, threshold = 400, 3, 3
	end
	points = affine ? data : [vcat(point,[1]) for point in data]
	@polyvar var[1:length(points[1])]
	startValuesEigen, startValuesVander, outputValues, err = [], [], [], Inf
	numEq, n = sum([entry[2] for entry in listOfDegrees]), maximum([entry[1] for entry in listOfDegrees])

	for entry in listOfDegrees
		EigenValueStart, EigenValueVander = comparisonOfMethods(entry[1], points, entry[2], threshold)
		append!(startValuesEigen, [[entry[2], EigenValueStart]])
		append!(startValuesVander, [[entry[2], EigenValueVander]])
	end

	startValueCombinations = quick ? [makeCombinations(startValuesEigen)] : [makeCombinations(startValuesEigen), makeCombinations(startValuesVander)]

	@polyvar w[1:binomial(n+length(points[1])-1,n),1:numEq]
	veronese = affineVeronese(n,var[1:length(data[1])])
	println("The search space has dimension: ", binomial(n+length(data[1]),n))

	for combinations in startValueCombinations
		for combination in combinations
			combination, zeroEntries = fillUpWithZeros(combination, n, numEq, length(var))
			intermediateValues = [[round(co/norm(comb), digits=14) for co in comb] for comb in combination]
			cl, intermediateValues = cleanUp(intermediateValues, listOfDegrees, length(var))
			cl || continue # ideal generators are not minimal
			result = [vector'*veronese for vector in intermediateValues]
			currentError = calculateMeanDistanceToVariety(points, result, var)
			currentError!="dim>0" || continue # the ideal has a higher dimension than anticipated.
			if currentError < err
				println("Ansatz without iterations takes the cake! Error: ", currentError)
				err = currentError
				outputValues = [comb for comb in intermediateValues]
			end
			currentError>=TOL || break

			for i in 1:epochs
				saverArray = sampsonDistance(points, numEq, n, w, intermediateValues)
				intermediateValues, _ = weightedGradientDescent(points, n, w, [value for value in intermediateValues], numEq, maxiter, saverArray, zeroEntries)
				result = [vector'*veronese for vector in intermediateValues]
				currentError = calculateMeanDistanceToVariety(points, result, var)
				currentError!="dim>0" || break
				if currentError < err
					println("Ansatz with ", i,  " iterations takes the cake! Error: ",currentError)
					err = currentError
					outputValues = [vector for vector in intermediateValues]
				end
			end
		end
	end
	time2 = time()
	println("The algorithm took: ",time2-time1,"s")
	_, outputValues = cleanUp(outputValues, listOfDegrees, length(var))
	return([[round(val/norm(value), digits=6) for val in value] for value in outputValues], err)
end

#=
"""
  This method learns polynomials for a given degree of the variety and
  given codimension by iterating through all possible combinations of
  degrees that make up the degree of the variety.
"""
function approximateVanishingIdealDegreeVariety(points, degVariety::Int, numEq::Int, quick=false)
	@polyvar var[1:length(points[1])]
	listsOfDegrees = findAllCompositions(numEq, degVariety)
	error = Inf
	finaloutput = []
	finalList = []
	for list in listsOfDegrees
		indicator = 0
		for entry in list
			if entry >= 6
				indicator = 1
			end
		end
		if indicator == 1
			continue
		end
		n = maximum(list)
		try
			output, currentError = approximateVanishingIdeal(points, list, quick)
			currentError = calculateMeanDistanceToVariety(points,[value'*affineVeronese(maximum(list),var) for value in output],var)
			if(currentError < error)
				error = currentError
				finaloutput = [vector for vector in output]
				finalList = list
			end
		catch e
			println("Fehler gefunden!")
			println(e)
		end
	end
	return(finaloutput, error,finalList)
end

"""
  This method is a statistical learning ansatz. it iterates through all
  degrees up to a given terminationDegree and calculates equations in this degree
"""
function inferringVanishingIdeal(points, terminationDegree, numEq)
	@polyvar var[1:length(points[1])]
	error = Inf
	outputValues = []
	errorsInDegrees = []
	finalDegree = 0
	for deg in 1:terminationDegree
		values, currentError, finalList = approximateVanishingIdealDegreeVariety(points, deg, numEq, true)
		println("Error of the variety with degree ",deg," has error: ",currentError)
		if(currentError < error)
			error = currentError
			outputValues = [vector for vector in values]
			finalDegree = deg
		end
		append!(errorsInDegrees,[currentError])
		display(errorsInDegrees)
	end
	return(outputValues, error, deg)
end
=#

end


text:
using PowerSimulationsDynamics
using Sundials

"""
Case 24:
This case study a 15-state grid following inverter against an infinite bus located at bus 1, with the inverter located at bus 2.
The perturbation increase the reference power (analogy for mechanical power) from 0.5 to 0.7.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test24.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

#time span
tspan = (0.0, 2.0);

#Define Fault using Callbacks
Pref_change = ControlReferenceChange(1.0, case_inv, PSID.P_ref_index, 0.7)

@testset "Test 24 Grid Following Inverter ImplicitModel" begin
    path = (joinpath(pwd(), "test-24"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(
            ImplicitModel,
            omib_sys, # system
            path,
            tspan,
            Pref_change,
        )

        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Solve problem in equilibrium
        execute!(sim, Sundials.IDA(), dtmax = 0.001)

        #Obtain frequency data
        series = get_state_series(sim, ("generator-102-1", :p_oc))
        t = series[1]
        p = series[2]

        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test24_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)
        @test LinearAlgebra.norm(eigs - test24_eigvals) < 1e-3
        @test sim.solution.retcode == :Success

        power = PSID.get_activepower_series(sim, "generator-102-1")
        rpower = PSID.get_reactivepower_series(sim, "generator-102-1")
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 24 Grid Following Inverter ImplicitModel" begin
    path = (joinpath(pwd(), "test-24"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            omib_sys, # system
            path,
            tspan,
            Pref_change,
        )

        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Solve problem in equilibrium
        execute!(sim, Rodas5(), dtmax = 0.001)

        #Obtain frequency data
        series = get_state_series(sim, ("generator-102-1", :p_oc))
        t = series[1]
        p = series[2]

        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test24_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)
        @test LinearAlgebra.norm(eigs - test24_eigvals) < 1e-3
        @test sim.solution.retcode == :Success

        power = PSID.get_activepower_series(sim, "generator-102-1")
        rpower = PSID.get_reactivepower_series(sim, "generator-102-1")
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
@testset "Unramified extension" begin
    Qx,x = FlintQQ["x"]
    f = Qx([1, 8, -40, -46, 110, 71, -113, -43, 54, 11, -12, -1, 1])
    L = number_field(f)[1]
    P = prime_decomposition(maximal_order(L),7)[1][1]
    lp, mp = Hecke.generic_completion(L,P)
    Qy,y = PolynomialRing(lp,"y")
    f, mf = ResidueField(lp)
    N = Hecke.unramified_extension(y^3+preimage(mf,(gen(f)))+4)[1]
    F, mF = ResidueField(N)
    @test order(F) == 7^6
    G, mG = automorphism_group(N)
    @test order(G) == 3
    @test mG(G[1]^2) == mG(G[1])^2
    b = rand(f)
    @test norm(Hecke.norm_equation(F, b)) == b
    for i = 1:3
        c = 1+uniformizer(lp)^i
        chk = norm(Hecke.norm_equation_unramified(N, c))*c^-1 -1
        @test iszero(chk) || Int(ee*valuation(chk)) >= precision(c)
    end
    n = Int(ceil(absolute_ramification_index(lp)//(Int(prime(lp))-1)))+1
    l = base_field(lp)
    ee = absolute_ramification_index(l)
    for i = n:n+2
        c = 1+uniformizer(l)^i
        chk = norm(Hecke.norm_equation(lp, c))*c^-1 -1
        @test iszero(chk) || Int(ee*valuation(chk)) >= precision(c)
    end

end


text:
using TypographicUnits
using Unitful.DefaultSymbols
using Test

const ABSOLUTE = pt, bp, pc, dd, cc, nd, nc, sp
const RELATIVE = em, ex, px, en

for x in ABSOLUTE, y in ABSOLUTE
    (3.4x + 2y//9)/4
end
for x in RELATIVE
    9*(9x - 3.14x)
end

for x in ABSOLUTE, y in RELATIVE
    y == en && continue
    5y//15 - 1x + 2.4y
end

const io = IOBuffer()
result() = String(take!(io))

show(io, (2pt + 4em + 6ex + 8px)/2)
@test result() == "1.0 pt + 2.0 em + 3.0 ex + 4.0 px"

show(io, 3.1pt - 4.3ex)
@test result() == "3.1 pt - 4.3 ex"

show(io, 3px + 4.0em)
@test result() == "4.0 em + 3.0 px"

show(io, 1cm + 1pt)
@test result() == "7481//254 pt"

show(IOContext(io, :compact => true), float(3cm + 4pt + 4em))
@test result() == "89.3583 pt+4.0 em"


text:
# BoundaryEntropyDissipation functional

#------------------------------------------------------------------------------
# Boundary dissipation
function calcBndryFunctional(mesh::AbstractDGMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::BoundaryEntropyDiss{Topt}) where {Tmsh, Tsol, Tres, Topt}

  val = zero(Topt)
  for i=1:mesh.numBC
    if i in func.bcnums
      functor_i = mesh.bndry_funcs[i]
      start_index = mesh.bndry_offsets[i]
      end_index = mesh.bndry_offsets[i+1] - 1
      idx_range = start_index:end_index
      bndry_facenums_i = sview(mesh.bndryfaces, idx_range)

      val += _calcBndryEntropyDiss(mesh, sbp, eqn, opts, func, functor_i, idx_range, bndry_facenums_i)
    end
  end

  return MPI.Allreduce(val, MPI.SUM, eqn.comm)
end

function calcBndryFunctional(mesh::AbstractDGMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::NegBoundaryEntropyDiss{Topt}) where {Tmsh, Tsol, Tres, Topt}

  return -calcBndryFunctional(mesh, sbp, eqn, opts, func.func)
end

function calcBndryFunctional_revm(mesh::AbstractDGMesh{Tmsh},
     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
     opts, func::BoundaryEntropyDiss, _val_bar::Number=1,
     ) where {Tmsh, Tsol, Tres, Tdim}

  for i=1:mesh.numBC
    if i in func.bcnums
      functor_i = mesh.bndry_funcs_revm[i]
      start_index = mesh.bndry_offsets[i]
      end_index = mesh.bndry_offsets[i+1] - 1
      idx_range = start_index:end_index
      bndry_facenums_i = sview(mesh.bndryfaces, idx_range)

      _calcBndryEntropyDiss_revm(mesh, sbp, eqn, opts, func, functor_i, idx_range, bndry_facenums_i, _val_bar)
    end
  end

  return nothing
end

function calcBndryFunctional_revm(mesh::AbstractDGMesh{Tmsh},
     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
     opts, func::NegBoundaryEntropyDiss, _val_bar::Number=1,
     ) where {Tmsh, Tsol, Tres, Tdim}

  calcBndryFunctional_revm(mesh, sbp, eqn, opts, func.func, -_val_bar)
end

function calcFunctionalDeriv(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol}, opts,
                           func::BoundaryEntropyDiss,
                           func_deriv_arr::Abstract3DArray,
                           _val_bar::Number=1) where {Tmsh, Tsol}

  for i=1:mesh.numBC
    if i in func.bcnums
      functor_i = mesh.bndry_funcs_revq[i]
      start_index = mesh.bndry_offsets[i]
      end_index = mesh.bndry_offsets[i+1] - 1
      idx_range = start_index:end_index
      bndry_facenums_i = sview(mesh.bndryfaces, idx_range)

      _calcBndryEntropyDiss_revq(mesh, sbp, eqn, opts, func, functor_i, idx_range, bndry_facenums_i, func_deriv_arr, _val_bar)
    end
  end

  return nothing
end

function calcFunctionalDeriv(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol}, opts,
                           func::NegBoundaryEntropyDiss,
                           func_deriv_arr::Abstract3DArray,
                           _val_bar::Number=1) where {Tmsh, Tsol}

  calcFunctionalDeriv(mesh, sbp, eqn, opts, func.func, func_deriv_arr, -_val_bar)
end

function _calcBndryEntropyDiss(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::BoundaryEntropyDiss{Topt},
            functor::BCType,
            global_nums::UnitRange, faces::AbstractArray{Boundary}) where {Tmsh, Tsol, Tres, Topt}

  qface = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  w_face_j = zeros(Tsol, mesh.numDofPerNode)
  q_bndry_j = zeros(Tsol, mesh.numDofPerNode)
  w_bndry_j = zeros(Tsol, mesh.numDofPerNode)
  delta_w = zeros(Tsol, mesh.numDofPerNode)
  q_avg = zeros(Tsol, mesh.numDofPerNode)
  aux_vars = Tres[0]
  kernel = eqn.params.entropy_lf_kernel
  flux = zeros(Tres, mesh.numDofPerNode)

  val = zero(Tres)
  for i=1:length(faces)  # loop over faces with this BC
    bndry_i = faces[i]
    global_facenum = global_nums[i]
    q_el = sview(eqn.q, :, :, bndry_i.element)
    boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, q_el, qface)

    val_i = zero(Tres)
    for j=1:mesh.numNodesPerFace

      # get components
      # convert to conservative variables if needed
      q_j = sview(qface, :, j)
      aux_vars[1] = calcPressure(eqn.params, q_j)
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      nrm_xy = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
      bndry_node = BoundaryNode(bndry_i, i, j)

      getDirichletState(functor, eqn.params, q_j, aux_vars, coords, nrm_xy,
                        q_bndry_j)

      # compute Lambda(delta_w)
      convertToIR(eqn.params, q_j, w_face_j)
      convertToIR(eqn.params, q_bndry_j, w_bndry_j)
      for k=1:mesh.numDofPerNode
        q_avg[k]   = 0.5*(q_j[k] + q_bndry_j[k])
        delta_w[k] = w_face_j[k] - w_bndry_j[k]
      end

      applyEntropyKernel(kernel, eqn.params, q_avg, delta_w, nrm_xy, flux)

      # compute delta_w*B*Lambda(delta_w)
      for k=1:mesh.numDofPerNode
        val_i += delta_w[k]*mesh.sbpface.wface[j]*flux[k]
      end
    end  # end loop j
    val += val_i
  end  # end loop i

  return -val
end

function _calcBndryEntropyDiss_revm(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::BoundaryEntropyDiss{Topt},
            functor::BCType_revm,
            global_nums::UnitRange, faces::AbstractArray{Boundary},
            val_bar::Number) where {Tmsh, Tsol, Tres, Topt}

  qface = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  w_face_j = zeros(Tsol, mesh.numDofPerNode)
  q_bndry_j = zeros(Tsol, mesh.numDofPerNode)
  w_bndry_j = zeros(Tsol, mesh.numDofPerNode)
  delta_w = zeros(Tsol, mesh.numDofPerNode)
  q_avg = zeros(Tsol, mesh.numDofPerNode)
  aux_vars = Tres[0]
  kernel = eqn.params.entropy_lf_kernel
  flux = zeros(Tres, mesh.numDofPerNode)

  q_avg_bar = zeros(Tres, mesh.numDofPerNode)
  delta_w_bar = zeros(Tres, mesh.numDofPerNode)
  flux_bar = zeros(Tres, mesh.numDofPerNode)
  w_face_j_bar = zeros(Tres, mesh.numDofPerNode)
  w_bndry_j_bar = zeros(Tres, mesh.numDofPerNode)
  q_bndry_j_bar = zeros(Tres, mesh.numDofPerNode)
  qface_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  A0inv = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode)

  val_bar = -val_bar

  for i=1:length(faces)  # loop over faces with this BC
    bndry_i = faces[i]
    global_facenum = global_nums[i]
    q_el = sview(eqn.q, :, :, bndry_i.element)
    boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, q_el, qface)

    for j=1:mesh.numNodesPerFace

      # get components
      # convert to conservative variables if needed
      q_j = sview(qface, :, j)
      q_j_bar = sview(qface_bar, :, j)
      aux_vars[1] = calcPressure(eqn.params, q_j)
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      coords_bar = sview(mesh.coords_bndry_bar, :, j, global_facenum)
      nrm_xy = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
      nrm_xy_bar = sview(mesh.nrm_bndry_bar, :, j, global_facenum)
      bndry_node = BoundaryNode(bndry_i, i, j)

      getDirichletState(functor, eqn.params, q_j, aux_vars, coords,
                        nrm_xy, q_bndry_j, bndry_node)

      # compute Lambda(delta_w)
      convertToIR(eqn.params, q_j, w_face_j)
      convertToIR(eqn.params, q_bndry_j, w_bndry_j)
      for k=1:mesh.numDofPerNode
        q_avg[k]   = 0.5*(q_j[k] + q_bndry_j[k])
        delta_w[k] = w_face_j[k] - w_bndry_j[k]
      end

      applyEntropyKernel(kernel, eqn.params, q_avg, delta_w, nrm_xy, flux)

      ## compute delta_w*B*Lambda(delta_w)
      #for k=1:mesh.numDofPerNode
      #  val += delta_w[k]*mesh.sbpface.wface[j]*flux[k]
      #end

      # reverse sweep
      fill!(delta_w_bar, 0); fill!(flux_bar, 0); fill!(q_avg_bar, 0)
      fill!(w_face_j_bar, 0); fill!(w_bndry_j_bar, 0); fill!(q_bndry_j_bar, 0)
      fill!(q_j_bar, 0)
      for k=1:mesh.numDofPerNode
        delta_w_bar[k] += mesh.sbpface.wface[j]*flux[k]*val_bar
        flux_bar[k]    += delta_w[k]*mesh.sbpface.wface[j]*val_bar
      end

      applyEntropyKernel_revm(kernel, eqn.params, q_avg, delta_w, nrm_xy,
                              nrm_xy_bar, flux, flux_bar)
      applyEntropyKernel_revq(kernel, eqn.params, q_avg, q_avg_bar, delta_w,
                              delta_w_bar, nrm_xy, flux, flux_bar)

      for k=1:mesh.numDofPerNode
        #w_face_j_bar[k] += delta_w_bar[k]
        w_bndry_j_bar[k] -= delta_w_bar[k]

        #q_j_bar[k]      += 0.5*q_avg_bar[k]
        q_bndry_j_bar[k] += 0.5*q_avg_bar[k]
      end

      # reverse of convertToIR
      #getIRA0inv(eqn.params, q_j, A0inv)
      #smallmatTvec_kernel!(A0inv, w_face_j_bar, q_j_bar, 1, 1)
      getIRA0inv(eqn.params, q_bndry_j, A0inv)
      smallmatTvec_kernel!(A0inv, w_bndry_j_bar, q_bndry_j_bar, 1, 1)

      getDirichletState_revm(functor, eqn.params, q_j, aux_vars,
                             coords, coords_bar, nrm_xy, nrm_xy_bar,
                             q_bndry_j_bar, bndry_node)

    end  # end loop j
  end  # end loop i

  return nothing
end

function _calcBndryEntropyDiss_revq(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::BoundaryEntropyDiss{Topt},
            functor::BCType_revq,
            global_nums::UnitRange, faces::AbstractArray{Boundary},
            func_deriv_arr::Abstract3DArray,
            val_bar::Number) where {Tmsh, Tsol, Tres, Topt}

  qface = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  w_face_j = zeros(Tsol, mesh.numDofPerNode)
  q_bndry_j = zeros(Tsol, mesh.numDofPerNode)
  w_bndry_j = zeros(Tsol, mesh.numDofPerNode)
  delta_w = zeros(Tsol, mesh.numDofPerNode)
  q_avg = zeros(Tsol, mesh.numDofPerNode)
  aux_vars = Tres[0]
  kernel = eqn.params.entropy_lf_kernel
  flux = zeros(Tres, mesh.numDofPerNode)

  q_avg_bar = zeros(Tres, mesh.numDofPerNode)
  delta_w_bar = zeros(Tres, mesh.numDofPerNode)
  flux_bar = zeros(Tres, mesh.numDofPerNode)
  w_face_j_bar = zeros(Tres, mesh.numDofPerNode)
  w_bndry_j_bar = zeros(Tres, mesh.numDofPerNode)
  q_bndry_j_bar = zeros(Tres, mesh.numDofPerNode)
  qface_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  A0inv = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode)

  val_bar = -val_bar
  for i=1:length(faces)  # loop over faces with this BC
    bndry_i = faces[i]
    global_facenum = global_nums[i]
    q_el = sview(eqn.q, :, :, bndry_i.element)
    boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, q_el, qface)

    for j=1:mesh.numNodesPerFace

      # get components
      # convert to conservative variables if needed
      q_j = sview(qface, :, j)
      q_j_bar = sview(qface_bar, :, j)
      aux_vars[1] = calcPressure(eqn.params, q_j)
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      nrm_xy = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
      bndry_node = BoundaryNode(bndry_i, i, j)

      getDirichletState(functor, eqn.params, q_j, aux_vars, coords, nrm_xy,
                        q_bndry_j, bndry_node)

      # compute Lambda(delta_w)
      convertToIR(eqn.params, q_j, w_face_j)
      convertToIR(eqn.params, q_bndry_j, w_bndry_j)
      for k=1:mesh.numDofPerNode
        q_avg[k]   = 0.5*(q_j[k] + q_bndry_j[k])
        delta_w[k] = w_face_j[k] - w_bndry_j[k]
      end

      applyEntropyKernel(kernel, eqn.params, q_avg, delta_w, nrm_xy, flux)

      ## compute delta_w*B*Lambda(delta_w)
      #for k=1:mesh.numDofPerNode
      #  val += delta_w[k]*mesh.sbpface.wface[j]*flux[k]
      #end

      # reverse sweep
      fill!(delta_w_bar, 0); fill!(flux_bar, 0); fill!(q_avg_bar, 0)
      fill!(w_face_j_bar, 0); fill!(w_bndry_j_bar, 0); fill!(q_bndry_j_bar, 0)
      fill!(q_j_bar, 0)
      for k=1:mesh.numDofPerNode
        delta_w_bar[k] += mesh.sbpface.wface[j]*flux[k]*val_bar
        flux_bar[k]    += delta_w[k]*mesh.sbpface.wface[j]*val_bar
      end

      applyEntropyKernel_revq(kernel, eqn.params, q_avg, q_avg_bar, delta_w,
                              delta_w_bar, nrm_xy, flux, flux_bar)

      for k=1:mesh.numDofPerNode
        w_face_j_bar[k]  += delta_w_bar[k]
        w_bndry_j_bar[k] -= delta_w_bar[k]

        q_j_bar[k]       += 0.5*q_avg_bar[k]
        q_bndry_j_bar[k] += 0.5*q_avg_bar[k]
      end

      # reverse of convertToIR
      getIRA0inv(eqn.params, q_j, A0inv)
      smallmatTvec_kernel!(A0inv, w_face_j_bar, q_j_bar, 1, 1)
      getIRA0inv(eqn.params, q_bndry_j, A0inv)
      smallmatTvec_kernel!(A0inv, w_bndry_j_bar, q_bndry_j_bar, 1, 1)

      getDirichletState_revq(functor, eqn.params, q_j, q_j_bar, aux_vars,
                             coords, nrm_xy, q_bndry_j_bar, bndry_node)

    end  # end loop j

    q_el_bar = sview(func_deriv_arr, :, :, bndry_i.element)
    boundaryFaceInterpolate_rev!(mesh.sbpface, bndry_i.face, q_el_bar, qface_bar)

  end  # end loop i

  return nothing
end


text:
function F_autoreg(x, parms)
  gene, P2_gene, RNA, P, P2 = x
  k1, k1r, k2, k3, k4, k4r, k5, k6 = parms

  repression_binding = k1 * gene * P2
  reverse_repression = k1r * P2_gene
  transcription = k2 * gene
  translation = k3 * RNA
  dimerization = k4 * P * (P - 1) / 2
  dissociation = k4r * P2
  RNA_degradation = k5 * RNA
  protein_degradation = k6 * P

  return [repression_binding, reverse_repression, transcription, translation, dimerization, dissociation, RNA_degradation, protein_degradation]
end

function gillespie_autoreg(;
  gene = 10, P2_gene = 0, RNA = 0, P = 0, P2 = 0,
  k1 = 1.0, k1r = 10.0, k2 = 0.01, k3 = 10.0, k4 = 1.0, k4r = 1.0, k5 = 0.1, k6 = 0.01)
  # initial conditions
  x0 = [gene, P2_gene, RNA, P, P2]
  # stoichiometries
  nu = [
    [-1  1  0  0 -1]; # gene + P2 --> P2_gene
    [ 1 -1  0  0  1]; # P2_gene --> gene + P2
    [ 0  0  1  0  0]; # gene --> gene + RNA
    [ 0  0  0  1  0]; # RNA --> RNA + P
    [ 0  0  0 -2  1]; # P + P --> P2
    [ 0  0  0  2 -1]; # P2 --> P + P
    [ 0  0 -1  0  0]; # RNA --> 0
    [ 0  0  0 -1  0]  # P --> 0
  ]
  # parameter vector
  parms = [k1, k1r, k2, k3, k4, k4r, k5, k6]
  return x0, nu, parms
end

x0, nu, parms = gillespie_autoreg()
F = F_autoreg