text:
abstract type TimeseriesFileFormat end
abstract type TimeseriesFormatPeriodAsColumn <: TimeseriesFileFormat end
abstract type TimeseriesFormatYMDPeriodAsColumn <: TimeseriesFormatPeriodAsColumn end
abstract type TimeseriesFormatDateTimePeriodAsColumn <: TimeseriesFormatPeriodAsColumn end
abstract type TimeseriesFormatPeriodAsHeader <: TimeseriesFileFormat end
abstract type TimeseriesFormatYMDPeriodAsHeader <: TimeseriesFormatPeriodAsHeader end
abstract type TimeseriesFormatComponentsAsColumnsNoTime <: TimeseriesFileFormat end

"""
    read_time_series(file_path::AbstractString, component_name=nothing)

Return a TimeArray from a CSV file.

Pass component_name when the file does not have the component name in a column header.
"""
function read_time_series(file_path::AbstractString, component_name=nothing; kwargs...)
    if !isfile(file_path)
        msg = "Timeseries file doesn't exist : $file_path"
        throw(DataFormatError(msg))
    end

    file = CSV.File(file_path)
    @debug "Read CSV data from $file_path."

    return read_time_series(get_timeseries_format(file), file, component_name; kwargs...)
end

"""Return the timeseries format used in the CSV file."""
function get_timeseries_format(file::CSV.File)
    columns = propertynames(file)
    has_ymd = :Year in columns && :Month in columns && :Day in columns
    has_period = :Period in columns
    has_datetime = :DateTime in columns

    if has_period
        if has_datetime
            format = TimeseriesFormatDateTimePeriodAsColumn
        else
            format = TimeseriesFormatYMDPeriodAsColumn
        end
    elseif has_ymd
        format = TimeseriesFormatYMDPeriodAsHeader
    elseif !has_datetime
        format = TimeseriesFormatComponentsAsColumnsNoTime
    else
        @assert(false, "Unknown timeseries format in $(file.name)")
    end

    if format in (TimeseriesFormatYMDPeriodAsColumn, TimeseriesFormatYMDPeriodAsHeader)
        if !has_ymd
            throw(DataFormatError("$(file.name) is missing required Year/Month/Day"))
        end
    end

    return format
end

"""Return the column names with values (components)."""
function get_value_columns(::Type{TimeseriesFormatYMDPeriodAsColumn},
                           file::CSV.File)
    return [x for x in propertynames(file) if !in(x, (:Year, :Month, :Day, :Period))]
end

function get_value_columns(::Type{TimeseriesFormatDateTimePeriodAsColumn},
                           file::CSV.File)
    return [x for x in propertynames(file) if !in(x, (:DateTime, :Period))]
end

"""Return the column names with values."""
function get_value_columns(::Type{TimeseriesFormatComponentsAsColumnsNoTime},
                           file::CSV.File)
    return propertynames(file)
end

"""Return the column names that specify the Period."""
function get_period_columns(::Type{TimeseriesFormatPeriodAsColumn},
                            file::CSV.File)
    return [:Period]
end

function get_period_columns(::Type{TimeseriesFormatYMDPeriodAsHeader},
                            file::CSV.File)
    return [x for x in propertynames(file) if !in(x, (:Year, :Month, :Day))]
end

"""Return a vector of dicts of unique timestamps and their counts."""
function get_unique_timestamps(::Type{T}, file::CSV.File) where T <: TimeseriesFileFormat
    timestamps = Vector{Dict{String, Any}}()
    new_timestamp = x -> Dict("timestamp" => x, "count" => 1)

    for i in 1:length(file)
        timestamp = get_timestamp(T, file, i)
        if i == 1
            push!(timestamps, new_timestamp(timestamp))
        else
            if timestamp == timestamps[end]["timestamp"]
                timestamps[end]["count"] += 1
            else
                push!(timestamps, new_timestamp(timestamp))
            end
        end
    end

    @assert length(timestamps) > 0
    for timestamp in timestamps[2:end]
        @assert timestamp["count"] == timestamps[1]["count"]
    end

    return timestamps
end

"""Return a Dates.DateTime for the row in the CSV file."""
function get_timestamp(::Type{TimeseriesFormatYMDPeriodAsColumn}, file::CSV.File,
                       row_index::Int)
    return Dates.DateTime(file.Year[row_index], file.Month[row_index], file.Day[row_index])
end

function get_timestamp(::Type{TimeseriesFormatDateTimePeriodAsColumn}, file::CSV.File,
                       row_index::Int)
    return Dates.DateTime(file.DateTime[row_index])
end

function get_timestamp(::Type{TimeseriesFormatYMDPeriodAsHeader}, file::CSV.File,
                       row_index::Int)
    return get_timestamp(TimeseriesFormatYMDPeriodAsColumn, file, row_index)
end

"""Return a TimeSeries.TimeArray representing the CSV file.

This version of the function only has component_name to match the interface. It is unused.
"""
function read_time_series(
                          ::Type{T},
                          file::CSV.File,
                          component_name=nothing;
                          kwargs...
                         ) where T <: TimeseriesFormatPeriodAsColumn
    timestamps = Vector{Dates.DateTime}()
    step = get_step_time(T, file, file.Period)

    # All timestamps must be sequential by step, so we can ignore the timestamps in the
    # file after the first one.
    # They were validated in get_step_time.
    first = get_timestamp(T, file, 1)
    push!(timestamps, first)
    for i in 2:length(file)
        timestamp = first + step * (i - 1)
        push!(timestamps, timestamp)
    end

    value_columns = get_value_columns(T, file)
    vals = [getproperty(file, x) for x in value_columns]

    return TimeSeries.TimeArray(timestamps, hcat(vals...), value_columns)
end

"""This version of the function supports the format where there is no column header for
a component, so the component_name must be passed in.
"""
function read_time_series(
                          ::Type{T},
                          file::CSV.File,
                          component_name::AbstractString;
                          kwargs...
                         ) where T <: TimeseriesFormatPeriodAsHeader
    timestamps = Vector{Dates.DateTime}()

    period_cols_as_symbols = get_period_columns(T, file)
    period = [parse(Int, string(x)) for x in period_cols_as_symbols]
    step = get_step_time(T, file, period)

    # All timestamps must be sequential by step, so we can ignore the timestamps in the
    # file after the first one.
    # They were validated in get_step_time.

    first = Dates.DateTime(Dates.today())
    count = 0
    for i in 1:length(file)
        if i == 1
            first = get_timestamp(T, file, 1)
        end
        for j in 1:length(period)
            timestamp = first + step * count
            count += 1
            push!(timestamps, timestamp)
        end
    end

    vals = Vector{Float64}()
    for i in 1:length(file)
        for period in period_cols_as_symbols
            val = getproperty(file, period)[i]
            push!(vals, val)
        end
    end

    return TimeSeries.TimeArray(timestamps, vals, Symbol.([component_name]))
end

"""This version of the function only has component_name to match the interface.
It is unused.

Set start_datetime as a keyword argument for the starting timestamp, otherwise the current
day is used.
"""
function read_time_series(
                          ::Type{T},
                          file::CSV.File,
                          component_name=nothing;
                          kwargs...
                         ) where T <: TimeseriesFormatComponentsAsColumnsNoTime
    timestamps = Vector{Dates.DateTime}()
    step = get_step_time(T, file)

    start = get(kwargs, :start_datetime, Dates.DateTime(Dates.today()))
    for i in 1:length(file)
        timestamp = start + step * (i - 1)
        push!(timestamps, timestamp)
    end

    value_columns = get_value_columns(T, file)
    vals = [getproperty(file, x) for x in value_columns]

    return TimeSeries.TimeArray(timestamps, hcat(vals...), value_columns)
end

"""Return a DateTime for the step between values as specified by the period in the file."""
function get_step_time(
                       ::Type{T},
                       file::CSV.File,
                       period::AbstractArray,
                      ) where T <: TimeseriesFileFormat
    timestamps = get_unique_timestamps(T, file)
    if T <: TimeseriesFormatPeriodAsColumn
        num_steps = timestamps[1]["count"]
    elseif T <: TimeseriesFormatPeriodAsHeader
        num_steps = period[end]
    else
        error("unsupported $T")
    end

    @debug timestamps, num_steps
    if length(timestamps) == 1
        # TODO: Not sure how to handle this. We could make specific functions for each type.
        # For any YMD format the lowest resolution is Day.
        # What is it for DateTime? We can't infer from one value.
        resolution = Dates.Day(1)
    else
        resolution = timestamps[2]["timestamp"] - timestamps[1]["timestamp"]
        if length(timestamps) > 2
            for i in 3:length(timestamps)
                diff = timestamps[i]["timestamp"] - timestamps[i - 1]["timestamp"]
                if diff != resolution
                    msg = "conflicting resolution=$resolution i=$i diff=$diff"
                    throw(DataFormatError(msg))
                end
            end
        end
    end

    return calculate_step_time(resolution, num_steps)
end

function get_step_time(
                       ::Type{T},
                       file::CSV.File,
                      ) where T <: TimeseriesFormatComponentsAsColumnsNoTime
    resolution = Dates.Day(1)
    num_steps = length(file)
    return calculate_step_time(resolution, num_steps)
end

function calculate_step_time(resolution::Dates.Period, num_steps::Int)
    # Seconds should be the lowest possible resolution.
    step = Dates.Second(resolution) / num_steps
    @debug "file has step time of $step" resolution num_steps
    return step
end


text:
## Uniform ##

@scalar_rule(
    uniformlogpdf(a, b, x),
    @setup(
        insupport = a <= x <= b,
        diff = b - a,
        c = insupport ? inv(diff) : inv(one(diff)),
    ),
    (c, -c, ZeroTangent()),
)

ext: jl


text:
function Base.show(io::IO, ::MIME"text/html", sim::Simulation)
    show_simulation_table(io, sim, backend = :html)
    println(io)
end

function Base.show(io::IO, ::MIME"text/plain", sim::Simulation)
    show_simulation_table(io, sim, backend = :auto)
    println(io)
end

function Base.show(io::IO, ::SimulationInputs)
    println(io, "SimulationInputs()")
end

function Base.show(io::IO, smr::SmallSignalOutput)
    val = smr.stable ? "is" : "is not"
    println(io, "The system $(val) small signal stable")
end

function Base.show(io::IO, pert::ControlReferenceChange)
    println(
        io,
        "Change of $(pert.signal) at time t = $(pert.time) of device $(PSY.get_name(pert.device)) to $(pert.ref_value) per unit",
    )
end

function Base.show(io::IO, ::MIME"text/html", res::SimulationResults)
    show_results_table(io, res, backend = :html)
    println(io)
end

function Base.show(io::IO, ::MIME"text/plain", res::SimulationResults)
    show_results_table(io, res, backend = :auto)
    println(io)
end

function show_results_table(io::IO, res::SimulationResults; kwargs...)
    header = ["Property", "Value"]
    sys = get_system(res)
    sol = get_solution(res)
    table = [
        "System Base Power [MVA]" string(PSY.get_base_power(sys))
        "System Base Frequency [Hz]" string(PSY.get_frequency(sys))
        "Time Span" string((sol.t[1], sol.t[end]))
        "Total Time Steps" string(length(sol.t))
        "Number of States" string(length(sol.u[1]))
        "Total solve time" string(res.time_log[:timed_solve_time])
    ]
    PrettyTables.pretty_table(
        io,
        table,
        header;
        title = "Simulation Results Summary",
        alignment = :l,
        kwargs...,
    )
end

function show_simulation_table(
    io::IO,
    sim::Simulation{T};
    kwargs...,
) where {T <: SimulationModel}
    header = ["Property", "Value"]
    val_multimachine = sim.multimachine ? "Yes" : "No"
    val_initialized = sim.initialized ? "Yes" : "No"
    val_model = T == ResidualModel ? "Residual Model" : "Mass Matrix Model"
    table = [
        "Status" sim.status
        "Simulation Type" val_model
        "Initialized?" val_initialized
        "Multimachine system?" val_multimachine
        "Time Span" string(sim.tspan)
        "Number of States" string(length(sim.x0_init))
        "Number of Perturbations" string(length(sim.perturbations))
    ]
    PrettyTables.pretty_table(
        io,
        table,
        header;
        title = "Simulation Summary",
        alignment = :l,
        kwargs...,
    )
end


text:
using Test, SafeTestsets

@safetestset "InpParser Tests" begin
    include("inp_parser/parser.jl")
end
@safetestset "TopOptProblems Tests" begin
    include("topopt_problems/problems.jl")
    include("topopt_problems/metadata.jl")
end
@safetestset "Functions" begin
    include("functions.jl")
end
# @safetestset "Solver" begin
#     include("fea/solvers.jl")
# end
@safetestset "Truss Problem Tests" begin
    include("truss_topopt_problems/test_problem.jl")
    include("truss_topopt_problems/test_fea.jl")
end
@safetestset "Examples" begin
    @safetestset "CSIMP" begin
        include("examples/csimp.jl")
    end
    @safetestset "Global Stress" begin
        include("examples/global_stress.jl")
    end
    @safetestset "More examples" begin
        include("examples/test_examples.jl")
    end
end

# @safetestset "CSIMP Tests" begin include("csimp.jl") end


text:
###############################################################################
# Stochiometry for a given reaction is a vector of pairs mapping species id to
# stochiometric coefficient.
###############################################################################

@inline @fastmath function evalrxrate(speciesvec::AbstractVector{T}, rxidx::S,
                              majump::MassActionJump{U,V,W,X})::R where {T,S,R,U <: AbstractVector{R},V,W,X}
    val = one(T)
    @inbounds for specstoch in majump.reactant_stoch[rxidx]
        specpop = speciesvec[specstoch[1]]
        val    *= specpop
        @inbounds for k = 2:specstoch[2]
            specpop -= one(specpop)
            val     *= specpop
        end
    end

    @inbounds return val * majump.scaled_rates[rxidx]
end

@inline @fastmath function executerx!(speciesvec::AbstractVector{T}, rxidx::S,
                                      majump::M) where {T,S,M <: AbstractMassActionJump}
    @inbounds net_stoch = majump.net_stoch[rxidx]
    @inbounds for specstoch in net_stoch
        speciesvec[specstoch[1]] += specstoch[2]
    end
    nothing
end

@inline @fastmath function executerx(speciesvec::SVector{T}, rxidx::S,
                                      majump::M) where {T,S,M <: AbstractMassActionJump}
    @inbounds net_stoch = majump.net_stoch[rxidx]
    @inbounds for specstoch in net_stoch
        speciesvec = setindex(speciesvec,speciesvec[specstoch[1]]+specstoch[2],specstoch[1])
    end
    speciesvec

    #=
    map(net_stoch) do stoch
        @inbounds speciesvec[stoch[1]] + stoch[2]
    end
    =#
end

function scalerates!(unscaled_rates::AbstractVector{U}, stochmat::AbstractVector{V}) where {U,S,T,W <: Pair{S,T}, V <: AbstractVector{W}}
    @inbounds for i in eachindex(unscaled_rates)
        coef = one(T)
        @inbounds for specstoch in stochmat[i]
            coef *= factorial(specstoch[2])
        end
        unscaled_rates[i] /= coef
    end
    nothing
end

function scalerates!(unscaled_rates::AbstractMatrix{U}, stochmat::AbstractVector{V}) where {U,S,T,W <: Pair{S,T}, V <: AbstractVector{W}}
    @inbounds for i in size(unscaled_rates, 1)
        coef = one(T)
        @inbounds for specstoch in stochmat[i]
            coef *= factorial(specstoch[2])
        end
        unscaled_rates[i,:] /= coef
    end
    nothing
end

function scalerate(unscaled_rate::U, stochmat::AbstractVector{Pair{S,T}}) where {U <: Number, S, T}
    coef = one(T)
    @inbounds for specstoch in stochmat
        coef *= factorial(specstoch[2])
    end
    unscaled_rate /= coef
end

###############################################################################
# dependency graph when MassActionJump uses pairs to represent (species,stoich)
###############################################################################

# map from species to reactions depending on that species
# uses a Vector instead of a Set as the latter requires isEqual,
# and by using an underlying Dict can be slower for small numbers
# of dependencies
function var_to_jumps_map(numspec, ma_jumps::AbstractMassActionJump)

    numrxs = get_num_majumps(ma_jumps)

    # map from a species to reactions that depend on it
    spec_to_dep_rxs = [Vector{Int}() for n = 1:numspec]
    for rx in 1:numrxs
        for (spec,stoch) in ma_jumps.reactant_stoch[rx]
            push!(spec_to_dep_rxs[spec], rx)
        end
    end

    foreach(s -> unique!(sort!(s)), spec_to_dep_rxs)
    spec_to_dep_rxs
end

"make a map from reactions to dependent species"
function jump_to_vars_map(majumps)
    [[s for (s, c) in majumps.net_stoch[i]] for i in 1:get_num_majumps(majumps)]
end

# dependency graph is a map from a reaction to a vector of reactions
# that should depend on species it changes
function make_dependency_graph(numspec, ma_jumps::AbstractMassActionJump)

    numrxs          = get_num_majumps(ma_jumps)
    spec_to_dep_rxs = var_to_jumps_map(numspec, ma_jumps)

    # create map from rx to reactions depending on it
    dep_graph = [Vector{Int}() for n = 1:numrxs]
    for rx in 1:numrxs

        # rx changes spec, hence rxs depending on spec depend on rx
        for (spec,stoch) in ma_jumps.net_stoch[rx]
            for dependent_rx in spec_to_dep_rxs[spec]
                push!(dep_graph[rx], dependent_rx)
            end
        end
    end

    add_self_dependencies!(dep_graph, dosort=false)
    foreach(deps -> unique!(sort!(deps)), dep_graph)
    dep_graph
end

# update dependency graph to make sure jumps depend on themselves
function add_self_dependencies!(dg; dosort=true)
    for (i,jump_deps) in enumerate(dg)
        if !any(y->isequal(y,i), jump_deps)
            push!(jump_deps, i)
            dosort && sort!(jump_deps)
        end
    end
end


text:
for f in (:run, :success)
    @eval begin
        """
            $($f)(container::DockerContainer, config::DockerConfig, user_cmd::Cmd; kwargs...)
        """
        function $f(container::DockerContainer,
                    config::DockerConfig,
                    cmd::Cmd;
                    kwargs...)
            container_cmd = construct_container_command(
                container,
                config,
                cmd,
            )
            docker_cmd = pipeline(
                container_cmd;
                config.stdin,
                config.stdout,
                config.stderr,
            )
            if config.verbose
                @info("Running sandboxed command", cmd.exec)
            end
            return $f(docker_cmd; kwargs...)
        end
    end
end


text:
function _make_janelia_gui()

    janelia_grid=Grid()
    janelia_label=Label("Janelia Parameters")
    janelia_grid[1,1]=janelia_label

    janelia_seed_thres=SpinButton(0.01:.01:1.0)
    setproperty!(janelia_seed_thres,:value,0.99)
    janelia_grid[1,2]=janelia_seed_thres
    janelia_grid[2,2]=Label("Seed Threshold")

    janelia_seed_iterations=SpinButton(1:1:10)
    setproperty!(janelia_seed_iterations,:value,1)
    janelia_grid[1,3]=janelia_seed_iterations
    janelia_grid[1,3]=Label("Seed Iterations")

    janelia_win=Window(janelia_grid)
    Gtk.showall(janelia_win)
    visible(janelia_win,false)

    j_widgets=janelia_widgets(janelia_win,janelia_seed_thres,janelia_seed_iterations)
end

function add_janelia_callbacks(w::janelia_widgets,handles::Tracker_Handles)

    signal_connect(jt_seed_thres_cb,w.jt_seed_thres_button,"value-changed",Void,(),false,(handles,))
    signal_connect(jt_seed_iterations_cb,w.jt_seed_iterations_button,"value-changed",Void,(),false,(handles,))

    nothing
end

function jt_seed_thres_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    thres=getproperty(han.janelia_widgets.jt_seed_thres_button,:value,Float64)

    change_JT_param(:paramSEED_THRESH,convert(Float32,thres))

    nothing
end

function jt_seed_iterations_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    iterations=getproperty(han.janelia_widgets.jt_seed_iterations_button,:value,Int64)

    change_JT_param(:paramSEED_ITERATIONS,convert(Int32,iterations))

    nothing
end


text:
using PyPlot
using SMC

function render_hmm!(hmm::HiddenMarkovModel)
    # renders a HMM in the current figure
    # example usage:
    # plt[:figure]()
    # render!(hmm)
    # plt[:savefig]("hmm.png")

    num_states = hmm.num_states
    num_obs = hmm.num_obs

    # three subplots left to right:
    # 1. the prior (column vector)
    # 2. the transition matrix (square matrix)
    # 3. the observation matrix (rectangle)
    # the height is always num_states
    width_ratios = [num_states, num_obs]
    gs = matplotlib[:gridspec][:GridSpec](1, 2, width_ratios=width_ratios)
    trans_ax = plt[:subplot](gs[1])
    obs_ax = plt[:subplot](gs[2])
    prior_mat = Array{Float64,2}(num_states,1)
    prior_mat[:,1] = hmm.initial_state_prior
    trans_ax[:imshow](hmm.transition_model, cmap="gray", origin="lower", vmin=0.0, vmax=1.0)
    obs_ax[:imshow](hmm.observation_model, cmap="gray", origin="lower", vmin=0.0, vmax=1.0)
    trans_ax[:set_title]("transition model")
    obs_ax[:set_title]("observation model")
    for ax in [trans_ax, obs_ax]
        ax[:set_yticks](0:num_states-1)
        ax[:set_yticklabels](1:num_states)
    end

    trans_ax[:set_ylabel]("previous state", fontsize=12)
    trans_ax[:set_xlabel]("next state", fontsize=12)
    trans_ax[:set_xticks](0:num_states-1)
    trans_ax[:set_xticklabels](1:num_states)
    trans_ax[:set_xlim]([-0.5, num_states - 0.5])
    trans_ax[:set_ylim]([-0.5, num_states- 0.5])
    trans_ax[:set_xticks]((1:num_states-1) - 0.5, minor=true)
    trans_ax[:set_yticks]((1:num_states-1) - 0.5, minor=true)

    obs_ax[:set_ylabel]("state", fontsize=12)
    obs_ax[:set_xlabel]("observation", fontsize=12)
    obs_ax[:set_xticks](0:num_obs-1)
    obs_ax[:set_xticklabels](1:num_obs)
    obs_ax[:set_xlim]([-0.5, num_obs - 0.5])
    obs_ax[:set_ylim]([-0.5, num_states- 0.5])
    obs_ax[:set_xticks]((1:num_obs-1) - 0.5, minor=true)
    obs_ax[:set_yticks]((1:num_states-1) - 0.5, minor=true)

    plt[:tight_layout]()
end

function render_hmm_states!(hmm::HiddenMarkovModel, states::Array{Int,1},
                            grid::Bool, xticks, yticks)
    if maximum(states) > hmm.num_states || minimum(states) < 1
        error("bad states")
    end
    mat = zeros(hmm.num_states, length(states))
    for i=1:length(states)
        mat[states[i], i] = 1.0
    end
    ax = plt[:gca]()

    # white will be 0.0, and black will be 1.0 (to save ink)
    ax[:imshow](1.0 - mat, cmap="gray", origin="lower", vmin=0.0, vmax=1.0)

    ax[:set_yticks](yticks-1)
    ax[:set_xticks](xticks-1)
    ax[:set_yticklabels](yticks)
    ax[:set_xticklabels](xticks)

    ax[:set_xlim]([-0.5, length(states) - 0.25])
    ax[:set_ylim]([-0.5, hmm.num_states - 0.25])

    if grid
        # set grid lines using minor ticks
        ax[:set_xticks]((1:length(states)-1) - 0.5, minor=true);
        ax[:set_yticks]((1:hmm.num_states-1) - 0.5, minor=true);
        ax[:grid](which="minor", color="orange", linewidth=1)
    end
end

function render_hmm_posterior_marginals!(hmm::HiddenMarkovModel,
                                         marginals::Array{Float64,2},
                                         grid::Bool, xticks, yticks)
    if size(marginals)[1] != hmm.num_states
        error("bad states")
    end
    num_steps = size(marginals)[2]
    ax = plt[:gca]()

    # white will be 0.0, and black will be 1.0 (to save ink)
    ax[:imshow](1.0 - marginals, cmap="gray", origin="lower", vmin=0.0, vmax=1.0)

    ax[:set_yticks](yticks-1)
    ax[:set_xticks](xticks-1)
    ax[:set_yticklabels](yticks)
    ax[:set_xticklabels](xticks)
    ax[:set_xlim]([-0.5, num_steps - 0.25])
    ax[:set_ylim]([-0.5, hmm.num_states - 0.25])

    # set grid lines using minor ticks
    if grid
        ax[:set_xticks]((1:num_steps-1) - 0.5, minor=true);
        ax[:set_yticks]((1:hmm.num_states-1) - 0.5, minor=true);
        ax[:grid](which="minor", color="orange", linewidth=1)
    end
end

function render_hmm_observations!(hmm::HiddenMarkovModel,
                                  observations::Array{Int,1})
    if maximum(observations) > hmm.num_obs || minimum(observations) < 1
        error("bad observations")
    end
    mat = zeros(hmm.num_obs, length(observations))
    for i=1:length(observations)
        mat[observations[i], i] = 1.0
    end
    ax = plt[:gca]()
    ax[:imshow](mat, cmap="gray", origin="lower", vmin=0.0, vmax=1.0)
    ax[:set_ylabel]("observation", fontsize=12)
    ax[:set_xlabel]("time step", fontsize=12)
    ax[:set_yticks](0:2:hmm.num_obs-1)
    ax[:set_yticklabels](1:2:hmm.num_obs)
    ax[:set_xticks](0:2:length(observations))
    ax[:set_xticklabels](1:2:length(observations))
    ax[:set_xlim]([-0.5, length(observations) - 0.5])
    ax[:set_ylim]([-0.5, hmm.num_obs - 0.5])

    # set grid lines using minor ticks
    ax[:set_xticks]((1:length(observations)-1) - 0.5, minor=true);
    ax[:set_yticks]((1:hmm.num_obs-1) - 0.5, minor=true);
end


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

@testset "正方格子的功能" begin
    qs = common_square_lattice(0.)
    println(qs)
end

@testset "三角格子的功能" begin
    ts = common_triangle_lattice(0.)
    println(ts)
end

@testset "kagome格子的功能" begin
    ks = upperband_kagome_lattice(0.)
    plt = draw_lines(ks.brillouin.edges)
    @savepng plt "test_latt_ks1"
end


text:
module TestDistributionsDependency

using Distributions
using Test

@testset "Distributions" begin

    ## reference values, taken from qnorm [R]
    qnorm = ((0.1, -1.281552),
         (0.5, 0.0),
         (0.9, 1.281552),
         (0.95, 1.644854),
         (0.99, 2.326348))

    for (alpha, ref) in qnorm
        v = quantile(Normal(), alpha)
        @test isapprox(v, ref, atol = 1e-6)
    end

end

end


text:
using jInv.Mesh
using jInv.LinearSolvers
using SparseArrays
using LinearAlgebra
using Test
println("===  Test Julia wrapper with symmetric (real) matrix ====");
domain = [0.0, 1.0, 0.0, 1.0];
n      = [33,17];
Mr     = getRegularMesh(domain,n)
G      = getNodalGradientMatrix(Mr);
m      = sparse(Diagonal(exp.(randn(size(G,1)))));
Ar     = G'*m*G;
Ar     = Ar + 1e-1*norm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
N      = size(Ar,2);
b      = Ar*rand(N,1);
bs     = Ar*rand(N,3);
Bs     = (b,sparse(b),bs,sparse(bs));

sNonSym = getJuliaSolver(sym = 0);
sSym    = getJuliaSolver(sym = 1);
sGenSym = getJuliaSolver(sym = 2);
sols    = (sNonSym,sSym,sGenSym)

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	for k=1:length(sols)
		print("sym=$(sols[k].sym), ")
		xtt1,  = solveLinearSystem(Ar,Bs[j],sols[k]);
		@test norm(Ar*xtt1-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-10
	end
	print("\n")
end
println("\n")

println("===  Test Julia wrapper with shifted (complex) Laplacian ====");
Ar     = G'*m*G +  (1+1im)*sparse(1.0I,size(Ar,2),size(Ar,2));
b      = Ar*rand(N,1);
bs     = Ar*rand(N,3);
Bs     = (b,sparse(b),bs,sparse(bs));

sNonSym = getJuliaSolver(sym = 0);
sSym    = getJuliaSolver(sym = 1);
sGenSym = getJuliaSolver(sym = 2);
sols    = (sNonSym,sSym,sGenSym)

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	for k=1:length(sols)
		print("sym=$(sols[k].sym), ")
		xtt2,  = solveLinearSystem(Ar,Bs[j],sols[k]);
		@test norm(Ar*xtt2-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-10
	end
	print("\n")
end
println("\n")

println("===  Test Julia Wrapper: nonsymmetric matrices ====");
n = 50
sNonSym  = getJuliaSolver(sym = 0);
sNonSym  = copySolver(sNonSym);
A = sprandn(n,n,5/n) + 10*sparse(1.0I,n,n)
B = randn(n)

X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-10

# multiple sparse rhs -
# Bs = sprandn(n,10,.3)
# X, = solveLinearSystem(A,Bs,sNonSym,0);
# @test norm(A*X - Bs,Inf)/norm(Bs,Inf) < 1e-10

X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-10

X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-10

X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-10

println("===  End Test Julia Wrapper ====");


text:
"Generate a noise image for experimentation."
function noisepic_3(A, physheight, physwidth)
    n_rows, n_cols = size(A)
    no = Matrix{Float64}(undef, n_rows, n_cols)
    n_rows = size(no)[1]
    n_cols_in_division = div(size(no)[2], 3)
    length_one_cycle = CUTOFF_23 * DU_2
    length_one_pixel = physheight / n_rows
    octaves = round(Int, log2(length_one_cycle / length_one_pixel))
    dimension_aligned(OU_27 + (-0.4 * physwidth, 0.0m),  OU_27 +  (-0.4 * physwidth, length_one_cycle ))
    n_cycles = physheight / length_one_cycle
    for row in 1:n_rows
        for col in 1:n_cols
            divno = 1 + col ÷ n_cols_in_division
            rate = 1.414 * divno * n_cycles / n_rows
            # Noise is a function generating noisy data,
            # but where the lowest frequency has wavelength 1.41.
            # I.e., one full cycle when x or y varies from 0 to 1.41.
            no[row, col] = noise(row * rate, col * rate, detail = octaves, persistence = 1.0) +
                noise(row * rate * 0.75, col * rate * 0.75, detail = octaves, persistence = 1.0) +
                noise(row * rate * 0.25, col * rate * 0.25, detail = octaves, persistence = 1.0)
        end
    end
    normalize_datarange(no)
end
"Lowpass filter coefficients, n is an even number"
function lowpass_coefficients_1(n)
    @assert iseven(n)
    # Single side, odd number
    m = div(n - 1,  2)
    ωc = π / 16
    rng = -m:1:-1
    # Filter coefficients h
    h = sin.(ωc * rng) ./ (rng * π)
    push!(h, ωc / π)
    append!(h, reverse(h[1:m]))
    h
end

"Highpass filter coefficients, n is an even number"
function highpass_coefficients_1(n)
    @assert iseven(n)
    # Single side, odd number
    m = div(n - 1,  2)
    ωc = π / 16
    rng = -m:1:-1
    # Filter coefficients h
    h = - sin.(ωc * rng) ./ (rng * π)
    push!(h, 1 - ωc / π)
    append!(h, reverse(h[1:m]))
    h
end

"Blackman window coeffients, n is an even number"
function blackman_coefficients_1(n)
    @assert iseven(n)
    m = div(n - 1,  2)
    freqs = (-m:1:m) * π / m
    0.42 .+ 0.5 * cos.( freqs) .+ 0.08 * cos.(2 * freqs)
end

"Finite impulse response filter, high pass, total length n (even)"
function fir_coefficients_1(n)
    lowpass_coefficients_1(n) .* blackman_coefficients_1(n)
end

function convolute_image_3(xs, ys, f_xy, n_xy, h, cutoff)
    # Output buffer
    M = Array{Float64}(undef, length(ys), length(xs)) # Image processing convention: a column has a horizontal line of pixels
    fir = SVector{NS_3 * 2 -1}(fir_coefficients_1(NS_3 * 2))
    @assert eltype(xs) <: Quantity{Float64}
    @assert eltype(ys) <: Quantity{Float64}
    @assert n_xy(xs[1], ys[1]) isa Float64
    @assert first(ys) < last(ys)
    # prepare a mutable but fast working buffer forward, and another backward.
    vxf = similar(SVector{NS_3}(fill(xs[1], NS_3)))
    vyf = similar(vxf)
    vxb = similar(vxf)
    vyb = similar(vxf)
    rowsy = length(ys)
    for (i::Int64, x::Quantity{Float64}) in enumerate(xs), (j::Int64, y::Quantity{Float64}) in enumerate(ys)
        vxf[1], vyf[1], vxb[1], vyb[1] = x, y, x, y
        # Find the streamline path
        # Put coordinates to sample in forward and aft buffers (no need to store them in memory really?)
        rk4_steps!(f_xy, vxf, vyf, h)
        rk4_steps!(f_xy, vxb, vyb, -h)
        vx = vcat(reverse(vxf), vxb[1:(end -1)])
        vy = vcat(reverse(vyf), vyb[1:(end -1)])
        # Find the intensity for our one pixel.
        pv  = convolute_pixel(fir, vx, vy, n_xy, h, cutoff)
        # Find the original indexes and update the image matrix
        ii = rowsy + 1 - j
        jj = i
        M[ii, jj] = pv / oneunit(pv)
    end
    M
end

"""
    convolute_pixel(wf, wb, vxf, vyf, vxb, vyb, nxy, h, cutoff))
where
    wf    Window forward: A vector defining a convolution (filter) window, or kernel, with the same length as buffers.
          The second value is used for weighting the second coordinate.
    wb    Window backward: Similar, but for the backwards direction. The second value is used for weighting the second coordinate.
          Hence, it will be identical to wf in many cases.
    vxf   Vector or buffer of x coordinates, forward direction
    vyf   Vector y coordinates, forward.
    vxb   Vector x, backward
    vyb   Vector y, backwards direction
    nxy   Function, noise field
    h     Step length, duration
    cutoff Velocity, substituted to calculate out-of bounds streamlines
"""
function convolute_pixel(wf, wb, vxf, vyf, vxb, vyb, nxy, h, cutoff)
    @assert length(vxf) == length(vyf)
    @assert length(vxb) == length(vyb)
    x0 = vxf[1]
    y0 = vyf[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    pv = 0.0 * prevcontrib
    for i = 2:length(vxf)
        prevcontrib = _convolve_contrib(prevcontrib, wf, nxy, vxf, vyf, i)
        pv += prevcontrib
    end
    x0 = vxb[1]
    y0 = vyb[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    for i = 2:length(vxb)
        prevcontrib = _convolve_contrib(prevcontrib, wb, nxy, vxb, vyb, i)
        pv += prevcontrib
    end
    pv
end
"""
    convolute_pixel(w, vx,  vy, nxy, h, cutoff)
where
    w     Window : A vector defining a convolution (filter) window, or kernel, with the same length as buffers.
    vx    Vector or buffer of x coordinates
    vy    Vector y coordinates
    nxy   Function, noise field
    h     Step length, duration
    cutoff Maximum magnitude
"""
function convolute_pixel(w, vx, vy, nxy, h, cutoff)
    @assert length(vx) == length(vy)
    @assert length(w) == length(vy)
    x0 = vx[1]
    y0 = vy[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    pv = 0.0 * prevcontrib
    for i = 2:length(vx)
        prevcontrib = _convolve_contrib(prevcontrib, w, nxy, vx, vy, i)
        pv += prevcontrib
    end
    pv
end
function _convolve_contrib(prev, w, nxy, vx, vy, i)
    @assert i > 1
    x = vx[i]
    y = vy[i]
    x0 = vx[i - 1]
    y0 = vy[i - 1]
    if isnan(x0) || isnan(y0) || isnan(x) || isnan(y)
        return prev
    end
    k = w[i]
    ds =  hypot(x - x0, y - y0)
    n = nxy(x, y)
    n * k * ds
end


text:
include("./src/SFFM.jl")
using Plots, LinearAlgebra, KernelDensity, StatsBase

## Define the model
T = [-2.0 1.0 1.0; 1.0 -2.0 1.0; 1.0 1.0 -2.0]
C = [1.0; -2.0; 0.0]
r = (
    r = function (x)
        [(1.1 .+ sin.(π*x)) (sqrt.(x .* (x .> 0)) .+ 1) ((x .> 0) .* x .+ 1)]
    end, # r = function (x); [1.0.+0.01*x 1.0.+0.01*x 1*ones(size(x))]; end,
    R = function (x)
        [(1.1 .* x .- cos.(π*x)./π) ((x .* (x .> 0)) .^ (3 / 2) .* 2 / 3 .+ 1 * x) (
            (x .> 0) .* x .^ 2 / 2.0 .+ 1 * x
        )]
    end, # R = function (x); [1*x.+0.01.*x.^2.0./2 1*x.+0.01.*x.^2.0./2 1*x]; end
) # [1*(x.<=2.0).-2.0*(x.>1.0) -2.0*(x.<=2.0).+(x.>2.0)] # [1.0 -2.0].*ones(size(x))#

# r = (
#     r = function (x)
#         [ones(size(x)) ones(size(x)) ones(size(x))]
#     end, # r = function (x); [1.0.+0.01*x 1.0.+0.01*x 1*ones(size(x))]; end,
#     R = function (x)
#         [x x x]
#     end, # R = function (x); [1*x.+0.01.*x.^2.0./2 1*x.+0.01.*x.^2.0./2 1*x]; end
# )

Bounds = [-10 10; -Inf Inf]
model = SFFM.Model(T = T, C = C, r = r, Bounds = Bounds)

# in out Y-level
y = 10

## Simulate the model
NSim = 50000
IC = (φ = ones(Int, NSim), X = zeros(NSim), Y = zeros(NSim))
# IC = (φ = 2 .*ones(Int, NSim), X = -10*ones(NSim), Y = zeros(NSim))
# IC = (
#     φ = sum(rand(NSim) .< [1 / 3 2 / 3 1], dims = 2),
#     X = 20.0 .* rand(NSim) .- 10,
#     Y = zeros(NSim),
# )
sims =
    SFFM.SimSFFM(model = model, StoppingTime = SFFM.FirstExitY(u = -Inf, v = y), InitCondition = IC)

## Define the mesh
Δ = 5
Nodes = collect(Bounds[1, 1]:Δ:Bounds[1, 2])
# Fil = Dict{String,BitArray{1}}(
#     "1+" => trues(length(Nodes) - 1),
#     "2+" => trues(length(Nodes) - 1),
#     "3+" => trues(length(Nodes) - 1),
#     "p2+" => trues(1),
#     "p3+" => trues(1),
#     "q1+" => trues(1),
#     "q3+" => trues(1),
# )
NBases = 1
Basis = "lagrange"
mesh = SFFM.MakeMesh(model = model, Nodes = Nodes, NBases = NBases, Basis=Basis)

## Construct all DG operators
All = SFFM.MakeAll(model = model, mesh = mesh, approxType = "projection")
Matrices = All.Matrices
MatricesR = All.MatricesR
B = All.B
R = All.R
D = All.D
DR = All.DR

## initial condition
initpm = [
    zeros(sum(model.C.<=0)) # LHS point mass
    zeros(sum(model.C.>=0)) # RHS point mass
]
initprobs = zeros(Float64,1,mesh.NIntervals,model.NPhases)
initprobs[1,(1+mesh.NIntervals)÷2+1,1] = 1
initdist = (pm = initpm, distribution = initprobs, x = Matrix(mesh.CellNodes[1,:]'+mesh.Δ'/2), type = "probability")
x0 = SFFM.Dist2Coeffs(model = model, mesh = mesh, Distn = initdist)
# p = SFFM.PlotSFM(model=model,mesh=mesh,Dist=initdist)

## approximations to exp(Dy)
h = 0.0001
@time yvals = SFFM.EulerDG(D = D["++"](s = 0), y = y, x0 = x0, h = h)

# convert to densities
# densities
simdensity = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="density")
DGdensity = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=yvals,type="density")

# convert to probabilities
# get estimates of probabilities
simprobs = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")
probs = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=yvals,type="probability")

## plots
# plot solutions
# densities
p = SFFM.PlotSFM(model=model,mesh=mesh,Dist=DGdensity,color=7)
# plot sims
p = SFFM.PlotSFM!(p;model=model,mesh=mesh,Dist=simdensity,color=4)

# probabilities
p = SFFM.PlotSFM(model=model,mesh=mesh,Dist=probs)
# plot sims
p = SFFM.PlotSFM!(p;model=model,mesh=mesh,Dist=simprobs,color=2)

## other analysis
# compute errors
derrpm = DGdensity.pm-simdensity.pm
derrdensity = DGdensity.distribution - simdensity.distribution
derrdist = (pm=derrpm,distribution=derrdensity,x=DGdensity.x,type="density")

# plot
p = SFFM.PlotSFM(model=model,mesh=mesh,Dist=derrdist)

# display point mass data
pmdata = [
    "." (Nodes[1] * ones(sum(model.C .<= 0)))' (Nodes[end] * ones(sum(model.C .>= 0)))'
    "sim" simprobs.pm'
    "pm" probs.pm'
]
display(pmdata)

# probabilities
# compute errors
perrpm = probs.pm-simprobs.pm
perrdensity = probs.distribution - simprobs.distribution
perrdist = (pm=perrpm,distribution=perrdensity,x=probs.x,type="probability")

p = SFFM.PlotSFM!(p;model=model,mesh=mesh,Dist=perrdist)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using Mustache
using Compat

function _genmodels(plot::Plot)
    bkplot = Bokehjs.Plot()
    doc = Bokehjs.uuid4()
    obs = Dict{String, BkAny}[]

    cdss = Bokehjs.ColumnDataSource[]
    renderers = Bokehjs.PlotObject[]
    legends = Tuple[]
    for datacolumn in plot.datacolumns
        cds = Bokehjs.ColumnDataSource(datacolumn.data)
        push!(cdss, cds)
        pushdict!(obs, cds, doc)
        glyph = datacolumn.glyph
        pushdict!(obs, glyph, doc)
        glyphrenderer = Bokehjs.GlyphRenderer(cds, glyph, nothing, glyph)
        push!(renderers, glyphrenderer)
        pushdict!(obs, glyphrenderer, doc)
        if datacolumn.legend != nothing
            push!(legends, (datacolumn.legend, glyphrenderer))
        end
    end

    if length(legends) > 0
        legend = Bokehjs.Legend(bkplot, legends, plot.legendsgo)
        push!(renderers, legend)
        pushdict!(obs, legend, doc)
    end

    dr1x = Bokehjs.DataRange1d(cdss, String["x"])
    dr1y = Bokehjs.DataRange1d(cdss, String["y"])
    pushdict!(obs, dr1x, doc)
    pushdict!(obs, dr1y, doc)

    # TODO: currently auto and log aren't implemented properly
    axis_types = @compat Dict{Symbol, Tuple}(
    :auto =>     (:BasicTicker, :BasicTickFormatter),
    :linear =>   (:BasicTicker, :BasicTickFormatter),
    :log =>      (:BasicTicker, :BasicTickFormatter),
    :datetime => (:DatetimeTicker, :DatetimeTickFormatter),
    )
    xticker_type, xtickform_type = axis_types[plot.x_axis_type]
    yticker_type, ytickform_type = axis_types[plot.y_axis_type]

    ticker0 = Bokehjs.Ticker(xticker_type)
    ticker1 = Bokehjs.Ticker(yticker_type)
    pushdict!(obs, ticker0, doc)
    pushdict!(obs, ticker1, doc)

    tickform0 = Bokehjs.TickFormatter(xtickform_type)
    tickform1 = Bokehjs.TickFormatter(ytickform_type)
    pushdict!(obs, tickform0, doc)
    pushdict!(obs, tickform1, doc)

    axis0 = Bokehjs.LinearAxis(0, tickform0, ticker0, bkplot)
    axis1 = Bokehjs.LinearAxis(1, tickform1, ticker1, bkplot)
    pushdict!(obs, axis0, doc)
    pushdict!(obs, axis1, doc)
    grid0 = Bokehjs.Grid(0, bkplot, ticker0)
    grid1 = Bokehjs.Grid(1, bkplot, ticker1)
    pushdict!(obs, grid0, doc)
    pushdict!(obs, grid1, doc)

    tools = Bokehjs.PlotObject[]
    if in(:pan, plot.tools)
        pantool = Bokehjs.Metatool("PanTool", bkplot, String["width", "height"])
        pushdict!(obs, pantool, doc)
        push!(tools, pantool)
    end
    if in(:wheelzoom, plot.tools)
        wheelzoomtool = Bokehjs.Metatool("WheelZoomTool", bkplot, String["width", "height"])
        pushdict!(obs, wheelzoomtool, doc)
        push!(tools, wheelzoomtool)
    end
    if in(:boxzoom, plot.tools)
        boxzoomtool = Bokehjs.Metatool("BoxZoomTool", bkplot)
        pushdict!(obs, boxzoomtool, doc)
        push!(tools, boxzoomtool)
    end
    if in(:resize, plot.tools)
        resizetool = Bokehjs.Metatool("ResizeTool", bkplot)
        pushdict!(obs, resizetool, doc)
        push!(tools, resizetool)
    end
    if in(:reset, plot.tools)
        resettool = Bokehjs.Metatool("ResetTool", bkplot)
        pushdict!(obs, resettool, doc)
        push!(tools, resettool)
    end

    push!(renderers, axis0)
    push!(renderers, axis1)
    push!(renderers, grid0)
    push!(renderers, grid1)

    axes = @compat Dict{Symbol, Vector}(
        :above => [],
        :below => [axis0],
        :left => [axis1],
        :right => []
    )

    bkplot = Bokehjs.Plot(bkplot,
                dr1x,
                dr1y,
                renderers,
                axes,
                tools,
                plot.title,
                plot.height,
                plot.width)
    pushdict!(obs, bkplot, doc)

    plotcontext = Bokehjs.PlotContext(bkplot)
    pushdict!(obs, plotcontext, doc)

    indent = DEBUG ? 2 : 0
    method_exists(json, (Dict, Int)) ? (json(obs, indent), plotcontext):
                                       (json(obs), plotcontext)
end

function _obdict(ob::Bokehjs.PlotObject, doc::Bokehjs.UUID)
    d = @compat Dict{String, BkAny}()
    d["id"] = ob.uuid
    extra_attrs = typeof(ob).names
    d["type"] = in(:_type_name, extra_attrs) ? ob._type_name : typeof(ob)
    attrs = @compat Dict{String, Any}()
    attrs["id"] = d["id"]
    attrs["doc"] = doc
    special = [:_type_name]
    for name in extra_attrs[2:end]
        in(name, special) && continue
        ob.(name) == Bokehjs.omit && continue
        key = string(name)
        # key = begingswith(key, "_") ? key[2:end] : key
        attrs[key] = ob.(name)
    end
    d["attributes"] = attrs
    return d
end

pushdict!(obs::Any, ob::Bokehjs.PlotObject, doc::Bokehjs.UUID) = push!(obs, _obdict(ob, doc))

_get_resources_dir() = Pkg.dir("Bokeh", "templates")

function _gettemplate(template::String, path::Union(String, Nothing)=nothing)
    path = path == nothing ? _get_resources_dir() : path
    fname = joinpath(path, template)
    open(readall, fname, "r")
end

function _bokehjs_paths(minified::Bool=true)
    dir = joinpath(dirname(@__FILE__), "..", "deps", "bokehjs")
    jspath = joinpath(dir, minified ? "bokeh.min.js" : "bokeh.js")
    csspath = joinpath(dir, minified ? "bokeh.min.css" : "bokeh.css")
    (jspath, csspath)
end

function _render_jscss(jspath::String, csspath::String, buildin::Bool)
    if buildin
        css = open(readall, csspath, "r")
        js = open(readall, jspath, "r")
        return "<style>\n$css\n</style>\n<script type=\"text/javascript\">\n$js\n</script>\n"
    else
        return "<link rel=\"stylesheet\" href=\"$csspath\" type=\"text/css\" />\n"*
               "<script type=\"text/javascript\" src=\"$jspath\"></script>\n"
    end
end

function _rendertemplate(models::String, plotcon::Bokehjs.PlotContext, isijulia::Bool)
    base = isijulia ? _gettemplate("ijulia.html") : _gettemplate("standalone.html")
    main = _gettemplate("main.html")
    body = _gettemplate("body.html")
    jspath, csspath = _bokehjs_paths(!DEBUG)
    builtin = isijulia || INCLUDE_JS
    jscss = _render_jscss(jspath, csspath, builtin)
    if DEBUG
        open(replace(PLOTFILE, ".html", "") * ".json", "w") do f
            print(f, models)
        end
    end
    context = @compat Dict{String, String}(
        "model_id" => string(plotcon.uuid),
        "all_models" => models,
        "div_id" => string(Bokehjs.uuid4())
    )
    main = render(main, context)
    body = render(body, context)
    maincontext = @compat Dict{String, String}(
        "jscss" => jscss,
        "main" => main,
        "body" => body,
    )
    result = render(base, maincontext)
end

function renderplot(plot::Plot, isijulia::Bool=false)
    modelsjson, plotcontext = _genmodels(plot)
    _rendertemplate(modelsjson, plotcontext, isijulia)
end

renderplot() = renderplot(CURPLOT)

function genplot(p::Plot, filename::NullString=nothing)
    filename = filename == nothing ? p.filename : filename
    html = renderplot(p)
    if ispath(filename) && FILE_WARNINGS
        if WARN_FILE != filename
            println("$(filename) already exists, overwriting")
            global WARN_FILE = filename
        end
    end
    open(filename, "w") do f
        print(f, html)
    end
end

genplot() = genplot(CURPLOT)
genplot(filename::String) = genplot(CURPLOT, filename)

max_stars_repo_name: Laeeth/dbokeh

max_issues_repo_name: Laeeth/dbokeh

max_forks_repo_name: Laeeth/dbokeh


text:
# a convenient way of interacting with a simulation
# maintained by @zsunberg

"""
    sim(polfunc::Function, mdp::MDP; [<keyword arguments>])
    sim(polfunc::Function, pomdp::POMDP; [<keyword arguments>])

Alternative way of running a simulation with a function specifying how to calculate the action at each timestep.

# Usage

    sim(mdp) do s
        # code that calculates action `a` based on `s` - this is the policy
        # you can also do other things like display something
        return a
    end

for an MDP or

    sim(pomdp) do o
        # code that calculates 'a' based on observation `o`
        # optionally you could save 'o' in a global variable or do a belief update
        return a
    end

or with a POMDP

    sim(pomdp, updater) do b
        # code that calculates 'a' based on belief `b`
        # `b` is calculated by `updater`
        return a
    end

for a POMDP and a belief updater.

# Keyword Arguments

## All Versions

- `initialstate`: the initial state for the simulation
- `simulator`: keyword argument to specify any simulator to run the simulation. If nothing is specified for the simulator, a HistoryRecorder will be used as the simulator, with all keyword arguments forwarded to it, e.g.
  ```
  sim(mdp, max_steps=100, show_progress=true) do s
      # ...
  end
  ```
  will limit the simulation to 100 steps.

## POMDP version

- `initialobs`: this will control the initial observation given to the policy function. If this is not defined, `rand(initialobs(m, s))` will be used if it is available. If it is not, `missing` will be used.

## POMDP and updater version

- `initialbelief`: `initialize_belief(updater, initialbelief)` is the first belief that will be given to the policy function.
"""
function sim end

function sim(polfunc::Function, mdp::MDP;
             initialstate=nothing,
             simulator=nothing,
             kwargs...
            )

    kwargd = Dict(kwargs)
    if initialstate==nothing && statetype(mdp) != Nothing
        initialstate = default_init_state(mdp)
    end
    if simulator==nothing
        simulator = HistoryRecorder(;kwargd...)
    end
    policy = FunctionPolicy(polfunc)
    simulate(simulator, mdp, policy, initialstate)
end

function sim(polfunc::Function, pomdp::POMDP;
             initialstate=nothing,
             simulator=nothing,
             initialobs=nothing,
             kwargs...
            )

    kwargd = Dict(kwargs)
    if initialstate==nothing && statetype(pomdp) != Nothing
        initialstate = default_init_state(pomdp)
    end
    if simulator==nothing
        simulator = HistoryRecorder(;kwargd...)
    end
    updater = PreviousObservationUpdater()
    if initialobs==nothing && obstype(pomdp) != Nothing
        initialobs = default_init_obs(pomdp, initialstate)
    end
    policy = FunctionPolicy(polfunc)
    simulate(simulator, pomdp, policy, updater, initialobs, initialstate)
end

function sim(polfunc::Function, pomdp::POMDP, updater::Updater;
             is=nothing,
             simulator=nothing,
             initialbelief=initialstate(pomdp),
             kwargs...
            )

    kwargd = Dict(kwargs)
    if is==nothing && statetype(pomdp) != Nothing
        is = default_init_state(pomdp)
    end
    if simulator==nothing
        simulator = HistoryRecorder(;kwargd...)
    end
    policy = FunctionPolicy(polfunc)
    simulate(simulator, pomdp, policy, updater, initialbelief, is)
end

function default_init_obs(p::POMDP, s)
    try
        return rand(initialobs(p, s))
    catch ex
        if ex isa MethodError
            return missing
        else
            rethrow(ex)
        end
    end
end

@generated function default_init_state(p::Union{MDP,POMDP})
    if implemented(initialstate, Tuple{p, typeof(Random.GLOBAL_RNG)})
        return :(rand(Random.GLOBAL_RNG, initialstate(p)))
    else
        return quote
            error("""
                  Error in sim(::$(typeof(p))): No initial state specified.

                  Please supply it as a keyword argument or provide POMDPs.initialstate(::$(typeof(p))).

                  """)
        end
    end
end


text:
export pcap_lookupdev

#----------
# lookup default device
#----------
function pcap_lookupdev()
    dev = ccall((:pcap_lookupdev, "libpcap"), Ptr{UInt8}, ())
    if dev == C_NULL
        return Union()
    end
    bytestring(dev)
end # function pcap_lookupdev


text:
# https://adventofcode.com/2020/day/8

# The boot code is represented as a text file with one instruction per line of text.
# Each instruction consists of an operation (`acc`, `jmp`, or `nop`) and an argument
# (a signed number like `+4` or `-20`).
#
#  * `acc` increases or decreases a single global value called the accumulator by the
#    value given in the argument. For example, `acc +7` would increase the accumulator
#    by 7. The accumulator starts at 0. After an acc instruction, the instruction
#    immediately below it is executed next.
#  * `jmp` jumps to a new instruction relative to itself. The next instruction to execute
#    is found using the argument as an offset from the jmp instruction; for example,
#    `jmp +2` would skip the next instruction, `jmp +1` would continue to the instruction
#    immediately below it, and `jmp -20` would cause the instruction 20 lines above to be
#    executed next.
#  * `nop` stands for No OPeration - it does nothing. The instruction immediately below
#    it is executed next.

using Printf

example = "
nop +0
acc +1
jmp +4
acc +3
jmp -3
acc -99
acc +1
jmp -4
acc +6
"

const commands = Dict("nop" => :nop, "acc" => :acc, "jmp" => :jmp)

struct Command
    name::Symbol
    value::Integer

    Command(name, value) = new(name, value)

    function Command(string::AbstractString)
        name, value = split(strip(string), r"\s+")
        new(commands[name], parse(Int, value))
    end
end

cmd(obj::Command) = obj.name
val(obj::Command) = obj.value

function parser(string::AbstractString)
    program = Command[]
    for row in split(string, '\n', keepempty=false)
        push!(program, Command(row))
    end
    return program
end

function run_code(program::Vector{Command}, max_loops=1000, verbose=false)
    exec_count = zeros(Int, length(program))
    status = false
    accumulator = 0
    offset = 1

    verbose && print("       <start>     =>  ")
    while true
        verbose && println("$(accumulator)")

        if offset == length(program) + 1
            status = true
            break
        elseif offset > length(program)
            verbose && println("<error>")
            break
        end
        exec_count[offset] += 1
        command = program[offset]
        verbose && @printf("%3.0f (%d): %s %+4.0f  =>  ", offset, exec_count[offset], cmd(command), val(command))

        if exec_count[offset] > max_loops
            verbose && println("<break>")
            break
        end

        if cmd(command) == :nop
            offset += 1
            continue
        end

        if cmd(command) == :acc
            accumulator += val(command)
            offset += 1
        elseif cmd(command) == :jmp
            offset += val(command)
        else
            verbose && println("<error>")
            break
        end
    end
    verbose && println()

    return accumulator, status
end

@assert run_code(parser(example), 1) == (5, false)

code_that_finishes = "
nop +0
acc +10
jmp +3
nop +0
acc +100
acc -5
"

@assert run_code(parser(code_that_finishes), 1) == (5, true)

# Run your copy of the boot code. Immediately before any instruction is executed a second
# time, what value is in the accumulator?

function part1(code)
    program = parser(code)
    accumulator, _ = run_code(program, 1)
    return accumulator
end

# Fix the program so that it terminates normally by changing exactly one `jmp` (to `nop`)
# or `nop` (to `jmp`). What is the value of the accumulator after the program terminates?

function swap(command::Command)
    if cmd(command) == :acc
        error("invalid command: $(command)")
    end
    name = cmd(command) == :nop ? :jmp : :nop
    return Command(name, val(command))
end

function part2(code, verbose=false)
    program_orig = parser(code)
    command_names = reverse([cmd(command) for command in program_orig])
    start_search = 1
    success = false
    accumulator = nothing

    while !success
        program = copy(program_orig)
        idx = findfirst(x -> x in (:nop, :jmp), command_names[start_search:end])
        idx += start_search - 1
        # because we search starting from the back
        program[end - idx + 1] = swap(program_orig[end - idx + 1])
        accumulator, success = run_code(program, 2)
        start_search = idx + 1

        verbose && @printf("Progress: %.2f%%\n", idx / length(program) * 100)
    end

    return accumulator
end

@assert part2(example) == 8

test = read("data/day-08.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 1832
@assert result2 == 662

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
function onehot(l, labels)
    i = something(findfirst(isequal(l), labels), 0)
    i > 0 || error("Value $l is not in labels")
    x = zeros(Int8, length(labels))'
    x[i] = 1
    x
end

function onehot(l::AbstractArray, labels)
    m = Matrix{Int8}(undef, length(labels), length(l))
    for (i, elem) in enumerate(l)
        m[:, i] = onehot(elem, labels)
    end
    m
end

function onehot(l)
    u = unique(l)
    onehot(l, minimum(u):maximum(u))
end

function onecold(l::AbstractMatrix, labels)
    m = Matrix{eltype(labels)}(undef, 1, size(l, 2))
    for (i, col) in enumerate(eachcol(l))
        m[1, i] = labels[argmax(col)]
    end
    m
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

@testset "Universal Generating Function" begin

    @testset "Optional Reduction" begin
        ugfᵍᵉⁿ = UGF(:power, [0.0u"MW",0.0u"MW",2.0u"MW"], [0.1,0.2,0.7])
        @test isapprox(ugfᵍᵉⁿ.val, [0.0u"MW",2.0u"MW"], rtol=1e-6)
        @test isapprox(ugfᵍᵉⁿ.prb, [0.3,0.7], rtol=1e-6)

        ugfᵍᵉⁿ = UGF(:power, [0.0u"MW",0.0u"MW",2.0u"MW"], [0.1,0.2,0.7], rdc=false)
        @test isapprox(ugfᵍᵉⁿ.val, [0.0u"MW",0.0u"MW",2.0u"MW"], rtol=1e-6)
        @test isapprox(ugfᵍᵉⁿ.prb, [0.1, 0.2, 0.7], rtol=1e-6)
    end

    @testset "STD → UGF" begin
        stdᵍᵉⁿ = solvedSTD(prob  = [0.1,0.2,0.7],
                           power = [0.0u"MW",0.0u"MW",2.0u"MW"])
        ugfᵍᵉⁿ = UGF(:power, stdᵍᵉⁿ)
        @test isapprox(ugfᵍᵉⁿ.val, [0.0u"MW",2.0u"MW"], rtol=1e-6)
        @test isapprox(ugfᵍᵉⁿ.prb, [0.3,0.7], rtol=1e-6)
    end
end


text:
#
# Helper
#
@inline _extract_val(::Val{T}) where {T} = T

#! format: off
# turn formatting for for the following functions
# due to the if with returns inside (formatter puts a return upfront the if)
function _split_signature(sig::Expr)
    if sig.head == :where
        where_exprs = sig.args[2:end]
        call_expr = sig.args[1]
    elseif sig.head == :call
        where_exprs = []
        call_expr = sig
    else
        error("Incorrect syntax in $ex. Expected a :where or :call expression.")
    end
    fname = call_expr.args[1]
    if isa(call_expr.args[2], Expr) && call_expr.args[2].head == :parameters
        # we have keyword arguments
        callargs = call_expr.args[3:end]
        kwargs_list = call_expr.args[2].args
    else
        callargs = call_expr.args[2:end]
        kwargs_list = []
    end
    argnames = map(callargs) do arg
        if isa(arg, Expr)
            return arg.args[1]
        else
            return arg
        end
    end
    argtypes = map(callargs) do arg
        if isa(arg, Expr)
            return arg.args[2]
        else
            return Any
        end
    end

    kwargs_call = map(kwargs_list) do kwarg
        if kwarg.head === :...
            return kwarg
        else
            if isa(kwarg.args[1], Symbol)
                kwargname = kwarg.args[1]
            else
                kwargname = kwarg.args[1].args[1]
            end
            return :($kwargname = $kwargname)
        end
    end

    return (;
        fname = fname,
        where_exprs = where_exprs,
        callargs = callargs,
        kwargs_list = kwargs_list,
        argnames = argnames,
        argtypes = argtypes,
        kwargs_call = kwargs_call,
        fname__parent = Symbol(string(fname) * "__parent"),
        fname__transparent = Symbol(string(fname) * "__transparent"),
        fname__intransparent = Symbol(string(fname) * "__intransparent"),
    )
end
#! format: on
function _split_function(ex::Expr)
    if ex.head == :function
        sig = ex.args[1]
        body = ex.args[2]
    else
        error("Incorrect syntax in $ex. Expected :function.")
    end

    return (; body = body, _split_signature(sig)...)
end

#
# Type
#
"""
    AbstractDecoratorManifold{𝔽} <: Manifold{𝔽}

An `AbstractDecoratorManifold` indicates that to some extent a manifold subtype
decorates another [`Manifold`](@ref) in the sense that it either

* it extends the functionality of a manifold with further features
* it defines a new manifold that internally uses functions from the decorated manifold

with the main intent that several or most functions of [`Manifold`](@ref) are transparently
passed through to the manifold that is decorated. This way a function implemented for a
decorator acts transparent on all other decorators, i.e. they just pass them through. If
the decorator the function is implemented for is not among the decorators, an error is
issued. By default all base manifold functions, for example [`exp`](@ref) and [`log`](@ref)
are transparent for all decorators.

Transparency of functions with respect to decorators can be specified using the macros
[`@decorator_transparent_fallback`](@ref), [`@decorator_transparent_function`](@ref) and
[`@decorator_transparent_signature`](@ref).
"""
abstract type AbstractDecoratorManifold{𝔽} <: Manifold{𝔽} end

#
# Macros
#
"""
    @decorator_transparent_fallback(ex)
    @decorator_transparent_fallback(fallback_case = :intransparent, ex)

This macro introduces an additional implementation for a certain additional case.
This can especially be used if for an already transparent function and an abstract
intermediate type a change in the default is required.
For implementing a concrete type, neither this nor any other trick is necessary. One
just implements the function as before. Note that a decorator that [`is_default_decorator`](@ref)
still dispatches to the transparent case.

* `:transparent` states, that the function is transparently passed on to the manifold that
  is decorated by the [`AbstractDecoratorManifold`](@ref) `M`, which is determined using
  the function [`decorated_manifold`](@ref).
* `:intransparent` states that an implementation for this decorator is required, and if
  none of the types provides one, an error is issued. Since this macro provides such an
  implementation, this is the default.
* `:parent` states, that this function passes on to the supertype instead of to the
  decorated manifold.

Inline definitions are not supported. The function signature however may contain
keyword arguments and a where clause. It does not allow for parameters with default values.

# Examples

```julia
@decorator_transparent_fallback function log!(M::AbstractGroupManifold, X, p, q)
    log!(decorated_manifold(M), X, p, Q)
end
@decorator_transparent_fallback :transparent function log!(M::AbstractGroupManifold, X, p, q)
    log!(decorated_manifold(M), X, p, Q)
end
```
"""
macro decorator_transparent_fallback(ex)
    return esc(quote
        @decorator_transparent_fallback :intransparent ($ex)
    end)
end
macro decorator_transparent_fallback(fallback_case, input_ex)
    ex = macroexpand(__module__, input_ex)
    parts = _split_function(ex)
    callargs = parts[:callargs]
    where_exprs = parts[:where_exprs]
    fname_fallback = Symbol(string(parts.fname) * "__" * string(fallback_case)[2:end])
    return esc(
        quote
            function ($(fname_fallback))(
                $(callargs...);
                $(parts[:kwargs_list]...),
            ) where {$(where_exprs...)}
                return ($(parts[:body]))
            end
        end,
    )
end

"""
    @decorator_transparent_function(ex)
    @decorator_transparent_function(fallback_case = :intransparent, ex)

Introduce the function specified by `ex` to act transparently with respect to
[`AbstractDecoratorManifold`](@ref)s. This introduces the possibility to modify the kind of
transparency the implementation is done for. This optional first argument, the `Symbol`
within `fallback_case`. This macro can be used to define a function and introduce it as
transparent to other decorators. Note that a decorator that [`is_default_decorator`](@ref)
still dispatches to the transparent case.

The cases of transparency are

* `:transparent` states, that the function is transparently passed on to the manifold that
  is decorated by the [`AbstractDecoratorManifold`](@ref) `M`, which is determined using
  the function [`decorated_manifold`](@ref).
* `:intransparent` states that an implementation for this decorator is required, and if
  none of the types provides one, an error is issued. Since this macro provides such an
  implementation, this is the default.
* `:parent` states, that this function passes on to the supertype instead of to the
  decorated manifold. Passing is performed using the `invoke` function where the type of
  manifold is replaced by its supertype.

Innkoline-definitions are not yet covered – the function signature however may contain
keyword arguments and a where clause.

# Examples

```julia
@decorator_transparent_function log!(M::AbstractDecoratorManifold, X, p, q)
    log!(decorated_manifold(M), X, p, Q)
end
@decorator_transparent_function :parent log!(M::AbstractDecoratorManifold, X, p, q)
    log!(decorated_manifold(M), X, p, Q)
end
```
"""
macro decorator_transparent_function(ex)
    return esc(quote
        @decorator_transparent_function :intransparent ($ex)
    end)
end
macro decorator_transparent_function(fallback_case, input_ex)
    ex = macroexpand(__module__, input_ex)
    parts = _split_function(ex)
    kwargs_list = parts[:kwargs_list]
    callargs = parts[:callargs]
    fname = parts[:fname]
    where_exprs = parts[:where_exprs]
    body = parts[:body]
    argnames = parts[:argnames]
    argtypes = parts[:argtypes]
    kwargs_call = parts[:kwargs_call]
    fname_fallback = Symbol(string(parts.fname) * "__" * string(fallback_case)[2:end])

    return esc(
        quote
            function ($fname)(
                $(argnames[1])::AbstractDecoratorManifold,
                $(callargs[2:end]...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                transparency = ManifoldsBase._acts_transparently($fname, $(argnames...))
                return if transparency === Val(:parent)
                    return ($(parts.fname__parent))($(argnames...); $(kwargs_call...))
                elseif transparency === Val(:transparent)
                    return ($(parts.fname__transparent))($(argnames...); $(kwargs_call...))
                elseif transparency === Val(:intransparent)
                    return ($(parts.fname__intransparent))(
                        $(argnames...);
                        $(kwargs_call...),
                    )
                else
                    error("incorrect transparency: $transparency")
                end
            end
            function ($(parts[:fname__transparent]))(
                $(argnames[1])::AbstractDecoratorManifold,
                $(callargs[2:end]...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return ($fname)(
                    ManifoldsBase.decorated_manifold($(argnames[1])),
                    $(argnames[2:end]...);
                    $(kwargs_call...),
                )
            end
            function ($(parts[:fname__intransparent]))(
                $(argnames[1])::AbstractDecoratorManifold,
                $(callargs[2:end]...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                error_msg = ManifoldsBase.manifold_function_not_implemented_message(
                    $(argnames[1]),
                    $fname,
                    $(argnames[2:end]...),
                )
                return error(error_msg)
            end
            function ($fname)(
                $(argnames[1])::Manifold,
                $(callargs[2:end]...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return error(
                    string(
                        ManifoldsBase.manifold_function_not_implemented_message(
                            $(argnames[1]),
                            $fname,
                            $(argnames[2:end]...),
                        ),
                        " Usually this is implemented for a ",
                        $(argtypes[1]),
                        ". Maybe you missed to implement this function for a default?",
                    ),
                )
            end
            function ($(parts[:fname__parent]))(
                $(argnames[1])::AbstractDecoratorManifold,
                $(callargs[2:end]...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return invoke(
                    $fname,
                    Tuple{supertype($(argtypes[1])),$(argtypes[2:end]...)},
                    $(argnames...);
                    $(kwargs_call...),
                )
            end
            function ($fname_fallback)(
                $(callargs[1]),
                $(callargs[2:end]...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return ($body)
            end
            function decorator_transparent_dispatch(
                ::typeof($fname),
                $(callargs...),
            ) where {$(where_exprs...)}
                return Val($fallback_case)
            end
        end,
    )
end
#! format: off
# due to the if with returns inside (formatter puts a return upfront the if)
"""
    @decorator_transparent_signature(ex)

Introduces a given function to be transparent with respect to all decorators.
The function is adressed by its signature in `ex`.

Supports standard, keyword arguments and `where` clauses. Doesn't support parameters with
default values. It introduces a dispatch on several transparency modes

The cases of transparency are

* `:transparent` states, that the function is transparently passed on to the manifold that
  is decorated by the [`AbstractDecoratorManifold`](@ref) `M`, which is determined using
  the function [`decorated_manifold`](@ref). This is the default.
* `:intransparent` states that an implementation for this decorator is required, and if
  none of the types provides one, an error is issued.
* `:parent` states, that this function passes on to the supertype instead of to the
  decorated manifold.

Inline definitions are not supported. The function signature however may contain
keyword arguments and a where clause.

The dispatch kind can later still be set to something different, see [`decorator_transparent_dispatch`](@ref)

# Examples:

```julia
@decorator_transparent_signature log!(M::AbstractDecoratorManifold, X, p, q)
@decorator_transparent_signature log!(M::TD, X, p, q) where {TD<:AbstractDecoratorManifold}
@decorator_transparent_signature isapprox(M::AbstractDecoratorManifold, p, q; kwargs...)
```
"""
macro decorator_transparent_signature(ex)
    parts = _split_signature(ex)
    kwargs_list = parts[:kwargs_list]
    callargs = parts[:callargs]
    fname = parts[:fname]
    where_exprs = parts[:where_exprs]
    argnames = parts[:argnames]
    argtypes = parts[:argtypes]
    kwargs_call = parts[:kwargs_call]
    #! format: off
    return esc(
        quote
            function ($fname)($(callargs...); $(kwargs_list...)) where {$(where_exprs...)}
                transparency = ManifoldsBase._acts_transparently($fname, $(argnames...))
                if transparency === Val(:parent)
                    return ($(parts.fname__parent))($(argnames...); $(kwargs_call...))
                elseif transparency === Val(:transparent)
                    return ($(parts.fname__transparent))($(argnames...); $(kwargs_call...))
                elseif transparency === Val(:intransparent)
                    return ($(parts.fname__intransparent))(
                        $(argnames...);
                        $(kwargs_call...),
                    )
                else
                    error("incorrect transparency: $transparency")
                end
            end
            function ($(parts[:fname__transparent]))(
                $(callargs...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return ($fname)(
                    ManifoldsBase.decorated_manifold($(argnames[1])),
                    $(argnames[2:end]...);
                    $(kwargs_call...),
                )
            end
            function ($(parts[:fname__intransparent]))(
                $(callargs...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                error_msg = ManifoldsBase.manifold_function_not_implemented_message(
                    $(argnames[1]),
                    $fname,
                    $(argnames[2:end]...),
                )
                return error(error_msg)
            end
            function ($(parts[:fname__parent]))(
                $(callargs...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return invoke(
                    $fname,
                    Tuple{supertype($(argtypes[1])),$(argtypes[2:end]...)},
                    $(argnames...);
                    $(kwargs_call...),
                )
            end
        end,
    )
end
#! format: on

#
# Functions
#

"""
    is_default_decorator(M) -> Bool

For any manifold that is a subtype of [`AbstractDecoratorManifold`](@ref), this function
indicates whether a certain manifold `M` acts as a default decorator.

This yields that _all_ functions are passed through to the decorated [`Manifold`](@ref)
if `M` is indicated as default. This overwrites all [`is_decorator_transparent`](@ref)
values.

This yields the following advantange: For a manifold one usually implicitly assumes for
example a metric. To avoid reimplementation of this metric when introducing a second metric,
the first metric can be set to be the default, i.e. its implementaion is already given by
the undecorated case.

Value returned by this function is determined by [`default_decorator_dispatch`](@ref),
which returns a `Val`-wrapped boolean for type stability of certain functions.
"""
is_default_decorator(M::Manifold) = _extract_val(default_decorator_dispatch(M))

"""
    default_decorator_dispatch(M) -> Val

Return whether by default to dispatch the the inner manifold of
a decorator (`Val(true)`) or not (`Val(false`). For more details see
[`is_decorator_transparent`](@ref).
"""
default_decorator_dispatch(M::Manifold) = Val(false)

"""
    is_decorator_transparent(f, M::Manifold, args...) -> Bool

Given a [`Manifold`](@ref) `M` and a function `f(M, args...)`, indicate, whether an
[`AbstractDecoratorManifold`](@ref) acts transparently for `f`. This means, it
just passes through down to the internally stored manifold.
Transparency is only defined for decorator manifolds and by default all decorators are transparent.
A function that is affected by the decorator indicates this by returning `false`. To change
this behaviour, see [`decorator_transparent_dispatch`](@ref).

If a decorator manifold is not in general transparent, it might still pass down
for the case that a decorator is the default decorator, see [`is_default_decorator`](@ref).
"""
function is_decorator_transparent(f, M::Manifold, args...)
    return decorator_transparent_dispatch(f, M, args...) === Val(:transparent)
end

"""
    decorator_transparent_dispatch(f, M::Manifold, args...) -> Val

Given a [`Manifold`](@ref) `M` and a function `f(M,args...)`, indicate, whether a
function is `Val(:transparent)` or `Val(:intransparent)` for the (decorated)
[`Manifold`](@ref) `M`. Another possibility is, that for `M` and given `args...`
the function `f` should invoke `M`s `Val(:parent)` implementation, see
[`@decorator_transparent_function`](@ref) for details.
"""
decorator_transparent_dispatch(f, M::Manifold, args...) = Val(:transparent)

function _acts_transparently(f, M::Manifold, args...)
    return _val_or(
        default_decorator_dispatch(M),
        decorator_transparent_dispatch(f, M, args...),
    )
end

_val_or(::Val{true}, ::Val{T}) where {T} = Val(:transparent)
_val_or(::Val{false}, val::Val) = val

#
# Functions overwritten with decorators
#

function base_manifold(M::AbstractDecoratorManifold, depth::Val{N} = Val(-1)) where {N}
    N == 0 && return M
    N < 0 && return base_manifold(decorated_manifold(M), depth)
    return base_manifold(decorated_manifold(M), Val(N - 1))
end

@decorator_transparent_signature check_manifold_point(
    M::AbstractDecoratorManifold,
    p;
    kwargs...,
)

@decorator_transparent_signature check_tangent_vector(
    M::AbstractDecoratorManifold,
    p,
    X;
    kwargs...,
)

"""
    decorated_manifold(M::AbstractDecoratorManifold)

Return the manifold decorated by the decorator `M`. Defaults to `M.manifold`.
"""
decorated_manifold(M::Manifold) = M.manifold

@decorator_transparent_signature distance(M::AbstractDecoratorManifold, p, q)

@decorator_transparent_signature embed(M::AbstractDecoratorManifold, p, X)
@decorator_transparent_signature embed(M::AbstractDecoratorManifold, p)

@decorator_transparent_signature embed!(M::AbstractDecoratorManifold, q, p)
@decorator_transparent_signature embed!(M::AbstractDecoratorManifold, Y, p, X)

@decorator_transparent_signature exp(M::AbstractDecoratorManifold, p, X)

@decorator_transparent_signature exp!(M::AbstractDecoratorManifold, q, p, X)

@decorator_transparent_signature injectivity_radius(M::AbstractDecoratorManifold)
@decorator_transparent_signature injectivity_radius(M::AbstractDecoratorManifold, p)
@decorator_transparent_signature injectivity_radius(
    M::AbstractDecoratorManifold,
    m::AbstractRetractionMethod,
)
@decorator_transparent_signature injectivity_radius(
    M::AbstractDecoratorManifold,
    m::ExponentialRetraction,
)
@decorator_transparent_signature injectivity_radius(
    M::AbstractDecoratorManifold,
    p,
    m::AbstractRetractionMethod,
)
@decorator_transparent_signature injectivity_radius(
    M::AbstractDecoratorManifold,
    p,
    m::ExponentialRetraction,
)

@decorator_transparent_signature inner(M::AbstractDecoratorManifold, p, X, Y)

@decorator_transparent_signature inverse_retract(
    M::AbstractDecoratorManifold,
    p,
    q,
    m::AbstractInverseRetractionMethod,
)

@decorator_transparent_signature inverse_retract(
    M::AbstractDecoratorManifold,
    p,
    q,
    m::LogarithmicInverseRetraction,
)

@decorator_transparent_signature inverse_retract!(
    M::AbstractDecoratorManifold,
    X,
    p,
    q,
    m::AbstractInverseRetractionMethod,
)

@decorator_transparent_signature inverse_retract!(
    M::AbstractDecoratorManifold,
    X,
    p,
    q,
    m::LogarithmicInverseRetraction,
)

@decorator_transparent_signature isapprox(M::AbstractDecoratorManifold, p, q; kwargs...)
@decorator_transparent_signature isapprox(M::AbstractDecoratorManifold, p, X, Y; kwargs...)

@decorator_transparent_signature log(M::AbstractDecoratorManifold, p, q)
@decorator_transparent_signature log!(M::AbstractDecoratorManifold, X, p, q)

@decorator_transparent_signature manifold_dimension(M::AbstractDecoratorManifold)

@decorator_transparent_signature mid_point(M::AbstractDecoratorManifold, p1, p2)
@decorator_transparent_signature mid_point!(M::AbstractDecoratorManifold, q, p1, p2)

@decorator_transparent_signature number_system(M::AbstractDecoratorManifold)

@decorator_transparent_signature project(M::AbstractDecoratorManifold, p)
@decorator_transparent_signature project!(M::AbstractDecoratorManifold, q, p)

@decorator_transparent_signature project(M::AbstractDecoratorManifold, p, X)
@decorator_transparent_signature project!(M::AbstractDecoratorManifold, Y, p, X)

@decorator_transparent_signature representation_size(M::AbstractDecoratorManifold)

@decorator_transparent_signature retract(
    M::AbstractDecoratorManifold,
    p,
    X,
    m::AbstractRetractionMethod,
)

@decorator_transparent_signature retract(
    M::AbstractDecoratorManifold,
    p,
    X,
    m::ExponentialRetraction,
)

@decorator_transparent_signature retract!(
    M::AbstractDecoratorManifold,
    q,
    p,
    X,
    m::AbstractRetractionMethod,
)

@decorator_transparent_signature retract!(
    M::AbstractDecoratorManifold,
    q,
    p,
    X,
    m::ExponentialRetraction,
)

@decorator_transparent_signature vector_transport_along(
    M::AbstractDecoratorManifold,
    p,
    X,
    c,
)
@decorator_transparent_signature vector_transport_along!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    c,
)
@decorator_transparent_signature vector_transport_along!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::AbstractVectorTransportMethod,
)
@decorator_transparent_signature vector_transport_along!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::PoleLadderTransport,
)
@decorator_transparent_signature vector_transport_along!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::SchildsLadderTransport,
)

@decorator_transparent_signature vector_transport_direction(
    M::AbstractDecoratorManifold,
    p,
    X,
    d,
)
@decorator_transparent_signature vector_transport_direction!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    d,
)

@decorator_transparent_signature vector_transport_to(
    M::AbstractDecoratorManifold,
    p,
    X,
    q,
    m::AbstractVectorTransportMethod,
)
@decorator_transparent_signature vector_transport_to!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    q,
    m::AbstractVectorTransportMethod,
)
@decorator_transparent_signature vector_transport_to!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    q,
    m::ProjectionTransport,
)
@decorator_transparent_signature vector_transport_to!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    q,
    m::PoleLadderTransport,
)
@decorator_transparent_signature vector_transport_to!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    q,
    m::ScaledVectorTransport,
)
@decorator_transparent_signature vector_transport_to!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    X,
    q,
    m::SchildsLadderTransport,
)

@decorator_transparent_signature zero_tangent_vector(M::AbstractDecoratorManifold, p)
@decorator_transparent_signature zero_tangent_vector!(M::AbstractDecoratorManifold, X, p)


text:
@testset "Load Localizations" begin

    @testset "Nikon Text" begin
        expected = Localization[Localization(1,  "647",  8394.1,  7259.0, 0.0, 10.97216,    1,  1),
                                Localization(2,  "647", 22137.7, 13796.5, 0.0, 12.56947, 2114,  1),
                                Localization(3,  "647", 36718.1, 28336.0, 0.0, 13.24367, 3160,  1),
                                Localization(4,  "647", 13799.5, 26063.8, 0.0, 13.50408, 7444,  1),
                                Localization(5,  "647",  5691.9, 26318.1, 0.0,  8.73353, 11009, 2),
                                Localization(6,  "488",  7934.0, 10576.6, 0.0,  4.75511, 11011, 3),
                                Localization(7,  "488", 11612.8, 20987.3, 0.0, 12.5259,  11135, 1),
                                Localization(8,  "488", 12917.9, 19160.3, 0.0, 16.88895, 14050, 1),
                                Localization(9,  "488", 38189.6, 22873.5, 0.0, 14.11827, 18538, 1),
                                Localization(10, "488", 24023.7, 16153.2, 0.0,  9.34327, 22014, 2)]

        localizations = loadlocalizations("nikontestdata.bin.txt", nikonelementstext)

        @test all(expected .== localizations)
    end

end


text:
# # Themes

"""
Namespace for all themes that are provided with the package.

Themes defined here should be 0-arity functions that return the absolute path
to the folder containing the theme definition.

Provided themes are currently:

  - [`default`](# "`Themes.default`")
"""
module Themes

using Pkg.Artifacts

"""
The "default" theme used by this package.
"""
default() = artifact"publish_theme_default"

end

# ## Theme Loader

"""
    loadtheme(tree, env)

Searches the given `env` configuration for a Publish.jl theme and loads it into
the provided `tree`. Since `tree` is immutable the updated version is returned
along with a new `env` containing any merged data from the theme.

`loadtheme` looks for a `publish.theme` key in the `env` and tries to call the
function defined by the value of the key.

The themes included with this package are listed in the [`Themes`](#) module.
"""
function loadtheme(tree::FileTree, env::AbstractDict)
    theme = env["publish"]["theme"]
    file = joinpath(path(tree), theme, "Theme.toml")
    if !isfile(file)
        bind = binding(theme)
        dir = resolve_theme(Docs.defined(bind) ? Docs.resolve(bind) : Themes.default)
        file = joinpath(dir, "Theme.toml")
        isfile(file) || error("'$theme' theme does not exist.")
    end
    dict = TOML.parsefile(string(file))
    env = rmerge(env, Dict("publish" => dict))
    tree′ = FileTree(dirname(file))
    tree′ = FileTrees.load(tree′; lazy=LAZY[]) do file
        loadfile(env, joinpath(basename(tree′), path(file)))
    end
    return merge(tree, rename(tree′, basename(tree))), env
end

resolve_theme(f::Function) = f()
resolve_theme(s::AbstractString) = s


text:
struct CartPoleParams{T}
    gravity::T
    masscart::T
    masspole::T
    totalmass::T
    halflength::T
    polemasslength::T
    forcemag::T
    tau::T
    thetathreshold::T
    xthreshold::T
    maxsteps::Int64
end

mutable struct CartPole{T} <: AbstractEnv
    params::CartPoleParams{T}
    actionspace::DiscreteSpace
    observationspace::BoxSpace{T}
    state::Array{T, 1}
    action::Int
    done::Bool
    t::Int64
end

function CartPole(; T = Float64, gravity = T(9.8), masscart = T(1.),
                  masspole = T(.1), halflength = T(.5), forcemag = T(10.),
                  maxsteps = 200)
    params = CartPoleParams(gravity, masscart, masspole, masscart + masspole,
                            halflength, masspole * halflength, forcemag,
                            T(.02), T(2 * 12 * π /360), T(2.4), maxsteps)
    high = [2 * params.xthreshold, T(1e38),
            2 * params.thetathreshold, T(1e38)]
    cp = CartPole(params, DiscreteSpace(2, 1), BoxSpace(-high, high),
                  zeros(T, 4), 2, false, 0)
    reset!(cp)
    cp
end

actionspace(env::CartPole) = env.actionspace

function reset!(env::CartPole{T}) where T <: Number
    env.state[:] = T(.1) * rand(T, 4) .- T(.05)
    env.t = 0
    env.action = 2
    env.done = false
    (observation=env.state,)
end

getstate(env::CartPole) = (observation=env.state, isdone=env.done)

function interact!(env::CartPole{T}, a) where T <: Number
    env.action = a
    env.t += 1
    force = a == 2 ? env.params.forcemag : -env.params.forcemag
    x, xdot, theta, thetadot = env.state
    costheta = cos(theta)
    sintheta = sin(theta)
    tmp = (force + env.params.polemasslength * thetadot^2 * sintheta) /
        env.params.totalmass
    thetaacc = (env.params.gravity * sintheta - costheta * tmp) /
        (env.params.halflength *
            (4/3 - env.params.masspole * costheta^2/env.params.totalmass))
    xacc = tmp - env.params.polemasslength * thetaacc * costheta /
        env.params.totalmass
    env.state[1] += env.params.tau * xdot
    env.state[2] += env.params.tau * xacc
    env.state[3] += env.params.tau * thetadot
    env.state[4] += env.params.tau * thetaacc
    env.done = abs(env.state[1]) > env.params.xthreshold ||
               abs(env.state[3]) > env.params.thetathreshold ||
               env.t >= env.params.maxsteps
    (observation=env.state, reward=1., isdone=env.done)
end

function plotendofepisode(x, y, d)
    if d
        setmarkercolorind(7)
        setmarkertype(-1)
        setmarkersize(6)
        polymarker([x], [y])
    end
    return nothing
end
function plotenv(env::CartPole)
    s, a, d = env.state, env.action, env.done
    x, xdot, theta, thetadot = s
    l = 2 * env.params.halflength
    clearws()
    setviewport(0, 1, 0, 1)
    xthreshold = env.params.xthreshold
    setwindow(-xthreshold, xthreshold, -.1, l + .1)
    fillarea([x-.5, x-.5, x+.5, x+.5], [-.05, 0, 0, -.05])
    setlinecolorind(4)
    setlinewidth(3)
    polyline([x, x + l * sin(theta)], [0, l * cos(theta)])
    setlinecolorind(2)
    drawarrow(x + (a == 1) - .5, -.025, x + 1.4 * (a==1) - .7, -.025)
    plotendofepisode(xthreshold - .2, l, d)
    updatews()
end


text:
using Xpress
using LinearAlgebra
using MathOptInterface

const MOI = MathOptInterface

function callback_simple_model()
    model = Xpress.Optimizer(
        PRESOLVE = 0,
        CUTSTRATEGY = 0,
        HEURSTRATEGY = 0,
        SYMMETRY = 0,
        OUTPUTLOG = 0
    )

    MOI.Utilities.loadfromstring!(model, """
        variables: x, y
        maxobjective: y
        c1: x in Integer()
        c2: y in Integer()
        c3: x in Interval(0.0, 2.5)
        c4: y in Interval(0.0, 2.5)
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    return model, x, y
end

model, x, y = callback_simple_model()

data = Matrix(1.0I, 3, 3)

function foo(cb::Xpress.CallbackData)

    cb.data[1] = 98

    cols = Xpress.getintattrib(cb.model,Xpress.Lib.XPRS_COLS)
    rows = Xpress.getintattrib(cb.model,Xpress.Lib.XPRS_ROWS)
    Xpress.getdblattrib(cb.model, Xpress.Lib.XPRS_BESTBOUND)

    ans_variable_primal = Vector{Float64}(undef,Int(cols))
    ans_linear_primal = Vector{Float64}(undef,Int(cols))

    Xpress.Lib.XPRSgetlpsol(cb.model,
        ans_variable_primal,
        ans_linear_primal,
        C_NULL, C_NULL)

    return
end

func_ptr, data_ptr = Xpress.set_callback_preintsol!(model.inner, foo, data)

@test data[1] == 1
MOI.optimize!(model)
@test data[1] == 98

@test typeof(data_ptr) <: Any
@test typeof(func_ptr) <: Ptr{Cvoid}


text:
# tests from sandra kubler, ryan mcdonald, joakim nivre 09 "dependency
# parsing" (https://doi.org/10.2200/S00169ED1V01Y200901HLT002)

@testset "Kübler et al 09" begin
    @test projective_only(ArcEager())

    @testset "Figure 3.7" begin
        tree = test_sentence("economicnews.conll")
        oracle = StaticOracle(ArcEager(), arc=typed)

        gold_transitions = [Shift(), LeftArc("ATT"), Shift(), LeftArc("SBJ"),
                            RightArc("PRED"), Shift(), LeftArc("ATT"),
                            RightArc("OBJ"), RightArc("ATT"), Shift(),
                            LeftArc("ATT"), RightArc("PC"), Reduce(),
                            Reduce(), Reduce(), RightArc("PU")]

        cfg = initconfig(ArcEager(), tree)

        for t in gold_transitions
            @test !isfinal(cfg)
            @test oracle.oracle_function(cfg, tree, typed) == t
            cfg = t(cfg)
        end
        @test isfinal(cfg)
        result = DependencyTree(tokens(cfg))
        @test tree == result

        oracle = StaticOracle(ArcEager(), arc=typed)
        pairs = xys(oracle, tree)
        @test collect(last.(pairs)) == gold_transitions
        @test collect(xys(oracle, [tree])) == collect(pairs)
    end
end