text:
"""
    JGraph

Maintain the graph state comprising of nodes, edges, layout, animation ordering etc.

This will be a part of the Javis [`Object`](@ref) metadata, when a new graph is created.

# Fields
- `adjacency_list`: A light internal representation of the graph structure. Can be initialized with a known graph type.
    - For undirected graphs the underlying graph type used is `SimpleGraph` from [Graphs.jl]() and for directed graphs it is `SimpleDiGraph`.
- `width::Int`: The width of the graph on the canvas.
- `height::Int`: The height of the graph on the canvas.
- `mode::Symbol`: The animation of the graph can be done in two ways.
    - `static`: A lightweight animation which does not try to animate every detail during adding and deletion of nodes.
    - `dynamic`: The graph layout change is animated on addition of a new node. Can be computationally heavy depending on the size of graph.
- `layout::Symbol`: The graph layout to be used. Can be one of :-
    - `:none`
    - `:spring`
    - `:spectral`
- `get_node_attribute`: A function that enables fetching properties defined for nodes in the input graph data structure.
    - Required only when a node property like `cost` needs to be mapped to a drawing property like `radius`.
- `get_edge_attribute`: Similar to `get_node_attribute` but for edge properties.
- `ordering`: Store the relative ordering used to add nodes and edges to a graph using [`GraphNode`](@ref) and [`GraphEdge`](@ref)
    - If input graph is of a known type, defaults to a simple BFS ordering starting at the root node.
- `node_property_limits`: The minima and maxima calculated on the node properties in the input graph.
    - This is internally created and updated when [`updateGraph`](@ref) or the final render function is called.
    - This is skipped for node properties of non-numeric types.
    - Used to scale drawing property values within sensible limits.
- `edge_property_limits`: The minima and maxima calculated on the edge properties in the input graph.
    - Similar to `node_attribute_fn`.
"""
struct JGraph <: AbstractJavisGraph
    width::Int
    height::Int
    mode::Symbol
    layout::Symbol
    graph::ReferenceGraph
    object::AbstractObject
    opts::Dict{Symbol, Any}
    ordering::Vector{AbstractGraphElement}
end

GRAPHS = Vector{AbstractJavisGraph}()
CURRENT_GRAPH = Array{AbstractJavisGraph, 1}()

"""
    JGraph(directed::Bool, width::Int, height::Int)

Create an empty graph on the canvas.
"""
JGraph(directed::Bool, width::Int, height::Int, frames=:same; layout::Symbol=:spring) =
    directed ? JGraph(ReferenceGraph(SimpleDiGraph()), width, height, frames; layout=layout) :
    JGraph(ReferenceGraph(SimpleGraph()), width, height, frames; layout=layout)

JGraph(directed::Bool, width::Int, height::Int, start_pos::Point, frames=:same; layout::Symbol=:spring) =
    directed ? JGraph(ReferenceGraph(SimpleDiGraph()), width, height, frames, start_pos; layout=layout) :
    JGraph(ReferenceGraph(SimpleGraph()), width, height, frames, start_pos; layout=layout)

"""
    JGraph(graph, width::Int, height::Int; <keyword arguments>)

Creates a Javis object for the graph and assigns its `Metadata` field to the object created by this struct.

# Arguments
- `graph`: A known data structure storing information about nodes, edges, properties, etc.
    - Graph types of type `AbstractGraph` from the [Graphs.jl]() package are supported.
    - Using this eliminates the requirement to create nodes and edges separately.
- `width::Int`: Size of the graph along the horizontal direction.
- `height::Int`: Size of the graph along the vertical direction.

# Keywords
- `mode`: `:static` or `:dynamic`. Default is `:static`.
- `layout`: `:none`, `:spring` or `:spectral`. Default is `:spring`.
- `get_node_attribute::Function`: A `Function` with a signature `(graph, node::Int, attr::Any)::Any`
    - `graph` is the object containing node properties
    - Returns a value corresponding to the type of the node property `attr`.
        - Must be either a Julia primitive type or `String`.
- `get_edge_attribute::Function`: A `Function` with a signature `(graph, from_node::Int, to_node::Int, attr::Any)::Any`
    - Similar to `get_node_attribute`.

# Examples
```julia
using Javis
function ground()
    background("white")
    sethue("black")
end

video = Video(300, 300)
Background(1:100, ground)
# A star graph
graph = [[2, 3, 4, 5, 6], [], [], [], [], []]
ga = @Object(1:100, JGraph(false, 100, 100), O)
render(video; pathname="graph_animation.gif")
```

# Implementation
To be filled in ...

"""
function JGraph(
    graph::ReferenceGraph,
    width::Int,
    height::Int,
    frames=:same,
    start_pos::Point=O;
    mode::Symbol = :static,
    layout::Symbol = :spring
)
    # Check available layouts
    if !(layout in [:none, :spring, :spectral])
        layout = :spring
        @warn "Unknown layout '$(layout)', defaulting to a spring layout"
    end

    styles = OrderedDict{Symbol, Function}()
    styles[:global_property_limits] = _global_property_limits
    if mode == :static
        styles[:global_layout] = _global_layout
    elseif mode == :dynamic
    end
    object = Object(frames, get_draw(:graph), start_pos; _graph_idx = length(GRAPHS)+1)
    opts = Dict{Symbol, Any}()
    opts[:styles] = styles
    object.opts[:_style_opts_cache] = Dict{Symbol, Any}()
    jgraph = JGraph(
        width,
        height,
        mode,
        layout,
        graph,
        object,
        opts,
        Vector{AbstractGraphElement}(),
    )
    push!(GRAPHS, jgraph)
    if isempty(CURRENT_GRAPH)
        push!(CURRENT_GRAPH, jgraph)
    else
        CURRENT_GRAPH[1] = jgraph
    end
    return jgraph
end

function _global_property_limits(video, object, frames; kwargs...)
    g = GRAPHS[object.opts[:_graph_idx]]
    for el in g.ordering
        get_property = typeof(el) <: AbstractGraphVertex ? g.graph.get_vertex_property : g.graph.get_edge_property
        limits = typeof(el) <: AbstractGraphVertex ? g.graph.vertex_property_limits : g.graph.edge_property_limits
        for (k, _) in el.style_property_map
            val = get_property(k)
            if !(typeof(val) <: Real)
                throw("Cannot calculate limits. Property $(k) of $(typeof(el)) is not of `Real` type")
            end
            if !(k in keys(limits))
                limits[k] = (val, val)
            end
            limits[k] = (min(limits[k][1], val), max(limits[k][2], val))
        end
    end
end

function Graphs.edges(g::JGraph, in_order=false)
    e = []
    if !in_order
        return g.ordering[map(e -> get_prop(g.graph.adjacency_graph, src(e), dst(e)), edges(g.graph.adjacency_graph))]
    end
    for el in g.ordering
        if el isa AbstractGraphEdge
            push!(e, el)
        end
    end
    return e
end

function Graphs.vertices(g::JGraph, in_order=false)
    if !in_order
        return g.ordering[map(v -> get_prop(g.graph.adjacency_graph, v), vertices(g.graph.adjacency_graph))]
    end
    v = []
    for el in g.ordering
        if el isa AbstractGraphVertex
            push!(v, el)
        end
    end
    return v
end

function Graphs.neighbors(g::JGraph, v::Integer; strict=false)
    n = copy(neighbors(g.graph.adjacency_graph, v))
    if strict
        filter!(x -> x ≠ v, n)
    end
    return map(x -> g.ordering[x], get_prop.([g.graph.adjacency_graph], n))
end

function _global_layout(video, object, frame; kwargs...)
    g = GRAPHS[object.opts[:_graph_idx]]
    if frame == first(get_frames(object))
        layout_x = []
        layout_y = []
        if g.layout == :none
            object.opts[:layout_coords] = Vector{Point}()
            for v in vertices(g)
                push!(object.opts[:layout_coords], v.object.opts[:position])
            end
        else
            if nv(g.graph.adjacency_graph) <= 1
                layout_x, layout_y = [O], [O]
            elseif g.layout == :spring
                layout_x, layout_y = spring_layout(g.graph.adjacency_graph.graph)
            elseif g.layout == :spectral
                # Check special property layout_weight is defined on edges and collect weights
                weights = map(
                    (e) -> get(e.object.opts, :weight, 1),
                    edges(g, true),
                )
                # ToDo: Why does just g.graph.adjacency_graph not work? (is_directed not implemented error)
                layout_x, layout_y = spectral_layout(g.graph.adjacency_graph.graph, weights)
            end
            # Coordinates are between -1 and 1
            coords = map((p) -> Point(p), collect(zip(layout_x, layout_y)))
            # Scale to graph dimensions
            coords = coords .* [(g.width/2, g.height/2)]
            object.opts[:layout_coords] = coords
        end
    end
    # Now assign positions back to all nodes
    for (idx, v) in  enumerate(vertices(g))
        v.object.opts[:position] = object.opts[:layout_coords][idx] + object.start_pos
    end
    # Define keyword arguments for edges defining endpoint position
    for e in edges(g)
        from_vertex = g.ordering[get_prop(g.graph.adjacency_graph, e.from_vertex.vertex_id)]
        to_vertex = g.ordering[get_prop(g.graph.adjacency_graph, e.to_vertex.vertex_id)]
        # ToDo: Do we need to set these?
        e.object.opts[:p1] = from_vertex.object.opts[:position]
        e.object.opts[:p2] = to_vertex.object.opts[:position]
        e.object.opts[:from_vertex_bbx] = get(from_vertex.object.opts, :bounding_box, (O, O) .+ e.object.opts[:p1])
        e.object.opts[:to_vertex_bbx] = get(to_vertex.object.opts, :bounding_box, (O, O) .+ e.object.opts[:p2])
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ cbdf4e42-2985-11eb-2041-dbd79e434686
using Plots, OTSIN, StatsBase, NamedArrays

# ╔═╡ d94c4784-29bb-11eb-1050-517f3dd80fa0
md"""
Illustration of the optimal transportation on a seed dispersal dataset `M_SD_004`.

```

Spindalis portoricensis
Euphonia musica
Nesospingus sp1 M_SD_004
Loxigilla portoricensis
Vireo altiloquous
Melanerpes portoricensis
Columba squamosa
Margarops fuscatus
Tyrannus dominicensis
Zenaida asiatica
Coereba flaveola
Dendroica tigrina
Vireo flavirostris
Dendroica caerulescens
Icterus dominicensis
Myarchus antillarum
Tiaris bicolor
Todus mexicanus
Tyrannus caudifasciatus
Vireo latimeri
Phoradendron sp1 M_SD_004	10	55	3	9	1	0	0	0	0	0	1	0	1	0	0	0	0	0	0	0
Schefflera morototoni	56	1	4	1	3	0	5	1	4	0	1	1	0	0	0	0	0	0	1	0
Cecropia schreberiana	29	0	6	6	4	3	0	1	0	0	2	3	0	1	0	0	0	0	0	0
Ficus sp1 M_SD_004	34	0	4	0	1	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0
Guarea guidonia	0	0	5	17	14	0	0	1	0	0	0	0	1	0	0	0	0	0	0	0
Cordia sulcata	28	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Alchornea latifolia	3	0	0	14	3	0	0	0	0	5	0	0	0	0	0	0	0	1	0	1
Miconia serrulata	6	2	12	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Clusia rosea	0	0	1	0	11	3	0	0	1	0	0	0	0	0	1	1	0	0	0	0
Anthurium scandens	1	16	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Miconia affinis	2	0	10	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Inga vera	0	0	8	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Musa acuminata	1	0	3	4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Inga laurina	4	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Dendropemon bicolor	0	0	3	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Myrcine coriacea	0	0	2	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Syzigium jambos	0	0	1	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Andira inermis	1	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Coffea arabica	0	0	0	1	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Miconia racemosa	0	0	2	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Palicourea guianensis	0	0	1	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Tilandsia sp1 M_SD_004	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Buchenavia capitata	0	0	0	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Casearia arborea	1	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Citrus sinensis	1	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Dendropanax arboreus	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Henriettea fascicularis	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Momordica charantia	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Palicourea crocea	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Philodendron angustatum	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Piper sp1 M_SD_004	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
Solanum rugosum	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0
Unidentified sp1 M_SD_004	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Zanthoxylum martinicensis	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
```

"""

# ╔═╡ 8dae9448-29bc-11eb-18b0-fb55e9fb1be3
plant_species = ["Phoradendron sp1 M_SD_004",
	"Schefflera morototoni",
	"Cecropia schreberiana",
	"Ficus sp1 M_SD_004",
	"Guarea guidonia",
	"Cordia sulcata",
	"Alchornea latifolia",
	"Miconia serrulata",
	"Clusia rosea",
	"Anthurium scandens",
	"Miconia affinis",
"Inga vera",
"Musa acuminata",
"Inga laurina",
"Dendropemon bicolor",
"Myrcine coriacea",
"Syzigium jambos",
"Andira inermis",
"Coffea arabica",
"Miconia racemosa",
"Palicourea guianensis",
"Tilandsia sp1 M_SD_004",
"Buchenavia capitata",
"Casearia arborea",
"Citrus sinensis",
"Dendropanax arboreus",
"Henriettea fascicularis",
"Momordica charantia",
"Palicourea crocea",
"Philodendron angustatum",
"Piper sp1 M_SD_004",
"Solanum rugosum",
"Unidentified sp1 M_SD_004",
"Zanthoxylum martinicensis"]

# ╔═╡ 67d66742-2b32-11eb-1a58-370b8ca72b53
length(plant_species)

# ╔═╡ 37e622b0-29bc-11eb-2ede-a9430c312543
bird_species = ["Spindalis portoricensis",
	"Euphonia musica",
	"Nesospingus sp1 M_SD_004",
	"Loxigilla portoricensis",
	"Vireo altiloquous",
	"Melanerpes portoricensis",
	"Columba squamosa",
	"Margarops fuscatus",
	"Tyrannus dominicensis",
	"Zenaida asiatica",
	"Coereba flaveola",
	"Dendroica tigrina",
	"Vireo flavirostris",
	"Dendroica caerulescens",
	"Icterus dominicensis",
	"Myarchus antillarum",
	"Tiaris bicolor",
	"Todus mexicanus",
	"Tyrannus caudifasciatus",
	"Vireo latimeri"]

# ╔═╡ 60219008-2b32-11eb-3279-37c1b9489162
length(bird_species)

# ╔═╡ ed1c921c-2986-11eb-28ad-e1f7d3489193
Y = [10 55 3 9 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0;
56 1 4 1 3 0 5 1 4 0 1 1 0 0 0 0 0 0 1 0;
29 0 6 6 4 3 0 1 0 0 2 3 0 1 0 0 0 0 0 0;
34 0 4 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 5 17 14 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0;
28 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
3 0 0 14 3 0 0 0 0 5 0 0 0 0 0 0 0 1 0 1;
6 2 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 11 3 0 0 1 0 0 0 0 0 1 1 0 0 0 0;
1 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
2 0 10 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 8 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 3 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
4 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 2 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

# ╔═╡ 06c42ff4-29be-11eb-23d3-bffa5ba9a554
n_interactions = sum(Y)

# ╔═╡ 08688ab4-298a-11eb-0ee0-0f00e8ee6248
n, m = size(Y)

# ╔═╡ 3ec1637c-2987-11eb-27c3-7d291ca7306e
a, b = marginals(Y)

# ╔═╡ 4859100e-2987-11eb-3f77-b5825e36a555
P = normalize(Y)

# ╔═╡ 398ed8c4-29bb-11eb-3089-afafbc8345af
plot_Y = marginalsheatmap(Y, cmap=cgrad([:white, :blue]), showaxis=true, xticks=(5:5:m), yticks=5:5:n, xlabel="bird species", ylabel="plant species")

# ╔═╡ 50b53a72-2987-11eb-27b5-a53dc0b26e5d
M = fitM(P, fix_a=true, fix_b=true, γ=1e-4)

# ╔═╡ 88a8f162-2987-11eb-2ed8-d5689f6e2895
plot_M = heatmap(M, color=cgrad([:blue,:white,:red]), xlabel="bird species", ylabel="plant species", title="fitted M")

# ╔═╡ 66caee3c-2988-11eb-188a-ef54b7e65c39
Q = optimaltransport(M, a, b)

# ╔═╡ c72c0ccc-29bb-11eb-1819-f361cf9caa6f
Qpu = optimaltransport(M, ones(n)/n, b)

# ╔═╡ 8de8b87c-29be-11eb-3022-83ee7e3f562b
Qbo = optimaltransport(M, a)

# ╔═╡ 368d73e2-2989-11eb-3dc5-2f89432bff64
plot_Qpu = marginalsheatmap(Qpu, cmap=cgrad([:white, :purple]), showaxis=true, xticks=(5:5:m), yticks=5:5:n, xlabel="bird species", ylabel="plant species")

# ╔═╡ 98e014b6-29be-11eb-0713-1de8527dc3d2
plot_Qbo = marginalsheatmap(Qbo, cmap=cgrad([:white, :purple]), showaxis=true, xticks=(5:5:m), yticks=5:5:n, xlabel="bird species", ylabel="plant species")

# ╔═╡ 265d6474-29bf-11eb-2e20-a1180d19e9e8
p_total = plot(plot_Y, plot_M, plot_Qpu, plot_Qbo, size=(1200, 1000))

# ╔═╡ 7597e448-29c3-11eb-2262-951017f11586
savefig(p_total, "../figures/SD_sim.svg")

# ╔═╡ 7762f2c6-2988-11eb-33b1-a3817c1041b7
mean(abs.(Y .- Q))

# ╔═╡ 26ca4e9e-2989-11eb-06ae-fb1ec5ab97a8
kldivergence(P, normalize(Q))

# ╔═╡ Cell order:
# ╟─d94c4784-29bb-11eb-1050-517f3dd80fa0
# ╠═8dae9448-29bc-11eb-18b0-fb55e9fb1be3
# ╠═67d66742-2b32-11eb-1a58-370b8ca72b53
# ╠═37e622b0-29bc-11eb-2ede-a9430c312543
# ╠═60219008-2b32-11eb-3279-37c1b9489162
# ╠═cbdf4e42-2985-11eb-2041-dbd79e434686
# ╠═ed1c921c-2986-11eb-28ad-e1f7d3489193
# ╠═06c42ff4-29be-11eb-23d3-bffa5ba9a554
# ╠═08688ab4-298a-11eb-0ee0-0f00e8ee6248
# ╠═3ec1637c-2987-11eb-27c3-7d291ca7306e
# ╠═4859100e-2987-11eb-3f77-b5825e36a555
# ╠═398ed8c4-29bb-11eb-3089-afafbc8345af
# ╠═50b53a72-2987-11eb-27b5-a53dc0b26e5d
# ╠═88a8f162-2987-11eb-2ed8-d5689f6e2895
# ╠═66caee3c-2988-11eb-188a-ef54b7e65c39
# ╠═c72c0ccc-29bb-11eb-1819-f361cf9caa6f
# ╠═8de8b87c-29be-11eb-3022-83ee7e3f562b
# ╠═368d73e2-2989-11eb-3dc5-2f89432bff64
# ╠═98e014b6-29be-11eb-0713-1de8527dc3d2
# ╠═265d6474-29bf-11eb-2e20-a1180d19e9e8
# ╠═7597e448-29c3-11eb-2262-951017f11586
# ╠═7762f2c6-2988-11eb-33b1-a3817c1041b7
# ╠═26ca4e9e-2989-11eb-06ae-fb1ec5ab97a8

ext: jl

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
# module Solvers
#
# using Parameters
# using MathOptInterface, JuMP
# using GraphUtils
# using DataStructures
#
# using CRCBS
# using ..TaskGraphs

"""
    AbstractPCTAPFSolver

Abstract type of which all PC-TAPF solvers must be concrete subtypes. All
concrete solvers must implement the following interface for solving PC-TAPF
problems:
- `solution, cost = solve!(solver,problem_def)`
- `check_runtime(solver)` should trigger an interrupt + early return if the
    allowable runtime has been exceeded

Also, we need a good Logger type for keeping track of thing like runtime,
iterations, optimality gap (including upper and lower bound), etc.
"""
abstract type AbstractPCTAPFSolver end

abstract type SearchTrait end
struct Prioritized <: SearchTrait end
struct NonPrioritized <: SearchTrait end
function search_trait end

# get_primary_cost(::NonPrioritized,model,cost) = cost[1]
# get_primary_cost(::Prioritized,model,cost) = cost[2]
# get_primary_cost(solver,cost) = get_primary_cost(search_trait(solver),solver,cost)

primary_cost(solver,cost) = cost[1]
primary_cost_type(solver) = Float64

export AStarSC

"""
    AStarSC

Low-level path planner that employs Slack-and-Collision-Aware A*.
Fields:
- logger
- replan : if true, planner will replan with an empty conflict table following
    timeout.
"""
@with_kw struct AStarSC{C} <: AbstractAStarPlanner
    logger::SolverLogger{C} = SolverLogger{C}(iteration_limit=1000000, runtime_limit=1000)
    replan::Bool            = false
end
search_trait(solver::AStarSC) = NonPrioritized()
function CRCBS.logger_step_a_star!(solver::AStarSC, env, base_path, s, q_cost)
    increment_iteration_count!(solver)
    # if solver.DEBUG
    #     push!(solver.search_history, s)
    # end
    if check_iterations(solver)
        # if solver.DEBUG
        #     # Dump env to JLD2 environment
        #     filename = joinpath(DEBUG_PATH,string("A_star_dump_",get_debug_file_id(),".jld2"))
        #     mkpath(DEBUG_PATH)
        #     @log_info(-1,verbosity(solver),"Dumping A* env to $filename")
        #     agent_id = env.agent_idx
        #     history = map(s->(s.vtx,s.t),solver.search_history)
        #     start   = (get_final_state(base_path).vtx,get_final_state(base_path).t)
        #     goal    = (env.goal.vtx,env.goal.t)
        #     idx = isa(solver.AStarSC) ? 2 : 1
        #     paths   = get_cost_model(env).cost_models[idx].model.table.paths
        #     state_constraints = map(c->(c.a,(c.v.s.vtx,c.v.sp.vtx),c.t),collect(env.constraints.state_constraints))
        #     action_constraints = map(c->(c.a,(c.v.s.vtx,c.v.sp.vtx),c.t),collect(env.constraints.action_constraints))
        #     @save filename agent_id history start goal paths state_constraints action_constraints
        #     # @assert false "making a bogus assertion to hack my way out of this block"
        # end
        # throw(SolverAstarMaxOutException(string("# MAX OUT: A* limit of ",solver.LIMIT_A_star_iterations," exceeded.")))
    end
    @log_info(2,verbosity(solver),"A* iter $(iterations(solver)): s = $(string(s)), q_cost = $q_cost")
end

export DefaultAStarSC
DefaultAStarSC() = AStarSC{NTuple{5,Float64}}()

export
    PrioritizedAStarSC

"""
    PrioritizedAStarSC

Low-level proritized path planner that employs Slack-and-Collision-Aware A*.
"""
@with_kw struct PrioritizedAStarSC{C}
    logger::SolverLogger{C} = SolverLogger{C}()
    replan::Bool            = false
end
search_trait(solver::PrioritizedAStarSC) = Prioritized()
primary_cost(::PrioritizedAStarSC,cost::NTuple{5,Float64}) = cost[2]

"""
    construct_heuristic_model(solver,env_graph;kwargs...)

Construct the heuristic model to be used by solver.
"""
function construct_heuristic_model(trait::NonPrioritized,solver,env_graph,
        # ph = PerfectHeuristic(get_dist_matrix(env_graph))
        ph = EnvDistanceHeuristic(),
        )
    construct_composite_heuristic(ph,NullHeuristic(),ph,ph,NullHeuristic())
end
function construct_heuristic_model(trait::Prioritized,args...)
    h = construct_heuristic_model(NonPrioritized(),args...)
    construct_composite_heuristic(
        h.cost_models[2],
        h.cost_models[1],
        h.cost_models[3:end]...
    )
end
function construct_heuristic_model(solver, args...)
    construct_heuristic_model(search_trait(solver),solver,args...)
end

"""
    construct_cost_model(solver::AStarSC, args...;kwargs...)

Defines the cost model used by Slack- and Collision-aware A*.
This particular setting of cost model is crucial for good performance of A_star,
because it encourages depth first search. If we were to replace terms 3-5 with
SumOfTravelTime(), we would get worst-case exponentially slow breadth-first
search!
"""
function construct_cost_model(trait::NonPrioritized,
        solver, sched, cache, problem_spec, env_graph, primary_objective=SumOfMakeSpans();
        extra_T::Int=400)
    N = length(get_robot_ICs(sched))
    @assert N > 0 "num_robots should be > 0. We need at least one robot!"
    cost_model = construct_composite_cost_model(
        typeof(primary_objective)(sched,cache),
        HardConflictCost(env_graph,makespan(sched)+extra_T, N),
        SumOfTravelDistance(),
        FullCostModel(sum,NullCost()), # SumOfTravelTime(),
        FullCostModel(sum,TransformCostModel(c->-1*c,TravelTime()))
    )
    heuristic_model = construct_heuristic_model(trait,solver,env_graph)
    cost_model, heuristic_model
end
function construct_cost_model(trait::SearchTrait,
        solver, env::SearchEnv,
        primary_objective=get_problem_spec(env).cost_function
        ;
        kwargs...)
    construct_cost_model(trait,solver,
        get_schedule(env), get_cache(env), get_problem_spec(env), get_graph(env),
        primary_objective;kwargs...)
end
function construct_cost_model(trait::Prioritized,args...;kwargs...)
    c, h = construct_cost_model(NonPrioritized(),args...;kwargs...)
    # switch the first two elements of the cost and heuristic models
    cost_model = construct_composite_cost_model(
        c.cost_models[2],
        c.cost_models[1],
        c.cost_models[3:end]...
    )
    heuristic = construct_composite_heuristic(
        h.cost_models[2],
        h.cost_models[1],
        h.cost_models[3:end]...
    )
    cost_model, heuristic
end
function construct_cost_model(solver, args...;kwargs...)
    construct_cost_model(search_trait(solver),solver,args...;kwargs...)
end

export
    update_route_plan!

"""
    update_route_plan!()
"""
function update_route_plan!(solver,pc_mapf::AbstractPC_MAPF,env,v,path,cost,node,
    agent_id=get_id(get_default_robot_id(node)))
    set_solution_path!(env, path, agent_id)
    set_path_cost!(env, cost, agent_id)
    update_env!(solver,env,v,path,agent_id)
    set_cost!(env,aggregate_costs(get_cost_model(env),get_path_costs(env)))
    return env
end
function update_route_plan!(solver,pc_mapf::C_PC_MAPF,env,v,meta_path,meta_cost,schedule_node)
    # add to solution
    paths = MetaAgentCBS.split_path(meta_path)
    # @show length(paths)
    # @show length(meta_env.envs)
    # @show length(meta_cost.independent_costs)
    # @show length(schedule_node.instructions)
    for (new_path, sub_node) in zip(paths, sub_nodes(schedule_node))
        agent_id = get_id(get_robot_id(sub_node))
        path = new_path
        update_route_plan!(solver,PC_MAPF(get_env(pc_mapf)),env,v,path,get_cost(path),sub_node,agent_id)
        @log_info(3,verbosity(solver),"agent_id = ", agent_id)
        @log_info(3,verbosity(solver),string("agent_path = ", convert_to_vertex_lists(path)))
        @log_info(3,verbosity(solver),string("cost = ", get_cost(new_path)))
    end
    return env
end

export
    plan_path!

path_finder(solver::AStarSC,args...) = CRCBS.a_star!(solver,args...)

"""
    plan_path!

    Computes next path specified by the project schedule and updates the
    solution in the ConstraintTreeNode::node accordingly.
"""
function plan_path!(solver::AStarSC, pc_mapf::AbstractPC_MAPF, env::SearchEnv, node::N,
        schedule_node::T, v::Int;
        kwargs...) where {N<:ConstraintTreeNode,T}

    cache = get_cache(env)
    node_id = get_vtx_id(get_schedule(env),v)
    # println("hi")
    reset_solver!(solver)
    cbs_env = build_env(solver, pc_mapf, env, node, VtxID(v)) #schedule_node, v)
    base_path = get_base_path(solver,env,cbs_env)
    ### PATH PLANNING ###
    # solver.DEBUG ? validate(base_path,v) : nothing
    path, cost = path_finder(solver, cbs_env, base_path)
    @assert get_cost(path) == cost
    if cost == get_infeasible_cost(cbs_env)
        if solver.replan == true
            @log_info(-1,verbosity(solver),"A*: replanning without conflict cost", string(schedule_node))
            reset_solver!(solver)
            cost_model, _ = construct_cost_model(solver, env;
                primary_objective=get_problem_spec(env).cost_function)
            cbs_env = build_env(solver, pc_mapf, env, node, VtxID(v);cost_model=cost_model)
            base_path = get_base_path(solver,env,cbs_env)
            path, cost = path_finder(solver, cbs_env, base_path)
        end
        if cost == get_infeasible_cost(cbs_env)
            # @log_info(-1,verbosity(solver),"A*: r1eturned infeasible path fosssr 1 node ", string(schedule_node))
            return false
        end
    end
    # solver.DEBUG ? validate(path,v,cbs_env) : nothing
    #####################
    @log_info(2,verbosity(solver),string("A* iterations = ",iterations(solver)))
    # Make sure every robot sticks around for the entire time horizon
    if is_terminal_node(get_graph(get_schedule(env)),v)
        @log_info(2,verbosity(solver),"ISPS: length(path) = ",length(path),
        ". Extending terminal node ", string(schedule_node),
        " to ",makespan(get_schedule(env)))
        # extend_path!(cbs_env,path,maximum(cache.tF))
        extend_path!(cbs_env,path,makespan(get_schedule(env)))
        @log_info(2,verbosity(solver),"ISPS: length(path) = ",length(path),
        ". maximum(cache.tF) = ",makespan(get_schedule(env)))
        # solver.DEBUG ? validate(path,v) : nothing
    end
    # add to solution
    update_route_plan!(solver,pc_mapf,env,v,path,cost,schedule_node)
    @log_info(2,verbosity(solver),"ISPS: after update_route_plan! maximum(cache.tF) = ",
        makespan(get_schedule(env)))
    if get_cost(node) >= best_cost(solver)
        @log_info(0,verbosity(solver),"ISPS: get_cost(node) (",get_cost(node),
            ") >= best_cost(solver)",best_cost(solver)," ... Exiting early")
        return false
    end

    return true
end

export
    ISPS

"""
    ISPS

Path planner that employs Incremental Slack-Prioritized Search.
"""
@with_kw struct ISPS{L,C} <: BiLevelPlanner
    low_level_planner::L = DefaultAStarSC()
    logger::SolverLogger{C} = SolverLogger{cost_type(low_level_planner)}(
        iteration_limit = 2
    )
    tighten_paths::Bool = true
end
ISPS(planner) = ISPS(low_level_planner=planner)
construct_cost_model(solver::ISPS,args...;kwargs...) = construct_cost_model(low_level(solver),args...;kwargs...)
search_trait(solver::ISPS) = search_trait(low_level(solver))
primary_cost(solver::ISPS,cost) = primary_cost(low_level(solver),cost)
primary_cost_type(solver::ISPS) = primary_cost_type(low_level(solver))
function CRCBS.set_best_cost!(solver::ISPS,cost)
    set_best_cost!(get_logger(solver),cost)
    set_best_cost!(low_level(solver),cost)
end
function CRCBS.hard_reset_solver!(solver::ISPS)
    hard_reset_solver!(get_logger(solver))
    hard_reset_solver!(low_level(solver))
end
for op in [:set_deadline!,:set_runtime_limit!,:set_verbosity!]
    eval(quote
        CRCBS.$op(solver::ISPS,args...) = begin
            $op(get_logger(solver),args...)
            $op(low_level(solver),args...)
        end
    end)
end

export
    plan_next_path!,
    compute_route_plan!

"""
    `plan_next_path`
"""
function plan_next_path!(solver::ISPS, pc_mapf::AbstractPC_MAPF, env::SearchEnv, node::N
        ) where {N<:ConstraintTreeNode}

    valid_flag = true
    if ~isempty(get_cache(env).node_queue)
        v,priority = dequeue_pair!(get_cache(env).node_queue)
        node_id = get_vtx_id(get_schedule(env),v)
        schedule_node = get_node_from_id(get_schedule(env),node_id)
        if get_path_spec(get_schedule(env), v).plan_path == true
            try
                @log_info(2,verbosity(solver),sprint(show,env))
                @log_info(2,verbosity(solver),
                """
                ISPS:
                    schedule_node: $(string(schedule_node))
                    get_tF(sched,v): $(get_tF(env,v))
                    maximum(get_cache(env).tF): $(makespan(get_schedule(env)))
                """)
                valid_flag = plan_path!(low_level(solver),pc_mapf,env,node,schedule_node,v)
                @log_info(2,verbosity(solver),"""
                ISPS:
                    routes:
                """,
                sprint_indexed_list_array(
                    convert_to_vertex_lists(get_route_plan(env));leftaligned=true),
                """
                    cache.tF: $(get_tF(env,v))
                    get_tF(sched,v): $(get_tF(env,v))
                """)
            catch e
                if isa(e, SolverException)
                    handle_solver_exception(solver,e)
                    @log_info(-1,verbosity(solver),"Exiting ISPS with failed status")
                    valid_flag = false
                    # return valid_flag
                    return false
                else
                    rethrow(e)
                end
            end
        else
            update_planning_cache!(solver,env,v,-1)
        end
    end
    return valid_flag
end

"""
    compute_route_plan!

Computes all paths specified by the project schedule and updates the solution
in the ConstraintTreeNode::node accordingly.
"""
function compute_route_plan!(solver::ISPS, pc_mapf::AbstractPC_MAPF, node::N, env::SearchEnv=node.solution;
        kwargs...
        ) where {N<:ConstraintTreeNode}

    while length(get_cache(env).node_queue) > 0
        status = plan_next_path!(solver,pc_mapf,env,node)
        # @log_info(-1,verbosity(solver),"ISPS status = ",status)
        # if status == false
        #     @log_info(-1,verbosity(solver),"Exiting ISPS with failed status")
        #     return false
        # end
        status ? nothing : return false
        status = tighten_gaps!(solver,pc_mapf,env,node)
        status ? nothing : return false
        enforce_time_limit!(solver)
    end
    return true
end

function CRCBS.low_level_search!(solver::ISPS, pc_mapf::AbstractPC_MAPF,
        node::N=initialize_root_node(solver,pc_mapf),
        idxs::Vector{Int}=Vector{Int}()
    ) where {N<:ConstraintTreeNode}

    reset_solver!(solver)
    valid_flag = true
    search_env = node.solution
    for i in 1:iteration_limit(solver)
        increment_iteration_count!(solver)
        # reset solution
        reset_schedule_times!(get_schedule(search_env),get_schedule(get_env(pc_mapf)))
        reset_cache!(
            get_cache(search_env),
            get_schedule(search_env),
            # get_cache(get_env(pc_mapf)).t0,
            # get_cache(get_env(pc_mapf)).tF,
            # map(v->get_t0(search_env,v),vertices(get_schedule(get_env(pc_mapf)))),
            # map(v->get_tF(search_env,v),vertices(get_schedule(get_env(pc_mapf)))),
            )
        reset_route_plan!(node,get_route_plan(get_env(pc_mapf)))
        valid_flag = compute_route_plan!(solver, pc_mapf, node)
        if valid_flag == false
            @log_info(0,verbosity(solver),"ISPS: failed on ",i,"th repair iteration.")
            break
        end
    end
    return valid_flag
end
function CRCBS.low_level_search!(solver::CBSSolver, pc_mapf::AbstractPC_MAPF,node::ConstraintTreeNode,args...;kwargs...)
    low_level_search!(low_level(solver),pc_mapf,node,args...)
end

"""
    solve!(solver::ISPS,mapf::M) where {M<:PC_MAPF}

Solve a PC-MAPF problem with the ISPS algorithm.
"""

function CRCBS.solve!(solver::ISPS,pc_mapf::PC_MAPF)
    node = initialize_root_node(solver,pc_mapf)
    valid_flag = compute_route_plan!(solver,pc_mapf,node)#,node)
    return node.solution, get_cost(node.solution)
end

construct_cost_model(solver::CBSSolver,args...;kwargs...) = construct_cost_model(low_level(solver),args...;kwargs...)
search_trait(solver::CBSSolver) = search_trait(low_level(solver))
primary_cost(solver::CBSSolver,cost) = primary_cost(low_level(solver),cost)
primary_cost_type(solver::CBSSolver) = primary_cost_type(low_level(solver))
CRCBS.solve!(solver::CBSSolver,pc_mapf::PC_MAPF) = CRCBS.cbs!(solver,pc_mapf)

export
    TaskGraphsMILPSolver

"""
    TaskGraphsMILPSolver

Wrapper for MILP solver for assignment problem.
"""
@with_kw struct TaskGraphsMILPSolver{M,C}
    milp::M = SparseAdjacencyMILP()
    logger::SolverLogger{C} = SolverLogger{Int}()
end
TaskGraphsMILPSolver(milp) = TaskGraphsMILPSolver(milp,SolverLogger{Int}())
get_assignment_matrix(solver::TaskGraphsMILPSolver) = get_assignment_matrix(solver.milp)

function formulate_milp(solver::TaskGraphsMILPSolver,args...;kwargs...)
    milp = formulate_milp(solver.milp,args...;
        kwargs...)
    set_optimizer_attributes(milp,default_optimizer_attributes()...)
    set_time_limit_sec(milp, max(0,time_to_deadline(solver)))
    milp
end

export
    NBSSolver,
    assignment_solver,
    route_planner

"""
    NBSSolver{A,P}

A hierarchical PC-TAPF solver with an assignment level and a path-planning
level.
The solver works by alternating between assignment and path-planning until the
optimality gap between the lower bound (from task assignment) and the lower
bound (from path planning) disappears.
The input to the assignment problem is the full PC-TAPF problem specification.
The output of the assignment problem is a valid `OperatingSchedule`--that is,
an operating schedule wherein all assignments have been made in a legal way.
The input to the route planner is the PC-TAPF problem spec along with the
`OperatingSchedule` that comes from the assignment solution.
"""
@with_kw struct NBSSolver{A,P,C} <: AbstractPCTAPFSolver
    assignment_model::A     = TaskGraphsMILPSolver()
    path_planner    ::P     = CBSSolver(ISPS())
    logger          ::SolverLogger{C} = SolverLogger{primary_cost_type(path_planner)}()
end
NBSSolver(a,b) = NBSSolver(assignment_model=a,path_planner=b)
assignment_solver(solver::NBSSolver) = solver.assignment_model
route_planner(solver::NBSSolver) = solver.path_planner
construct_cost_model(solver::NBSSolver,args...;kwargs...) = construct_cost_model(route_planner(solver),args...;kwargs...)
search_trait(solver::NBSSolver) = search_trait(route_planner(solver))
primary_cost(solver::NBSSolver,cost) = primary_cost(route_planner(solver),cost)
primary_cost_type(solver::NBSSolver) = primary_cost_type(route_planner(solver))
function CRCBS.set_best_cost!(solver::NBSSolver,cost)
    set_best_cost!(get_logger(solver),cost)
    set_best_cost!(route_planner(solver),cost)
end
function CRCBS.hard_reset_solver!(solver::NBSSolver)
    hard_reset_solver!(get_logger(solver.logger))
    hard_reset_solver!(route_planner(solver))
end
function CRCBS.reset_solver!(solver::NBSSolver)
    reset_solver!(get_logger(solver))
    reset_solver!(route_planner(solver))
end
for op in [:set_deadline!,:set_runtime_limit!,:set_verbosity!]
    eval(quote
        CRCBS.$op(solver::NBSSolver,args...) = begin
            $op(get_logger(solver),args...)
            $op(assignment_solver(solver),args...)
            $op(route_planner(solver),args...)
        end
    end)
end

"""
    solve!(solver, base_env::SearchEnv;kwargs...) where {A,P}

Use the planner defined by `solver` to solve the PC-TAPF problem encoded by
`base_env`. For solvers of type `NBSSolver`, the algorithm involves
repeatedly solving an assignment problem followed by a route-planning problem.
Within the generic `solve!` method it is possible to initialize an assignment
problem (the type is not constrained) and then modify it via
`update_assignment_problem!` prior to each new call to
`solve_assignment_problem!`. This is the approach taken for various MILP-based
assignment solvers. It is also possible to reconstruct the assignment problem
from scratch within each call to `solve_assignment_problem!`.

Arguments:
- solver <: AbstractPCTAPFSolver
- base_env::SearchEnv : a PC-TAPF problem

Outputs:
- best_env : a `SearchEnv` data structure that encodes a solution to the problem
- cost : the cost of the solution encoded by `best_env`
"""
function CRCBS.solve!(solver::NBSSolver, prob::E;kwargs...) where {E<:AbstractPC_TAPF}
    best_env = SearchEnv(get_env(prob),route_plan=initialize_route_plan(get_env(prob)))
    set_cost!(best_env,get_infeasible_cost(best_env))
    set_best_cost!(solver, primary_cost(solver, get_cost(best_env)))
    assignment_problem = formulate_assignment_problem(assignment_solver(solver),
        prob;kwargs...)
    try
        while optimality_gap(solver) > 0
            enforce_time_limit!(solver)
            update_assignment_problem!(assignment_solver(solver),
                assignment_problem, prob)
            schedule, l_bound = solve_assignment_problem!(
                assignment_solver(solver),
                assignment_problem,
                prob)
            set_lower_bound!(solver,l_bound)
            if optimality_gap(solver) > 0
                env, cost = plan_route!(route_planner(solver),
                    schedule, prob;kwargs...)
                if cost < best_cost(solver)
                    best_env = env
                    set_best_cost!(solver,cost)
                end
            end
            increment_iteration_count!(solver)
            if check_iterations(solver)
                @log_info(1,verbosity(solver),
                    "NBS: Reached $(iteration_limit(solver))-iteration limit.")
                break
            end
        end
    catch e
        isa(e, SolverException) ? handle_solver_exception(solver,e) : rethrow(e)
    end
    return best_env, best_cost(solver)
end
function CRCBS.solve!(solver, env_graph, schedule::OperatingSchedule, problem_spec::ProblemSpec;
    kwargs...)
    base_env = construct_search_env(solver,schedule,problem_spec,env_graph;kwargs...)
    solve!(solver, base_env;kwargs...)
end
function CRCBS.solve!(solver, env_graph, project_spec::ProjectSpec, problem_spec::ProblemSpec,
        robot_ICs;kwargs...)
    schedule = construct_partial_project_schedule(
        project_spec,problem_spec,map(i->robot_ICs[i], 1:problem_spec.N))
    solve!(solver, env_graph, schedule, problem_spec;kwargs...)
end
CRCBS.solve!(solver::NBSSolver, env::SearchEnv;kwargs...) = solve!(solver,PC_TAPF(env);kwargs...)

export
    formulate_assignment_problem

"""
    formulate_assignment_problem(solver,prob;

Returns an assignment problem instance that can be updated (as opposed to being
reconstructed from scratch) on each call to `update_assignment_problem!` prior
to being resolved.
"""
function formulate_assignment_problem(solver,prob;
        kwargs...)
    formulate_milp(solver,get_env(prob);
        cost_model=get_problem_spec(get_env(prob)).cost_function,
        kwargs...)
end

function formulate_milp(milp_model::TaskGraphsMILP,env::SearchEnv;kwargs...)
    t0,tF = get_node_start_and_end_times(env)
    formulate_milp(milp_model,get_schedule(env),get_problem_spec(env);
        t0_=t0,
        tF_=tF,
        kwargs...
        )
end
function formulate_milp(model::AssignmentMILP,env::SearchEnv;kwargs...)
    formulate_milp(model,get_schedule(env),get_problem_spec(env);kwargs...)
end

export
    update_assignment_problem!

"""
    update_assignment_problem!(solver, assignment_problem)

A helper method for updating an instance of an assignment problem. In the case
    of MILP-based models, this method simply excludes all previous solutions by
    adding new constraints on the assignment/adjacency matrix.
"""
function update_assignment_problem!(solver, model::T, base_prob) where {T<:TaskGraphsMILP}
    exclude_solutions!(model) # exclude most recent solution in order to get next best solution
    # Trying this
    set_time_limit_sec(model, max(0,time_to_deadline(solver)))
    model
end

export
    solve_assignment_problem!

"""
    solve_assignment_problem!(solver,model,prob)

Solve the "assignment problem"--i.e., the relaxation of the full PC-TAPF problem
wherein we ignore collisions--using the algorithm encoded by solver.
"""
function solve_assignment_problem!(solver::TaskGraphsMILPSolver, model, prob)
    l_bound = lower_bound(solver)
    set_time_limit_sec(model, max(0,time_to_deadline(solver)))
    optimize!(model)
    if primal_status(model) != MOI.FEASIBLE_POINT
        throw(SolverException("Assignment problem is infeasible -- in `solve_assignment_problem!()`"))
    end
    set_lower_bound!(solver, Int(round(value(objective_bound(model)))) )
    set_best_cost!(solver, Int(round(value(objective_function(model)))) )
    if termination_status(model) == MOI.OPTIMAL
        @assert lower_bound(solver) <= best_cost(solver) "lower_bound($(solver_type(solver))) = $(value(objective_bound(model))) -> $(lower_bound(solver)) but should be equal to best_cost($(solver_type(solver))) = $(value(objective_function(model))) -> $(best_cost(solver))"
    end
    sched = deepcopy(get_schedule(get_env(prob)))
    update_project_schedule!(solver, model, sched, get_problem_spec(get_env(prob)))
    sched, lower_bound(solver)
end

function CRCBS.solve!(solver,pcta::PC_TA)
    prob = formulate_assignment_problem(solver,pcta)
    sched, l_bound = solve_assignment_problem!(solver,prob,pcta)
end

export
    plan_route!

"""
    plan_route!

Compute a route plan that corresponds to the OperatingSchedule.
Arguments:
- solver
- schedule
- search_env

Outputs:
- A `SearchEnv` the contains a valid solution
"""
function plan_route!(
        solver,
        schedule,
        prob;
        kwargs...)

    env = construct_search_env(solver, schedule, get_env(prob);kwargs...)
    pc_mapf = construct_routing_problem(prob,env)
    solution, cost = solve!(solver, pc_mapf)
    solution, primary_cost(solver,get_cost(solution))
end

CRCBS.solver_type(::NBSSolver)          = "NBSSolver"
CRCBS.solver_type(::CBSSolver)          = "CBSSolver"
CRCBS.solver_type(::PIBTPlanner)        = "PIBTPlanner"
CRCBS.solver_type(::ISPS)               = "ISPS"
CRCBS.solver_type(::AStarSC)            = "AStarSC"
CRCBS.solver_type(::PrioritizedAStarSC) = "PrioritizedAStarSC"
CRCBS.solver_type(::AStar)              = "AStar"
CRCBS.solver_type(::TaskGraphsMILPSolver{M,C}) where {M,C} = "TaskGraphsMILPSolver{$M}"

# end


text:
"""
    TextConfig

    TextConfig(; kw...)
    TextConfig(face, namepixels, namepos, timepixels, timepos, fcolor, bcolor)

Text configuration for printing timestep and grid name on the image.

# Arguments / Keywords

- `font`: A `FreeTypeAbstraction.FTFont`, or a `String` with the font name to look for.
    The `FTFont` may load more quickly. See also [`set_default_font`](@ref).
- `namepixels` and `timepixels`: the pixel size of the font.
- `timepos` and `namepos`: tuples that set the label positions, in `Int` pixels.
- `fcolor` and `bcolor`: the foreground and background colors, as `ARGB32`.
"""
struct TextConfig{F,NPi,NPo,TPi,TPo,FC,BC}
    face::F
    namepixels::NPi
    namepos::NPo
    timepixels::TPi
    timepos::TPo
    fcolor::FC
    bcolor::BC
end
function TextConfig(;
    font=autofont(), namepixels=12, timepixels=12,
    namepos=(3timepixels + namepixels, timepixels),
    timepos=(2timepixels, timepixels),
    fcolor=ARGB32(1.0), bcolor=ZEROCOL,
)
    if font isa FreeTypeAbstraction.FTFont
	    face = font
    elseif font isa AbstractString
        face = FreeTypeAbstraction.findfont(font)
        face isa Nothing && _fontnotfounderror(font)
    else
        _fontwrongtype(font)
    end
    TextConfig(face, namepixels, namepos, timepixels, timepos, fcolor, bcolor)
end

# isbits(FreeTypeAbstraction.FTFont) == false,
# hence isassigned can tell whether the cache has been initialized
const _default_font_ref = Ref{FreeTypeAbstraction.FTFont}()

"""
Set default font.

See also [`TextConfig`](@ref)

# Examples
Using `findfont` to get the font by name:
```julia
using DynamicGrids
using FreeTypeAbstraction: findfont

font = findfont("Times")
DynamicGrids.set_default_font(font)
```
Or giving the font path directly:
```julia
using DynamicGrids
using FreeTypeAbstraction: FTFont

font = FTFont("/usr/share/fonts/truetype/Adobe-Times-Regular.otb")
DynamicGrids.set_default_font(font)
```
"""
function set_default_font(font)
	_default_font_ref[] = font
end

function autofont()
    if isassigned(_default_font_ref)
        return _default_font_ref[]
    else
        names = if Sys.islinux()
            ("cantarell", "sans-serif", "Bookman")
        else
            ("arial", "sans-serif")
        end
        for name in names
            face = FreeTypeAbstraction.findfont(name)
            if face isa FreeTypeAbstraction.FTFont
                set_default_font(face)
                return face
            end
        end
        _nodefaultfonterror(names)
    end
end

# Render time `name` and `t` as text onto the image, following config settings.
function _rendertime! end

function _rendertext!(img, config::TextConfig, name, t)
    _rendername!(img, config::TextConfig, name)
    _rendertime!(img, config::TextConfig, t)
    img
end
_rendertext!(img, config::Nothing, name, t) = nothing

# Render `name` as text on the image following config settings.
function _rendername!(img, config::TextConfig, name)
    renderstring!(img, name, config.face, config.namepixels, config.namepos...;
        fcolor=config.fcolor, bcolor=config.bcolor
    )
    img
end
_rendername!(img, config::TextConfig, name::Nothing) = img
_rendername!(img, config::Nothing, name) = img
_rendername!(img, config::Nothing, name::Nothing) = img

# Render time `t` as text on the image following config settings.
function _rendertime!(img, config::TextConfig, t)
    renderstring!(img, string(t), config.face, config.timepixels, config.timepos...;
        fcolor=config.fcolor, bcolor=config.bcolor
    )
    img
end
_rendertime!(img, config::Nothing, t) = img
_rendertime!(img, config::TextConfig, t::Nothing) = img
_rendertime!(img, config::Nothing, t::Nothing) = img

@noinline _fontwrongtype(font) = throw(ArgumentError(
    """
    font must be either a String or a FreeTypeAbstraction.FTFont,
    got `$font` which is a `$(typeof(font))`.
    """
))

@noinline _fontnotfounderror(font) = throw(ArgumentError(
    """
    Font "$font" wasn't found in this system.
    Specify an existing font with the `font` keyword,
    or call DynamicGrids.set_default_font first,
    or use `text=nothing` to display no text.
    """
))

@noinline _nodefaultfonterror(font) = error(
    """
    Your system does not contain any of the default fonts
    $font.
    Use DynamicGrids.set_default_font first.
    """
)


text:
export Point, Point1, Point2, Point3, Point1f, Point2f, Point3f, Point1B,
       Point2B, Point3B
export coordinates, distance, distance², isCCW, midpoint, nan, ϵ_Point

"""
    Point{Dim, T}

A point in `Dim`-dimensional space with coordinates of type `T`.

The coordinates of the point are given with respect to the canonical
Euclidean basis.

## Examples

```julia
# 2D points
A = Point(0.0, 1.0) # double precision as expected
B = Point(0f0, 1f0) # single precision as expected
C = Point2(0, 1)    # explicitly ask for double precision
D = Point2f(0, 1)   # explicitly ask for single precision
E = Point2B(0, 1)   # explicitly ask for BigFloat

# 3D points
F = Point(1.0, 2.0, 3.0) # double precision as expected
G = Point(1f0, 2f0, 3f0) # single precision as expected
H = Point3(1, 2, 3)      # explicitly ask for double precision
I = Point3f(1, 2, 3)     # explicitly ask for single precision
J = Point3B(1, 2, 3)     # explicitly ask for BigFloat
```

### Notes

- Type aliases are `Point1`, `Point2`, `Point3`, `Point1f`, `Point2f`, `Point3f`,
    `Point1B`, `Point2B`, `Point3B`.
"""
struct Point{Dim,T} <:AbstractVector{T}
    coords::SVector{Dim,T}
    Point{Dim,T}(coords::SVector{Dim,T}) where {Dim,T} = new{Dim,T}(coords)
end

# constructors
Point{Dim,T}(coords...) where {Dim,T} = Point{Dim,T}(SVector{Dim,T}(coords...))
Point(coords::SVector{Dim,T}) where {Dim,T} = Point{Dim,T}(coords)
Point(coords::AbstractVector{T}) where {T} = Point{length(coords),T}(coords)
Point(coords...) = Point(SVector(coords...))

# conversions
Base.convert(::Type{Point{2,T}}, P::Point) where {T} = Point{2,T}(P[1], P[2])
SVector(P::Point{Dim,T}) where {Dim,T} = P.coords

# type aliases
const Point1  = Point{1,Float64}
const Point2  = Point{2,Float64}
const Point3  = Point{3,Float64}
const Point1f = Point{1,Float32}
const Point2f = Point{2,Float32}
const Point3f = Point{3,Float32}
const Point1B = Point{1,BigFloat}
const Point2B = Point{2,BigFloat}
const Point3B = Point{3,BigFloat}

# abstract array interface
Base.size(P::Point) = Base.size(P.coords)
Base.getindex(P::Point, i::Int) = Base.getindex(P.coords, i)
Base.IndexStyle(P::Point) = Base.IndexStyle(typeof(P.coords))

Base.zero(::Type{Point{Dim,T}}) where {Dim,T} = Point{Dim,T}(@SVector zeros(T, Dim))
nan(::Type{Point{Dim,T}}) where {Dim,T} = Point{Dim,T}(@SVector fill(T(NaN), Dim))

const ϵ_Point = 1e-5 # Points separated by 1e-5 cm = 0.1 micron are treated the same.

"""
    coordinates(p::Point)

Return the coordinates of the point with respect to the
canonical Euclidean basis.
"""
coordinates(p::Point) = p.coords

"""
    -(A::Point, B::Point)

Return the [`Vec`](@ref) displacement from point `B` to point `A`.
"""
Base.:-(A::Point, B::Point) = coordinates(A) - coordinates(B)

"""
    +(A::Point, v::Vec)

Return the point at the end of the vector `v` placed
at a reference (or start) point `A`.
"""
Base.:+(A::Point, v::Vec) = Point(coordinates(A) + v)

"""
    -(A::Point, v::Vec)

Return the point at the end of the vector `-v` placed
at a reference (or start) point `A`.
"""
Base.:-(A::Point, v::Vec) = Point(coordinates(A) - v)

"""
    isapprox(A::Point, B::Point)

Return if two points are approximately equal (less than ϵ_Point apart).
"""
Base.isapprox(A::Point{Dim,T}, B::Point{Dim,T}) where {Dim,T} = norm²(A-B) < T(ϵ_Point^2)

"""
    distance(A::Point, B::Point)

Return the distance from point `A` to point `B`.
"""
distance(A::Point, B::Point) = norm(B - A)

"""
    distance²(A::Point, B::Point)

Return the distance from point `A` to point `B` squared.
"""
distance²(A::Point, B::Point) = norm²(B - A)

"""
    midpoint(A::Point, B::Point)

Return the midpoint of the line segment from point `A` to point `B`.
"""
midpoint(A::Point, B::Point) = Point((coordinates(A) + coordinates(B))/2)

"""
    isCCW(A::Point{2}, B::Point{2}, C::Point{2})

If the triplet of 2-dimensional points is counter-clockwise oriented from A to B to C.
"""
isCCW(A::Point{2}, B::Point{2}, C::Point{2}) = 0 ≤ (B - A) × (C - A)

function Base.show(io::IO, point::Point)
    print(io, point.coords.data)
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#export simulate_and_get_coverage_hist_2,simulate_oceansislands,plot_oceansislands_simulation,no_islands_vs_no_reads_simulation,no_oceans_vs_no_reads_simulation,avg_length_oceans_vs_no_reads_simulation,avg_length_islands_vs_no_reads_simulation,simulate_and_get_coverage_hist_3,no_reads_vs_no_islands

using OffsetArrays
using ProgressMeter
using PyPlot
using LsqFit
#include("HistogramReadCoverages.jl")

Base.@kwdef struct Para
    genome_length::Int = 10_000_000
    read_length::Int = 100
    read_number::Int
end

function simulate_and_get_coverage_hist_2(p::Para)

    # simulation
    g = fill(zero(Int32), p.genome_length)
    for i in 1:p.read_number
        pos = rand(1:(p.genome_length - p.read_length))
        g[pos] += 1
        g[pos + p.read_length] -= 1
    end

#     @show g

    max_coverage_for_hist = 1
    hists = OffsetArray(map(x->Int[], 1:max_coverage_for_hist+1), 0:max_coverage_for_hist)

    # compute coverages
    g0 = 0 + g[1]
    p0 = 1
    for pos in 2:p.genome_length
        if (g0==0) && (g[pos] >0)
            increment!(hists[0], pos - p0)
            g0 += g[pos]
            p0 = pos
        elseif (g0 >0) && (g[pos] == - g0)
            increment!(hists[1], pos - p0)
            g0 += g[pos]
            p0 = pos
        elseif g[pos]!=0
            g0+=g[pos]
        end
    end
    hists
end

function simulate_oceansislands(n,genome_length,read_length,no_reads)
    paras = fill(Para(genome_length = genome_length,read_length=read_length,read_number=no_reads), n);
    res = @showprogress map(simulate_and_get_coverage_hist_2, paras);
    h = reduce(vec_hist_add!, res) ./ n
    return h
end

function no_islands_vs_no_reads_simulation(genome_length,read_length,number_reads,width,n)
    counter=[]
    for i in 1:width:number_reads
        n = 30
        paras = fill(Para(genome_length,read_length, i), n);
        h = mapreduce(simulate_and_get_coverage_hist_2, vec_hist_add!, paras)
        h= h./ n
        push!(counter,sum(h[1]))
    end
    return counter
end

function no_oceans_vs_no_reads_simulation(genome_length,read_length,number_reads,width,n)
    counter=[]
    for i in 1:width:number_reads
        n = 30
        paras = fill(Para(genome_length,read_length, i), n);
        h = mapreduce(simulate_and_get_coverage_hist_2, vec_hist_add!, paras)
        h= h./ n
        push!(counter,sum(h[0]))
    end
    return counter
end

function avg_length_oceans_vs_no_reads_simulation(genome_length,read_length,number_reads,width,n)
    counter=[]
    for i in 1:width:number_reads
        n = 30
        paras = fill(Para(genome_length,read_length, i), n);
        h = mapreduce(simulate_and_get_coverage_hist_2, vec_hist_add!, paras)
        h= h./ n
        val=sum(h[0].*[i for i in 1:length(h[0])])/sum(h[0])
        push!(counter,val)
    end
    return counter
end

function avg_length_islands_vs_no_reads_simulation(genome_length,read_length,number_reads,width,n)
    counter=[]
    for i in 1:width:number_reads
        n = 30
        paras = fill(Para(genome_length,read_length, i), n);
        h = mapreduce(simulate_and_get_coverage_hist_2, vec_hist_add!, paras)
        h= h./ n
        val=sum(h[1].*[i for i in 1:length(h[1])])/sum(h[1])
        push!(counter,val)
    end
    return counter
end

function simulate_and_get_coverage_hist_3(p::Para)
    # simulation
    g = fill(zero(Int32), p.genome_length)
    g1 = fill(zero(Int32), p.genome_length)
    for i in 1:p.read_number
        pos = rand(1:(p.genome_length - p.read_length))
        g[pos] += 1
        g[pos + p.read_length] -= 1
        g1[pos] += 1
    end
    hists = zeros(1:p.read_number)
    reads=0
    g0 = 0
    for pos in 1:p.genome_length
        if g[pos] != 0 || g1[pos]>0
            g0 += g[pos]
            reads+=g1[pos]
            if g0==0
                hists[reads]+=1
                reads=0
            end
        end
    end
    hists
end

function no_reads_vs_no_islands(n,genome_length,read_length,read_number)
    p=Para(genome_length=genome_length,read_length=read_length,read_number=read_number)
    a=simulate_and_get_coverage_hist_3(p);
    for j in 1:n-1
        a=a+simulate_and_get_coverage_hist_3(p);
    end
    a=a./n;
    return a
end


text:
"""
    freqdim(x)
    freqdim!(x, val)

Which spatial dimension (1, 2, or 3) corresponds to phase acquisition.
"""
@defprop FrequencyDimension{:freqdim}::Int

"""
    phasedim(x)
    phasedim!(x, val)

Which spatial dimension (1, 2, or 3) corresponds to phase acquisition.
"""
@defprop PhaseDimension{:phasedim}::Int

"""
    slice_start(x)
    slice_start!(x, val)

Which slice corresponds to the first slice acquired during MRI acquisition (i.e. not padded slices).
"""
@defprop SliceStart{:slice_start}::Int

"""
    slice_end(x)
    slice_end!(x, val)

Which slice corresponds to the last slice acquired during MRI acquisition (i.e. not padded slices).
"""
@defprop SliceEnd{:slice_end}::Int

"""
    slicedim(x)
    slicedim!(x, val)

Which dimension slices where acquired at throughout MRI acquisition.
"""
@defprop SliceDimension{:slicedim}::Int

"""
    slice_duration(x)
    slice_duration!(x, val)

Time to acquire one slice
"""
@defprop SliceDuration{:slice_duration}::Float64

"""
    phase_encoding_direction(x) -> EncodingDirection

The phase encoding direction is defined as the direction along which phase was
modulated which may result in visible distortions. Note that this is not the
same as the DICOM term "in_plane_phase_encoding_direction" which can have "ROW" or
"COL" values. This parameter is REQUIRED if corresponding fieldmap data is present
or when using multiple runs with different phase encoding directions (which can
be later used for field inhomogeneity correction).
"""
@defprop PhaseEncodingDirection{:phase_encoding_direction}::EncodingDirection begin
    @getproperty self -> EncodingDirection(phasedim(self))
    @setproperty! (self, val) -> phasedim!(self, val)
end

"""
    slice_encoding_direction(x) -> EncodingDirection

Values ending in "*neg" indicate that the contents of slice_timing are defined
in reverse order (the first entry corresponds to the slice with the largest index,
and the final entry corresponds to slice index zero. When present, the axis
defined by `slice_encoding` needs to be consistent with the `slicedim` field in the
NIfTI header. When absent, the entries in slice_timing must be in the order of increasing
slice index as defined by the NIfTI header.
"""
@defprop SliceEncodingDirection{:slice_encoding_direction}::EncodingDirection begin
    @getproperty self -> EncodingDirection(slicedim(self))
    @setproperty! (self, val) -> slicedim!(self, val)
end

"""
    EncodingDirectionMetadata

Metadata structure for general MRI sequence information.

## Supported Properties
$(description_list(freqdim, phasedim, slicedim, slice_start, slice_end, slice_duration,slice_encoding_direction,phase_encoding_direction))

## Examples

```jldoctest
julia> using NeuroCore

julia> using NeuroCore.NeuroMetadata

julia> m = EncodingDirectionMetadata(1, 2, 3, 4, 5, 6)
EncodingDirectionMetadata(1, 2, 3, 4, 5, 6.0)

julia> m.slice_encoding_direction
kpos::EncodingDirection = 3

julia> m.slice_encoding_direction == slice_encoding_direction(m)
true

julia> m.phase_encoding_direction
jpos::EncodingDirection = 2

julia> m.phase_encoding_direction == phase_encoding_direction(m)
true

julia> m.freqdim
1

julia> m.freqdim == freqdim(m)
true

julia> m.phasedim
2

julia> m.phasedim == phasedim(m)
true

julia> m.slicedim
3

julia> m.slicedim == slicedim(m)
true

julia> m.slice_start
4

julia> m.slice_start == slice_start(m)
true

julia> m.slice_end
5

julia> m.slice_end == slice_end(m)
true

julia> m.slice_duration
6.0

julia> m.slice_duration == slice_duration(m)
true
```
"""
struct EncodingDirectionMetadata
    freqdim::Int
    phasedim::Int
    slicedim::Int
    slice_start::Int
    slice_end::Int
    slice_duration::Float64
end

@properties EncodingDirectionMetadata begin
    freqdim(self) => :freqdim
    phasedim(self) => :phasedim
    slicedim(self) => :slicedim
    slice_start(self) => :slice_start
    slice_end(self) => :slice_end
    slice_duration(self) => :slice_duration
    slice_encoding_direction(self) -> slice_encoding_direction(self)
    phase_encoding_direction(self) -> phase_encoding_direction(self)
end


text:
using Base.Test
using StringParserPEG

include("../examples/calc1.jl")
@test transformed == 9

include("../examples/calc2.jl")
@test ast == 9

@test StringParserPEG.grammargrammar == Grammar(StringParserPEG.grammargrammar_string)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 17deace8-9071-11eb-3ff4-e70cf05a1bcb
begin
using Images, TestImages
using ImageQualityIndexes

img = testimage("cameraman") .|> float64
noisy_img = img .+ 0.1 .* randn(size(img))
assess_ssim(noisy_img, img) # 0.24112
assess_psnr(noisy_img, img) # 19.9697

kernel = ones(3, 3)./9 # mean filter
denoised_img = imfilter(noisy_img, kernel)
assess_psnr(denoised_img, img) # 28.4249
assess_ssim(denoised_img, img) # 0.6390
assess_msssim(denoised_img, img) # 0.8460

img = testimage("fabio");
colorfulness(img) # 68.5530
end

# ╔═╡ Cell order:
# ╠═17deace8-9071-11eb-3ff4-e70cf05a1bcb

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
@time @testset "Reachability.options" begin include("unit_options.jl") end
@time @testset "Reachability.solve" begin include("unit_solve.jl") end


text:
# Pkg.test runs with --check_bounds=1, forcing all bounds checks.
# This is incompatible with CUDAnative (see JuliaGPU/CUDAnative.jl#98)
if Base.JLOptions().check_bounds == 1
  file = @__FILE__
  run(```
    $(Base.julia_cmd())
    --color=$(Base.have_color ? "yes" : "no")
    --compiled-modules=$(Bool(Base.JLOptions().use_compiled_modules) ? "yes" : "no")
    --startup-file=$(Base.JLOptions().startupfile != 2 ? "yes" : "no")
    --code-coverage=$(["none", "user", "all"][1+Base.JLOptions().code_coverage])
    $(file)
    ```)
  exit()
end

using Flux, Test, Random
using Random

Random.seed!(0)

# So we can use the system CuArrays
insert!(LOAD_PATH, 2, "@v#.#")

@testset "Flux" begin

@info "Testing Basics"

include("utils.jl")
include("onehot.jl")
include("optimise.jl")
include("data.jl")

@info "Testing Layers"

include("layers/basic.jl")
include("layers/normalisation.jl")
include("layers/stateless.jl")
include("layers/conv.jl")

@info "Running Gradient Checks"

include("tracker.jl")

if Base.find_package("CuArrays") != nothing
  include("cuda/cuda.jl")
end

end


text:
include("utils.jl")

using HybridArrays, Random

Random.seed!(42)

@testset "Power manifold" begin

    Ms = Sphere(2)
    Ms1 = PowerManifold(Ms, 5)
    @test manifold_dimension(Ms1) == 10
    @test injectivity_radius(Ms1) == π
    Ms2 = PowerManifold(Ms, 5, 7)
    @test manifold_dimension(Ms2) == 70
    Mr = Manifolds.Rotations(3)
    Mr1 = PowerManifold(Mr, 5)
    Mrn1 = PowerManifold(Mr, Manifolds.NestedPowerRepresentation(), 5)
    @test manifold_dimension(Mr1) == 15
    Mr2 = PowerManifold(Mr, 5, 7)
    @test manifold_dimension(Mr2) == 105

    @test is_default_metric(Ms1, PowerMetric()) == Val(true)
    types_s1 = [Array{Float64,2},
                HybridArray{Tuple{3,StaticArrays.Dynamic()}, Float64, 2}]
    types_s2 = [Array{Float64,3},
                HybridArray{Tuple{3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Float64, 3}]

    types_r1 = [Array{Float64,3},
                HybridArray{Tuple{3,3,StaticArrays.Dynamic()}, Float64, 3}]

    types_rn1 = [Vector{Matrix{Float64}},
                 Vector{MMatrix{3,3,Float64}}]
    types_r2 = [Array{Float64,4},
                HybridArray{Tuple{3,3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Float64, 4}]

    retraction_methods = [Manifolds.PowerRetraction(ManifoldsBase.ExponentialRetraction())]
    inverse_retraction_methods = [Manifolds.InversePowerRetraction(ManifoldsBase.LogarithmicInverseRetraction())]

    sphere_dist = Manifolds.uniform_distribution(Ms, @SVector [1.0, 0.0, 0.0])
    power_s1_pt_dist = Manifolds.PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5))
    power_s2_pt_dist = Manifolds.PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7))
    sphere_tv_dist = Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0)
    power_s1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms1), rand(power_s1_pt_dist), sphere_tv_dist)
    power_s2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms2), rand(power_s2_pt_dist), sphere_tv_dist)

    id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0)
    power_r1_pt_dist = Manifolds.PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5))
    power_rn1_pt_dist = Manifolds.PowerPointDistribution(Mrn1, rotations_dist, [randn(Float64, 3, 3) for i in 1:5])
    power_r2_pt_dist = Manifolds.PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7))
    rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0)
    power_r1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr1), rand(power_r1_pt_dist), rotations_tv_dist)
    power_rn1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mrn1), rand(power_rn1_pt_dist), rotations_tv_dist)
    power_r2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr2), rand(power_r2_pt_dist), rotations_tv_dist)

    trim(s::String) = s[1:min(length(s), 20)]

    basis_types = (ArbitraryOrthonormalBasis(),
        ProjectedOrthonormalBasis(:svd)
    )
    for T in types_s1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_s1_pt_dist)) for _ in 1:3]
            @test injectivity_radius(Ms1,pts1[1]) == π
            basis_diag = DiagonalizingOrthonormalBasis(log(Ms1, pts1[1], pts1[2]))
            basis_arb = get_basis(Ms1, pts1[1], ArbitraryOrthonormalBasis())
            test_manifold(Ms1,
                pts1;
                test_reverse_diff = true,
                test_musical_isomorphisms = true,
                test_injectivity_radius = false,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_s1_pt_dist],
                tvector_distributions = [power_s1_tv_dist],
                basis_types_to_from = (basis_diag, basis_arb, basis_types...),
                rand_tvector_atol_multiplier = 6.0,
                retraction_atol_multiplier = 12.0,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end
    for T in types_s2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_s2_pt_dist)) for _ in 1:3]
            test_manifold(
                Ms2,
                pts2;
                test_reverse_diff = true,
                test_musical_isomorphisms = true,
                test_injectivity_radius = false,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_s2_pt_dist],
                tvector_distributions = [power_s2_tv_dist],
                rand_tvector_atol_multiplier = 6.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end

    for T in types_r1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_r1_pt_dist)) for _ in 1:3]
            test_manifold(
                Mr1,
                pts1;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_r1_pt_dist],
                tvector_distributions = [power_r1_tv_dist],
                basis_types_to_from = basis_types,
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0
            )
        end
    end

    for T in types_rn1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_rn1_pt_dist)) for _ in 1:3]
            test_manifold(
                Mrn1,
                pts1;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_rn1_pt_dist],
                tvector_distributions = [power_rn1_tv_dist],
                basis_types_to_from = basis_types,
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0
            )
        end
    end
    for T in types_r2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_r2_pt_dist)) for _ in 1:3]
            test_manifold(
                Mr2,
                pts2;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_r2_pt_dist],
                tvector_distributions = [power_r2_tv_dist],
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end

    @testset "Power manifold of Circle" begin
        pts_t = [[0.0, 1.0, 2.0], [1.0, 1.0, 2.4], [0.0, 2.0, 1.0]]
        MT = PowerManifold(Circle(), 3)
        @test representation_size(MT) == (3,)
        test_manifold(
            MT,
            pts_t;
            test_reverse_diff = false,
            test_forward_diff = false,
            test_injectivity_radius = false,
            test_musical_isomorphisms = true,
            retraction_methods = retraction_methods,
            inverse_retraction_methods = inverse_retraction_methods,
            rand_tvector_atol_multiplier = 5.0,
            retraction_atol_multiplier = 12,
            is_tangent_atol_multiplier = 12.0,
        )
    end

end


text:
@testset "hcat" begin
    simple = joinpath(pwd(), "data", "simple.fasta")
    msa = read(simple, FASTA, generatemapping=true)
    setannotresidue!(msa, "ONE", "example", "ab")
    setannotresidue!(msa, "TWO", "example", "cd")
    setannotresidue!(msa, "ONE", "OnlyONE", "xx")
    setannotresidue!(msa, "TWO", "OnlyTWO", "yy")
    msa_2 = copy(msa)
    setannotcolumn!(msa_2, "example", "HE")

    same_ref = Residue[
        'A' 'R' 'A' 'R'
        'R' 'A' 'R' 'A'
    ]

    diff_ref = Residue[
        'A' 'R' 'R' 'A'
        'R' 'A' 'A' 'R'
    ]

    annot_same = hcat(msa, msa_2)
    annot_diff = hcat(msa_2, msa[[2, 1], :])
    msa_same = hcat(
        MultipleSequenceAlignment(msa),
        MultipleSequenceAlignment(msa))
    msa_diff = hcat(
        MultipleSequenceAlignment(msa),
        MultipleSequenceAlignment(msa)[[2, 1], :])

    @testset "Same sequence names" begin
        for concatenated_msa in (annot_same, msa_same)
            @test size(concatenated_msa) == (2, 4)
            @test concatenated_msa == same_ref
            @test sequencenames(concatenated_msa) == ["ONE", "TWO"]
            @test columnnames(concatenated_msa) == ["1_1", "1_2", "2_1", "2_2"]
            @test getcolumnmapping(concatenated_msa) == [1, 2, 1, 2]
            if concatenated_msa isa AnnotatedMultipleSequenceAlignment
                @test getsequencemapping(concatenated_msa, "ONE") == [1, 2, 1, 2]
                @test getsequencemapping(concatenated_msa, "TWO") == [1, 2, 1, 2]
                @test getannotresidue(concatenated_msa, "ONE", "example") == "abab"
                @test getannotresidue(concatenated_msa, "TWO", "example") == "cdcd"
                @test getannotresidue(concatenated_msa, "ONE", "OnlyONE") == "xxxx"
                @test getannotresidue(concatenated_msa, "TWO", "OnlyTWO") == "yyyy"
                @test getannotcolumn(concatenated_msa, "example") == "  HE"
            end
            @test gethcatmapping(concatenated_msa) == [1, 1, 2, 2]
        end
    end

    @testset "Different sequence names" begin
        for concatenated_msa in (annot_diff, msa_diff)
            @test size(concatenated_msa) == (2, 4)
            @test concatenated_msa == diff_ref
            @test sequencenames(concatenated_msa) == ["ONE_&_TWO", "TWO_&_ONE"]
            @test columnnames(concatenated_msa) == ["1_1", "1_2", "2_1", "2_2"]
            @test getcolumnmapping(concatenated_msa) == [1, 2, 1, 2]
            if concatenated_msa isa AnnotatedMultipleSequenceAlignment
                @test getsequencemapping(concatenated_msa, "ONE_&_TWO") == [1, 2, 1, 2]
                @test getsequencemapping(concatenated_msa, "TWO_&_ONE") == [1, 2, 1, 2]
                @test getannotresidue(concatenated_msa, "ONE_&_TWO", "example") == "abcd"
                @test getannotresidue(concatenated_msa, "TWO_&_ONE", "example") == "cdab"
                @test getannotresidue(concatenated_msa, "ONE_&_TWO", "OnlyONE") == "xx  "
                @test getannotresidue(concatenated_msa, "ONE_&_TWO", "OnlyTWO") == "  yy"
                @test getannotresidue(concatenated_msa, "TWO_&_ONE", "OnlyONE") == "  xx"
                @test getannotresidue(concatenated_msa, "TWO_&_ONE", "OnlyTWO") == "yy  "
                @test getannotcolumn(concatenated_msa, "example") == "HE  "
            end
            @test gethcatmapping(concatenated_msa) == [1, 1, 2, 2]
        end
    end

    @testset "Multiple concatenations and annot column/residue" begin
        concatenated_msa = hcat(msa, msa_2, msa, msa_2, msa)
        @test getannotcolumn(concatenated_msa, "example") == "  HE  HE  "
    end

    @testset "IO" begin
        path = tempdir()
        tmp_file = joinpath(path, ".tmp.stockholm")
        try
            write(tmp_file, annot_diff, Stockholm)
            out_msa = read(tmp_file, Stockholm)
            @test columnnames(out_msa) == columnnames(annot_diff)
            @test out_msa == annot_diff
        finally
            if isfile(tmp_file)
                rm(tmp_file)
            end
        end
    end

    @testset "Inception" begin
        concatenated_in = hcat(msa, msa_2)
        concatenated_diff_a = hcat(msa[[2, 1], :], msa_2)
        concatenated_diff_b = hcat(msa_2, msa[[2, 1], :])

        @testset "concatenated concatenated" begin
            concatenated_out = hcat(concatenated_in, concatenated_in)
            concat_ab = hcat(concatenated_diff_a, concatenated_diff_b)

            @test size(concatenated_out) == (2, 8)
            @test sequencenames(concatenated_out) == ["ONE", "TWO"]
            @test columnnames(concatenated_out) == [
                "1_1", "1_2", "2_1", "2_2", "3_1", "3_2", "4_1", "4_2"
                ]
            @test getcolumnmapping(concatenated_out) == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concatenated_out, "ONE") == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concatenated_out, "TWO") == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concatenated_out, "ONE") == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concatenated_out, "TWO") == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getannotresidue(concatenated_out, "ONE", "example") == "abababab"
            @test getannotresidue(concatenated_out, "TWO", "example") == "cdcdcdcd"
            @test getannotresidue(concatenated_out, "ONE", "OnlyONE") == "xxxxxxxx"
            @test getannotresidue(concatenated_out, "TWO", "OnlyTWO") == "yyyyyyyy"
            @test getannotcolumn(concatenated_out, "example") == "  HE  HE"
            @test gethcatmapping(concatenated_out) == [1, 1, 2, 2, 3, 3, 4, 4]

            @test size(concat_ab) == (2, 8)
            @test sequencenames(concat_ab) == [
                "TWO_&_ONE_&_ONE_&_TWO", "ONE_&_TWO_&_TWO_&_ONE"]
            @test columnnames(concat_ab) == [
                "1_1", "1_2", "2_1", "2_2", "3_1", "3_2", "4_1", "4_2"]
            @test getcolumnmapping(concat_ab) == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concat_ab,
                "TWO_&_ONE_&_ONE_&_TWO") == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concat_ab,
                "ONE_&_TWO_&_TWO_&_ONE") == [1, 2, 1, 2, 1, 2, 1, 2]
            @test getannotresidue(concat_ab,
                "TWO_&_ONE_&_ONE_&_TWO", "example") == "cdababcd"
            @test getannotresidue(concat_ab,
                "ONE_&_TWO_&_TWO_&_ONE", "example") == "abcdcdab"
            @test getannotresidue(concat_ab,
                "TWO_&_ONE_&_ONE_&_TWO", "OnlyONE") == "  xxxx  "
            @test getannotresidue(concat_ab,
                "TWO_&_ONE_&_ONE_&_TWO", "OnlyTWO") == "yy    yy"
            @test getannotcolumn(concat_ab, "example") == "  HEHE  "
            @test gethcatmapping(concat_ab) == [1, 1, 2, 2, 3, 3, 4, 4]
        end

        @testset "concatenated non_concatenated" begin
            concatenated_msas = [
                hcat(concatenated_in, msa),
                hcat(msa, concatenated_in)
            ]
            for (i, concatenated_out) in enumerate(concatenated_msas)
                @test size(concatenated_out) == (2, 6)
                @test sequencenames(concatenated_out) == ["ONE", "TWO"]
                @test columnnames(concatenated_out) == ["1_1", "1_2", "2_1", "2_2", "3_1", "3_2"]
                @test getcolumnmapping(concatenated_out) == [1, 2, 1, 2, 1, 2]
                @test getsequencemapping(concatenated_out, "ONE") == [1, 2, 1, 2, 1, 2]
                @test getsequencemapping(concatenated_out, "TWO") == [1, 2, 1, 2, 1, 2]
                @test getsequencemapping(concatenated_out, "ONE") == [1, 2, 1, 2, 1, 2]
                @test getsequencemapping(concatenated_out, "TWO") == [1, 2, 1, 2, 1, 2]
                @test getannotresidue(concatenated_out, "ONE", "example") == "ababab"
                @test getannotresidue(concatenated_out, "TWO", "example") == "cdcdcd"
                @test getannotresidue(concatenated_out, "ONE", "OnlyONE") == "xxxxxx"
                @test getannotresidue(concatenated_out, "TWO", "OnlyTWO") == "yyyyyy"
                if i == 1
                    @test getannotcolumn(concatenated_out, "example") == "  HE  "
                else
                    @test getannotcolumn(concatenated_out, "example") == "    HE"
                end
                @test gethcatmapping(concatenated_out) == [1, 1, 2, 2, 3, 3]
            end

            concat_a = hcat(concatenated_diff_a, msa)

            @test size(concat_a) == (2, 6)
            @test sequencenames(concat_a) == [
                "TWO_&_ONE_&_ONE", "ONE_&_TWO_&_TWO"]
            @test columnnames(concat_a) == [
                "1_1", "1_2", "2_1", "2_2", "3_1", "3_2"]
            @test getcolumnmapping(concat_a) == [1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concat_a, "TWO_&_ONE_&_ONE") == [1, 2, 1, 2, 1, 2]
            @test getsequencemapping(concat_a, "ONE_&_TWO_&_TWO") == [1, 2, 1, 2, 1, 2]
            @test getannotresidue(concat_a, "TWO_&_ONE_&_ONE", "example") == "cdabab"
            @test getannotresidue(concat_a, "ONE_&_TWO_&_TWO", "example") == "abcdcd"
            @test getannotresidue(concat_a,
                "TWO_&_ONE_&_ONE", "OnlyONE") == "  xxxx"
            @test getannotresidue(concat_a,
                "TWO_&_ONE_&_ONE", "OnlyTWO") == "yy    "
            @test getannotcolumn(concat_a, "example") == "  HE  "
            @test gethcatmapping(concat_a) == [1, 1, 2, 2, 3, 3]
        end
    end
end


text:
using Distributions

@testset "gains" begin
    _,vis, amp, lcamp, cphase = load_data()

    st = scantable(vis)
    stamp = scantable(amp)
    stcp = scantable(cphase)
    stlca = scantable(lcamp)

    gcache = GainCache(st)

    tel = stations(vis)
    gamp_prior = NamedTuple{Tuple(tel)}(ntuple(_->LogNormal(0.0, 0.1), length(tel)))
    gph_prior = NamedTuple{Tuple(tel)}(ntuple(_->Normal(0.0, 0.1), length(tel)))

    gamp = GainPrior(gamp_prior, st)
    gpha = GainPrior(gph_prior, st)

    rga = rand(gamp)
    rph = rand(gpha)
    @inferred logdensityof(gamp, rga)
    @inferred logdensityof(gpha, rph)
end


text:
# ---
# title: 1028. Recover a Tree From Preorder Traversal
# id: problem1028
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/recover-a-tree-from-preorder-traversal/description/>
# hidden: true
# ---
#
# We run a preorder depth-first search (DFS) on the `root` of a binary tree.
#
# At each node in this traversal, we output `D` dashes (where `D` is the depth
# of this node), then we output the value of this node.  If the depth of a node
# is `D`, the depth of its immediate child is `D + 1`.  The depth of the `root`
# node is `0`.
#
# If a node has only one child, that child is guaranteed to be **the left
# child**.
#
# Given the output `S` of this traversal, recover the tree and return _its_
# `root`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/04/08/recover-a-tree-from-
# preorder-traversal.png)
#
#
#
#     Input: S = "1-2--3--4-5--6--7"
#     Output: [1,2,5,3,4,6,7]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/04/11/screen-
# shot-2019-04-10-at-114101-pm.png)
#
#
#
#     Input: S = "1-2--3---4-5--6---7"
#     Output: [1,2,5,3,null,6,null,4,null,7]
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2019/04/11/screen-
# shot-2019-04-10-at-114955-pm.png)
#
#
#
#     Input: S = "1-401--349---90--88"
#     Output: [1,401,null,349,88,90]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the original tree is in the range `[1, 1000]`.
#   * `1 <= Node.val <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Interpolations, Test

@testset "==" begin

    # issue #333
    knots = ([1,1.2],)
    vals = [1.0, 2.0]
    scheme = Gridded(Linear())
    itp = interpolate(knots, vals , scheme)
    @test itp == itp
    @test itp == deepcopy(itp)
    @test itp != interpolate(([1,1.3],), vals, scheme)
    @test itp != interpolate(knots, [1.0, 3.0], scheme)
    @test itp != interpolate(knots, vals, Gridded(Constant()))

    scheme = BSpline(Quadratic(Reflect(OnCell())))
    vals = [1.0, 2.0]
    itp = interpolate(vals, scheme)
    @test itp == itp
    @test itp == deepcopy(itp)
    @test itp != interpolate(vals, BSpline(Linear()))
    @test itp != interpolate(randn(2), scheme)

    vals = rand(Float64, 2,3)
    itp1 = interpolate(vals, BSpline(Quadratic(Flat(OnGrid()))))
    itp2 = interpolate(Float64.(vals), BSpline(Quadratic(Flat(OnGrid()))))
    @test itp1 == itp2
end

@testset "Core" begin
    A = reshape([0], 1, 1, 1, 1, 1)
    wis = ntuple(d->Interpolations.WeightedAdjIndex(1, (1,)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedAdjIndex(1, (1.0,)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,), (1,)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,), (1.0,)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,1), (1,0)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,1), (1.0,0.0)), ndims(A))
    @test @inferred(A[wis...]) === 0.0

    wi = Interpolations.WeightedAdjIndex(2, (0.2, 0.8))
    @test wi*2 === Interpolations.WeightedAdjIndex(2, (0.4, 1.6))
    @test wi/2 === Interpolations.WeightedAdjIndex(2, (0.1, 0.4))
    @test Interpolations.indextuple(wi) == (2, 3)
    wi = Interpolations.WeightedArbIndex((2, -1), (0.2, 0.8))
    @test wi*2 === Interpolations.WeightedArbIndex((2, -1), (0.4, 1.6))
    @test wi/2 === Interpolations.WeightedArbIndex((2, -1), (0.1, 0.4))
    @test Interpolations.indextuple(wi) == (2, -1)
end

@testset "Hash" begin
    # Issue 339: Base.hash accesses out-of-bound indexes due to AbstractArray interface
    # Relevant to Distributed computing
    xr = collect(0.0:1.0:5.0)
    yr = collect(1.0:1.0:6.0)
    itp = interpolate((xr,),yr,Gridded(Linear()))
    # using Distributed
    # addprocs(2)
    # @everywhere begin
    #    using Interpolations
    #    callitp(itp, x) = itp(x)
    # end
    # r = remotecall(callitp, 2, itp, 1.2)
    # fetch(r) # 2.2
    @test hash(itp) != 0
    etp = LinearInterpolation([2, 3], [4, 5])
    @test hash(etp.itp) != 0
    @test hash(etp) != 0
end


text:
println("# Test Point Estimation using VI")

println(" * test data..")

Z = hcat(vec([1 1 2 2]),
	vec([1 1 2 2]),
	vec([1 2 2 2]),
	vec([2 2 1 1]))

println(" * compute psm")

psm = compute_psm(Z)

println(" * test variation_of_information_lb..")

cls = vec([1 1 1 1])

@test_approx_eq_eps BNP.variation_of_information_lb(cls, psm)[1] 0.9207175 1e-6

println(" * find point estimate")

Z = hcat(vec([1 1 2 2]),
	vec([1 1 2 2]),
	vec([1 2 2 2]),
	vec([1 1 2 2]))

psm = compute_psm(Z)

(idx, value) = point_estimate(psm, maxk = 2)

@test length(unique(idx)) == 2

println("# FINISHED")


text:
module ClimatePlots

using Dates
using CFTime
using AxisArrays
using Reexport
@reexport using ClimateBase
using PyCall
using PyPlot
using Statistics
import Base.show

const mpl = PyNULL()
const plt = PyNULL()
const cmocean = PyNULL()
const cartopy = PyNULL()

function __init__()
    copy!(mpl, pyimport_conda("matplotlib", "matplotlib", "conda-forge"))
    copy!(plt, pyimport_conda("matplotlib.pyplot", "matplotlib", "conda-forge"))
    copy!(cmocean, pyimport_conda("cmocean", "cmocean", "conda-forge"))
    copy!(cartopy, pyimport_conda("cartopy", "cartopy", "conda-forge"))
end

include("types.jl")
include("maps.jl")
include("plots.jl")
include("maps_definition.jl")
include("utils.jl")

export plot, mapclimgrid, hist
export contourf, contour, pcolormesh
export BBox
export PlotInfo

end # module


text:
export IndexDefinition, IndexType

baremodule IndexType
    const HASH = 1
    const JSON = 2
end

abstract type AbstractIndex end

"""
    IndexDefinition(args)

Object used to create a specific index in Redis.

# Fields
- `args::Vector{Union{String,Number}}`: List of arguments used to define a searchable index.
"""
mutable struct IndexDefinition <: AbstractIndex
    args::Vector{Union{String,Number}}
end

"""
    IndexDefinition(;
        prefix=[],
        filter="",
        language_field="",
        language="",
        score_field="",
        score=1.0,
        payload_field="",
        index_type=IndexType.HASH
    ) -> IndexDefinition

IndexDefinition is used to define a index definition for automatic indexing on
Hash or Json update.

# KeyWords
- `prefix::AbstractArray{AbstractString}`: Tells the index which keys it should index. Each
prefix should finish with a `:`.
- `filter::AbstractString`: Filter expression with the full RediSearch aggregation
expression language.
- `language_field::AbstractString`: If set indicates the document attribute that should be
used as the document language.
- `language::AbstractString`: If set indicates the default language for documents in the
index. Default to English.
- `score_field::AbstractString`: If set indicates the document attribute that should be used
as the document's rank based on the user's ranking. Ranking must be between 0.0 and 1.0.
If not set the default score is 1
- `score::Float64`: If set indicates the default score for documents in the index.
Default score is 1.0.
- `payload_field::AbstractString`: If set indicates the document attribute that should be
used as a binary safe payload string to the document, that can be evaluated at query time by
a custom scoring function, or retrieved to the client.
- `index_type::Number`: Currently supports HASH (default) and JSON.
"""
function IndexDefinition(;
    prefix::AbstractArray=String[],
    filter::AbstractString="",
    language_field::AbstractString="",
    language::AbstractString="",
    score_field::AbstractString="",
    score::Float64=1.0,
    payload_field::AbstractString="",
    index_type::Number=IndexType.HASH,
)
    args = Vector{Union{AbstractString,Number}}()

    if index_type == IndexType.HASH
        append!(args, ["ON", "HASH"])
    elseif index_type == IndexType.JSON
        append!(args, ["ON", "JSON"])
    elseif !isnothing(index_type)
        throw(ErrorException("Cannot use index type " * index_type))
    end

    if length(prefix) > 0
        prefix_args = vcat(["PREFIX", length(prefix)], [p for p in prefix])
        append!(args, prefix_args)
    end

    if !isempty(filter)
        append!(args, ["FILTER", filter])
    end

    if !isempty(language_field)
        append!(args, ["LANGUAGE_FIELD", language_field])
    end

    if !isempty(language)
        append!(args, ["LANGUAGE", language])
    end

    if !isempty(score_field)
        append!(args, ["SCORE_FIELD", score_field])
    end

    if score != 1.0
        append!(args, ["SCORE", score])
    end

    if !isempty(payload_field)
        append!(args, ["PAYLOAD_FIELD", payload_field])
    end

    return IndexDefinition(args)
end


text:
module Model
    using Agents: ABM, AgentBasedModel, GraphSpace, add_agent_pos!, allagents, nearby_agents, nagents
    using LightGraphs: SimpleGraph, nv
    using Statistics: mean

    include("agent.jl")

    function build_model(; G::SimpleGraph, cost_model::Symbol, r::Float64, c::Float64 = 1.0, C_rate::Float64 = 0.5)::ABM
        # "r" is the PGG enhancement factor.
        model = ABM(Agent.Player, GraphSpace(G), properties = Dict(:cost_model => cost_model, :r => r, :c => c))

        # モデル上にエージェントを配置する
        foreach(id -> add_agent_pos!(Agent.Player(id, C_rate), model), 1:nv(G))

        return model
    end

    function model_step!(model::AgentBasedModel)
        # 全部Cもしくは全部Dの場合は処理をスキップする
        average_strategy = mean([agent.strategy for agent in allagents(model)])
        (average_strategy == 0 || average_strategy == 1) && return

        calc_payoffs!(model)
        set_next_strategies!(model)
        update_strategies!(model)
    end

    function calc_payoffs!(model::AgentBasedModel)
        foreach(allagents(model)) do agent
            # The incomes of a defector and a cooperator in one group are given by P_D = c * r * n_C / (k_x + 1) and P_C = P_D - c
            # C individuals with kx neighbours contribute a cost c/(k_x + 1) per game, such that the individual contribution of each C equals c independently of the number of social ties.
            neighbors = nearby_agents(agent, model)
            k_x = length(neighbors)
            n_C = length([n for n in neighbors if n.strategy])
            _c = model.cost_model == :fixed_cost ? model.c : model.c / (k_x + 1)

            agent.payoff += _c * model.r * n_C / (k_x + 1)
            agent.strategy && (agent.payoff -= _c)
            foreach(neighbors) do neighbor
                neighbor.payoff += _c * model.r * n_C / (k_x + 1)
                neighbor.strategy && (neighbor.payoff -= _c)
            end
        end
    end

    function set_next_strategies!(model::AgentBasedModel)
        for agent in allagents(model)
            # When a site x is updated, a neighbor y is drawn at random among all kx neighbors
            neighbors::Vector{Agent.Player} = collect(nearby_agents(agent, model))
            neighbor::Agent.Player = neighbors[rand(1:end)]

#             max_payoff::Float64 = 0.0
#             for n in neighbors
#                 if n.payoff > max_payoff
#                     max_payoff = n.payoff
#                     neighbor = n
#                 end
#             end

            # only if Py > Px the strategy of chosen neighbor y replaces that of x
            # with probability given by (Py - Px) / M
            Px = agent.payoff
            Py = neighbor.payoff
            agent.next_strategy = agent.strategy
#             println("agent: $(agent.pos), $(agent.strategy), $(agent.payoff)")
#             println("neighbor: $(neighbor.pos), $(neighbor.strategy), $(neighbor.payoff)")
#             println(["$(n.pos), $(n.strategy), $(n.payoff)" for n in neighbors])

            (Py < Px) && continue
            M = max([n.payoff for n in neighbors]...) - Px
#             println("Px = $(Px), Py = $(Py), M = $(M), (Py - Px) / M = $((Py - Px) / M)")

            ((Py - Px) / M < rand()) && continue
            agent.next_strategy = neighbor.strategy
#             println("$(agent.next_strategy)\n")
        end
    end

    update_strategies!(model::AgentBasedModel) = foreach(allagents(model)) do agent
            agent.strategy = agent.next_strategy
            agent.payoff = 0f0
        end
end


text:
using Documenter, ComorbidityIndexes

makedocs(
	 sitename = "Comorbidity Indexes",
   modules=[ComorbidityIndexes],
      format = Documenter.HTML(
        canonical = "https://timbp.github.io/ComorbidityIndexes.jl/stable/",
        edit_link = "main",
    ),
)

deploydocs(
    repo = "github.com/timbp/ComorbidityIndexes.jl.git",
    devbranch = "main",
)


text:
function inline_pass()
    #= pass =#
end

function inline_ellipsis()
end

function indexing()::Int64
    sum = 0
    a::Vector{Int64} = []
    for i = 0:9
        push!(a, i)
        sum += a[i+1]
    end
    return sum
end

function infer_bool(code::Int64)::Bool
    return code ∈ [1, 2, 4]
end

function show()
    a1 = 10
    b1 = 15
    b2 = 15
    @assert(b1 == 15)
    @assert(b2 == 15)
    b9 = 2
    b10 = 2
    @assert(b9 == b10)
    a2::Float64 = 2.1
    @assert(a2 == 2.1)
    res = []
    for i = 0:9
        push!(res, i)
    end
    for i = 0:2:9
        push!(res, i)
    end
    for i = 1:9
        push!(res, i)
    end
    @assert(res == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 2, 4, 6, 8, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    a3 = -(a1)
    @assert(a3 == -10)
    a4 = a3 + a1
    @assert(a4 == 0)
    t1 = a1 > 5 ? (10) : (5)
    @assert(t1 == 10)
    sum1 = indexing()
    @assert(sum1 == 45)
    a5 = [1, 2, 3]
    println(length(a5))
    a9::Vector{String} = ["a", "b", "c", "d"]
    @assert(length(a9) == 4)
    @assert(a9 == ["a", "b", "c", "d"])
    a6 = Set([1, 2, 3, 4])
    i = (x for x in a6)
    @assert(length(a6) == 4)
    @assert(next(i, nothing) == 1)
    @assert(next(i, nothing) == 2)
    @assert(next(i, nothing) == 3)
    @assert(next(i, nothing) == 4)
    @assert(next(i, nothing) === nothing)
    a7 = Dict("a" => 1, "b" => 2)
    @assert(length(a7) == 2)
    @assert(a7["a"] == 1)
    @assert(a7["b"] == 2)
    a8 = true
    a9 = false
    @assert(a8 == true)
    @assert(a9 == false)
    if true
        a1 += 1
    end
    @assert(a1 == 11)
    @assert(Bool(1))
    @assert(infer_bool(1))
    @assert(1 != 2 != 3)
    if a8
        println("true")
    elseif a4 > 0
        println("never get here")
    end
    if a1 == 11
        println("false")
    else
        println("true")
    end
    if 1 != nothing
        println("World is sane")
    end
    println(true ? ("True") : ("False"))
    if true
        println("true")
    end
    inline_pass()
    s = "1    2"
    @assert(s == "1    2")
    _escape_quotes = " foo \"bar\" baz "
    @assert(findfirst("bbc", "aaabbccc") != Nothing)
    2
    _c1, _c2 = (1, 3)
    (_c3, _, _c4) = [1, 2, 3]
    @assert(_c1 == 1)
    @assert(_c2 == 3)
    @assert(_c3 == 1)
    @assert(_c4 == 3)
end

if abspath(PROGRAM_FILE) == @__FILE__
    show()
end


text:
using BinDeps
using Libdl

const pathsep = Sys.iswindows() ? "\\" : "/"

# get the complete directory path an library filename

function libdirpath(libname::String, prefix_lib::Bool=true)
    if !startswith(libname, "lib") && prefix_lib
       libname = string("lib", libname)
    end

    libstr = Libdl.find_library(libname)
    length(libstr) === 0 && throw(ErrorException("lib$libname not found"))
    libstr = dlpath(libstr)
    return abspath(libstr)
end

function package_directory(pkgName::String)
    pkgdir = Base.find_package(pkgName)
    return abspath(joinpath(split(pkgdir, pkgName)[1], pkgName))
end

function library_filepath(libsdir::String, filenames::Vector{String}, libname::String)
    libfile = filenames[ findfirst([startswith(x,libname) for x in filenames]) ]
    return joinpath( libsdir, libfile )
end

const ArbLibsDir = abspath(normpath(joinpath(@__DIR__,"..","..","local","lib")))

libnames = readdir(ArbLibsDir)

haslibarb   = any(startswith(x,"libarb") for x in libnames)
haslibflint = any(startswith(x,"libflint") for x in libnames)

if !(haslibarb)
   throw(ErrorException(string("the Arb C library was not found in ",ArbLibsDir)))
elseif !(haslibflint)
   throw(ErrorException(string("the Flint C library was not found in ",ArbLibsDir)))
end

const LibArb   = library_filepath( ArbLibsDir, libnames, "libarb"   )
const LibFlint = library_filepath( ArbLibsDir, libnames, "libflint" )

const LibArbHandle   = dlopen(LibArb)
const LibFlintHandle = dlopen(LibFlint)

# @ccall(@libarb(library_function), ReturnType, (arg types), args)

macro libarb(libraryfunction)
    (:($libraryfunction), LibArb)
end

macro libflint(libraryfunction)
    (:($libraryfunction), LibFlint)
end

macro arblib(libraryfunction)
    :(dlsym(LibArbHandle, $libraryfunction))
end

macro flintlib(libraryfunction)
    :(dlsym(LibFlintHandle, $libraryfunction))
end