text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ._ProtoBuf_Top_.opentelemetry

const SeverityNumber = (;
    [
        Symbol("SEVERITY_NUMBER_UNSPECIFIED") => Int32(0),
        Symbol("SEVERITY_NUMBER_TRACE") => Int32(1),
        Symbol("SEVERITY_NUMBER_TRACE2") => Int32(2),
        Symbol("SEVERITY_NUMBER_TRACE3") => Int32(3),
        Symbol("SEVERITY_NUMBER_TRACE4") => Int32(4),
        Symbol("SEVERITY_NUMBER_DEBUG") => Int32(5),
        Symbol("SEVERITY_NUMBER_DEBUG2") => Int32(6),
        Symbol("SEVERITY_NUMBER_DEBUG3") => Int32(7),
        Symbol("SEVERITY_NUMBER_DEBUG4") => Int32(8),
        Symbol("SEVERITY_NUMBER_INFO") => Int32(9),
        Symbol("SEVERITY_NUMBER_INFO2") => Int32(10),
        Symbol("SEVERITY_NUMBER_INFO3") => Int32(11),
        Symbol("SEVERITY_NUMBER_INFO4") => Int32(12),
        Symbol("SEVERITY_NUMBER_WARN") => Int32(13),
        Symbol("SEVERITY_NUMBER_WARN2") => Int32(14),
        Symbol("SEVERITY_NUMBER_WARN3") => Int32(15),
        Symbol("SEVERITY_NUMBER_WARN4") => Int32(16),
        Symbol("SEVERITY_NUMBER_ERROR") => Int32(17),
        Symbol("SEVERITY_NUMBER_ERROR2") => Int32(18),
        Symbol("SEVERITY_NUMBER_ERROR3") => Int32(19),
        Symbol("SEVERITY_NUMBER_ERROR4") => Int32(20),
        Symbol("SEVERITY_NUMBER_FATAL") => Int32(21),
        Symbol("SEVERITY_NUMBER_FATAL2") => Int32(22),
        Symbol("SEVERITY_NUMBER_FATAL3") => Int32(23),
        Symbol("SEVERITY_NUMBER_FATAL4") => Int32(24),
    ]...
)

const LogRecordFlags = (;
    [
        Symbol("LOG_RECORD_FLAG_UNSPECIFIED") => Int32(0),
        Symbol("LOG_RECORD_FLAG_TRACE_FLAGS_MASK") => Int32(255),
    ]...
)

mutable struct LogRecord <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function LogRecord(; kwargs...)
        obj = new(meta(LogRecord), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct LogRecord
const __meta_LogRecord = Ref{ProtoMeta}()
function meta(::Type{LogRecord})
    ProtoBuf.metalock() do
        if !isassigned(__meta_LogRecord)
            __meta_LogRecord[] = target = ProtoMeta(LogRecord)
            fnum = Int[1, 11, 2, 3, 4, 5, 6, 7, 8, 9, 10]
            wtype = Dict(
                :time_unix_nano => :fixed64,
                :observed_time_unix_nano => :fixed64,
                :flags => :fixed32,
            )
            allflds = Pair{Symbol,Union{Type,String}}[
                :time_unix_nano=>UInt64,
                :observed_time_unix_nano=>UInt64,
                :severity_number=>Int32,
                :severity_text=>AbstractString,
                :name=>AbstractString,
                :body=>opentelemetry.proto.common.v1.AnyValue,
                :attributes=>Base.Vector{opentelemetry.proto.common.v1.KeyValue},
                :dropped_attributes_count=>UInt32,
                :flags=>UInt32,
                :trace_id=>Vector{UInt8},
                :span_id=>Vector{UInt8},
            ]
            meta(
                target,
                LogRecord,
                allflds,
                ProtoBuf.DEF_REQ,
                fnum,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                wtype,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_LogRecord[]
    end
end
function Base.getproperty(obj::LogRecord, name::Symbol)
    if name === :time_unix_nano
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :observed_time_unix_nano
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :severity_number
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :severity_text
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :body
        return (obj.__protobuf_jl_internal_values[name])::opentelemetry.proto.common.v1.AnyValue
    elseif name === :attributes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            opentelemetry.proto.common.v1.KeyValue,
        }
    elseif name === :dropped_attributes_count
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :flags
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :trace_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :span_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct InstrumentationLibraryLogs <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function InstrumentationLibraryLogs(; kwargs...)
        obj = new(meta(InstrumentationLibraryLogs), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct InstrumentationLibraryLogs
const __meta_InstrumentationLibraryLogs = Ref{ProtoMeta}()
function meta(::Type{InstrumentationLibraryLogs})
    ProtoBuf.metalock() do
        if !isassigned(__meta_InstrumentationLibraryLogs)
            __meta_InstrumentationLibraryLogs[] =
                target = ProtoMeta(InstrumentationLibraryLogs)
            allflds = Pair{Symbol,Union{Type,String}}[
                :instrumentation_library=>opentelemetry.proto.common.v1.InstrumentationLibrary,
                :log_records=>Base.Vector{LogRecord},
                :schema_url=>AbstractString,
            ]
            meta(
                target,
                InstrumentationLibraryLogs,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_InstrumentationLibraryLogs[]
    end
end
function Base.getproperty(obj::InstrumentationLibraryLogs, name::Symbol)
    if name === :instrumentation_library
        return (obj.__protobuf_jl_internal_values[name])::opentelemetry.proto.common.v1.InstrumentationLibrary
    elseif name === :log_records
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{LogRecord}
    elseif name === :schema_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct ResourceLogs <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ResourceLogs(; kwargs...)
        obj = new(meta(ResourceLogs), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ResourceLogs
const __meta_ResourceLogs = Ref{ProtoMeta}()
function meta(::Type{ResourceLogs})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ResourceLogs)
            __meta_ResourceLogs[] = target = ProtoMeta(ResourceLogs)
            allflds = Pair{Symbol,Union{Type,String}}[
                :resource=>opentelemetry.proto.resource.v1.Resource,
                :instrumentation_library_logs=>Base.Vector{InstrumentationLibraryLogs},
                :schema_url=>AbstractString,
            ]
            meta(
                target,
                ResourceLogs,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ResourceLogs[]
    end
end
function Base.getproperty(obj::ResourceLogs, name::Symbol)
    if name === :resource
        return (obj.__protobuf_jl_internal_values[name])::opentelemetry.proto.resource.v1.Resource
    elseif name === :instrumentation_library_logs
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            InstrumentationLibraryLogs,
        }
    elseif name === :schema_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct LogsData <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function LogsData(; kwargs...)
        obj = new(meta(LogsData), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct LogsData
const __meta_LogsData = Ref{ProtoMeta}()
function meta(::Type{LogsData})
    ProtoBuf.metalock() do
        if !isassigned(__meta_LogsData)
            __meta_LogsData[] = target = ProtoMeta(LogsData)
            allflds =
                Pair{Symbol,Union{Type,String}}[:resource_logs=>Base.Vector{ResourceLogs}]
            meta(
                target,
                LogsData,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_LogsData[]
    end
end
function Base.getproperty(obj::LogsData, name::Symbol)
    if name === :resource_logs
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ResourceLogs}
    else
        getfield(obj, name)
    end
end

export SeverityNumber,
    LogRecordFlags, LogsData, ResourceLogs, InstrumentationLibraryLogs, LogRecord


text:
#=
Simplest case : input only.
No recurrence or inhibition here
=#

push!(LOAD_PATH, abspath(@__DIR__,".."))

using LinearAlgebra,Statistics,StatsBase
using Plots,NamedColors ; theme(:dark)
using ProgressMeter
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using Distributions
using Random ; Random.seed!(0)
using FileIO
##

const Ne = 500
const Ni = 0
const Ntot = Ne+Ni

const Nas = 3
const p_as = 0.20

const Δt_as = 0.2
const Δt_blank = 0.1
const t_as_delay = 0.3
const Ttot = 4.0

# input to all neurons
const lowrate = 1E3
# additional input to assembly neurons
const highrate = 3E3
##

as_idxs = map(_->findall(rand(Ne) .< p_as),1:Nas)

inputfun=S.pattern_functor(Δt_as,Ttot,lowrate,highrate,Ntot,
  as_idxs;t_pattern_delay=t_as_delay,Δt_pattern_blank=Δt_blank)

## let's visualize it !
_  = let times = range(-1,Ttot+3;length=500)
  ys=Vector{Float64}[]
  for t in times
    y = inputfun.(t,1:Ne)
    push!(ys,y)
  end
  ys = hcat(ys...)
  heatmap(times,1:Ntot,ys)
end

idxs_sort = S.order_by_pattern_idxs(as_idxs,Ne)

_  = let times = range(-1,Ttot+3;length=500)
  ys=Vector{Float64}[]
  for t in times
    y = inputfun.(t,1:Ne)
    push!(ys,y)
  end
  ys = hcat(ys...)[idxs_sort,:]
  heatmap(times,1:Ntot,ys)
end

# we also need the upper limit
inputfun_upper = S.pattern_functor_upperlimit(lowrate,highrate,Ne,as_idxs)

_  = let times = range(-1,Ttot+10;length=500)
  ys=Vector{Float64}[]
  for t in times
    y = inputfun_upper.(t,1:Ne)
    push!(ys,y)
  end
  ys = hcat(ys...)[idxs_sort,:]
  heatmap(times,1:Ntot,ys;title="upper limit for all times")
end

##
# let's say kernel and v_reversal are same as E population
# excitatory neurons are Auguste's
const taueplus = 6E-3
const taueminus = 1E-3
const eifslope = 2.0
const τe = 20E-3
const Cap = 300.0
const vthexp = -50.0
const v_rev_e = 0.0
const vth_e = 20.0
const v_reset_e = -60.0
const v_rest_e = -60.0
const τrefr  = 5E-3
nt_e = let sker = S.SKExpDiff(taueplus,taueminus)
  sgen = S.SpikeGenEIF(vthexp,eifslope)
  S.NTLIFConductance(sker,sgen,τe,Cap,
   vth_e,v_reset_e,v_rest_e,τrefr,v_rev_e)
end
ps_e = S.PSLIFConductance(nt_e,Ne)

const w_in = 3.0
nt_in = let sker = S.SKExpDiff(taueplus,taueminus)
  sgen = S.SGPoissonFExact(inputfun,inputfun_upper)
  S.NTInputConductance(sgen,sker,v_rev_e)
end
ps_in = S.PSInputPoissonConductanceExact(nt_in,w_in,Ntot)

##

pop_e = S.Population(ps_e,(S.FakeConnection(),ps_in))

## network
const dt = 0.1E-3
ntw = S.RecurrentNetwork(dt,pop_e)

# record spiketimes and internal potential
krec = 1
n_e_rec = Ne
t_wup = 0.0
rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot;idx_save=collect(1:n_e_rec),t_warmup=t_wup)
rec_spikes_e = S.RecSpikes(ps_e,500.0,Ttot;idx_save=collect(1:Ne),t_warmup=t_wup)

## Run

times = (0:ntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_state_e,rec_spikes_e])
S.reset!.([ps_e,ps_in])
# initial conditions
ps_e.state_now .= v_reset_e #rand(Uniform(v_reset_e,v_reset_e+0.3),Ne)

@time begin
  @showprogress 5.0 "network simulation " for (k,t) in enumerate(times)
    rec_state_e(t,k,ntw)
    rec_spikes_e(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end
end

S.add_fake_spikes!(1.0vth_e,rec_state_e,rec_spikes_e)

## visualize activity

_ = let neu = 1,
  plt=plot()
  plot!(plt,rec_state_e.times,rec_state_e.state_now[neu,:];
    linewidth=2,leg=false,
    xlims=(0,Ttot))
end

##
idxs_sort = S.order_by_pattern_idxs(as_idxs,Ne)
_ = let Traster = Ttot
  tmp = S.raster_png(1E-3,rec_spikes_e;spike_height=3,Ttot=Traster,Nneurons=Ne,
   reorder=idxs_sort)
  save("/tmp/tmp.png",tmp)
end

##
_ = let (_times,binned) =S.binned_spikecount(0.01,rec_spikes_e)
  heatmap(_times,1:size(binned,1),binned[idxs_sort,:])
end

ext: jl


text:
using Documenter
using VLOptionsModeling

makedocs(
    sitename = "Model",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [VLOptionsModeling],
    pages = [
        "VLOptionsModeling" => "index.md",
    ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/varnerlab/VLOptionsModeling.jl.git",
    devurl = "stable",
    devbranch = "main",
)


text:
function ConvexHull2D(arg0::Vector{Vector2D}, arg1::jdouble)
    return ConvexHull2D((Vector{Vector2D}, jdouble), arg0, arg1)
end

function create_region(obj::ConvexHull2D)
    return jcall(obj, "createRegion", Region, ())
end

function get_line_segments(obj::ConvexHull2D)
    return jcall(obj, "getLineSegments", Vector{Segment}, ())
end

function get_vertices(obj::ConvexHull2D)
    return jcall(obj, "getVertices", Vector{Vector2D}, ())
end


text:
"""
    IndexablePoint2D <: AbstractPoint2D

A special 2D point type, compatible with `VoronoiDelaunay`, that carries an
identifying `_index` field.

By default, `_index` is -1.

See also: `setindex!`, `getindex`
"""
mutable struct IndexablePoint2D <: AbstractPoint2D
    _x::Float64
    _y::Float64
    _index::Int
end
IndexablePoint2D(x::Float64, y::Float64) = IndexablePoint2D(x, y, -1)
getx(p::IndexablePoint2D) = p._x
gety(p::IndexablePoint2D) = p._y
Base.:+(p1::IndexablePoint2D, p2::IndexablePoint2D) = IndexablePoint2D(getx(p1)+getx(p2), gety(p1)+gety(p2), -1)
Base.:/(p1::IndexablePoint2D, n::Real) = IndexablePoint2D(getx(p1)/n, gety(p1)/n, -1)

"""
    getindex(p::IndexablePoint2D)

Return `_index` field of `p`.
"""
Base.getindex(p::IndexablePoint2D) = p._index

"""
    setindex!(p::IndexablePoint2D)

Set the `_index` field of `p`.

# Example

```jldoctest
julia> p = IndexablePoint2D(103.4, 25.0);

julia> getindex(p)
-1

julia> setindex!(p, 90);

julia> getindex(p)
90
```
"""
Base.setindex!(p::IndexablePoint2D, v::Int) = setfield!(p, :_index, v)

"""
    same(e1::DelaunayEdge{IndexablePoint2D}, e2::DelaunayEdge{IndexablePoint2D})

Return `true` if edge `e1` is the same as `e2`, regardless of direction.
Otherwise, return `false`. This is similar to `==`.

To check if two edges are identical in the sense that no program could
distinguish them, use `===`.
"""
@inline function same(e1::DelaunayEdge{IndexablePoint2D},e2::DelaunayEdge{IndexablePoint2D})
    e1 === e2 && return true

    # if they're reversed
    e1a, e1b = geta(e1)._index, getb(e1)._index
    e2a, e2b = geta(e2)._index, getb(e2)._index
    (e1a == e2b) && (e1b == e2a) && return true

    return false
end

"""
    sameunique!

In-place `unique!` that considers both directions of edges the same (just like
[`same`](@ref)).
"""
function sameunique!(v::AbstractVector{DelaunayEdge{IndexablePoint2D}})
    deleteidxs = trues(length(v))
    for i in eachindex(v)
        @inbounds edgei = v[i]
        eia = getindex(geta(edgei))
        eib = getindex(getb(edgei))
        revdupe = false
        for j in eachindex(v)
            @inbounds edgej = v[j]
            eja = getindex(geta(edgej))
            ejb = getindex(getb(edgej))
            if (eia == ejb) && (eib == eja)
                revdupe = true
                break
            end
        end
        if !revdupe
            deleteidxs[i] = false
        end
    end

    deleteat!(v, deleteidxs)
end

# Improved performance of `delaunayedges()`
# see: https://github.com/JuliaGeometry/VoronoiDelaunay.jl/issues/47
function delaunayedges_fast(t::DelaunayTessellation2D{T}) where T <: AbstractPoint2D
    result = DelaunayEdge{T}[]
    @inbounds for ix in 2:t._last_trig_index
        tr = t._trigs[ix]
        isexternal(tr) && continue

        ix_na = tr._neighbour_a
        if (ix_na > ix) | isexternal(t._trigs[ix_na])
            push!(result, DelaunayEdge(getb(tr), getc(tr)))
        end
        ix_nb = tr._neighbour_b
        if (ix_nb > ix) | isexternal(t._trigs[ix_nb])
            push!(result, DelaunayEdge(geta(tr), getc(tr)))
        end
        ix_nc = tr._neighbour_c
        if (ix_nc > ix) | isexternal(t._trigs[ix_nc])
            push!(result, DelaunayEdge(geta(tr), getb(tr)))
        end
    end
    return result
end


text:
function execute(problem::PotentialFlowProblem{:analytical};kwargs...)

  # Parameters
  L = 2*π
  H = 1.0
  n = 20
  order = 2
  g = 9.81
  ξ = 0.1
  λ = L/2
  k = 2*π/L
  h = L/n
  ω = √(g*k*tanh(k*H))
  t₀ = 0.0
  tf = 8*π
  Δt = h/(10*λ*ω)
  θ = 0.5

  # Exact solution
  ϕₑ(x,t) = ω/k * ξ * (cosh(k*(x[2]))) / sinh(k*H) * sin(k*x[1] - ω*t)
  ηₑ(x,t) = ξ * cos(k*x[1] - ω*t)
  ϕₑ(t::Real) = x -> ϕₑ(x,t)
  ηₑ(t::Real) = x -> ηₑ(x,t)

  # Domain
  domain = (0,L,0,H)
  partition = (n,n)
  model = CartesianDiscreteModel(domain,partition;isperiodic=(true,false))

  # Boundaries
  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"bottom",[1,2,5])
  add_tag_from_tags!(labels,"free_surface",[3,4,6])
  bgface_to_mask = get_face_mask(labels,[3,4,6],1)
  model_Γ =BoundaryDiscreteModel(Polytope{1},model,bgface_to_mask)

  # Triangulation
  Ω = Triangulation(model)
  Γ = Triangulation(model_Γ)
  dΩ = Measure(Ω,2*order)
  dΓ = Measure(Γ,2*order)

  # FE spaces
  reffe = ReferenceFE(lagrangian,Float64,order)
  V = TestFESpace(model,reffe,conformity=:H1)
  V_Γ = TestFESpace(model_Γ,reffe,conformity=:H1)
  U = TransientTrialFESpace(V)
  U_Γ = TransientTrialFESpace(V_Γ)
  X = TransientMultiFieldFESpace([U,U_Γ])
  Y = MultiFieldFESpace([V,V_Γ])

  # Weak form
  α = 2/Δt
  m((ϕt,ηt),(w,v)) = ∫( 0.5*(α/g*(w*ϕt) + v*ϕt) - (w*ηt) )dΓ
  a((ϕ,η),(w,v)) = ∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( 0.5*(α*(w*η) + g*v*η) )dΓ
  b((w,v)) = ∫( 0.0*w )dΓ
  op = TransientConstantFEOperator(m,a,b,X,Y)

  # Solver
  ls = LUSolver()
  odes = ThetaMethod(ls,Δt,θ)
  solver = TransientFESolver(odes)

  # Initial solution
  x₀ = interpolate_everywhere([ϕₑ(0.0),ηₑ(0.0)],X(0.0))

  # Solution
  sol_t = solve(solver,op,x₀,t₀,tf)

  # Post-process
  l2_Ω(w) = √(∑(∫(w*w)dΩ))
  l2_Γ(v) = √(∑(∫(v*v)dΓ))
  E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ )
  E_pot(v) = g*0.5*∑( ∫(v*v)dΓ )
  Eₑ = 0.5*g*ξ^2*L

  folderName = "ϕFlow-results"
  if !isdir(folderName)
    mkdir(folderName)
  end
  filePath_Ω = folderName * "/fields_Ω"
  filePath_Γ = folderName * "/fields_Γ"
  pvd_Ω = paraview_collection(filePath_Ω, append=false)
  pvd_Γ = paraview_collection(filePath_Γ, append=false)
  for ((ϕn,ηn),tn) in sol_t
    E = E_kin(ϕn) + E_pot(ηn)
    error_ϕ = l2_Ω(ϕn-ϕₑ(tn))
    error_η = l2_Γ(ηn-ηₑ(tn))
    #println(E/Eₑ," ", error_ϕ," ",error_η)

    pvd_Ω[tn] = createvtk(Ω,filePath_Ω * "_$tn.vtu",cellfields = ["phi" => ϕn])
    pvd_Γ[tn] = createvtk(Γ,filePath_Γ * "_$tn.vtu",cellfields = ["eta" => ηn])
  end
  vtk_save(pvd_Ω)
  vtk_save(pvd_Γ)
end

include("PotentialFlowBeam.jl")


text:
function test_matrix_rows()
    @test rows([1 2; 3 4]) == Array[[1,2],[3,4]]
end

function test_matrix_cols()
    @test cols([1 2; 3 4]) == Array[[1,3],[2,4]]
end

function test_matrix_diags()
    @test diags([1 2; 3 4]) == Array[[3],[1,4],[2]]
end

function test_matrix_antidiags()
    @test antidiags([1 2; 3 4]) == Array[[1],[2,3],[4]]
end

function test_matrix_all()
    print(rpad("DataStructure.Matrix...", 50, ' '))

    test_matrix_rows()
    test_matrix_cols()
    test_matrix_diags()
    test_matrix_antidiags()

    println("PASS")
end


text:
#this will construct types from mod files related to channels or whatever

#=

Everything in the mod files that is defined as a PARAMETER, STATE, or DERIVATIVE ends up being a part of param
Additionally, the last two entries in param seem to always be v, and _g (conductance)?
Things in param can be shared with other Props (like voltage is equal at every Prop for each node)

ppvar contains other values, looks like usually related to ions

So as far as Julia goes, at each node there are state variables unique to one ion channel(e.g. sodium gate m), state variables shared among ion channels (e.g. voltage), constant paramters for each channel
=#

type ParsedHoc
    file::Array{ByteString,1}

end

function ParsedHoc(filename::ASCIIString)
    f = open(filename)
    a=readlines(f)
    close(f)
    ParsedHoc()
end

function parsehoc(ph::ParsedHoc)
    linenum=1
    cursection=""
    sections=["COMMENT","UNITS","NEURON","PARAMETER","STATE","ASSIGNED","BREAKPOINT","INITIAL",
              "DERIVATIVE","PROCEDURE","FUNCTION"]

    while linenum<length(ph.file)

        if cursection==""
            for i in sections
                if contains(ph.file[linenum],i)
                    cursection=i
                    break
                end
            end

        else
            if contains(ph.file[linenum],"}")
                cursection=""
            elseif contains(ph.file[linenum],"ENDCOMMENT")
                cursection=""
                #skip blank lines or lines with comments
            elseif cursection=="COMMENT"
            elseif cursection=="UNITS"
                #don't care for now
            elseif cursection=="NEURON"
                #word after READ should be collected to add to myvars dictionary
            elseif cursection=="PARAMETER"
                first=split(ph.file[linenum],"=")
                var=replace(first[1]," ", "")
                num=split(first[2], "(")
                num=replace(num[1], " ", "")
                parsehoc.param[var]=float(num)
            elseif cursection=="STATE"
                #get characters not equal to whitespace
            elseif cursection=="ASSIGNED"
                #get characters not in parenthesis
            elseif cursection=="BREAKPOINT"
                #second half of prop_calc method
            elseif cursection=="INITIAL"
                #prop_init method
            elseif cursection=="DERIVATIVE"
                #first half of prop_calc method
            elseif cursection=="PROCEDURE"
                #for this and function, collect whole thing, figure out what variables from where, and get rid of comments and delcarations
            elseif cursection=="FUNCTION"
            end

        end

    end

end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
# Population density benchmark from https://github.com/SWI-Prolog/bench
# Author:  David H. D. Warren
#
# query population and area database to find coun-
# tries of approximately equal population density

clauses = @julog [
    query([C1,D1,C2,D2]) <<=
        density(C1,D1) &
        density(C2,D2) &
        (D1 > D2) &
        is(T1, 20*D1) &
        is(T2, 21*D2) &
        (T1 < T2),

    density(C,D) <<=
        pop(C,P) &
        area(C,A) &
        is(D, (P*100)/A),

    # populations in 100000's
    pop(china,	8250) <<= true,
    pop(india,	5863) <<= true,
    pop(ussr,	2521) <<= true,
    pop(usa,	2119) <<= true,
    pop(indonesia,	1276) <<= true,
    pop(japan,	1097) <<= true,
    pop(brazil,	1042) <<= true,
    pop(bangladesh,	 750) <<= true,
    pop(pakistan,	 682) <<= true,
    pop(w_germany,	 620) <<= true,
    pop(nigeria,	 613) <<= true,
    pop(mexico,	 581) <<= true,
    pop(uk,		 559) <<= true,
    pop(italy,	 554) <<= true,
    pop(france,	 525) <<= true,
    pop(philippines, 415) <<= true,
    pop(thailand,	 410) <<= true,
    pop(turkey,	 383) <<= true,
    pop(egypt,	 364) <<= true,
    pop(spain,	 352) <<= true,
    pop(poland,	 337) <<= true,
    pop(s_korea,	 335) <<= true,
    pop(iran,	 320) <<= true,
    pop(ethiopia,	 272) <<= true,
    pop(argentina,	 251) <<= true,

    # areas in 1000's of square miles
    area(china,     3380) <<= true,
    area(india,     1139) <<= true,
    area(ussr,      8708) <<= true,
    area(usa,       3609) <<= true,
    area(indonesia,  570) <<= true,
    area(japan,      148) <<= true,
    area(brazil,    3288) <<= true,
    area(bangladesh,  55) <<= true,
    area(pakistan,   311) <<= true,
    area(w_germany,   96) <<= true,
    area(nigeria,    373) <<= true,
    area(mexico,     764) <<= true,
    area(uk,          86) <<= true,
    area(italy,      116) <<= true,
    area(france,     213) <<= true,
    area(philippines, 90) <<= true,
    area(thailand,   200) <<= true,
    area(turkey,     296) <<= true,
    area(egypt,      386) <<= true,
    area(spain,      190) <<= true,
    area(poland,     121) <<= true,
    area(s_korea,     37) <<= true,
    area(iran,       628) <<= true,
    area(ethiopia,   350) <<= true,
    area(argentina, 1080) <<= true,
]

goal = @julog(query([C1,D1,C2,D2]))
println("Population (single run):")
@time resolve(goal, clauses)

function benchmark(n)
    for i = 1:n
        resolve(@julog(query([C1,D1,C2,D2])), clauses)
    end
end

println("Population (100 runs):")
@time benchmark(100)


text:
@testset "helpers" begin
    # unvcat: 1D
    arrs = Any[i * ones(i) for i=1:4]
    A = 2 * vcat(arrs...)
    for i=1:4
        @test 2 * arrs[i] == unvcat(A, i, arrs...)
    end

    # unvcat: 2D
    arrs = Any[i * ones(i, 3) for i=1:4]
    A = 2 * vcat(arrs...)
    for i=1:4
        @test 2 * arrs[i] == unvcat(A, i, arrs...)
    end

    # unvcat: 4D
    arrs = Any[i * ones(i, 3, 4, 5) for i=1:4]
    A = 2 * vcat(arrs...)
    for i=1:4
        @test 2 * arrs[i] == unvcat(A, i, arrs...)
    end

    # unvcat: 2D
    arrs = Any[i * ones(3, i) for i=1:4]
    A = 2 * hcat(arrs...)
    for i=1:4
        @test 2 * arrs[i] == unhcat(A, i, arrs...)
    end

    # unvcat: 4D
    arrs = Any[i * ones(3, i, 5, 6) for i=1:4]
    A = 2 * hcat(arrs...)
    for i=1:4
        @test 2 * arrs[i] == unhcat(A, i, arrs...)
    end
end


text:
s = BenchmarkGroup()
suite["matching"] = s

Random.seed!(17)
n = 100
p = 0.2
w = (w=sprand(n,n,p); w=w+w'; w-Diagonal(w))
g = Graph(w)
s["match n=$n p=$p"] = @benchmarkable minimum_weight_perfect_matching($g, $(EdgeMap(g, w)))

n = 300
p = 0.1
w = (w=sprand(n,n,p); w=w+w'; w-Diagonal(w))
g = Graph(w)
s["match n=$n p=$p"] = @benchmarkable minimum_weight_perfect_matching($g, $(EdgeMap(g, w)))


text:
using ..SearchProblem
using ..SearchProblem: reset!
using Random

const SP_ACTIONS_MEANINGS = Symbol[:UP, :DOWN, :LEFT, :RIGHT]  # In this order to generate children
const SP_ACTIONS_DIRECTIONS = Tuple{Int, Int}[(-1, 0), (1, 0), (0, -1), (0, 1)]

mutable struct SlidingPuzzle <: SearchProblem.AbstractSearchProblem{Matrix{UInt8}, Symbol}
    side_range::AbstractArray{Int}
    manhat_range::AbstractArray{Int}
    inverse::Bool
    side::Int
    puzzle::Array{UInt8,2}
    blank_loc::Tuple{Int,Int}
    puzzle_manhat::Int
    goal_puzzle::Array{UInt8,2}
    heuritic_cache::Dict{Array{UInt8,2},Float64}
    lk::ReentrantLock
    rng::MersenneTwister
    function SlidingPuzzle(side_range::AbstractArray{Int}, manhat_range::AbstractArray{Int}; inverse::Bool=false)
        sp = new(side_range, manhat_range, inverse)
        sp.rng = MersenneTwister()
        sp.heuritic_cache = Dict{Array{UInt8,2},Float64}()
        sp.lk = ReentrantLock()
        reset!(sp)
        return sp
    end
end

function max_possible_manhat(n::Int)::Int
    s = 0
    k = 1
    while n - k >= 0
        s += n - k
        k += 2
    end
    return 4n * s
end

function initial_puzzle!(puzzle::Array{UInt8,2})
    side = size(puzzle)[1]
    for column in 1:side
        @simd for row in 1:side
            @inbounds puzzle[row, column] = (row - 1) * side + column
        end
    end
    puzzle[side, side] = 0
    return puzzle
end

function initial_puzzle(side::Int)::Matrix{UInt8}
    puzzle = Array{UInt8}(undef, side, side)
    return initial_puzzle!(puzzle)
end

@inline is_valid_blank_location(side::Int, blank_location::Tuple{Int,Int})::Bool = all((1, 1) .<= blank_location .<= (side, side))

function puzzle_legal_actions(side::Int, blank_location::Tuple{Int,Int})
    return Iterators.filter(1:4) do action_id
        @inbounds direction::Tuple{Int, Int} = SP_ACTIONS_DIRECTIONS[action_id]
        return is_valid_blank_location(side, blank_location .+ direction)
    end
end

@inline l1_norm(p1, p2) = sum(abs.(p1 .- p2))

function find_tile_location(puzzle::Array{UInt8, 2}, tile::UInt8)::Tuple{Int, Int}
    side = size(puzzle)[1]
    for c in 1:side
        for r in 1:side
            @inbounds puzzle[r, c] == tile  &&  return r, c
        end
    end
    return -1, -1
end

function compute_manhat(puzzle::Array{UInt8,2}, to_puzzle::Array{UInt8,2}, inverse=false)::Union{Int, Float64}
    side = size(puzzle)[1]
    manhat = 0.0
    for to_c::Int in 1:side
        for to_r::Int in 1:side
            @inbounds value = to_puzzle[to_r, to_c]
            value == 0 && continue
            r, c = find_tile_location(puzzle, value)
            manhat += inverse ? l1_norm((r, c), (to_r, to_c)) / Float64(value) : l1_norm((r, c), (to_r, to_c))
        end
    end
    return manhat
end

function compute_next_puzzle(puzzle::Array{UInt8,2}, blank_loc::Tuple{Int,Int}, action::Int)::Tuple{Matrix{UInt8}, Tuple{Int, Int}}
    new_puzzle::Array{UInt8, 2} = copy(puzzle)
    blank_r, blank_c = blank_loc
    @inbounds new_blank_r, new_blank_c = blank_loc .+ SP_ACTIONS_DIRECTIONS[action]
    @inbounds new_puzzle[blank_r, blank_c] = new_puzzle[new_blank_r, new_blank_c]
    @inbounds new_puzzle[new_blank_r, new_blank_c] = 0
    new_blank_loc = (new_blank_r, new_blank_c)
    return new_puzzle, new_blank_loc
end

function shuffle_puzzle_hillclimb!(sp::SlidingPuzzle, target_manhat::Int, ϵ=0.1)::Int
    cache = sp.inverse ? Dict{Array{UInt8,2},Float64}() : sp.heuritic_cache
    manhat = haskey(cache, sp.puzzle) ? cache[sp.puzzle] : compute_manhat(sp.puzzle, sp.goal_puzzle)
    cache[sp.puzzle] = manhat
    while manhat < target_manhat
        legal_actions = puzzle_legal_actions(sp.side, sp.blank_loc)
        action = rand(sp.rng, collect(legal_actions))
        new_puzzle, new_blank_loc = compute_next_puzzle(sp.puzzle, sp.blank_loc, action)
        new_manhat = haskey(cache, new_puzzle) ? cache[new_puzzle] : compute_manhat(new_puzzle, sp.goal_puzzle)
        cache[new_puzzle] = new_manhat
        if new_manhat > manhat || rand(sp.rng) < ϵ
            sp.puzzle = new_puzzle
            sp.blank_loc = new_blank_loc
            manhat = new_manhat
        end
    end
    return manhat
end

function SearchProblem.info(sp::SlidingPuzzle)
    return Dict(
        :puzzle_manhat => sp.puzzle_manhat,
        :inverse_cost => sp.inverse
    )
end

@inline function SearchProblem.obs(sp::SlidingPuzzle)::Vector{Float64}
    return Float64[sp.puzzle_manhat / maximum(sp.manhat_range)]
end

function SearchProblem.reset!(sp::SlidingPuzzle)
    empty!(sp.heuritic_cache)
    sp.side = rand(sp.rng, sp.side_range)
    sp.goal_puzzle = initial_puzzle(sp.side)
    sp.puzzle = copy(sp.goal_puzzle)
    sp.blank_loc = (sp.side, sp.side)
    target_manhat = rand(sp.rng, sp.manhat_range)
    sp.puzzle_manhat = shuffle_puzzle_hillclimb!(sp, target_manhat)
    return nothing
end

@inline SearchProblem.start_state(sp::SlidingPuzzle)::Matrix{UInt8} = sp.puzzle

function findmovedtile(from_puzzle::Array{UInt8,2}, to_puzzle::Array{UInt8,2})::UInt8
    side = size(from_puzzle)[1]
    for c in 1:side
        @inbounds for r in 1:side
            if from_puzzle[r, c] != to_puzzle[r, c]
                return max(from_puzzle[r, c], to_puzzle[r, c])  # one of them is zero, the other is the moved tile
            end
        end
    end
end

@inline function SearchProblem.cost(sp::SlidingPuzzle, state::Array{UInt8,2}, action::Symbol, next_state::Array{UInt8,2})::Float64
    if sp.inverse
        moved_tile = findmovedtile(state, next_state)
        return 1.0 / moved_tile
    else
        return 1.0
    end
end

@inline SearchProblem.goal_test(sp::SlidingPuzzle, state::Array{UInt8,2})::Bool = state == sp.goal_puzzle

@inline function SearchProblem.heuristic(sp::SlidingPuzzle, state::Array{UInt8,2})::Float64
    lock(sp.lk) do
        if haskey(sp.heuritic_cache, state)
            return sp.heuritic_cache[state]
        end
    end

    h = compute_manhat(state, sp.goal_puzzle, sp.inverse)
    lock(sp.lk) do
        sp.heuritic_cache[state] = h
    end
    return h
end

function SearchProblem.successors(sp::SlidingPuzzle, state::Array{UInt8,2})
    blank_loc = find_tile_location(state, 0x00)
    return Iterators.map(puzzle_legal_actions(sp.side, blank_loc)) do action_id
        @inbounds action::Symbol = SP_ACTIONS_MEANINGS[action_id]
        next_puzzle, _ = compute_next_puzzle(state, blank_loc, action_id)
        return next_puzzle, action
    end
end


text:
ker(x::Real,y::Real; sigma::Float64=2.0)::Float64 = exp( -sigma*((x-y)^2) )

# Assuming equally weighted particles
function mmd!(val::Vector{Float64}, a::Vector{Float64}, b::Vector{Float64})
  len = length(a)
  @assert len == length(b)
  reci_len = 1.0/len
  val = 0.0
  @inbounds @fastmath for i in 1:len
     @simd for j in 1:len
       val -= ker(a[i], b[j])
     end
  end
  @fastmath val *= 2.0
  @inbounds @fastmath for i in 1:len
    @simd for j in 1:len
      val += ker(a[i], a[j])
      val += ker(b[i], b[j])
    end
  end
  @fastmath val *= reci_len
  @fastmath val *= reci_len
  return val
end


text:
function pascal(n)
   (n<=0) && error("Pascal trinalge can not have zero or negative rows")
   r=Vector{Int}(undef,n)
   pr=Vector{Int}(undef,n)
   pr[1]=r[1]=1
   println(@view pr[1])
   for i=2:n
      r[1]=r[i]=1
      for j=2:i-1
         r[j]=pr[j-1]+pr[j]
      end
      println(join(view(r,1:i), " "))
      r,pr=pr,r
   end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module PlotUtils

using ColorSchemes
using Dates
using Reexport
using Printf
@reexport using Colors
import Base: getindex
import Random: MersenneTwister

export ColorGradient,
    ColorPalette,
    cgrad,
    palette,
    color_list,
    cvec,
    rgb_string,
    rgba_string,
    invisible,
    get_color_palette,
    isdark,
    plot_color,
    adapted_grid,
    default_cgrad,
    zscale

include("color_utils.jl")
include("colors.jl")
include("colorschemes.jl")
include("adapted_grid.jl")
include("intervals.jl")

export optimize_ticks, optimize_datetime_ticks

include("ticks.jl")

const _default_colorscheme = generate_colorscheme()

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end # module


text:
simpleloop = quote
    while true
        x = 1/2*x
        y = 2y
    end
end

petter(n) = quote
    x, y = 1, 0
    while true
        x = x + y^$n
        y = y + 1
    end
end

cohencu = quote
    while n<=N
        n = n+1
        x = x+y
        y = y+z
        z = z+6
    end
end

freire1 = quote
    while x>r
        x = x-r
        r = r+1
    end
end

freire2 = quote
    while x-s > 0
        x = x-s
        s = s+6*r+3
        r = r+1
    end
end

psolv1s = quote
    k = 0
    t1, t2, a, b, c, d = 1, 1, 1, 1, 1, 1
    while true
        t1 = t2
        t2 = a
        a = 5*(k+2)*t2+6*(k^2+3k + 2)*t1
        b = 2*b
        c = 3*(k+2)*c
        d = (k+2)*d
        k = k+1
    end
end

psolv2s = quote
    t1 = 1; t2 = 1;
    s1 = 1; s2 = 2;
    a = 3; b = 1;
    c = 1; d = 3;
    e = 2; f = 5
    n = 1
    while true
        a = 3*(n+3/2)*a
        s1 = s2; s2 = b
        b = 5*(n+3/2)*s2- 3/2*(1+2*n)*(3+2*n)*s1
        t1 = t2; t2 = d
        d = 4*(4 + n)*t2 - 3*(3 + n)*(4 + n)*t1
        e = (n + 4)*e
        f = 2*f
        c = -3*c + 2
        n = n+1
    end
end

intcbrt = quote
    while true
        x = x-s
        s = s + 6r + 3
        r = r + 1
    end
end

cubes = quote
    c, k, m, n = 0, 1, 6, 0
    while true
        c = c + k
        k = k + m
        m = m + 6
        n = n + 1
    end
end

eucliddiv = quote
    while true
        r = r - y
        q = q + 1
    end
end


text:
using ZulipReminderBot

cfg = DotEnv.config()
include("configuration.jl")
setupbot!(token = cfg["REMINDER_BOT_TOKEN"],
          host = cfg["REMINDER_BOT_HOST"],
          port = cfg["REMINDER_BOT_PORT"],
          email = cfg["REMINDER_BOT_EMAIL"],
          apikey = cfg["REMINDER_BOT_API_KEY"],
          ep = cfg["REMINDER_BOT_ZULIPCHAT"])

# const db = SQLite.DB(cfg["REMINDER_BOT_DB"])

ZulipReminderBot.run(conn)


text:
using QuasiGeostrophy
using Test
import QuasiGeostrophy:compute

include(pwd() * "/test/test_utils.jl")

@wrapper a=1 b=2

@testset "Impero operator matching" begin

    println("for $a = ", a.data, " and $b = ", b.data)
    for unary_operator in unary_operators
        b_name, b_symbol = Meta.parse.(unary_operator)
        @eval c = $b_symbol(a)
        println("The value of $c is ", compute(c))
        @test @eval ( compute(c) == $b_symbol(a.data) )
    end

    for binary_operator in binary_operators
        b_name, b_symbol = Meta.parse.(binary_operator)
        @eval c = $b_symbol(a, b)
        println("The value of $c is ", compute(c))
        @test @eval compute(c) == $b_symbol(a.data, b.data)
    end

end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
struct SimpleTimeDependentSystem
    de::Function
    function SimpleTimeDependentSystem(de::Function)
        return new(de)
    end
end

function stiff_de!(xx::Vector{T}, x::Vector{T}, s::SimpleTimeDependentSystem, t::Float64=0.0) where {T}
    xx .= t
    return nothing
end

my_stiff_struct = SimpleTimeDependentSystem(stiff_de!)

x0 = zeros(Float64, 1) .+ 0.0
rr = makeRadauIntegrator(my_stiff_struct, x0, 1.0e-16, 2)
rr.rule.s = 2
t_final = 1.0
update_h!(rr, t_final)
h, x_final, t_f = solveRadau(rr, x0, 0.0)

@testset "time dependent" begin
    @test x_final[1] ≈ 0.5
end


text:
using BenchmarkTools
using ConstraintSolver, JuMP, MathOptInterface
using GLPK, JSON, Cbc

const CS = ConstraintSolver
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const SUITE = BenchmarkGroup()

dir = pkgdir(ConstraintSolver)
include(joinpath(dir, "benchmark/sudoku/benchmark.jl"))

SUITE["sudoku"] = BenchmarkGroup(["alldifferent"])
sudoku_grids = from_file(joinpath(dir, "benchmark/sudoku/data/top95.txt"))
for i in 1:5:95
    SUITE["sudoku"]["top95_$i"] = @benchmarkable solve_sudoku($sudoku_grids[$i]) seconds = 2
end

include(joinpath(dir, "benchmark/eternity/benchmark.jl"))

SUITE["eternity"] = BenchmarkGroup(["alldifferent", "table", "equal"])
# compiling run
solve_eternity("eternity_6x5"; height = 6, width = 5)
SUITE["eternity"]["6x5"] =
    @benchmarkable solve_eternity("eternity_6x5"; height = 6, width = 5) seconds = 30
SUITE["eternity"]["6x5_ABS"] =
    @benchmarkable solve_eternity("eternity_6x5"; height = 6, width = 5, branch_strategy=:ABS) seconds = 30
SUITE["eternity"]["5x5_opt"] =
    @benchmarkable solve_eternity("eternity_5x5"; height = 5, width = 5, optimize = true) seconds =
        120
SUITE["eternity"]["5x5_opt_ind"] = @benchmarkable solve_eternity(
    "eternity_5x5";
    height = 5,
    width = 5,
    optimize = true,
    indicator = true,
) seconds = 120
SUITE["eternity"]["5x5_opt_rei"] = @benchmarkable solve_eternity(
    "eternity_5x5";
    height = 5,
    width = 5,
    optimize = true,
    reified = true,
) seconds = 120
SUITE["eternity"]["5x5_all"] =
    @benchmarkable solve_eternity("eternity_5x5"; all_solutions = true) seconds = 30

include(joinpath(dir, "benchmark/lp/benchmark.jl"))

SUITE["lp"] = BenchmarkGroup(["objective", "less_than"])
# compiling run
solve_lp()
SUITE["lp"]["issue_83"] = @benchmarkable solve_lp() seconds = 2

include(joinpath(dir, "benchmark/killer_sudoku/benchmark.jl"))

SUITE["killer_sudoku"] = BenchmarkGroup(["alldifferent", "equal"])
# compiling run
solve_killer_sudoku("niallsudoku_5500")
SUITE["killer_sudoku"]["niall_5500"] =
    @benchmarkable solve_killer_sudoku("niallsudoku_5500") seconds = 5
SUITE["killer_sudoku"]["niall_5501"] =
    @benchmarkable solve_killer_sudoku("niallsudoku_5501") seconds = 5
SUITE["killer_sudoku"]["niall_5500_ABS"] =
    @benchmarkable solve_killer_sudoku("niallsudoku_5500"; branch_strategy=:ABS) seconds = 5
SUITE["killer_sudoku"]["niall_5501_ABS"] =
    @benchmarkable solve_killer_sudoku("niallsudoku_5501"; branch_strategy=:ABS) seconds = 5
SUITE["killer_sudoku"]["niall_5500_special"] =
    @benchmarkable solve_killer_sudoku("niallsudoku_5500"; special = true) seconds = 15
SUITE["killer_sudoku"]["niall_5501_special"] =
    @benchmarkable solve_killer_sudoku("niallsudoku_5501"; special = true) seconds = 15

include(joinpath(dir, "benchmark/graph_color/benchmark.jl"))

SUITE["graph_coloring"] = BenchmarkGroup(["notequal", "equal", "svc"])
# compiling run
solve_us_graph_coloring()
SUITE["graph_coloring"]["US_8+equal"] =
    @benchmarkable solve_us_graph_coloring(; num_colors = 8, equality = true) seconds = 5
SUITE["graph_coloring"]["US_50colors+equal"] =
    @benchmarkable solve_us_graph_coloring(; num_colors = 50, equality = true) seconds = 5
SUITE["graph_coloring"]["US_8+equal_ABS"] =
    @benchmarkable solve_us_graph_coloring(; num_colors = 8, equality = true, branch_strategy=:ABS) seconds = 5
SUITE["graph_coloring"]["US_50colors+equal_ABS"] =
    @benchmarkable solve_us_graph_coloring(; num_colors = 50, equality = true, branch_strategy=:ABS) seconds = 5
SUITE["graph_coloring"]["US_50colors"] =
    @benchmarkable solve_us_graph_coloring(; num_colors = 50, equality = false) seconds = 5
SUITE["graph_coloring"]["queen7_7"] =
    @benchmarkable color_graph(joinpath(dir, "benchmark/graph_color/data/queen7_7.col"), 7) seconds =
        10
SUITE["graph_coloring"]["le450_5d"] =
    @benchmarkable color_graph(joinpath(dir, "benchmark/graph_color/data/le450_5d.col"), 5) seconds =
        30

include(joinpath(dir, "benchmark/scheduling/benchmark.jl"))

SUITE["scheduling"] = BenchmarkGroup(["cumulative", "equal", "less_than"])
# compiling run
furniture_moving()
SUITE["scheduling"]["furniture_moving"] =
    @benchmarkable furniture_moving() seconds = 5

# Problem instance
organize_day_problem = Dict(

    #task id     1      2      3      4
    :tasks => ["Work","Mail","Shop","Bank"],

    # duration of the four tasks
    :durations => [4,1,2,1],

    # precedences
    # [A,B] : task A must be completed before task B
    :precedences => [
                        4 3;
                        2 1
                    ],
    # Time limits
    :start_time => 9,
    :end_time => 17
)

SUITE["scheduling"]["organize_day"] =
    @benchmarkable organize_day(organize_day_problem) seconds = 5

include(joinpath(dir, "benchmark/small/benchmark.jl"))

SUITE["small"] = BenchmarkGroup(["all_different_except", "reified", "indicator", "less_than"])
# compiling run
all_different_except_0(5)
SUITE["scheduling"]["all_different_except_0_len10"] =
    @benchmarkable all_different_except_0(10) seconds = 10
SUITE["scheduling"]["all_different_except_0_len12"] =
    @benchmarkable all_different_except_0(12) seconds = 15


text:
using SymbolicUtils

"""
A function that creates a tuple of CartesianIndices of unit length and `N` dimensions, one pointing along each dimension.
"""
function unit_indices(N::Int) #create unit CartesianIndex for each dimension
    out = Vector{CartesianIndex{N}}(undef, N)
    null = zeros(Int64, N)
    for i in 1:N
        unit_i = copy(null)
        unit_i[i] = 1
        out[i] = CartesianIndex(Tuple(unit_i))
    end
    Tuple(out)
end

function cartesian_to_linear(I::CartesianIndex, s)  #Not sure if there is a builtin that does this - convert cartesian index to linear index of an array of size s
    out = I[1]
    for i in 1:length(s)-1
        out += (I[i+1]-1)*prod(s[1:i])
    end
    return out
end

# Counts the Differential operators for given variable x. This is used to determine
# the order of a PDE.
function count_differentials(term, x::Symbolics.Symbolic)
    S = Symbolics
    SU = SymbolicUtils
    if !S.istree(term)
        return 0
    else
        op = SU.operation(term)
        count_children = sum(map(arg -> count_differentials(arg, x), SU.arguments(term)))
        if op isa Differential && op.x === x
            return 1 + count_children
        end
        return count_children
    end
end

# return list of differential orders in the equation
function differential_order(eq, x::Symbolics.Symbolic)
    S = Symbolics
    SU = SymbolicUtils
    orders = Set{Int}()
    if S.istree(eq)
        op = SU.operation(eq)
        if op isa Differential
            push!(orders, count_differentials(eq, x))
        else
            for o in map(ch -> differential_order(ch, x), SU.arguments(eq))
                union!(orders, o)
            end
        end
    end
    return filter(!iszero, orders)
end

add_dims(A::AbstractArray, n::Int; dims::Int = 1) = cat(ndims(A) + n, A, dims = dims)

""
perpindex(A, dim::Integer) = A[setdiff(1:length(A), dim)]

"""
the size of A perpendicular to dim
"""
perpsize(A::AbstractArray{T,N}, dim::Integer) where {T,N} = size(A)[setdiff(1:N, dim)]


text:
# AbstractDiscretizer

abstract AbstractDiscretizer{N,D}
    # N indicates the decoded, or natural type
    # D indicates the encoded, or discrete type

encoded_type{N,D}(::AbstractDiscretizer{N,D}) = D
decoded_type{N,D}(::AbstractDiscretizer{N,D}) = N

# DiscretizationAlgorithm

abstract DiscretizationAlgorithm

function get_discretization_counts(disc::AbstractDiscretizer, data::AbstractArray)
    counts = zeros(Int, nlabels(disc))
    for v in data
        counts[encode(disc, v)] += 1
    end
    counts
end

function get_histogram_plot_arrays{R<:Real, I<:Integer}(binedges::Vector{R}, counts::AbstractVector{I})
    n = length(binedges)
    n == length(counts)+1 || error("binedges must have exactly one more entry than counts!")

    arr_x = Array{R}(n+1)
    arr_x[1] = binedges[1]
    arr_x[2:n+1] = binedges

    arr_y = Array{I}(n+1)
    arr_y[1] = zero(I)
    arr_y[2:n] = counts
    arr_y[end] = zero(I)

    (arr_x, arr_y)
end
function get_histogram_plot_arrays{R<:Real, F<:Real}(binedges::Vector{R}, pdfs::AbstractVector{F})
    n = length(binedges)
    n == length(pdfs)+1 || error("binedges must have exactly one more entry than counts!")

    arr_x = Array{R}(n+1)
    arr_x[1] = binedges[1]
    arr_x[2:n+1] = binedges

    arr_y = Array{F}(n+1)
    arr_y[1] = zero(F)
    arr_y[2:n] = pdfs
    arr_y[end] = zero(F)

    (arr_x, arr_y)
end


text:
@testset "FactorData Tests" begin

  f = PriorData(Z = Normal(0.0, 10.0))
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"Z\":{\"mu\":0.0,\"sigma\":10.0,\"_type\":\"IncrementalInference.PackedNormal\"}},\"multihypo\":[],\"certainhypo\":[1],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  f = PriorPose2Data(Z = FullNormal([0.0, 0.0, 0.0], diagm([0.01, 0.01, 0.01])))
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"Z\":{\"mu\":[0.0,0.0,0.0],\"cov\":[0.01,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.01],\"_type\":\"IncrementalInference.PackedFullNormal\"}},\"multihypo\":[],\"certainhypo\":[1],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  f = PriorPoint2Data(Z = FullNormal([0.0, 0.0], diagm([0.01, 0.01])))
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"Z\":{\"mu\":[0.0,0.0],\"cov\":[0.01,0.0,0.0,0.01],\"_type\":\"IncrementalInference.PackedFullNormal\"}},\"multihypo\":[],\"certainhypo\":[1],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  f = LinearRelativeData(Z = Normal(5.0, 10.0))
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"Z\":{\"mu\":5.0,\"sigma\":10.0,\"_type\":\"IncrementalInference.PackedNormal\"}},\"multihypo\":[],\"certainhypo\":[1,2],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  f = Pose2Pose2Data(Z = FullNormal([0.0, 0.0, 0.0], diagm([0.01, 0.01, 0.01])))
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"Z\":{\"mu\":[0.0,0.0,0.0],\"cov\":[0.01,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.01],\"_type\":\"IncrementalInference.PackedFullNormal\"}},\"multihypo\":[],\"certainhypo\":[1,2],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  # Pose2AprilTag4CornersData
  f = Pose2AprilTag4CornersData(
    0,
    zeros(8),
    zeros(9),
    K=[300.0, 0.0, 0.0, 0.0, 300.0, 0.0, 180.0, 120.0, 1.0],
    taglength=0.25
  )
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"corners\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\"homography\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\"K\":[300.0,0.0,0.0,0.0,300.0,0.0,180.0,120.0,1.0],\"taglength\":0.25,\"id\":0,\"_type\":\"/application/JuliaLang/PackedPose2AprilTag4Corners\"},\"multihypo\":[],\"certainhypo\":[1,2],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  # Pose2Point2BearingRangeData
  f = Pose2Point2BearingRangeData(
    bearing=Normal(0, 5),
    range=Normal(20, 1)
  )
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"bearstr\":{\"mu\":0.0,\"sigma\":5.0,\"_type\":\"IncrementalInference.PackedNormal\"},\"rangstr\":{\"mu\":20.0,\"sigma\":1.0,\"_type\":\"IncrementalInference.PackedNormal\"}},\"multihypo\":[],\"certainhypo\":[1,2],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  # Point2Point2RangeData
  f = Point2Point2RangeData(
    range=Normal(20, 1)
  )
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"Z\":{\"mu\":20.0,\"sigma\":1.0,\"_type\":\"IncrementalInference.PackedNormal\"}},\"multihypo\":[],\"certainhypo\":[1,2],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  # MixtureData
  f = MixtureData(
    LinearRelativeData,
    (hypo1=Normal(0, 2), hypo2=Uniform(30, 55)),
    [0.4, 0.6],
    2)
  @test JSON.json(f) == "{\"eliminated\":false,\"potentialused\":false,\"edgeIDs\":[],\"fnc\":{\"N\":2,\"F_\":\"PackedLinearRelative\",\"S\":[\"hypo1\",\"hypo2\"],\"components\":[{\"mu\":0.0,\"sigma\":2.0,\"_type\":\"IncrementalInference.PackedNormal\"},{\"a\":30.0,\"b\":55.0,\"_type\":\"IncrementalInference.PackedUniform\"}],\"diversity\":{\"p\":[0.4,0.6],\"_type\":\"IncrementalInference.PackedCategorical\"}},\"multihypo\":[],\"certainhypo\":[],\"nullhypo\":0.0,\"solveInProgress\":0,\"inflation\":3.0}"

  # ScatterAlignPose2 and implicityly ManifoldKernelDensity
  Random.seed!(42)
  pts1 = [randn(2) for i in 1:50]
  pts2 = [randn(2) for i in 1:50]
  sap = ScatterAlignPose2Data("Point2",pts1,pts2)
  fid = open(joinpath(@__DIR__,"testdata","sap_test.json"))
  ref = readline(fid)
  close(fid)
  str = JSON.json(sap);
  @warn "Suppress test, Random.seed! or reference string not working the same as local tests which passed."
  # @test ref == str

end


text:
# Cyber Problem
# POMDP, Robust POMDP, and Robust ρPOMDP formulations

struct CyberPOMDP <: POMDP{Vector{Int}, Vector{Int}, Vector{Int}}
    discount::Float64
    b0::Vector{Float64}
end

struct CyberRPOMDP <: RPOMDP{Vector{Int}, Vector{Int}, Vector{Int}}
    discount::Float64
    b0::Vector{Float64}
end

struct CyberIPOMDP <: IPOMDP{Vector{Int}, Vector{Int}, Vector{Int}}
    discount::Float64
    inforeward::Vector{Vector{Float64}}
    b0::Vector{Float64}
end

struct CyberRIPOMDP <: RIPOMDP{Vector{Int}, Vector{Int}, Vector{Int}}
    discount::Float64
    inforeward::Vector{Vector{Float64}}
    b0::Vector{Float64}
end

# Discount
const discount_cyber = 0.75

# State space
const nmoe = 3 # number of measures of effectiveness (MOE)
const nlevels = 3 # number of levels at which the MOE is measured
const lmin = 1 # min level
const lmax = 3 # max level
const nS = nlevels ^ nmoe # number of states
S = Vector{Vector{Int}}(nS) # states, as [moe_1_level, ...]
ind = 1
for i = 1:nlevels, j = 1:nlevels, k = 1:nlevels
  S[ind] = [i,j,k]
  ind += 1
end
const states_cyber = copy(S)

# Action space
const na = 2 # number of assessment assets
const nA = 6 # total number of actions
A = Vector{Vector{Int}}() # actions, as [asset_1_assignment, ...]
for i = 1:nmoe, j = 1:nmoe
  (i != j) && (push!(A, [i,j]))
end
const actions_cyber = copy(A)

# Observation space
const nZ = nlevels ^ na
Z = Vector{Vector{Int}}(nZ)
ind = 1
for i = 1:nlevels, j = 1:nlevels
  Z[ind] = [i,j]
  ind += 1
end
const observations_cyber = copy(Z)

# Belief reward (PWLC vectors: simple)
br = Vector{Vector{Float64}}(nS)
br[1] = vcat(1.0, fill(-1/(nS-1), nS - 1))
br[nS] = vcat(fill(-1/(nS-1), nS - 1), 1.0)
for i = 2:(nS-1)
  br[i] = vcat(fill(-1/(nS-1), i - 1), 1.0, fill(-1/(nS-1), nS - i))
end
const inforeward_cyber = copy(br)

# Initial state and belief
e11 = zeros(nS); e11[1] = 1.0;
const e1 = copy(e11)        # intial state
const b0_cyber = copy(e11)  # initial belief

CyberPOMDP() = CyberPOMDP(discount_cyber, b0_cyber)

CyberRPOMDP(disc::Float64) = CyberRPOMDP(disc, b0_cyber)
CyberRPOMDP() = CyberRPOMDP(discount_cyber, b0_cyber)

CyberIPOMDP(alphas::Vector{Vector{Float64}}) = CyberIPOMDP(discount_cyber, alphas, b0_cyber)
CyberIPOMDP(disc::Float64) = CyberIPOMDP(disc, inforeward_cyber, b0_cyber)
CyberIPOMDP() = CyberIPOMDP(discount_cyber, inforeward_cyber, b0_cyber)

CyberRIPOMDP(disc::Float64, alphas::Vector{Vector{Float64}}) = CyberRIPOMDP(disc, alphas, b0_cyber)
CyberRIPOMDP(disc::Float64) = CyberRIPOMDP(disc, inforeward_cyber, b0_cyber)
CyberRIPOMDP() = CyberRIPOMDP(discount_cyber, inforeward_cyber, b0_cyber)

states(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = states_cyber
actions(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = actions_cyber
observations(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = observations_cyber

n_states(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = nS
n_actions(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = nA
n_observations(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = nZ

state_index(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}, s::Vector{Int}) = find(x -> x == s, states_cyber)[1]
action_index(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}, a::Vector{Int}) = find(x -> x == a, actions_cyber)[1]
observation_index(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}, z::Vector{Int}) = find(x -> x == z, observations_cyber)[1]
obs_index(prob::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}, z::Vector{Int}) = observation_index(prob, z)

initial_state_distribution(::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = SparseCat(states_cyber, e1)
initial_belief(prob::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = prob.b0
initial_belief_distribution(prob::Union{CyberPOMDP, CyberRPOMDP, CyberIPOMDP, CyberRIPOMDP}) = SparseCat(states_cyber, initial_belief(prob))

# Transitions
const pϵ_cyber = 1e-6 # near-saturation bound
const delt = 0.01 # imprecision in transition distribution
const pd = 0.15 # prob of decrease
const ps = 0.55 # prob of stay
const pi = 0.3 # prob of increase
const psd = ps + pd # prob of decline or stay at lower border
const psi = ps + pi # prob of stay or improve at upper border

# transition function: probabilty of going from state s to state j
function psj(s::Vector{Int},j::Vector{Int})
  prob = 1.0
  count = 0
  for i = 1:nmoe
    level_s = s[i]
    level_j = j[i]
    if level_s == lmin # current state is lowest level
      if level_j == level_s
        prob = prob * psd
      elseif level_j == level_s + 1
        prob = prob * pi
      else # zero probability of skipping levels
        return prob = 0.0
      end
    elseif level_s == lmax # current state is highest level
      if level_j == level_s
        prob = prob * psi
      elseif level_j == level_s - 1
        prob = prob * pd
      else # zero probability of skipping levels
        return prob = 0.0
      end
    else # current state is an intermediate level
      if level_j == level_s
        prob = prob * ps
        count += 1
      elseif level_j == level_s + 1
        prob = prob * pi
      elseif level_j == level_s - 1
        prob = prob * pd
      else # zero probability of skipping levels
        return prob = 0.0
      end
    end
  end
  prob
end

# robust transition function: probabilty of going from state s to state j
function psj(s::Array{Int,1}, j::Array{Int,1}, par, delt)
  # add imprecision (one-sided)
  delt > 0 ? (pd,ps,pi) = min.(par .+ delt, 1 - pϵ_cyber) : (pd,ps,pi) = max.(par .+ delt, 0 + pϵ_cyber)
  psd = min(ps + pd, 1 - pϵ_cyber) # prob of decline or stay at lower border
  psi = min(ps + pi, 1 - pϵ_cyber) # prob of stay or improve at upper border
  prob = 1.0
  for i = 1:nmoe
    level_s = s[i]
    level_j = j[i]
    if level_s == lmin # current state is lowest level
      if level_j == level_s
        prob = prob * psd
      elseif level_j == level_s + 1
        prob = prob * pi
      else # zero probability of skipping levels
        return prob = 0.0
      end
    elseif level_s == lmax # current state is highest level
      if level_j == level_s
        prob = prob * psi
      elseif level_j == level_s - 1
        prob = prob * pd
      else # zero probability of skipping levels
        return prob = 0.0
      end
    else # current state is an intermediate level
      if level_j == level_s
        prob = prob * ps
      elseif level_j == level_s + 1
        prob = prob * pi
      elseif level_j == level_s - 1
        prob = prob * pd
      else # zero probability of skipping levels
        return prob = 0.0
      end
    end
  end
  prob
end

# calculate trasnsition matrix
function calcTArray(S::Vector{Vector{Int}}, A::Vector{Vector{Int}}, par, delt)
  ns = length(S)
  na = length(A)
  T = zeros(ns, na, ns)
  Tl = zeros(ns, na, ns)
  Tu = zeros(ns, na, ns)
  for (si,s) in enumerate(S), (ji,j) in enumerate(S)
      T[si,:,ji] = psj(s,j)
      tl = psj(s,j,par,-delt)
      tu = psj(s,j,par,delt)
      # tl = max(tl, 0.0 + pϵ_cyber)
      # tu = min(tu, 1.0 - pϵ_cyber)
      (tu < tl) && (tu = tl + pϵ_cyber / 2)
      Tl[si,:,ji] = tl
      Tu[si,:,ji] = tu
  end
  for (si,s) in enumerate(S), (ai,a) in enumerate(A)
    T[si,ai,:] = T[si,ai,:] ./ sum(T[si,ai,:])
  end
  T, Tl, Tu
end

# Nominal transition function array tarray_cyber[s,a,s'] = Pr(s' | s ,a)
const tarray_cyber, tlarray_cyber, tuarray_cyber = calcTArray(states_cyber, actions_cyber, (pd,ps,pi), delt)

# calculate trasnsition distribution
function calcTDist(ns::Int, na::Int, tarray::Array{Float64,3})
  T = Array{SparseCat}(ns, na)
  for si = 1:ns, ai = 1:na
      T[si,ai] = SparseCat(states_cyber, tarray[si,ai,:])
  end
  T
end

# Nominal transition function distributions
const tdist_cyber = calcTDist(nS, nA, tarray_cyber)

# Nominal transitions
transition(prob::Union{CyberIPOMDP, CyberPOMDP}, s::Vector{Int}, a::Vector{Int}) = tdist_cyber[state_index(prob, s), action_index(prob, a)]
transition(prob::Union{CyberIPOMDP, CyberPOMDP}, s::Vector{Int}, a::Vector{Int}, sp::Vector{Int}) = tarray_cyber[state_index(prob, s), action_index(prob, a), state_index(prob, sp)]

# Robust transitions
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PtLower = Plower(t|s,a) = t_lower[s,a,sp]
function transition(prob::Union{CyberRPOMDP, CyberRIPOMDP}, s::Vector{Int}, a::Vector{Int})
    si, ai = state_index(prob, s), action_index(prob, a)
    tlarray_cyber[si,ai,:], tuarray_cyber[si,ai,:]
end
function transition(prob::Union{CyberRPOMDP, CyberRIPOMDP}, s::Vector{Int}, a::Vector{Int}, sp::Vector{Int})
    si, ai, spi = state_index(prob, s), action_index(prob, a), state_index(prob, sp)
    tlarray_cyber[si,ai,spi], tuarray_cyber[si,ai,spi]
end

# Observation function
const acc1 = 0.8 # likelihood of correct observation error for sensor 1
const acc2 = 0.9 # likelihood of correct observation error for sensor 2
const sensor_accuracy = [acc1, acc2]
const delo1 = pϵ_cyber # sensor 1 immprecision
const delo2 = 0.4 # sensor 2 imprecision
const sensor_imprecision = [delo1, delo2]

# observation function
function o(a::Vector{Int}, sp::Vector{Int}, z::Vector{Int}, accvec::Vector{Float64})
  pz = 1.0
  nsensor = length(a)
  for i = 1:nsensor
    assingment_sensor = a[i]
    level_sp = sp[assingment_sensor]
    level_z = z[i]
    pacc = accvec[i]
    if level_sp == lmin # current state is lowest level
      if level_z == level_sp
        pz = pz * pacc
      elseif level_z == level_sp + 1
        pz = pz * (1 - pacc)
      else # zero probability of error > 1
        return pz = 0.0
      end
    elseif level_sp == lmax # current state is highest level
      if level_z == level_sp
        pz = pz * pacc
      elseif level_z == level_sp - 1
        pz = pz * (1 - pacc)
      else # zero probability of error > 1
        return pz = 0.0
      end
    else # current state is an intermediate level
      if level_z == level_sp
        pz = pz * pacc
      elseif level_z == level_sp + 1
        pz = pz * (1 - pacc) / 2
      elseif level_z == level_sp - 1
        pz = pz * (1 - pacc) / 2
      else # zero probability of error > 1
        return pz = 0.0
      end
    end
  end
  return pz
end

# calculate trasnsition matrix
function calcOArray(S::Vector{Vector{Int}}, A::Vector{Vector{Int}}, Z::Vector{Vector{Int}})
  O = zeros(nA,nS,nZ)
  Ol = zeros(nA,nS,nZ)
  Ou = zeros(nA,nS,nZ)
  for (spi,sp) in enumerate(S), (zi,z) in enumerate(Z), (ai,a) in enumerate(A)
    O[ai, spi, zi] = o(a, sp, z, sensor_accuracy)
    neg = o(a, sp, z, sensor_accuracy - sensor_imprecision)
    pos = o(a, sp, z, sensor_accuracy + sensor_imprecision)
    ol = min(neg, pos)
    ou = max(neg, pos)
    (ol < 0.0) && (ol = max(ol, 0.0 + pϵ_cyber))
    (ou > 1.0) && (ou = min(ou, 1.0 + pϵ_cyber))
    (ou < ol) && (ou = ol + pϵ_cyber / 2)
    Ol[ai, spi, zi] = ol
    Ou[ai, spi, zi] = ou
  end
  for (spi,sp) in enumerate(S), (ai,a) in enumerate(A)
    O[ai, spi,:] = O[ai,spi,:] ./ sum(O[ai,spi,:])
  end
  O, Ol, Ou
end

const oarray_cyber, olarray_cyber, ouarray_cyber = calcOArray(states_cyber, actions_cyber, observations_cyber)

# calculate observation distribution
function calcODist(ns::Int, na::Int, oarray::Array{Float64,3})
  O = Array{SparseCat}(na, ns)
  for ai = 1:na, si = 1:ns
      O[ai,si] = SparseCat(observations_cyber, oarray[ai,si,:])
  end
  O
end

# Nominal transition function distributions
const odist_cyber = calcODist(nS, nA, oarray_cyber)

# Nominal transitions
observation(prob::Union{CyberIPOMDP, CyberPOMDP}, a::Vector{Int}, sp::Vector{Int}) = odist_cyber[action_index(prob, a), state_index(prob, sp)]
observation(prob::Union{CyberIPOMDP, CyberPOMDP}, a::Vector{Int}, sp::Vector{Int}, z::Vector{Int}) = oarray_cyber[action_index(prob, a), state_index(prob, sp), observation_index(prob, z)]

# Robust observations
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PzLower = Plower(z|s,a,t) = o_lower[a,sp,z]
function observation(prob::Union{CyberRPOMDP, CyberRIPOMDP}, a::Vector{Int}, sp::Vector{Int})
    ai, spi = action_index(prob, a), state_index(prob, sp)
    olarray_cyber[ai,spi,:], ouarray_cyber[ai,spi,:]
end
function observation(prob::Union{CyberRPOMDP, CyberRIPOMDP}, a::Vector{Int}, sp::Vector{Int}, z::Vector{Int})
    ai, spi, zi = action_index(prob, a), state_index(prob, sp), observation_index(prob, z)
    olarray_cyber[ai,spi,zi], ouarray_cyber[ai,spi,zi]
end

function dynamics(prob::Union{CyberPOMDP,CyberIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        pt = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pz = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        p[tind, zind, sind, aind] = pz * pt
    end
    p
end

function dynamics(prob::Union{CyberRPOMDP,CyberRIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        ptl, ptu = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pzl, pzu = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        plower[tind, zind, sind, aind] = pzl * ptl
        pupper[tind, zind, sind, aind] = pzu * ptu
    end
    plower, pupper
end

function reward(prob::Union{CyberPOMDP, CyberRPOMDP}, s::Vector{Int}, a::Vector{Int})
  u = [0.5, 0.4, 0.1] # utility of each MOE
  r = 0.0
  for i = 1:length(a)
    pacc = sensor_accuracy[i]
    r += u[a[i]] * pacc
  end
  r
end
reward(prob::Union{CyberPOMDP,CyberRPOMDP}, b::Vector{Float64}, s::Vector{Int}, a::Vector{Int}, sp::Vector{Int}) = reward(prob,s,a)

function reward(prob::Union{CyberIPOMDP, CyberRIPOMDP}, b::Vector{Float64}, a::Vector{Int})
    rmax = -Inf
    for α in prob.inforeward
        rmax = max(rmax, dot(α, b))
    end
    rmax
end
reward(prob::Union{CyberIPOMDP,CyberRIPOMDP}, b::Vector{Float64}, s::Vector{Int}, a::Vector{Int}, sp::Vector{Int}) = reward(prob,b,a)

function rewardalpha(prob::Union{CyberIPOMDP, CyberRIPOMDP}, b::Vector{Float64}, a::Vector{Int})
    rmax = -Inf
    ralpha = nothing
    for α in prob.inforeward
        rnext = dot(α, b)
        if rnext > rmax
            rmax = rnext
            ralpha = copy(α)
        end
    end
    ralpha
end

discount(p::Union{CyberPOMDP,CyberRPOMDP,CyberIPOMDP,CyberRIPOMDP}) = p.discount

# Simulation Functions
function initial_state(prob::Union{CyberPOMDP,CyberIPOMDP,CyberRPOMDP,CyberRIPOMDP}, rng::AbstractRNG)
    d = initial_state_distribution(prob)
    return rand(rng, d)
end

function generate_sor(prob::Union{CyberPOMDP,CyberIPOMDP}, b, s, a, rng::AbstractRNG)
    sp = rand(rng, transition(prob, s, a))
    o = rand(rng, observation(prob, a, sp))
    r = reward(prob, b, s, a, sp)
    sp, o, r
end

function generate_sor(prob::Union{CyberRPOMDP,CyberRIPOMDP}, b, s, a, rng::AbstractRNG)
    tdist = SparseCat(states(prob), psample(transition(prob, s, a)...))
    sp = rand(rng, tdist)
    odist = SparseCat(observations(prob), psample(observation(prob, a, sp)...))
    o = rand(rng, odist)
    r = reward(prob, b, s, a, sp)
    sp, o, r
end


text:
using MakieLayout
using AbstractPlotting
using Test

include("debugrect.jl")

@testset "GridLayout 1" begin
    bbox = BBox(0, 1000, 0, 1000)
    layout = GridLayout(bbox = bbox, alignmode = Outside(0))
    dr = layout[1, 1] = DebugRect()

    @test dr.layoutnodes.computedbbox[] == bbox

    dr.topprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 900)
    dr.bottomprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 100, 900)
    dr.leftprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 1000, 100, 900)
    dr.rightprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 900, 100, 900)

    dr2 = layout[1, 2] = DebugRect()
    @test layout.nrows == 1 && layout.ncols == 2
    colgap!(layout, 1, Fixed(0))

    @test dr.layoutnodes.computedbbox[].widths == dr2.layoutnodes.computedbbox[].widths == Float32[400.0, 800.0]
end


text:
# Functions to parse strings to intervals

"""
    parse{T}(DecoratedInterval{T}, s::AbstractString)

Parse a string of the form `"[a, b]_dec"` as a `DecoratedInterval`
with decoration `dec`.
"""
function parse(::Type{DecoratedInterval{T}}, s::AbstractString) where T
    m = match(r"(\[.*\])(\_.*)?", s)

    if m == nothing  # matched

        m = match(r"(.*\?[a-z0-9]*)(\_.*)?", s)

        if m == nothing
            throw(ArgumentError("Unable to process string $x as decorated interval"))
        end

    end

    interval_string, decoration_string = m.captures
    interval = parse(Interval{T}, interval_string)

    # type unstable:
    if decoration_string == nothing
        decoration_string = "_com"
    end

    decoration_symbol = Symbol(decoration_string[2:end])
    decoration = getfield(IntervalArithmetic, decoration_symbol)

    return DecoratedInterval(interval, decoration)

end

"""
    parse{T}(Interval{T}, s::AbstractString)

Parse a string as an interval. Formats allowed include:
- "1"
- "[1]"
- "[3.5, 7.2]"
- "[-0x1.3p-1, 2/3]"  # use numerical expressions
"""
function parse(::Type{Interval{T}}, s::AbstractString) where T

    # Check version!
    if !(occursin("[", s))  # string like "3.1"

        m = match(r"(.*)±(.*)", s)
        if m != nothing
            a = parse(T, strip(m.captures[1]))
            b = parse(T, strip(m.captures[2]))

            return a ± b

        a = parse(T, s, RoundDown)
        b = parse(T, s, RoundUp)

        return Interval(a, b)

        end

        if m == nothing

            m = match(r"(\-?\d*\.?\d*)\?\?([ud]?)", s)   # match with string of form "34??"

            if m!= nothing
                if m.captures[2] == ""    # strings of the form "10??"
                    return entireinterval(T)
                end
                if m.captures[2] == "u"   # strings of the form "10??u"
                    lo = parse(Float64, m.captures[1])
                    hi = Inf
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end
                if m.captures[2] == "d"   # strings of the form "10??d"
                    lo = -Inf
                    hi = parse(Float64, m.captures[1])
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end
            end

            m = match(r"(\-?\d*)\.?(\d*)\?(\d*)([ud]?)(e-?\d*)?", s) # match with strings like "3.56?1u" or "3.56?1e2"

            if m != nothing  # matched
                if m.captures[3] != "" && m.captures[4] == "" && m.captures[5] == nothing # string of form "3.4?1" or "10?2"
                    d = length(m.captures[2])
                    x = parse(Float64, m.captures[3] * "e-$d")
                    n = parse(Float64, m.captures[1]*"."*m.captures[2])
                    lo = n - x
                    hi = n + x
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "" && m.captures[5] == nothing  # strings of the form "3.46?"
                    d = length(m.captures[2])
                    x = parse(Float64, "0.5" * "e-$d")
                    n = parse(Float64, m.captures[1]*"." * m.captures[2])
                    lo = n - x
                    hi = n + x
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "" && m.captures[5] != nothing  # strings of the form "3.46?e2"
                    d = length(m.captures[2])
                    x = parse(Float64, "0.5" * "e-$d")
                    n = parse(Float64, m.captures[1]*"." * m.captures[2])
                    lo = parse(Float64, string(n-x) * m.captures[5])
                    hi = parse(Float64, string(n+x) * m.captures[5])
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "u" && m.captures[5] == nothing # strings of the form "3.4?u"
                    d = length(m.captures[2])
                    x = parse(Float64, "0.5" * "e-$d")
                    n = parse(Float64, m.captures[1]*"."*m.captures[2])
                    lo = n
                    hi = n+x
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "d" && m.captures[5] == nothing # strings of the form "3.4?d"
                        d = length(m.captures[2])
                        x = parse(Float64, "0.5" * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = n - x
                        hi = n
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                end

                if m.captures[3] != ""
                    if m.captures[4] == "u" && m.captures[5] != nothing    #strings of the form "3.46?1u"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = parse(Float64, string(n) * m.captures[5])
                        hi = parse(Float64, string(n+x) * m.captures[5])
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "u" && m.captures[5] == nothing    #strings of the form "3.46?1u"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = n
                        hi = n+x
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "d" && m.captures[5] != nothing  #strings of the form "3.46?1d"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = parse(Float64, string(n-x) * m.captures[5])
                        hi = parse(Float64, string(n) * m.captures[5])
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "d" && m.captures[5] == nothing  #strings of the form "3.46?1d"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = n-x
                        hi = n
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "" && m.captures[5] != nothing    # strings of the form "3.56?1e2"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = parse(Float64, string(n-x)*m.captures[5])
                        hi = parse(Float64, string(n+x)*m.captures[5])
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end
                end
            end

            m = match(r"(-?\d*\.?\d*)", s)  # match strings of form "1" and "2.4"

            if m != nothing
                lo = parse(Float64, m.captures[1])
                hi = lo
            end

            if m == nothing
                throw(ArgumentError("Unable to process string $s as interval"))
            end

            interval = eval(make_interval(T, lo, [hi]))
            return interval
        end
    end

    # match string of form [a, b]_dec:
    m = match(r"\[(.*),(.*)\]", s)

    if m != nothing  # matched

        m.captures[1] = strip(m.captures[1], [' '])
        m.captures[2] = strip(m.captures[2], [' '])
        lo, hi = m.captures

        if m.captures[2] == "+infinity" || m.captures[2] == ""
            hi = "Inf"
        end

        if m.captures[1] == "-infinity" || m.captures[1] == ""
            lo = "-Inf"
        end

    end

    if m == nothing

        m = match(r"\[(.*)\]", s)  # string like "[1]"

        if m == nothing
            throw(ArgumentError("Unable to process string $s as interval"))
        end

        m.captures[1] = strip(m.captures[1], [' '])

        if m.captures[1] == "Empty" || m.captures[1] == ""
            return emptyinterval(T)
        end
        if m.captures[1] == "entire"
            return entireinterval(T)
        end
        if m.captures[1] == "nai" || m.captures[1] == "Nai"
            return nai(T)
        end

        lo = m.captures[1]
        hi = lo

    end

    expr1 = Meta.parse(lo)
    expr2 = Meta.parse(hi)
    if lo == "-Inf"
        expr1 = parse(Float64, lo)
    end
    if hi == "Inf"
        expr2 = parse(Float64, hi)
    end

    interval = eval(make_interval(T, expr1, [expr2]))

    return interval

end


text:
module Slices
export Slice, SliceHom

using AutoHashEquals
using ...GAT
using ..FreeDiagrams, ..Limits, ..CSets
using ...Theories: Category
import ...Theories: dom, codom, compose, id
import ..Limits: limit, colimit, universal
import ..FinSets: force, pushout_complement

"""
The data of the object of a slice category (say, some category C sliced over an
object X in Ob(C)) is the data of a homomorphism in Hom(A,X) for some ob A.
"""
@auto_hash_equals struct Slice{Hom}
  slice::Hom
end

"""
The data of the morphism of a slice category (call it h, and suppose a category
C is sliced over an object X in Ob(C)) between objects f and g is a homomorphism
in the underlying category that makes the following triangle commute.

   h
A --> B
f ↘ ↙ g
   X
"""
@auto_hash_equals struct SliceHom{Hom, Dom<:Slice, Codom<:Slice}
  dom::Dom
  codom::Codom
  f::Hom
end

function SliceHom(d::Dom, cd::Codom, f::Hom; strict::Bool=true) where {Dom,Codom,Hom}
  !strict || codom(d) == codom(cd) || error("$d and $cd not in same category")
  !strict || dom(d) == dom(f) || error("dom $d does not match $f")
  !strict || dom(cd) == codom(f) || error("codom $cd does not match $f")
  return SliceHom{Hom,Dom,Codom}(d, cd, f)
end

dom(s::Slice) = dom(s.slice)
codom(s::Slice) = codom(s.slice)
force(s::Slice)  = Slice(force(s.slice))
force(s::SliceHom) = SliceHom(
  force(dom(s)), force(codom(s)), force(s.f))

struct SliceLimit{Hom, Ob <: Slice{Hom}, Diagram,
                  Cone <: Multispan{Ob}} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  underlying::AbstractLimit
end

struct SliceColimit{Hom, Ob <: Slice{Hom}, Diagram,
                    Cocone <: Multicospan{Ob}} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  underlying::AbstractColimit
end

@instance Category{Slice, SliceHom} begin
  dom(f::SliceHom) = f.dom
  codom(f::SliceHom) = f.codom
  id(A::Slice) = SliceHom(A, A, id(dom(A.slice)))
  function compose(f::SliceHom, g::SliceHom; strict::Bool=false)
    !strict || codom(f) == dom(g) ||
      error("Domain mismatch in composition: $(codom(f)) != $(dom(g))")
    SliceHom(dom(f), codom(g), compose(f.f, g.f))
  end
end

"""
Get the underlying diagram of a slice category diagram which has the object
being sliced over explicitly, and arrows are ACSetTransformations

Encode the base of slice as the first object in the diagram
"""
function slice_diagram(f::FreeDiagram)::FreeDiagram
  obs = vcat([codom(f[:ob][1])], [dom(x.slice) for x in f[:ob]])
  homs = [(h.slice,i+1,1) for (i, h) in enumerate(f[:ob])]
  append!(homs, [(h.f, i+1, j+1) for (i,j,h) in zip(f[:src],f[:tgt],f[:hom])])
  FreeDiagram(obs,homs)
end

"""
Convert a limit problem in the slice category to a limit problem of the
underlying category.
"""
function limit(::Type{Tuple{S,H}}, diagram::FreeDiagram) where {S<:Slice, H<:SliceHom}
  lim = limit(slice_diagram(diagram))
  new_apex = Slice(first(legs(lim.cone)))
  new_cone_legs = [SliceHom(new_apex, ob, leg) for (ob, leg)
                   in zip(diagram[:ob],legs(lim.cone)[2:end])]
  return SliceLimit(diagram, Multispan(new_apex, new_cone_legs), lim)
end

colimit(s::Multispan{<:Slice}) = colimit(FreeDiagram{Slice,SliceHom}(s))

function colimit(::Type{Tuple{S,H}}, diagram::FreeDiagram) where {S<:Slice, H<:SliceHom}
  nparts(diagram, :V) > 0 ||
    error("Requires nonempty diagram in order to know what is sliced over")

  # discard all the slice info in the colimit diagram - it's irrelevant
  colim = colimit(map(diagram, ob = x -> dom(x.slice), hom = h -> h.f))

  # compute new apex using the universal property of colimits
  X = codom(first(diagram[:ob]))
  csp = Multicospan(X, map(x -> x.slice, diagram[:ob]))
  new_apex = Slice(universal(colim, csp))
  new_cocone_legs = [SliceHom(o, new_apex, l)
                     for (o, l) in zip(diagram[:ob],legs(colim))]
  return SliceColimit(diagram, Multicospan(new_apex, new_cocone_legs), colim)
end

function universal(lim::SliceLimit, sp::Multispan)
  # Use the universal property of the underlying category.
  apx = apex(sp)
  newspan = vcat([apx.slice],[x.f for x in sp])
  u = universal(lim.underlying, Multispan(dom(apx), newspan))
  apx2 = Slice(first(legs(lim.underlying.cone)))
  return SliceHom(apx, apx2, u)
end

"""    pushout_complement(f::SliceHom, g::SliceHom)
Compute a pushout complement in a slice category by using the pushout complement
in the underlying category.

     f
  B <-- A ---⌝
  | ↘ ↙      |
 g|  X       | f′
  ↓ ↗  ↖ cx  |
  D <--- C <--
      g′

"""
function pushout_complement(f::SliceHom, g::SliceHom)
    f′, g′ = pushout_complement(ComposablePair(f.f, g.f))
    D = codom(g)
    C = Slice(compose(g′, D.slice))
    return SliceHom(dom(f), C, f′) => SliceHom(C, D, g′)
end

end # module

max_forks_repo_name: epatters/CompCat


text:
"""
    Origin{T} <: Node{T}

    Global reference frame

    id: always 0
    name: always :origin
    state: State; defaults to zero values
    shape: Shape; defaults to EmptyShape
"""
mutable struct Origin{T} <: Node{T}
    id::Int64
    name::Symbol
    state::State{T}
    shape::Shape{T}

    function Origin{T}(;
        name::Symbol=:origin,
        state=State{T}(),
        shape::Shape=EmptyShape()) where T
        new{T}(getGlobalID(), name, state, shape)
    end

    function Origin(;
        name::Symbol=:origin,
        state=State{Float64}(),
        shape::Shape=EmptyShape())
        Origin{Float64}(; name, state, shape)
    end

    function Origin(body::Body{T}) where T
        new{T}(body.id, body.name, body.state, body.shape)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, origin::Origin)
    summary(io, origin)
    println(io, "")
    println(io, " id:      "*string(origin.id))
    println(io, " name:    "*string(origin.name))
end


text:
using BFloat
using Test

b1 = BFloat16(3.14f0)
b2 = BFloat16(1.2f0)

@testset "BFloat.jl" begin
    # TODO: Figure out which tests define FP numbers
    @test b1 + b2 == BFloat16(Float32(b1) + Float32(b2))

    # Conversions

    # Promotions

    # Functions of
end