text:
const _graph_funcs = Dict{Symbol,Any}(
    :spectral => spectral_graph,
    :sfdp => sfdp_graph,
    :circular => circular_graph,
    :shell => shell_graph,
    :spring => spring_graph,
    :stress => by_axis_local_stress_graph,
    :tree => tree_graph,
    :buchheim => buchheim_graph,
    :arcdiagram => arc_diagram,
    :chorddiagram => chord_diagram
)

const _graph_inputs = Dict{Symbol,Any}(
    :spectral => :adjmat,
    :sfdp => :adjmat,
    :circular => :adjmat,
    :shell => :adjmat,
    :stress => :adjmat,
    :spring => :adjmat,
    :tree => :sourcedestiny,
    :buchheim => :adjlist,
    :arcdiagram => :sourcedestiny,
    :chorddiagram => :sourcedestiny
)

function prepare_graph_inputs(method::Symbol, inputs...)
    input_type = get(_graph_inputs, method, :sourcedestiny)
    if input_type == :adjmat
        (get_adjacency_matrix(inputs...),)
    elseif input_type == :sourcedestiny
        get_source_destiny_weight(inputs...)
    elseif input_type == :adjlist
        (get_adjacency_list(inputs...),)
    end
end

# -----------------------------------------------------

function get_source_destiny_weight(mat::AbstractArray{T,2}) where T
    nrow, ncol = size(mat)    # rows are sources and columns are destinies

    nosymmetric = !issymmetric(mat) # plots only triu for symmetric matrices
    nosparse = !issparse(mat) # doesn't plot zeros from a sparse matrix

    L = length(mat)

    source  = Array{Int}(undef, L)
    destiny = Array{Int}(undef, L)
    weights  = Array{T}(undef, L)

    idx = 1
    for i in 1:nrow, j in 1:ncol
        value = mat[i, j]
        if !isnan(value) && ( nosparse || value != zero(T) ) # TODO: deal with Nullable

            if i < j
                source[idx]  = i
                destiny[idx] = j
                weights[idx]  = value
                idx += 1
            elseif nosymmetric && (i > j)
                source[idx]  = i
                destiny[idx] = j
                weights[idx]  = value
                idx += 1
            end

        end
    end

    resize!(source, idx-1), resize!(destiny, idx-1), resize!(weights, idx-1)
end

function get_source_destiny_weight(source::AbstractVector, destiny::AbstractVector)
    if length(source) != length(destiny)
        throw(ArgumentError("Source and destiny must have the same length."))
    end
    source, destiny, Float64[ 1.0 for i in source ]
end

function get_source_destiny_weight(source::AbstractVector, destiny::AbstractVector, weights::AbstractVector)
    if !(length(source) == length(destiny) == length(weights))
        throw(ArgumentError("Source, destiny and weights must have the same length."))
    end
    source, destiny, weights
end

function get_source_destiny_weight(adjlist::AbstractVector{V}) where V<:AbstractVector{Int}
    source = Int[]
    destiny = Int[]
    for (i,l) in enumerate(adjlist)
        for j in l
            push!(source, i)
            push!(destiny, j)
        end
    end
    get_source_destiny_weight(source, destiny)
end

# -----------------------------------------------------

function get_adjacency_matrix(mat::AbstractMatrix)
    mat
end

function get_adjacency_matrix(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector)
    n = max(maximum(source), maximum(destiny))
    Matrix(sparse(source, destiny, weights, n, n))
end

function get_adjacency_matrix(adjlist::AbstractVector{V}) where V<:AbstractVector{Int}
    s,d,w = get_source_destiny_weight(adjlist)
    get_adjacency_matrix(s, d, w)
end

# -----------------------------------------------------

function get_adjacency_list(mat::AbstractMatrix)
    get_adjacency_list(get_source_destiny_weight(mat))
end

function get_adjacency_list(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector)
    n = max(maximum(source), maximum(destiny))
    adjlist = [Int[] for i=1:n]
    for (s,d) in zip(source,destiny)
        push!(adjlist[s], d)
    end
    adjlist
end

function get_adjacency_list(adjlist::AbstractVector{V}) where V<:AbstractVector{Int}
    adjlist
end

# -----------------------------------------------------

function make_symmetric(A::AbstractMatrix)
    A = copy(A)
    for i=1:size(A,1), j=(i+1):size(A,2)
        A[i,j] = A[j,i] = A[i,j]+A[j,i]
    end
    A
end

function compute_laplacian(adjmat::AbstractMatrix, node_weights::AbstractVector)
    n, m = size(adjmat)
    # @show size(adjmat), size(node_weights)
    @assert n == m == length(node_weights)

    # scale the edge values by the product of node_weights, so that "heavier" nodes also form
    # stronger connections
    adjmat = adjmat .* sqrt(node_weights * node_weights')

    # D is a diagonal matrix with the degrees (total weights for that node) on the diagonal
    deg = vec(sum(adjmat; dims=1)) - diag(adjmat)
    D = diagm(0 => deg)

    # Laplacian (L = D - adjmat)
    L = eltype(adjmat)[i == j ? deg[i] : -adjmat[i,j] for i=1:n,j=1:n]

    L, D
end

import LightGraphs

# TODO: so much wasteful conversion... do better
function estimate_distance(adjmat::AbstractMatrix)
    source, destiny, weights = get_source_destiny_weight(sparse(adjmat))

    g = LightGraphs.Graph(adjmat)
    dists = convert(Matrix{Float64}, hcat(map(i->LightGraphs.dijkstra_shortest_paths(g, i).dists, LightGraphs.vertices(g))...))
    tot = 0.0; cnt = 0
    for (i,d) in enumerate(dists)
        if d < 1e10
            tot += d
            cnt += 1
        end
    end
    avg = cnt > 0 ? tot / cnt : 1.0
    for (i,d) in enumerate(dists)
        if d > 1e10
            dists[i] = 3avg
        end
    end
    dists
end

function get_source_destiny_weight(g::LightGraphs.AbstractGraph)
    source = Vector{Int}()
    destiny =  Vector{Int}()
    sizehint!(source, LightGraphs.nv(g))
    sizehint!(destiny, LightGraphs.nv(g))
    for e in LightGraphs.edges(g)
      push!(source, LightGraphs.src(e))
      push!(destiny, LightGraphs.dst(e))
    end
    get_source_destiny_weight(source, destiny)
end

function get_adjacency_matrix(g::LightGraphs.AbstractGraph)
    adjacency_matrix(g)
end

function get_adjacency_matrix(source::AbstractVector{Int}, destiny::AbstractVector{Int})
    get_adjacency_matrix(source, destiny, ones(length(source)))
end

function get_adjacency_list(g::LightGraphs.AbstractGraph)
    g.fadjlist
end

# -----------------------------------------------------

# a graphplot takes in either an (N x N) adjacency matrix
#   note: you may want to pass node weights to markersize or marker_z
# A graph has N nodes where adj_mat[i,j] is the strength of edge i --> j.  (adj_mat[i,j]==0 implies no edge)

# NOTE: this is for undirected graphs... adjmat should be symmetric and non-negative

const graph_aliases = Dict(:curvature_scalar => [:curvaturescalar,:curvature],
                           :node_weights => [:nodeweights],
                           :nodeshape => [:node_shape,:markershape],
                           :nodesize => [:node_size,:markersize],
                           :nodecolor => [:marker_color,:markercolor],
                           :node_z => [:marker_z],
                           :nodestrokealpha => [:markerstrokealpha],
                           :nodealpha => [:markeralpha],
                           :nodestrokewidth => [:markerstrokewidth],
                           :nodestrokealpha => [:markerstrokealpha],
                           :nodestrokecolor => [:markerstrokecolor],
                           :nodestrokestyle => [:markerstrokestyle],
                           :shorten => [:shorten_edge],
                           :axis_buffer => [:axisbuffer],
                           :edgewidth => [:edge_width,:ew],
                           :edgelabel => [:edge_label,:el],
                           :edgelabel_offset => [:edgelabeloffset,:elo],
                           :self_edge_size => [:selfedgesize,:ses],
                           :edge_label_box => [:edgelabelbox,:edgelabel_box,:elb],
)

"""
    graphplot(g; kwargs...)

Visualize the graph `g`, where `g` represents a graph via a matrix or a
`LightGraphs.graph`.
## Keyword arguments
```
dim = 2
free_dims = nothing
T = Float64
curves = true
curvature_scalar = 0.05
root = :top
node_weights = nothing
names = []
fontsize = 7
nodeshape = :hexagon
nodesize = 0.1
node_z = nothing
nodecolor = 1
nodestrokealpha = 1
nodealpha = 1
nodestrokewidth = 1
nodestrokecolor = :black
nodestrokestyle = :solid
nodestroke_z = nothing
x = nothing
y = nothing
z = nothing
method = :stress
func = get(_graph_funcs, method, by_axis_local_stress_graph)
shorten = 0.0
axis_buffer = 0.2
layout_kw = Dict{Symbol,Any}()
edgewidth = (s,d,w)->1
edgelabel = nothing
edgelabel_offset = 0.0
self_edge_size = 0.1
edge_label_box = true
edge_z = nothing
edgecolor = :black
edgestyle = :solid
```

See the [documentation]( http://docs.juliaplots.org/latest/graphrecipes/introduction/ ) for
more details.
"""
@userplot GraphPlot

@recipe function f(g::GraphPlot;
                   dim = 2,
                   free_dims = nothing,
                   T = Float64,
                   curves = true,
                   curvature_scalar = 0.05,
                   root = :top,
                   node_weights = nothing,
                   names = [],
                   fontsize = 7,
                   nodeshape = :hexagon,
                   nodesize = 0.1,
                   node_z = nothing,
                   nodecolor = 1,
                   nodestrokealpha = 1,
                   nodealpha = 1,
                   nodestrokewidth = 1,
                   nodestrokecolor = :black,
                   nodestrokestyle = :solid,
                   nodestroke_z = nothing,
                   x = nothing,
                   y = nothing,
                   z = nothing,
                   method = :stress,
                   func = get(_graph_funcs, method, by_axis_local_stress_graph),
                   shorten = 0.0,
                   axis_buffer = 0.2,
                   layout_kw = Dict{Symbol,Any}(),
                   edgewidth = (s,d,w)->1,
                   edgelabel = nothing,
                   edgelabel_offset = 0.0,
                   self_edge_size = 0.1,
                   edge_label_box = true,
                   edge_z = nothing,
                   edgecolor = :black,
                   edgestyle = :solid,
                  )
    # Process the args so that they are a LightGraphs.Graph.
    if length(g.args) <= 1 && !(eltype(g.args[1]) <: AbstractArray) && !(g.args[1] isa LightGraphs.AbstractGraph) && method != :chorddiagram && method != :arcdiagram
        if !LinearAlgebra.issymmetric(g.args[1]) || any(diag(g.args[1]) .!= zeros(length(diag(g.args[1]))))
            g.args = (LightGraphs.DiGraph(g.args[1]),)
        elseif LinearAlgebra.issymmetric(g.args[1])
            g.args = (LightGraphs.Graph(g.args[1]),)
        end
    end

    # To process aliases that are unique to graphplot, find aliases that are in
    # plotattributes and replace the attributes with their aliases. Then delete the alias
    # names from the plotattributes dictionary.
    @process_aliases plotattributes graph_aliases
    for arg in keys(graph_aliases)
        remove_aliases!(arg, plotattributes, graph_aliases)
    end
    # The above process will remove all marker properties from the plotattributes
    # dictionary. To enusre consistency between markers and nodes, we replace all marker
    # properties with the corresponding node property.
    marker_node_collection = zip([:markershape,:markersize,:markercolor,
                                  :marker_z,:markerstrokealpha,:markeralpha,
                                  :markerstrokewidth,:markerstrokealpha,
                                  :markerstrokecolor,:markerstrokestyle],
                                 [nodeshape,nodesize,nodecolor,node_z,nodestrokealpha,
                                  nodealpha,nodestrokewidth,nodestrokealpha,
                                  nodestrokecolor,nodestrokestyle])
    for (markerproperty, nodeproperty) in marker_node_collection
        plotattributes[markerproperty] = nodeproperty
    end

    @assert dim in (2, 3)
    _3d = dim == 3
    adj_mat = get_adjacency_matrix(g.args...)
    isdirected = (g.args[1] isa DiGraph || !issymmetric(adj_mat)) &&
                 !in(method, (:tree, :buchheim)) && !(get(plotattributes, :arrow, true) == false)
    if isdirected && (g.args[1] isa Matrix)
        g = GraphPlot((adjacency_matrix(DiGraph(g.args[1])),))
    end

    source, destiny, weights = get_source_destiny_weight(g.args...)
    if !(eltype(source) <: Integer)
        names = unique(sort(vcat(source,destiny)))
        source = Int[findfirst(names, si) for si in source]
        destiny = Int[findfirst(names, di) for di in destiny]
    end
    n = max(maximum(source), maximum(destiny))

    if isnothing(node_weights)
        node_weights = ones(n)
    end
    @assert length(node_weights) == n

    xyz = _3d ? (x,y,z) : (x,y)
    numnothing = count(isnothing, xyz)

    # do we want to compute coordinates?
    if numnothing > 0
        if isnothing(free_dims)
            # compute free_dims
            free_dims = findall(isnothing, xyz)
        end
        x, y, z = func(
            prepare_graph_inputs(method, source, destiny, weights)...;
            node_weights = node_weights,
            dim = dim,
            free_dims = free_dims,
            root = root,
            layout_kw...
        )
    end

    # reorient the points after root
    if root in (:left,:right)
        x,y = y,-x
    end
    if root == :left
        x,y = -x, y
    end
    if root == :bottom
        x,y = x, -y
    end

    # Since we do nodehapes manually, they only work with aspect_ratio=1.
    # TODO: rescale the nodeshapes based on the ranges of x,y,z.
    aspect_ratio --> 1
    if length(axis_buffer) == 1
        axis_buffer = fill(axis_buffer, dim)
    end

    # center and rescale to the widest of all dimensions
    if method == :arcdiagram
        xl, yl = arcdiagram_limits(x, source, destiny)
        xlims --> xl
        ylims --> yl
        ratio --> :equal
    elseif all(axis_buffer .< 0) # equal axes
        ahw = 1.2 * 0.5 * maximum(v -> maximum(v)-minimum(v), xyz)
        xcenter = mean(extrema(x))
        #xlims --> (xcenter-ahw, xcenter+ahw)
        ycenter = mean(extrema(y))
        #ylims --> (ycenter-ahw, ycenter+ahw)
        if _3d
            zcenter = mean(extrema(z))
            #zlims --> (zcenter-ahw, zcenter+ahw)
        end
    else
        xlims = ignorenan_extrema(x)
        if method != :chorddiagram && numnothing > 0
            x .-= mean(x)
            x /= (xlims[2] - xlims[1])
            y .-= mean(y)
            ylims = ignorenan_extrema(y)
            y /= (ylims[2] - ylims[1])
        end
        xlims --> extrema_plus_buffer(x, axis_buffer[1])
        ylims --> extrema_plus_buffer(y, axis_buffer[2])
        if _3d
            if method != :chorddiagram && numnothing > 0
                zlims = ignorenan_extrema(z)
                z .-= mean(z)
                z /= (zlims[2] - zlims[1])
            end
            zlims --> extrema_plus_buffer(z, axis_buffer[3])
        end
    end
    xyz = _3d ? (x,y,z) : (x,y)
    # Get the coordinates for the edges of the nodes.
    node_vec_vec_xy = []
    nodewidth = 0.0
    nodewidth_array = Vector{Float64}(undef, length(x))
    if !(nodeshape isa Array)
        nodeshape = repeat([nodeshape], length(x))
    end
    if !_3d
        for i in 1:length(x)
            node_number = i % length(nodeshape) == 0 ? length(nodeshape) : i % length(nodeshape)
            node_weight = isnothing(node_weights) ? 1 : (10 + 100node_weights[i]/sum(node_weights))/50
            xextent, yextent = if isempty(names)
                [x[i] .+ [-0.5nodesize*node_weight, 0.5nodesize*node_weight], y[i] .+ [-0.5nodesize*node_weight, 0.5nodesize*node_weight]]
            else
                annotation_extent(plotattributes,
                                  (x[i], y[i], names[ifelse(i % length(names) == 0, length(names),
                                                            i % length(names))], fontsize*nodesize*node_weight))
            end
            nodewidth = xextent[2] - xextent[1]
            nodewidth_array[i] = nodewidth
            if nodeshape[node_number] == :circle
                push!(node_vec_vec_xy, partialcircle(0, 2π, [x[i], y[i]],
                                                     80, nodewidth/2))
            elseif (nodeshape[node_number] == :rect) || (nodeshape[node_number] == :rectangle)
                push!(node_vec_vec_xy, [(xextent[1],yextent[1]),
                                        (xextent[2],yextent[1]),
                                        (xextent[2],yextent[2]),
                                        (xextent[1],yextent[2]),
                                        (xextent[1],yextent[1])])
            elseif nodeshape[node_number] == :hexagon
                push!(node_vec_vec_xy, partialcircle(0, 2π, [x[i], y[i]],
                                                     7, nodewidth/2))
            elseif nodeshape[node_number] == :ellipse
                nodeheight = (yextent[2] - yextent[1])
                push!(node_vec_vec_xy, partialellipse(0, 2π, [x[i], y[i]],
                                                      80, nodewidth/2, nodeheight/2))
            else
                error("Unknown nodeshape: $(nodeshape[node_number]). Choose from :circle, ellipse, :hexagon, :rect or :rectangle.")
            end
        end
    else
        @assert _3d # TODO Make 3d work.
    end
    # The node_perimter_info list contains the information needed to construct the
    # information in node_vec_vec_xy. For example, if (nodeshape[i]==:circle && !_3d),
    # then all of the information in node_vec_vec_xy[i] can be summarised with three
    # numbers describing the center and the radius of the circle.
    node_perimeter_info = []
    for i in 1:length(node_vec_vec_xy)
        if nodeshape[i] == :circle
            push!(node_perimeter_info,
                  GeometryTypes.Circle(Point((convert(T, x[i]),convert(T, y[i]))), nodewidth_array[i]/2))
        else
            push!(node_perimeter_info, node_vec_vec_xy[i])
        end
    end

    # generate a list of colors, one per segment
    segment_colors = get(plotattributes, :linecolor, nothing)
    edge_label_array = Vector{Tuple}()
    edge_label_box_vertices_array = Vector{Array}()
    if !isa(edgelabel, Dict) && !isnothing(edgelabel)
        tmp = Dict()
        if length(size(edgelabel)) < 2
            matrix_size = round(Int, sqrt(length(edgelabel)))
            edgelabel = reshape(edgelabel, matrix_size, matrix_size)
        end
        for i in 1:size(edgelabel)[1]
            for j in 1:size(edgelabel)[2]
                if islabel(edgelabel[i, j])
                    tmp[(i, j)] = edgelabel[i, j]
                end
            end
        end
        edgelabel = tmp
    end
    # If the edgelabel dictionary is full of length two tuples, then make all of the
    # tuples length three with last element 1. (i.e. a multigraph that has no extra
    # edges).
    if edgelabel isa Dict
        edgelabel = convert(Dict{Any,Any}, edgelabel)
        for key in keys(edgelabel)
            if length(key) == 2
                edgelabel[(key..., 1)] = edgelabel[key]
            end
        end
    end
    edge_has_been_seen = Dict()
    for edge in zip(source, destiny)
        edge_has_been_seen[edge] = 0
    end
    if length(curvature_scalar) == 1
        curvature_scalar = fill(curvature_scalar, size(g.args[1])[1], size(g.args[1])[1])
    end

    edges_list = (T[], T[], T[], T[])
    # TODO do a proper job of calculating nsegments.
    nsegments = if curves && (method in (:tree, :buchheim))
        4
    elseif method == :chorddiagram
        3
    elseif method == :arcdiagram
        30
    elseif curves
        50
    else
        2
    end

    for (i, (si, di, wi)) in enumerate(zip(source, destiny, weights))
        edge_has_been_seen[(si, di)] += 1
        xseg = Vector{Float64}()
        yseg = Vector{Float64}()
        zseg = Vector{Float64}()
        l_wg = Vector{Float64}()

        # add a line segment
        xsi, ysi, xdi, ydi = shorten_segment(x[si], y[si], x[di], y[di], shorten)
        θ = (edge_has_been_seen[(si, di)] - 1)*pi/8
        if isdirected && si != di && !_3d
            xpt, ypt = if method != :chorddiagram
                control_point(xsi, xdi,
                              ysi, ydi,
                              edge_has_been_seen[(si, di)]
                              *curvature_scalar[si, di]*sign(si - di))
            else
                (0.0, 0.0)
            end
            # For directed graphs, shorten the line segment so that the edge ends at
            # the perimeter of the destiny node.
            if isdirected
                _, _, xdi, ydi = nearest_intersection(xpt, ypt, x[di], y[di],
                                                      node_perimeter_info[di])
            end
        end
        if curves
            if method in (:tree, :buchheim)
                # for trees, shorten should be on one axis only
                # dist = sqrt((x[di]-x[si])^2 + (y[di]-y[si])^2) * shorten
                dist = shorten * (root in (:left,:bottom) ? 1 : -1)
                ishoriz = root in (:left,:right)
                xsi, xdi = (ishoriz ? (x[si]+dist,x[di]-dist) : (x[si],x[di]))
                ysi, ydi = (ishoriz ? (y[si],y[di]) : (y[si]+dist,y[di]-dist))
                xpts, ypts = directed_curve(xsi, xdi, ysi, ydi,
                            xview=get(plotattributes, :xlims, (0,1)), yview=get(plotattributes, :ylims, (0,1)), root=root)
                append!(xseg, xpts)
                append!(yseg, ypts)
                append!(l_wg, [ wi for i in 1:length(xpts) ] )
            elseif method == :arcdiagram
                r  = (xdi - xsi) / 2
                x₀ = (xdi + xsi) / 2
                θ = range(0, stop=π, length=30)
                xpts = x₀ .+ r .* cos.(θ)
                ypts = r .* sin.(θ) .+ ysi # ysi == ydi
                for x in xpts
                    push!(xseg, x)
                    push!(l_wg, wi)
                end
                # push!(xseg, NaN)
                for y in ypts
                    push!(yseg, y)
                end
                # push!(yseg, NaN)
            else
                xpt, ypt = if method != :chorddiagram
                    control_point(xsi, x[di],
                                  ysi, y[di],
                                  edge_has_been_seen[(si, di)]*curvature_scalar[si, di]*sign(si - di))
                else
                    (0.0, 0.0)
                end
                xpts = [xsi, xpt, xdi]
                ypts = [ysi, ypt, ydi]
                t = range(0, stop=1, length=3)
                A = hcat(xpts, ypts)
                itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)
                tfine = range(0, stop=1, length=nsegments)
                xpts, ypts = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
                if !isnothing(edgelabel) && haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                    q = control_point(xsi, x[di],
                                      ysi, y[di],
                                      (edgelabel_offset
                                      + edge_has_been_seen[(si, di)]*curvature_scalar[si, di])*sign(si - di))
                    push!(edge_label_array,
                          (q...,
                           string(edgelabel[(si, di, edge_has_been_seen[(si, di)])]), fontsize))
                    edge_label_box_vertices = (
                    annotation_extent(plotattributes,
                                      (q[1], q[2],
                                      edgelabel[(si, di, edge_has_been_seen[(si, di)])],
                                      0.05fontsize)))
                    if !any(isnan.(q))
                        push!(edge_label_box_vertices_array, edge_label_box_vertices)
                    end
                end
                if method != :chorddiagram && !_3d
                    append!(xseg, xpts)
                    append!(yseg, ypts)
                    push!(l_wg, wi)
                else
                    push!(xseg, xsi, xpt, xdi)
                    push!(yseg, ysi, ypt, ydi)
                    _3d && push!(zseg, z[si], z[si], z[di])
                    push!(l_wg, wi)
                end
            end
        else
            push!(xseg, xsi, xdi)
            push!(yseg, ysi, ydi)
            _3d && push!(zseg, z[si], z[di])
            if !isnothing(edgelabel) && haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                q = [(xsi + xdi)/2, (ysi + ydi)/2]
                push!(edge_label_array,
                      (q...,
                       string(edgelabel[(si, di, edge_has_been_seen[(si, di)])]), fontsize))
                edge_label_box_vertices = (
                annotation_extent(plotattributes,
                                  (q[1], q[2],
                                  edgelabel[(si, di, edge_has_been_seen[(si, di)])],
                                  0.05fontsize)))
                if !any(isnan.(q))
                    push!(edge_label_box_vertices_array, edge_label_box_vertices)
                end
            end
        end
        if si == di && !_3d
            inds = 1:n .!= si
            self_edge_angle = pi/8 + (edge_has_been_seen[(si, di)] - 1)*pi/8
            θ1 = unoccupied_angle(xsi, ysi, x[inds], y[inds]) - self_edge_angle/2
            θ2 = θ1 + self_edge_angle
            nodewidth = nodewidth_array[si]
            if nodeshape == :circle
                xpts = [xsi + nodewidth*cos(θ1)/2,
                        NaN, NaN, NaN,
                        xsi + nodewidth*cos(θ2)/2]
                xpts[2] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ1)
                xpts[3] = mean([xpts[1], xpts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*cos((θ1 + θ2)/2)
                xpts[4] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ2)
                ypts = [ysi + nodewidth*sin(θ1)/2,
                        NaN, NaN, NaN,
                        ysi + nodewidth*sin(θ2)/2]
                ypts[2] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ1)
                ypts[3] = mean([ypts[1], ypts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*sin((θ1 + θ2)/2)
                ypts[4] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ2)
                t = range(0, stop=1, length=5)
                A = hcat(xpts, ypts)
                itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)
                tfine = range(0, stop=1, length=nsegments)
                xpts, ypts = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
            else
                _, _,
                start_point1,
                start_point2 = nearest_intersection(xsi, ysi,
                                                    xsi + 2nodewidth*cos(θ1),
                                                    ysi + 2nodewidth*sin(θ1),
                                                    node_vec_vec_xy[si])
                _, _, end_point1,
                end_point2 = nearest_intersection(xsi + edge_has_been_seen[(si, di)]*(nodewidth + self_edge_size)*cos(θ2),
                                                  ysi + edge_has_been_seen[(si, di)]*(nodewidth + self_edge_size)*sin(θ2),
                                                  xsi,
                                                  ysi,
                                                  node_vec_vec_xy[si])
                xpts = [start_point1,
                        NaN, NaN, NaN,
                        end_point1]
                xpts[2] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ1)
                xpts[3] = mean([xpts[1], xpts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*cos((θ1 + θ2)/2)
                xpts[4] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ2)
                ypts = [start_point2,
                        NaN, NaN, NaN,
                        end_point2]
                ypts[2] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ1)
                ypts[3] = mean([ypts[1], ypts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*sin((θ1 + θ2)/2)
                ypts[4] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ2)
                t = range(0, stop=1, length=5)
                A = hcat(xpts, ypts)
                itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)
                tfine = range(0, stop=1, length=nsegments)
                xpts, ypts = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
            end
            append!(xseg, xpts)
            append!(yseg, ypts)
            mid_ind = div(length(xpts), 2)
            q = [xpts[mid_ind] + edgelabel_offset*cos((θ1 + θ2)/2),
                 ypts[mid_ind] + edgelabel_offset*sin((θ1 + θ2)/2)]
            if !isnothing(edgelabel) && haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                push!(edge_label_array,
                      (q...,
                       string(edgelabel[(si, di, edge_has_been_seen[(si, di)])]), fontsize))
                edge_label_box_vertices = annotation_extent(plotattributes, (q...,
                                                            edgelabel[(si, di, edge_has_been_seen[(si, di)])],
                                                            0.05fontsize))
                if !any(isnan.(q))
                    push!(edge_label_box_vertices_array, edge_label_box_vertices)
                end
            end
        end
        append!(edges_list[1], xseg[.!isnan.(xseg)])
        append!(edges_list[2], yseg[.!isnan.(yseg)])
        _3d && append!(edges_list[3], zseg[.!isnan.(zseg)])
        append!(edges_list[4], l_wg[.!isnan.(l_wg)])
    end

    if !_3d
        edges_list = (reshape(edges_list[1], nsegments,
                              round(Int, length(edges_list[1])/nsegments)),
                      reshape(edges_list[2], nsegments,
                              round(Int, length(edges_list[2])/nsegments)))
    else
        edges_list = (reshape(edges_list[1], 3,
                              round(Int, length(edges_list[1])/3)),
                      reshape(edges_list[2], 3,
                              round(Int, length(edges_list[2])/3)),
                      reshape(edges_list[3], 3,
                              round(Int, length(edges_list[3])/3)))
    end

    @series begin
        seriestype := if method in (:tree, :buchheim, :chorddiagram)
            :curves
        else
            if _3d
                # TODO make curves work
                if curves
                    :curves
                end
            else
                :path
            end
        end

        colorbar_entry := true

        edge_z = process_edge_attribute(edge_z, source, destiny, weights)
        edgewidth = process_edge_attribute(edgewidth, source, destiny, weights)
        edgecolor = process_edge_attribute(edgecolor, source, destiny, weights)
        edgestyle = process_edge_attribute(edgestyle, source, destiny, weights)

        !isnothing(edge_z) && (line_z := edge_z)
        linewidthattr = get(plotattributes, :linewidth, 1)
        linewidth := linewidthattr*edgewidth
        fillalpha := 1
        linecolor := edgecolor
        linestyle := get(plotattributes, :linestyle, edgestyle)
        markershape := :none
        markersize := 0
        markeralpha := 0
        markercolor := :black
        marker_z := nothing
        isdirected && (arrow --> :simple, :head, 0.3, 0.3)
        primary := false

        _3d ? edges_list[1:3] : edges_list[1:2]
    end
    # The boxes around edge labels are defined as another list of series that sits on top
    # of the series for the edges.
    edge_has_been_seen = Dict()
    for edge in zip(source, destiny)
        edge_has_been_seen[edge] = 0
    end
    if edge_label_box && !isnothing(edgelabel)
        index = 0
        for (i, (si, di, wi)) in enumerate(zip(source, destiny, weights))
            edge_has_been_seen[(si, di)] += 1
            if haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                index += 1
                @series begin
                    seriestype := :shape

                    colorbar_entry --> false
                    fillcolor --> get(plotattributes, :background_color, :white)
                    linewidth --> 0
                    linealpha --> 0
                    edge_label_box_vertices = edge_label_box_vertices_array[index]
                    ([edge_label_box_vertices[1][1], edge_label_box_vertices[1][2],
                      edge_label_box_vertices[1][2], edge_label_box_vertices[1][1],
                      edge_label_box_vertices[1][1]],
                     [edge_label_box_vertices[2][1], edge_label_box_vertices[2][1],
                      edge_label_box_vertices[2][2], edge_label_box_vertices[2][2],
                      edge_label_box_vertices[2][1]])
                end
            end
        end
    end

    framestyle := :none
    axis := nothing
    legend --> false
    if method == :chorddiagram
        seriestype := :scatter
        markersize := 0
        markeralpha := 0
        ratio --> :equal
        if length(names) == length(x)
            annotations := [(x[i], y[i], names[i]) for i in 1:length(x)]
        end
        @series begin
            seriestype := :shape
            N = length(x)
            angles = Vector{Float64}(undef, N)
            for i in 1:N
                if y[i] > 0
                    angles[i] = acos(x[i])
                else
                    angles[i] = 2pi - acos(x[i])
                end
            end
            δ = 0.4 * (angles[2] - angles[1])
            vec_vec_xy = [ arcshape(Θ-δ,Θ+δ) for Θ in angles ] # Shape
            [ [ xy[1] for xy in vec_xy ] for vec_xy in vec_vec_xy ], [ [ xy[2] for xy in vec_xy ] for vec_xy in vec_vec_xy ]
        end
    else
        if _3d
            # seriestype := :scatter3d
            # linewidth := 0
            # linealpha := 0
            # markercolor := nodecolor
            # series_annotations --> map(string,names)
            # markersize --> (10 .+ (100 .* node_weights) ./ sum(node_weights))
        else
            @series begin
                seriestype := :shape

                colorbar_entry := true
                fill_z --> node_z
                fillalpha := nodealpha
                fillcolor := nodecolor
                markersize := 0
                markeralpha := 0
                fillalpha := 1
                linewidth := nodestrokewidth
                linealpha := nodestrokealpha
                linecolor := nodestrokecolor
                linestyle := nodestrokestyle
                line_z := nodestroke_z

                nodeperimeters = (T[], T[])
                for (i, vec_xy) in enumerate(node_vec_vec_xy)
                        append!(nodeperimeters[1], [xy[1] for xy in vec_xy])
                        push!(nodeperimeters[1], NaN)

                        append!(nodeperimeters[2], [xy[2] for xy in vec_xy])
                        push!(nodeperimeters[2], NaN)
                end
                nodeperimeters

                # if _3d
                #     seriestype := :volume
                #     ([[xyz[1] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz],
                #      [[xyz[2] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz],
                #      [[xyz[3] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz])
                # end
            end

            if isempty(names)
                seriestype := :scatter

                colorbar_entry --> false
                markersize --> 0
                markeralpha --> 0
                !isnothing(edgelabel) && (annotation --> edge_label_array)
            else
                seriestype := :scatter

                colorbar_entry --> false
                markersize --> 0
                markeralpha --> 0
                annotations --> [edge_label_array ; [(x[i], y[i],
                                                    names[ifelse(i % length(names) == 0,
                                                                  length(names),
                                                                  i % length(names))],
                                                    fontsize) for i in 1:length(x)]]
            end
        end
    end
    xyz
end

@recipe function f(g::AbstractGraph)
    GraphPlot(get_source_destiny_weight(get_adjacency_list(g)))
end


text:
using SamplingTaskFarm
using Documenter

DocMeta.setdocmeta!(SamplingTaskFarm, :DocTestSetup, :(using SamplingTaskFarm); recursive=true)

makedocs(;
    modules=[SamplingTaskFarm],
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com> and contributors",
    repo="https://github.com/jagot/SamplingTaskFarm.jl/blob/{commit}{path}#{line}",
    sitename="SamplingTaskFarm.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jagot.github.io/SamplingTaskFarm.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jagot/SamplingTaskFarm.jl",
)

ext: jl


text:
struct MPPI{DT<:AbstractFloat,nu,Covar<:AbsMat{DT},Value,Env,Init,Obs,State}
    # MPPI parameters
    K::Int
    H::Int
    lambda::DT
    gamma::DT
    value::Value
    envs::Vector{Env} # one per thread
    initfn!::Init

    # internal
    noise::Array{DT,3}
    covar_ul::Covar
    meantrajectory::Matrix{DT}
    trajectorycosts::Vector{DT}
    obsbuffers::Vector{Obs}
    statebuffers::Vector{State}

    function MPPI{DT}(
        env_tconstructor,
        K::Integer,
        H::Integer,
        covar::AbstractMatrix{<:Real},
        lambda::Real,
        gamma::Real,
        value,
        initfn!,
    ) where {DT<:AbstractFloat}
        envs = [e for e in env_tconstructor(Threads.nthreads())]

        ssp = statespace(first(envs))
        asp = actionspace(first(envs))
        osp = observationspace(first(envs))

        if !(asp isa Shapes.AbstractVectorShape)
            throw(ArgumentError("actionspace(env) must be a Shape.AbstractVectorShape"))
        end
        K > 0 || error("K must be > 0. Got $K.")
        H > 0 || error("H must be > 0. Got $H.")
        0 < lambda <= 1 || throw(ArgumentError("lambda must be in interval (0, 1]"))
        0 < gamma <= 1 || throw(ArgumentError("gamma must be in interval (0, 1]"))

        covar_ul = convert(AbsMat{DT}, cholesky(covar).UL)
        meantrajectory = zeros(DT, asp, H)
        trajectorycosts = zeros(DT, K)
        noise = zeros(DT, asp, H, K)
        obsbuffers = [allocate(osp) for _ = 1:Threads.nthreads()]
        statebuffers = [allocate(ssp) for _ = 1:Threads.nthreads()]

        new{
            DT,
            length(asp),
            typeof(covar_ul),
            typeof(value),
            eltype(envs),
            typeof(initfn!),
            eltype(obsbuffers),
            eltype(statebuffers),
        }(
            K,
            H,
            lambda,
            gamma,
            value,
            envs,
            initfn!,
            noise,
            covar_ul,
            meantrajectory,
            trajectorycosts,
            obsbuffers,
            statebuffers
        )
    end
end

"""
    $(TYPEDEF)

    MPPI{DT<:AbstractFloat}(args...; kwargs...) -> MPPI
    MPPI(args...; kwargs...) -> MPPI

Construct an instance of  `MPPI` with `args` and `kwargs`, where `DT <: AbstractFloat` is
the element type used for pre-allocated buffers, which defaults to Float32.

In the following explanation of the `MPPI` constructor, we use the
following notation:
- `U::Matrix`: the canonical control vector ``(u_{1}, u_{2}, \\dots, u_{H})``, where
    `size(U) == (length(actionspace(env)), H)`.

# Keywords

- `env_tconstructor`: a function with signature `env_tconstructor(n)` that returns `n`
    instances of `T`, where `T <: AbstractEnvironment`.
- `H::Integer`: Length of sampled trajectories.
- `K::Integer`: Number of trajectories to sample.
- `covar::AbstractMatrix`: The covariance matrix for the Normal distribution from which
    control pertubations are sampled from.
- `gamma::Real`: Reward discount, applied as `gamma^(t - 1) * reward[t]`.
- `lambda::Real`: Temperature parameter for the exponential reweighting of sampled
    trajectories. In the limit that lambda approaches 0, `U` is set to the highest reward
    trajectory. Conversely, as `lambda` approaches infinity, `U` is computed as the
    unweighted-average of the samples trajectories.
- `value`: a function mapping observations to scalar rewards, with the signature
    `value(obs::AbstractVector) --> reward::Real`
- `initfn!`: A function with the signature `initfn!(U::Matrix)` used for
    re-initializing `U` after shifting it. Defaults to setting the last
    element of `U` to 0.
"""
function MPPI{DT}(;
    env_tconstructor,
    covar,
    lambda,
    K,
    H,
    gamma = 1,
    value = zerofn,
    initfn! = default_initfn!,
) where {DT<:AbstractFloat}
    MPPI{DT}(env_tconstructor, K, H, covar, lambda, gamma, value, initfn!)
end

MPPI(args...; kwargs...) = MPPI{Float32}(args...; kwargs...)

"""
    $(TYPEDSIGNATURES)

Resets the canonical control vector to zeros.
"""
LyceumBase.reset!(m::MPPI) = (fill!(m.meantrajectory, 0); m)

"""
    $(SIGNATURES)

Starting from the environment's `state`, perform one step of the MPPI algorithm and
store the resulting action in `action`. The trajectory sampling portion of MPPI is
done in parallel using `nthreads` threads.
"""
@propagate_inbounds function LyceumBase.getaction!(
    action::AbstractVector,
    state,
    m::MPPI{DT,nu};
    nthreads::Integer = Threads.nthreads(),
) where {DT,nu}
    @boundscheck begin
        length(action) == nu || throw(ArgumentError("Expected action vector of length $nu. Got: $(length(action))"))
        require_one_based_indexing(action)
    end

    nthreads = min(m.K, nthreads)
    step!(m, state, nthreads)
    @inbounds copyto!(action, uview(m.meantrajectory, :, 1))
    shiftcontrols!(m)
    return action
end

function LyceumBase.step!(m::MPPI{DT,nu}, s, nthreads) where {DT,nu}
    randn!(m.noise)
    lmul!(m.covar_ul, reshape(m.noise, (nu, :)))
    if nthreads == 1
        # short circuit
        threadstep!(m, s, 1:m.K)
    else
        kranges = splitrange(m.K, nthreads)
        @sync for i = 1:nthreads
            Threads.@spawn threadstep!(m, s, kranges[i])
        end
    end
    combinetrajectories!(m)
end

function threadstep!(m::MPPI, s, krange)
    tid = Threads.threadid()
    for k in krange
        perturbedrollout!(m, s, k, tid)
    end
end

function perturbedrollout!(m::MPPI{DT,nu}, state, k, tid) where {DT,nu}
    env = m.envs[tid]
    obsbuf = m.obsbuffers[tid]
    statebuf = m.statebuffers[tid]
    mean = m.meantrajectory
    noise = m.noise

    setstate!(env, state)
    discountedreward = zero(DT)
    discountfactor = one(DT)
    @uviews mean noise @inbounds for t = 1:m.H
        mean_t = SVector{nu,DT}(view(mean, :, t))
        noise_tk = SVector{nu,DT}(view(noise, :, t, k))
        action_t = mean_t + noise_tk
        setaction!(env, action_t)

        step!(env)

        getobservation!(obsbuf, env)
        getstate!(statebuf, env)
        reward = getreward(statebuf, action_t, obsbuf, env)

        discountedreward += reward * discountfactor
        discountfactor *= m.gamma
    end # env at t=H+1
    getobservation!(obsbuf, env)
    terminalvalue = convert(DT, m.value(obsbuf))
    @inbounds m.trajectorycosts[k] = -(discountedreward + terminalvalue * discountfactor)
    return m
end

function combinetrajectories!(m::MPPI{DT,nu}) where {DT,nu}
    costs = m.trajectorycosts

    beta = minimum(costs)
    eta = zero(DT)
    for k = 1:m.K
        @inbounds costs[k] = softcost = exp((beta - costs[k]) / m.lambda)
        eta += softcost
    end

    costs ./= eta

    for k = 1:m.K, t = 1:m.H, u = 1:nu
        @inbounds m.meantrajectory[u, t] += costs[k] * m.noise[u, t, k]
    end

    m
end

function shiftcontrols!(m::MPPI{DT,nu}) where {DT,nu}
    for t = 2:m.H, u = 1:nu
        @inbounds m.meantrajectory[u, t-1] = m.meantrajectory[u, t]
    end
    m.initfn!(m.meantrajectory)
    m
end

@inline function default_initfn!(meantraj)
    @uviews meantraj @inbounds begin
        lastcontrols = view(meantraj, :, size(meantraj, 2))
        fill!(lastcontrols, 0)
    end
end


text:
# ---
# title: 564. Find the Closest Palindrome
# id: problem564
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/find-the-closest-palindrome/description/>
# hidden: true
# ---
#
# Given an integer n, find the closest integer (not including itself), which is
# a palindrome.
#
# The 'closest' is defined as absolute difference minimized between two
# integers.
#
# **Example 1:**
#
#
#
#     Input: "123"
#     Output: "121"
#
#
# **Note:**
#
#   1. The input **n** is a positive integer represented by string, whose length will not exceed 18.
#   2. If there is a tie, return the smaller one as answer.
#
#
## @lc code=start
using LeetCode

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


text: println("Fuck You Github")

max_stars_repo_name: ozcanyarimdunya/FuckYouGithub

max_issues_repo_name: ozcanyarimdunya/FuckYouGithub

max_forks_repo_name: ozcanyarimdunya/FuckYouGithub


text:
function wc(io, counter=Dict{String,Int}())
    for line in eachline(io), word in split(line)
        lword = lowercase(word)
        counter[lword] = get(counter, lword, 0) + 1
    end
    counter
end

sortdict(counter) = sort!(collect(counter), by=x->x[2], rev=true)

wc(stdin) |> sortdict .|> kv->println(kv[1], " ", kv[2])

max_stars_repo_name: smazurchuk/countwords

max_issues_repo_name: smazurchuk/countwords

max_forks_repo_name: smazurchuk/countwords


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

using Markdown
using InteractiveUtils

# ╔═╡ 22b67cb0-6687-11ec-1868-bb216a9703f4
begin
	using CondaPkg, PlutoUI
	CondaPkg.add.(("matplotlib", "numpy")); CondaPkg.resolve()
    using PythonCall
    import MarkdownLiteral: @mdx
end

# ╔═╡ f150770d-cbcd-4eda-a0a5-b759a21f5b9b
@mdx """
# Fun with 🐍

In this notebook we will showing some useage examples for the new package `PythonCall.jl`. It behaves similarly to the older package `PyCall.jl`, but with a few key differences that make interacting with Python quite nice:

* automatic, project-specific environments for easy reproducibility
* package management via [micromamba](https://mamba.readthedocs.io/en/latest/user_guide/micromamba.html) for blazing fast Python package management 🔥
* easy plot display support

Let's try some things out!
"""

# ╔═╡ 800e21bc-28e3-496d-8871-4c81de6713ca
TableOfContents()

# ╔═╡ 30181f53-97e2-4d21-83d2-9394ef64aba2
@mdx """
## `pyeval`
To see how we can start using this package, let's start by creating some simple Python objects:
"""

# ╔═╡ b1f1706d-da3b-48ae-bed1-4f5fb9e95726
py_list = @pyeval "[1, 2, 3]"

# ╔═╡ fa184670-718c-4e68-8e23-fe6e72a651c5
py_dict = @pyeval "{'one': 1, 'two': 2, 'three': 3}"

# ╔═╡ 53b89008-8ac0-4f20-955d-1d875530f121
py_range = @pyeval "range(0, 10, 2)"

# ╔═╡ 81a0a61d-0ff8-4610-915b-3f32a2742f5d
@mdx """
The `@pyeval` macro behaves like Python's [`eval`](https://docs.python.org/3/library/functions.html#eval) and evaluates whatever Python expression (written as a string literal or `cmd`) is passed to it. The returned objects can participate in the usual Python methods now:
"""

# ╔═╡ 15b2d0ab-bc1a-435a-9557-99537d117d91
py_list.append(4); py_list

# ╔═╡ ffa07e57-7797-4000-8be4-a166895bcc54
py_dict["three"]

# ╔═╡ 5e0919d8-0a02-46fc-92fe-0a5e40d1708c
py_range[0] # Python's zero-based indexing is automatically understood

# ╔═╡ 0b3322b5-cd15-49d9-a446-35224a4b7637
@mdx """
!!! tip
	We recommend passing strings ("") instead of cmds (\`\`) to play nice with Pluto's `ExpressionExplorer`

	Related issue? [#961](https://github.com/fonsp/Pluto.jl/issues/961)
"""

# ╔═╡ 3359adff-8eb5-4911-a914-047ecf3663fe
@mdx """
`PythonCall.jl` also provides a function version of `@pyeval` if we would like to do things like string interpolation first:
"""

# ╔═╡ 6be1822a-d905-4de9-8345-4e160e53ea64
let
	four = 4
	pyeval("[1, 2, 3, $four]", Main)
end

# ╔═╡ ed066bda-07b4-4c52-9302-15c66cb4e1d8
@mdx """
More info about the usage for each version is available in the Live docs.
"""

# ╔═╡ a7e24a65-cfa8-4e03-9664-f8f424324917
@mdx """
To a certain extent, Julia's functions can also operate on these objects automatically:
"""

# ╔═╡ 1e94e185-af8c-4796-9358-6a2c15c7fd43
py_range |> sum

# ╔═╡ b810a02b-8329-4f51-8cfc-c6920ad1cf5d
py_range |> collect

# ╔═╡ 40cb6860-7ac6-40ca-8b1f-30360fd6aef7
@mdx """
As an alternative, `PythonCall` makes it very easy to convert these objects to their native Julia counterparts with the `@pyconvert`/`pyconvert` macro/function:
"""

# ╔═╡ 57686014-99e8-4733-8f1d-ba80306b9da9
(@pyconvert StepRange py_range) |> collect

# ╔═╡ 1360029e-76f3-468a-9b2f-d27482bbc525
pyconvert(Dict, py_dict)

# ╔═╡ fde4db08-e5e1-4e92-b1c7-9a41edb04476
@mdx """
Next, we will take a look at running Python statements.
"""

# ╔═╡ 71ba808a-ebde-43ad-b264-ccf8e676891a
@mdx """
## `pyexec`

We saw how to run simple Python expressions with `pyeval`. Now we will turn to executing Python statements, which can be used to perform more complex tasks like storing values in variables and defining functions. This is accomplished with `PythonCalls.jl's` `@pyexec`/`pyexec` macro/function, which behave's like Python's [`exec`](https://docs.python.org/3/library/functions.html#exec):
"""

# ╔═╡ 193e65ae-c344-4dda-b92b-e0b136eb7581
@pyexec (x=1, y=2) => "ans = x + y" => ans::Int

# ╔═╡ 0c8e14ef-782d-420b-8906-3a139bacf331
pyexec(@NamedTuple{ans}, "ans = 1 + 2", Main)

# ╔═╡ 1eccfc11-4ea1-4b79-a688-2644d6d1d0fd
@mdx """
## Combining `pyexec` and `pyeval`

We can now compose these ideas to start interacting with full blocks of Python code:
"""

# ╔═╡ 5b597d2a-2483-4b80-860b-839fa3ddeaec
begin
	@pyexec """
	global greetings
	def greetings(name):
		return f"Hi {name} 👋"
	"""
	greetings(name) = @pyeval("greetings")(name)
end

# ╔═╡ 5aa1ba97-d55c-4794-a839-e90effb84bbe
greetings("Pluto citizen")

# ╔═╡ 53dd2eb5-ee4c-4f4c-a9ca-6f234a376a78
@mdx """
!!! note
	By default, this is executed in local scope. To make it accessible to other parts of the notebook, it is required to explicitly state this with `global`.
"""

# ╔═╡ e9e6c4e0-6834-4b6b-ac20-ff722f9a5cd9
@mdx """
## Using packages
`PythonCall.jl` has a nice companion package named `CondaPkg.jl`, which we can use to easily install Python packages into an environment in the same directory as this notebook:
"""

# ╔═╡ 0bf45621-7776-403e-b3da-5311a5c30e20
@mdx """
```julia
CondaPkg.add.(("matplotlib", "numpy"))
CondaPkg.resolve()
```
"""

# ╔═╡ f5d2228c-b45e-4dce-99fb-668f6c50df30
@mdx """
We now use the `@py` macro to import and interact with these packages:
"""

# ╔═╡ a159c36c-83c8-460e-a7e2-e85c7df8d9da
@py begin
	import matplotlib.pyplot as plt
	import numpy as np
end

# ╔═╡ c83fdb79-b0d9-4830-b6cf-74d2a669ceed
let
	xs = np.random.rand(10, 4, 4)

	fig, axes = plt.subplots(2, 2, sharex=true, sharey=true)

	for (ax, x) ∈ zip(axes.flat, xs)
		ax.plot(x)
	end

	fig.tight_layout()

	plt.gcf()
end

# ╔═╡ 681ece7a-6800-4779-8118-28c3179cd43a
@mdx """
## Environment details

Where is all this stuff being downloaded/run? Let's see!
"""

# ╔═╡ 655674de-56c1-4386-8fda-aa5c95b6271f
@with_terminal CondaPkg.status() # Will show directly in notebook in next Pluto release

# ╔═╡ d407d19b-40ee-496b-b585-ff685232b48b
@with_terminal run(CondaPkg.MicroMamba.cmd(`list`)) # Will ask for a better API for this

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
MarkdownLiteral = "736d6165-7244-6769-4267-6b50796e6954"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"

[compat]
CondaPkg = "~0.2.6"
MarkdownLiteral = "~0.1.1"
PlutoUI = "~0.7.32"
PythonCall = "~0.5.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "4aff51293dbdbd268df314827b7f409ea57f5b70"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.5"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CondaPkg]]
deps = ["JSON3", "Markdown", "MicroMamba", "Pkg", "TOML"]
git-tree-sha1 = "49589b0a9284fd8ea84c64a42e96739d1fe06751"
uuid = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
version = "0.2.6"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MarkdownLiteral]]
deps = ["CommonMark", "HypertextLiteral"]
git-tree-sha1 = "0d3fa2dd374934b62ee16a4721fe68c418b92899"
uuid = "736d6165-7244-6769-4267-6b50796e6954"
version = "0.1.1"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MicroMamba]]
deps = ["Pkg", "Scratch"]
git-tree-sha1 = "e61632c723f4d0fafa5419a506f08baafe973391"
uuid = "0b3b1443-0f03-428d-bdfb-f27f9c1191ea"
version = "0.1.4"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.PythonCall]]
deps = ["CondaPkg", "Dates", "Libdl", "MacroTools", "Markdown", "Pkg", "Requires", "Serialization", "Tables", "UnsafePointers"]
git-tree-sha1 = "402a4d004a9e7c9656758f479f456205f7f6a3ce"
uuid = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
version = "0.5.1"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.UnsafePointers]]
git-tree-sha1 = "c81331b3b2e60a982be57c046ec91f599ede674a"
uuid = "e17b2a0c-0bdf-430a-bd0c-3a23cae4ff39"
version = "1.0.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─f150770d-cbcd-4eda-a0a5-b759a21f5b9b
# ╠═800e21bc-28e3-496d-8871-4c81de6713ca
# ╟─30181f53-97e2-4d21-83d2-9394ef64aba2
# ╠═b1f1706d-da3b-48ae-bed1-4f5fb9e95726
# ╠═fa184670-718c-4e68-8e23-fe6e72a651c5
# ╠═53b89008-8ac0-4f20-955d-1d875530f121
# ╟─81a0a61d-0ff8-4610-915b-3f32a2742f5d
# ╠═15b2d0ab-bc1a-435a-9557-99537d117d91
# ╠═ffa07e57-7797-4000-8be4-a166895bcc54
# ╠═5e0919d8-0a02-46fc-92fe-0a5e40d1708c
# ╟─0b3322b5-cd15-49d9-a446-35224a4b7637
# ╟─3359adff-8eb5-4911-a914-047ecf3663fe
# ╠═6be1822a-d905-4de9-8345-4e160e53ea64
# ╟─ed066bda-07b4-4c52-9302-15c66cb4e1d8
# ╟─a7e24a65-cfa8-4e03-9664-f8f424324917
# ╠═1e94e185-af8c-4796-9358-6a2c15c7fd43
# ╠═b810a02b-8329-4f51-8cfc-c6920ad1cf5d
# ╟─40cb6860-7ac6-40ca-8b1f-30360fd6aef7
# ╠═57686014-99e8-4733-8f1d-ba80306b9da9
# ╠═1360029e-76f3-468a-9b2f-d27482bbc525
# ╟─fde4db08-e5e1-4e92-b1c7-9a41edb04476
# ╟─71ba808a-ebde-43ad-b264-ccf8e676891a
# ╠═193e65ae-c344-4dda-b92b-e0b136eb7581
# ╠═0c8e14ef-782d-420b-8906-3a139bacf331
# ╟─1eccfc11-4ea1-4b79-a688-2644d6d1d0fd
# ╠═5b597d2a-2483-4b80-860b-839fa3ddeaec
# ╠═5aa1ba97-d55c-4794-a839-e90effb84bbe
# ╟─53dd2eb5-ee4c-4f4c-a9ca-6f234a376a78
# ╟─e9e6c4e0-6834-4b6b-ac20-ff722f9a5cd9
# ╟─0bf45621-7776-403e-b3da-5311a5c30e20
# ╟─f5d2228c-b45e-4dce-99fb-668f6c50df30
# ╠═a159c36c-83c8-460e-a7e2-e85c7df8d9da
# ╠═c83fdb79-b0d9-4830-b6cf-74d2a669ceed
# ╟─681ece7a-6800-4779-8118-28c3179cd43a
# ╠═655674de-56c1-4386-8fda-aa5c95b6271f
# ╠═d407d19b-40ee-496b-b585-ff685232b48b
# ╠═22b67cb0-6687-11ec-1868-bb216a9703f4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
export MultiStart, n_points, n_params, optimize

"""
	$(SIGNATURES)

Contains information to run a multistart optimization.

# Fields

- optFct:
    Function to be optimized. Signature:
    `solnV, fVal, exitFlag = optFct(guessV)`
- guessesV:
    Vector of guesses. Each guess is a Vector{Float64}
- historyFn:
    Full path where optimization history will be stored.
- fTol, fTolN:
    Stop if fVal improves by less than `fTol` over the last `fTolN` points
    This can terminate early.
"""
mutable struct MultiStart{F, U}
    optFct :: F
    updateRule :: U
    guessesV :: Vector{Vector{Float64}}
    historyFn :: String
    maxHours :: Float64
    fTol :: Float64
    fTolN :: Integer
end

"""
	$(SIGNATURES)

Number of points to be evaluated.
"""
n_points(m :: MultiStart{F, U}) where {F, U} = length(m.guessesV);

"""
	$(SIGNATURES)

Number of calibrated parameters in each point.
"""
n_params(m :: MultiStart{F, U}) where {F, U} =
    length(first(m.guessesV));

"""
	$(SIGNATURES)

Run the multistart optimization.
"""
function optimize(m :: MultiStart{F, U},
        io = stdout) where {F, U}

    log_msg(1, io, start_header(m));
    startTime = Dates.now();
    h = init_history(m.historyFn, m.guessesV);

    done = false;
    j = 0;
    while !done
        j += 1;
        newGuessV = new_guess(m.updateRule, h, j);
        log_msg(2, io, point_start_header(j));
        # startTime = Dates.now();
        newPoint = m.optFct(newGuessV);
        # endTime = Dates.now();
        log_msg(2, io, finish_point_header(j, newPoint));
        add_to_history!(h, j, newPoint);
        save_history(h);
        continue_run(m, h, j, startTime; io)  ||  (done = true);
    end
    bestPoint = best_point(h);
    log_msg(1, io, finish_header(bestPoint, startTime, n_points(m)));
    return bestPoint.solnV, bestPoint.fVal, bestPoint.exitFlag
end

"""
	$(SIGNATURES)

Decide whether to continue running at the end of iteration `j`.
Report reason for termination.
"""
function continue_run(m :: MultiStart{F,U}, h :: History,
            j, startTime; io = stdout) where {F,U}
    cont = true;
    if reached_max_time(startTime, m.maxHours)
        log_msg(3, io, "Max time reached.");
        cont = false;
    end
    if j >= n_points(m)
        log_msg(3, io, "All points evaluated.");
        cont = false;
    end
    if j > m.fTolN
        fDiff = fval_improvement(h, j, m.fTolN);
        if fDiff < m.fTol
            fDiffStr = round(fDiff; digits = 2);
            log_msg(3, io, "Difference in fVal below fTol: $fDiffStr < $(m.fTol)");
            cont = false;
        end
    end
    return cont
end

# function fval_improvement_too_small(m :: MultiStart{F,U}, h :: History,
#         j) where {F,U}
#     if j > m.
#         return fval_improvement(h, j, n) < m.fTol;
#     else
#         return false
#     end
# end

# Improvement during last `n` points relative to best point of previous history.
function fval_improvement(h :: History, j :: Integer, n :: Integer)
    @assert (j > n)  "Do not have past $n points with history of length $j";
    fValV = f_values(h);
    fPast = minimum(fValV[1 : (j-n)]);
    fNew = minimum(fValV[(j-n+1) : j]);
    fDiff = fPast - fNew;
    return fDiff
end

reached_max_time(startTime, maxHours :: Float64) =
    seconds_elapsed(startTime, Dates.now()) / 3600 >= maxHours;

seconds_elapsed(startTime, endTime) = Dates.value(endTime - startTime) / 1000;

# +++++ move to general purpose pkg
function split_duration(dSeconds :: Integer)
    dHours, hRem = divrem(dSeconds, 3600);
    dMin, mRem = divrem(hRem, 60);
    dSec = round(Int, mRem);
    return dHours, dMin, dSec
end

function format_duration(dSeconds :: Integer)
    dHours, dMin, dSec = split_duration(dSeconds);
    d = DateTime(2020,1,1,dHours,dMin,dSec);
    if dHours > 0
        return Dates.format(d, "HH:MM") * " hr";
    elseif dMin > 0
        return Dates.format(d, "MM:SS") * " min";
    else
        return Dates.format(d, "SS") * " sec";
    end
end

## -----------  Headers to display

start_header(m :: MultiStart{F, U}) where {F,U} =
    "Multistart with $(n_points(m)) points begins at $(format_time())";

point_start_header(j :: Integer) =
    "Starting point $j at $(format_time())";

finish_point_header(j :: Integer, p :: Point{F}) where F =
    "Point $j completed at $(format_time()) with [fVal $(round(p.fVal, digits = 3))]  and  [exitFlag $(p.exitFlag)]";

function finish_header(bestPoint :: Point{F}, startTime, nPoints) where F
    avgTime = round(Int, average_time(startTime, Dates.now(), nPoints));
    avgTimeStr = format_duration(avgTime);
    return "MultiStart completed at $(format_time()) \n  fVal $(round(bestPoint.fVal, digits = 3)),  exitFlag $(bestPoint.exitFlag) \n  Average time per point: $(avgTimeStr)."
end

format_time(d = Dates.now()) = Dates.format(d, "u-d, HH:MM");

average_time(dStart, dStop, n) = seconds_elapsed(dStart, dStop) / n;

# -------------


text:
"""
    cauchyendgame(H::AbstractHomotopy, J_H, ∂H∂t, x, R, algorithm, kwargs..., pathtrackingkwargs...)

Execute a cauchy endgame for `H` starting at radius `R`. This assumes that ``H(x,t)=0``.
It is based on "A Parallel Endgame " by Bates, Hauenstein and Sommese. [^1]

## Optional arguments
* `prediction_tolerance=1e-4`:
* `geometric_series_factor=0.5`:
* `endgame_tolerance=1e-8`:
* `samples_per_loop=8`:
* `max_winding_number=8`:
* `loopclosed_tolerance=1e-3`
*  optional arguments to [`pathtracking`](@ref)

[^1]: Bates, Daniel J., Jonathan D. Hauenstein, and Andrew J. Sommese. "A Parallel Endgame." Contemp. Math 556 (2011): 25-35.
"""
function cauchyendgame(H::AbstractHomotopy{T}, J_H, ∂H∂t, x, t::Float64, algorithm::APCA{true};
    prediction_tolerance=1e-8,
    geometric_series_factor=0.5,
    endgame_tolerance=1e-8,
    samples_per_loop::Int=8,
    max_winding_number::Int=8,
    loopclosed_tolerance=1e-3,
    tolerance_infinity=1e-8,
    pathtrackingkwargs...
) where {T}
    λ = geometric_series_factor
    R = t
    steps = [R]
    xs = [x]
    lastprediction = Nullable(x, false)
    endgame_started = false
    k = 1
    while R > eps(Float64)
        res = trackpath(H, J_H, ∂H∂t, last(xs), algorithm, R, λ*R; pathtrackingkwargs...)
        if nosuccess(res)
            return CauchyEndgameResult(get(lastprediction, last(xs)), :IllConditionedZone, k, R, λ, xs, steps, Nullable{ConvergentCluster{T}}())
        end
        R = λ*R
        push!(steps, R)
        push!(xs, res.result)

        # we need at least 4 points to start with our first heuristic
        if length(xs) < 4
            continue
        end

        if atinfinity(get(lastprediction, last(xs)), tolerance_infinity)
            return CauchyEndgameResult(get(lastprediction, last(xs)), :AtInfinity, k, R, λ, xs, steps, Nullable{ConvergentCluster{T}}())
        end

        if endgame_started || firstheuristic(R, λ, xs[end-3], xs[end-2], xs[end-1], xs[end])
            retcode, samples =
                loop(H, J_H, ∂H∂t, last(xs), R, algorithm, samples_per_loop, max_winding_number,
                    loopclosed_tolerance, endgame_tolerance, apply_heuristic=!endgame_started)
            if retcode == :Success
                endgame_started = true
                prediction = predict_with_cif(samples)
                if !isnull(lastprediction)
                    # check whether to predictions are close enough, so we can say that
                    # the endgame terminates
                    Δprediction = projectivenorm(prediction, get(lastprediction))
                    if Δprediction < prediction_tolerance
                        cluster = ConvergentCluster(samples[1:samples_per_loop:end], prediction, R)
                        return CauchyEndgameResult(prediction, :Success, k, R, λ, xs, steps, Nullable(cluster))
                    end
                end
                lastprediction = Nullable(prediction)
            end
        end

        k += 1
    end

    CauchyEndgameResult(get(lastprediction, last(xs)), :MachineEpsilon, k, R, λ, xs, steps, Nullable{ConvergentCluster{T}}())
end

atinfinity(x, tolerance_infinity) = norm(normalize(x)[1]) < tolerance_infinity

"""
    loop(H, J_H, ∂H∂t, x, radius, algorithm, samples_per_loop, max_winding_number, loopclosed_tolerance, apply_heuristic=true, pathtrackingkwargs...)

Tracks the implicit defined path z(t) around the `n`-gon with vertices
``r⋅exp(i2πk/n)`` where `n=samples_per_loop`.
"""
function loop(
    H::AbstractHomotopy{T},
    J_H,
    ∂H∂t,
    x,
    radius::Real,
    algorithm::APCA,
    samples_per_loop::Int,
    max_winding_number::Int,
    loopclosed_tolerance::Float64,
    endgame_tolerance::Float64;
    apply_heuristic=true,
    pathtrackingkwargs...
) where {T}
    unitroots = UnitRootsIterator(radius, samples_per_loop)
    samples = [x]
    start = first(unitroots)
    for (k, finish) in enumerate(Iterators.drop(unitroots, 1))
        pathresult = trackpath(H, J_H, ∂H∂t, samples[end], algorithm, start, finish; pathtrackingkwargs...)
        if pathresult.returncode != :Success
            return (:BadBranch, samples)
        end

        # apply the second heuristic to tell us whether it makes sense to start the endgame
        if apply_heuristic && k == samples_per_loop &&
            !secondheuristic(radius, samples, β=endgame_tolerance)
            return (:HeuristicFailed, samples)
        end

        # check whether the loop is closed. This can only happen if we finished one circle.
        if k % samples_per_loop == 0
            loopnorm = projectivenorm(first(samples), pathresult.result)
            if loopnorm < loopclosed_tolerance
                break
            end
        end

        if k > max_winding_number * samples_per_loop
            return (:WindingNumberTooHigh, samples)
        end

        push!(samples, pathresult.result)
        start = finish
    end
    (:Success, samples)
end

"""
    firstheuristic(t, x_R, x_λR, x_λ2R, x_λ3R; L = 0.75)

This is the first heuristic when to start the endgame.
It checks whether two consectuive appromixations of the first nonvanishing
monomial of the Puiseux series expansion of x(t) at 0 agree.
For more details see "A Parallel Endgame " by Bates, Hauenstein and Sommese [^1],
page 8 and 9.

[^1]: Bates, Daniel J., Jonathan D. Hauenstein, and Andrew J. Sommese. "A Parallel Endgame." Contemp. Math 556 (2011): 25-35.
"""
function firstheuristic(t, λ, x_R, x_λR, x_λ2R, x_λ3R; L = 0.75)
    # failsafe mechanism
    if t < 1e-8
        return true
    end
    v = rand(Complex128, length(x_R))
    g(x_R, x_λR, x_λ2R) = log(abs2((v⋅x_λR - v⋅x_λ2R) / (v⋅x_R - v⋅x_λR))) / log(λ)

    g_R = g(x_R, x_λR, x_λ2R)
    g_λR = g(x_λR, x_λ2R, x_λ3R)

    if g_R <= 0
        return false
    end

    L < g_R / g_λR < 1 / L
end

"""
    secondheuristic(t, samplepoints, β; K=0.5)

This is the second heuristic when to start the endgame.
It enforces that the values around the lopp do not differ radically.
This functions assumes that `samplepoints` has been collected on ``x(Re^{iθ})`` for ``θ∈[0,2π]``.
See "A Parallel Endgame " by Bates, Hauenstein and Sommese [^1], page 9 for more details.

[^1]: Bates, Daniel J., Jonathan D. Hauenstein, and Andrew J. Sommese. "A Parallel Endgame." Contemp. Math 556 (2011): 25-35.
"""
function secondheuristic(radius, samplepoints; β=1e-8, K=0.5)
    # failsafe
    if radius < 1e-14
        true
    end

    norms = norm.(samplepoints)
    m = minimum(norms)
    M = maximum(norms)
    (M - m < β) || (m / M > K)
end

"""
    projectivenorm(a, b)

Calculate |a-b|, but first bring `a` and `b` on the same patch.
Brings fist `a` and `b` on the same patch by finding diving `a` through it's maximum value
(w.r.t. to the absolute value) with index `i` and diving `b` through `b[i]`.
Then computes the norm of the differences.
"""
function projectivenorm(a::Vector{<:Number}, b::Vector{<:Number})
    _, i = findmax(abs2.(a))
    norm(a ./ a[i] .- b ./ b[i])
end

"""
    predict_with_cif(samplepoints)

Predict the value of a function based on the samplepoints using the cauchy integral formula.
This function assumes that the samplepoints are sampled uniformly from a circle around 0,
since then the integral (approximated with the trapezoidal rule) is simply the mean of
all samplepoints.
"""
predict_with_cif(samplepoints) = mean(samplepoints)

"""
    UnitRootsIterator(r, n)

Construct an infinite iterator which returns the `n`-th scalded unit roots, i.e.
the values ``r⋅exp(i2πk/n)`` for ``k=0,1,...``.
"""
struct UnitRootsIterator
    radius::Float64
    order::Float64
end
UnitRootsIterator(r::Real, order::Real) = UnitRootsIterator(float(r), float(order))
Base.start(::UnitRootsIterator) = 0
Base.next(loop::UnitRootsIterator, k::Int) = (loop.radius *  exp(im * 2π * k / loop.order), k + 1)
Base.done(::UnitRootsIterator, ::Int) = false
Base.iteratorsize(::UnitRootsIterator) = Base.IsInfinite()
Base.iteratoreltype(::UnitRootsIterator) = Base.HasEltype()
Base.eltype(::UnitRootsIterator) = Complex128


text:
using Documenter, NFFT

# Updating doctests:
#
# julia> using Documenter, NFFT
# julia> DocMeta.setdocmeta!(NFFT, :DocTestSetup, :(using NFFT); recursive=true)
# julia> doctest(NFFT, fix=true)

DocMeta.setdocmeta!(NFFT, :DocTestSetup, :(using NFFT); recursive=true)

makedocs(;
    doctest = true,
    #strict = :doctest,
    modules = [NFFT],
    sitename = "NFFT",
    authors = "Tobias Knopp and contributors",
    repo="https://github.com/JuliaMath/NFFT.jl/blob/{commit}{path}#{line}",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tknopp.github.io/NFFT.jl",
        assets=String[],
    ),
    pages = [
        "Home" => "index.md",
        "Background" => "background.md",
        "Overview" => "overview.md",
        "AbstractNFFTs" => "abstract.md",
        "Performance" => "performance.md",
        "Tools" => "tools.md",
        #"Implementation" => "implementation.md",
        "API" => "api.md",
    ]
)

deploydocs(repo   = "github.com/JuliaMath/NFFT.jl.git")


text:
using ForwardDiff
using Rotations
using StaticArrays
using LinearAlgebra

using ConstrainedDynamics
using ConstrainedDynamics: vrotate, Lmat, Vᵀmat, Lᵀmat, VLᵀmat, ∂g∂pos, ∂g∂vel, skew

function transfunc3pos(vars)
    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    pa = vars[15:17]
    pb = vars[18:20]

    (xb + vrotate(SVector{3}(pb...), qb)) - (xa + vrotate(SVector{3}(pa), qa))
end

function transfunc3vel(vars)
    Δt = 0.01

    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    va = vars[15:17]
    wa = vars[18:20]
    vb = vars[21:23]
    wb = vars[24:26]

    pa = vars[27:29]
    pb = vars[30:32]

    xa = xa + Δt * va
    xb = xb + Δt * vb
    qa = Quaternion(Δt / 2 * (qa * Quaternion(SVector(sqrt((2 / Δt)^2 - wa' * wa), wa...))))
    qb = Quaternion(Δt / 2 * (qb * Quaternion(SVector(sqrt((2 / Δt)^2 - wb' * wb), wb...))))

    (xb + vrotate(SVector{3}(pb...), qb)) - (xa + vrotate(SVector{3}(pa...), qa))
end

function transfunc2pos(vars)
    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    pa = vars[15:17]
    pb = vars[18:20]

    V1 = vars[21:23]
    V2 = vars[24:26]
    V12 = [V1';V2']

    V12 * vrotate(SVector{3}(((xb + vrotate(SVector{3}(pb...), qb)) - (xa + vrotate(SVector{3}(pa...), qa)))...), inv(qa))
end

function transfunc2vel(vars)
    Δt = 0.01

    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    va = vars[15:17]
    wa = vars[18:20]
    vb = vars[21:23]
    wb = vars[24:26]

    pa = vars[27:29]
    pb = vars[30:32]

    V1 = vars[33:35]
    V2 = vars[36:38]
    V12 = [V1';V2']

    xa = xa + Δt * va
    xb = xb + Δt * vb
    qa = Quaternion(Δt / 2 * (qa * Quaternion(SVector(sqrt((2 / Δt)^2 - wa' * wa), wa...))))
    qb = Quaternion(Δt / 2 * (qb * Quaternion(SVector(sqrt((2 / Δt)^2 - wb' * wb), wb...))))

    V12 * vrotate(SVector{3}(((xb + vrotate(SVector{3}(pb...), qb)) - (xa + vrotate(SVector{3}(pa...), qa)))...), inv(qa))
end

function transfunc1pos(vars)
    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    pa = vars[15:17]
    pb = vars[18:20]

    v = vars[21:23]

    v' * vrotate(SVector{3}(((xb + vrotate(SVector{3}(pb...), qb)) - (xa + vrotate(SVector{3}(pa...), qa)))...), inv(qa))
end

function transfunc1vel(vars)
    Δt = 0.01

    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    va = vars[15:17]
    wa = vars[18:20]
    vb = vars[21:23]
    wb = vars[24:26]

    pa = vars[27:29]
    pb = vars[30:32]

    v = vars[33:35]

    xa = xa + Δt * va
    xb = xb + Δt * vb
    qa = Quaternion(Δt / 2 * (qa * Quaternion(SVector(sqrt((2 / Δt)^2 - wa' * wa), wa...))))
    qb = Quaternion(Δt / 2 * (qb * Quaternion(SVector(sqrt((2 / Δt)^2 - wb' * wb), wb...))))

    v' * vrotate(SVector{3}(((xb + vrotate(SVector{3}(pb...), qb)) - (xa + vrotate(SVector{3}(pa...), qa)))...), inv(qa))
end

function rotfunc3pos(vars)
    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    offset = Quaternion(SVector(vars[15:18]...))

    VLᵀmat(offset) * Lᵀmat(qa) * qb
end

function rotfunc3vel(vars)
    Δt = 0.01

    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    va = vars[15:17]
    wa = vars[18:20]
    vb = vars[21:23]
    wb = vars[24:26]

    offset = Quaternion(SVector(vars[27:30]...))

    xa = xa + Δt * va
    xb = xb + Δt * vb
    qa = Quaternion(Δt / 2 * (qa * Quaternion(SVector(sqrt((2 / Δt)^2 - wa' * wa), wa...))))
    qb = Quaternion(Δt / 2 * (qb * Quaternion(SVector(sqrt((2 / Δt)^2 - wb' * wb), wb...))))

    VLᵀmat(offset) * Lᵀmat(qa) * qb
end

function rotfunc2pos(vars)
    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    offset = Quaternion(SVector(vars[15:18]...))

    V1 = vars[19:21]
    V2 = vars[22:24]
    V12 = [V1';V2']

    V12 * (VLᵀmat(offset) * Lᵀmat(qa) * qb)
end

function rotfunc2vel(vars)
    Δt = 0.01

    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    va = vars[15:17]
    wa = vars[18:20]
    vb = vars[21:23]
    wb = vars[24:26]

    offset = Quaternion(SVector(vars[27:30]...))

    V1 = vars[31:33]
    V2 = vars[34:36]
    V12 = [V1';V2']

    xa = xa + Δt * va
    xb = xb + Δt * vb
    qa = Quaternion(Δt / 2 * (qa * Quaternion(SVector(sqrt((2 / Δt)^2 - wa' * wa), wa...))))
    qb = Quaternion(Δt / 2 * (qb * Quaternion(SVector(sqrt((2 / Δt)^2 - wb' * wb), wb...))))

    V12 * (VLᵀmat(offset) * Lᵀmat(qa) * qb)
end

function rotfunc1pos(vars)
    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    offset = Quaternion(SVector(vars[15:18]...))

    V3 = vars[19:21]

    V3' * (VLᵀmat(offset) * Lᵀmat(qa) * qb)
end

function rotfunc1vel(vars)
    Δt = 0.01

    xa = vars[1:3]
    qa = Quaternion(SVector(vars[4:7]...))
    xb = vars[8:10]
    qb = Quaternion(SVector(vars[11:14]...))

    va = vars[15:17]
    wa = vars[18:20]
    vb = vars[21:23]
    wb = vars[24:26]

    offset = Quaternion(SVector(vars[27:30]...))

    V3 = vars[31:33]

    xa = xa + Δt * va
    xb = xb + Δt * vb
    qa = Quaternion(Δt / 2 * (qa * Quaternion(SVector(sqrt((2 / Δt)^2 - wa' * wa), wa...))))
    qb = Quaternion(Δt / 2 * (qb * Quaternion(SVector(sqrt((2 / Δt)^2 - wb' * wb), wb...))))

    V3' * (VLᵀmat(offset) * Lᵀmat(qa) * qb)
end

function transtest3()
    xa = rand(3)
    qa = Quaternion(rand(RotMatrix{3}))
    xb = rand(3)
    qb = Quaternion(rand(RotMatrix{3}))

    va = rand(3)
    wa = rand(3)
    vb = rand(3)
    wb = rand(3)

    pa = rand(3)
    pb = rand(3)

    origin = Origin{Float64}()
    link1 = Body(Box(1., 1., 1., 1.))
    link2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(OriginConnection(origin, link1))
    oc2 = EqualityConstraint(OriginConnection(origin, link2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Translational3(link1, link2, pa, pb))

    bot = Mechanism(origin, [link1;link2], [oc1;oc2;joint1])

    setPosition!(bot, link1, x = xa, q = qa)
    setPosition!(bot, link2, x = xb, q = qb)
    link1.s1 = SVector([va;wa]...)
    link2.s1 = SVector([vb;wb]...)

    res = ForwardDiff.jacobian(transfunc3pos, [xa;qa;xb;qb;pa;pb])
    X1 = res[1:3,1:3]
    Q1 = res[1:3,4:7] * Lmat(Quaternion(qa)) * Vᵀmat()
    X2 = res[1:3,8:10]
    Q2 = res[1:3,11:14] * Lmat(Quaternion(qb)) * Vᵀmat()

    n11 = norm(X1 - ∂g∂pos(joint1, 1, bot)[1:3,1:3])
    n12 = norm(Q1 - ∂g∂pos(joint1, 1, bot)[1:3,4:6])
    n21 = norm(X2 - ∂g∂pos(joint1, 2, bot)[1:3,1:3])
    n22 = norm(Q2 - ∂g∂pos(joint1, 2, bot)[1:3,4:6])

    res = ForwardDiff.jacobian(transfunc3vel, [xa;qa;xb;qb;va;wa;vb;wb;pa;pb])
    V1 = res[1:3,15:17]
    W1 = res[1:3,18:20]
    V2 = res[1:3,21:23]
    W2 = res[1:3,24:26]

    n31 = norm(V1 - ∂g∂vel(joint1, 1, bot)[1:3,1:3])
    n32 = norm(W1 - ∂g∂vel(joint1, 1, bot)[1:3,4:6])
    n41 = norm(V2 - ∂g∂vel(joint1, 2, bot)[1:3,1:3])
    n42 = norm(W2 - ∂g∂vel(joint1, 2, bot)[1:3,4:6])

    display((n11, n12, n21, n22))
    display((n31, n32, n41, n42))
    return
end

function transtest2()
    xa = rand(3)
    qa = Quaternion(rand(RotMatrix{3}))
    xb = rand(3)
    qb = Quaternion(rand(RotMatrix{3}))

    va = rand(3)
    wa = rand(3)
    vb = rand(3)
    wb = rand(3)

    pa = rand(3)
    pb = rand(3)

    v = rand(3)

    origin = Origin{Float64}()
    link1 = Body(Box(1., 1., 1., 1.))
    link2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(OriginConnection(origin, link1))
    oc2 = EqualityConstraint(OriginConnection(origin, link2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Translational2(link1, link2, pa, pb, v))
    V12 = joint1.constraints[1].V12

    bot = Mechanism(origin, [link1;link2], [oc1;oc2;joint1])

    setPosition!(bot, link1, x = xa, q = qa)
    setPosition!(bot, link2, x = xb, q = qb)
    link1.s1 = SVector([va;wa]...)
    link2.s1 = SVector([vb;wb]...)

    res = ForwardDiff.jacobian(transfunc2pos, [xa;qa;xb;qb;pa;pb;V12[1,:];V12[2,:]])
    X1 = res[1:2,1:3]
    Q1 = res[1:2,4:7] * Lmat(Quaternion(qa)) * Vᵀmat()
    X2 = res[1:2,8:10]
    Q2 = res[1:2,11:14] * Lmat(Quaternion(qb)) * Vᵀmat()

    n11 = norm(X1 - ∂g∂pos(joint1, 1, bot)[1:2,1:3])
    n12 = norm(Q1 - ∂g∂pos(joint1, 1, bot)[1:2,4:6])
    n21 = norm(X2 - ∂g∂pos(joint1, 2, bot)[1:2,1:3])
    n22 = norm(Q2 - ∂g∂pos(joint1, 2, bot)[1:2,4:6])

    res = ForwardDiff.jacobian(transfunc2vel, [xa;qa;xb;qb;va;wa;vb;wb;pa;pb;V12[1,:];V12[2,:]])
    V1 = res[1:2,15:17]
    W1 = res[1:2,18:20]
    V2 = res[1:2,21:23]
    W2 = res[1:2,24:26]

    n31 = norm(V1 - ∂g∂vel(joint1, 1, bot)[1:2,1:3])
    n32 = norm(W1 - ∂g∂vel(joint1, 1, bot)[1:2,4:6])
    n41 = norm(V2 - ∂g∂vel(joint1, 2, bot)[1:2,1:3])
    n42 = norm(W2 - ∂g∂vel(joint1, 2, bot)[1:2,4:6])

    display((n11, n12, n21, n22))
    display((n31, n32, n41, n42))
    return
end

function transtest1()
    xa = rand(3)
    qa = Quaternion(rand(RotMatrix{3}))
    xb = rand(3)
    qb = Quaternion(rand(RotMatrix{3}))

    va = rand(3)
    wa = rand(3)
    vb = rand(3)
    wb = rand(3)

    pa = rand(3)
    pb = rand(3)

    v = rand(3)

    origin = Origin{Float64}()
    link1 = Body(Box(1., 1., 1., 1.))
    link2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(OriginConnection(origin, link1))
    oc2 = EqualityConstraint(OriginConnection(origin, link2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Translational1(link1, link2, pa, pb, v))
    v = joint1.constraints[1].V3'

    bot = Mechanism(origin, [link1;link2], [oc1;oc2;joint1])

    setPosition!(bot, link1, x = xa, q = qa)
    setPosition!(bot, link2, x = xb, q = qb)
    link1.s1 = SVector([va;wa]...)
    link2.s1 = SVector([vb;wb]...)

    res = ForwardDiff.gradient(transfunc1pos, [xa;qa;xb;qb;pa;pb;v])
    X1 = res[1:3]'
    Q1 = res[4:7]' * Lmat(Quaternion(qa)) * Vᵀmat()
    X2 = res[8:10]'
    Q2 = res[11:14]' * Lmat(Quaternion(qb)) * Vᵀmat()

    n11 = norm(X1 - ∂g∂pos(joint1, 1, bot)[1:3]')
    n12 = norm(Q1 - ∂g∂pos(joint1, 1, bot)[4:6]')
    n21 = norm(X2 - ∂g∂pos(joint1, 2, bot)[1:3]')
    n22 = norm(Q2 - ∂g∂pos(joint1, 2, bot)[4:6]')

    res = ForwardDiff.gradient(transfunc1vel, [xa;qa;xb;qb;va;wa;vb;wb;pa;pb;v])
    V1 = res[15:17]'
    W1 = res[18:20]'
    V2 = res[21:23]'
    W2 = res[24:26]'

    n31 = norm(V1 - ∂g∂vel(joint1, 1, bot)[1:3]')
    n32 = norm(W1 - ∂g∂vel(joint1, 1, bot)[4:6]')
    n41 = norm(V2 - ∂g∂vel(joint1, 2, bot)[1:3]')
    n42 = norm(W2 - ∂g∂vel(joint1, 2, bot)[4:6]')

    display((n11, n12, n21, n22))
    display((n31, n32, n41, n42))
    return
end

function rottest3()
    xa = rand(3)
    qa = Quaternion(rand(RotMatrix{3}))
    xb = rand(3)
    qb = Quaternion(rand(RotMatrix{3}))

    va = rand(3)
    wa = rand(3)
    vb = rand(3)
    wb = rand(3)

    offset = Quaternion(rand(RotMatrix{3}))

    origin = Origin{Float64}()
    link1 = Body(Box(1., 1., 1., 1.))
    link2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(OriginConnection(origin, link1))
    oc2 = EqualityConstraint(OriginConnection(origin, link2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Rotational3(link1, link2, offset = offset))

    bot = Mechanism(origin, [link1;link2], [oc1;oc2;joint1])

    setPosition!(bot, link1, x = xa, q = qa)
    setPosition!(bot, link2, x = xb, q = qb)
    link1.s1 = SVector([va;wa]...)
    link2.s1 = SVector([vb;wb]...)

    res = ForwardDiff.jacobian(rotfunc3pos, [xa;qa;xb;qb;offset])
    X1 = res[1:3,1:3]
    Q1 = res[1:3,4:7] * Lmat(Quaternion(qa)) * Vᵀmat()
    X2 = res[1:3,8:10]
    Q2 = res[1:3,11:14] * Lmat(Quaternion(qb)) * Vᵀmat()

    n11 = norm(X1 - ∂g∂pos(joint1, 1, bot)[1:3,1:3])
    n12 = norm(Q1 - ∂g∂pos(joint1, 1, bot)[1:3,4:6])
    n21 = norm(X2 - ∂g∂pos(joint1, 2, bot)[1:3,1:3])
    n22 = norm(Q2 - ∂g∂pos(joint1, 2, bot)[1:3,4:6])

    res = ForwardDiff.jacobian(rotfunc3vel, [xa;qa;xb;qb;va;wa;vb;wb;offset])
    V1 = res[1:3,15:17]
    W1 = res[1:3,18:20]
    V2 = res[1:3,21:23]
    W2 = res[1:3,24:26]

    n31 = norm(V1 - ∂g∂vel(joint1, 1, bot)[1:3,1:3])
    n32 = norm(W1 - ∂g∂vel(joint1, 1, bot)[1:3,4:6])
    n41 = norm(V2 - ∂g∂vel(joint1, 2, bot)[1:3,1:3])
    n42 = norm(W2 - ∂g∂vel(joint1, 2, bot)[1:3,4:6])

    display((n11, n12, n21, n22))
    display((n31, n32, n41, n42))
    return
end

function rottest2()
    xa = rand(3)
    qa = Quaternion(rand(RotMatrix{3}))
    xb = rand(3)
    qb = Quaternion(rand(RotMatrix{3}))

    va = rand(3)
    wa = rand(3)
    vb = rand(3)
    wb = rand(3)

    offset = Quaternion(rand(RotMatrix{3}))

    v = rand(3)

    origin = Origin{Float64}()
    link1 = Body(Box(1., 1., 1., 1.))
    link2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(OriginConnection(origin, link1))
    oc2 = EqualityConstraint(OriginConnection(origin, link2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Rotational2(link1, link2, v, offset = offset))
    V12 = joint1.constraints[1].V12

    bot = Mechanism(origin, [link1;link2], [oc1;oc2;joint1])

    setPosition!(bot, link1, x = xa, q = qa)
    setPosition!(bot, link2, x = xb, q = qb)
    link1.s1 = SVector([va;wa]...)
    link2.s1 = SVector([vb;wb]...)

    res = ForwardDiff.jacobian(rotfunc2pos, [xa;qa;xb;qb;offset;V12[1,:];V12[2,:]])
    X1 = res[1:2,1:3]
    Q1 = res[1:2,4:7] * Lmat(Quaternion(qa)) * Vᵀmat()
    X2 = res[1:2,8:10]
    Q2 = res[1:2,11:14] * Lmat(Quaternion(qb)) * Vᵀmat()

    n11 = norm(X1 - ∂g∂pos(joint1, 1, bot)[1:2,1:3])
    n12 = norm(Q1 - ∂g∂pos(joint1, 1, bot)[1:2,4:6])
    n21 = norm(X2 - ∂g∂pos(joint1, 2, bot)[1:2,1:3])
    n22 = norm(Q2 - ∂g∂pos(joint1, 2, bot)[1:2,4:6])

    res = ForwardDiff.jacobian(rotfunc2vel, [xa;qa;xb;qb;va;wa;vb;wb;offset;V12[1,:];V12[2,:]])
    V1 = res[1:2,15:17]
    W1 = res[1:2,18:20]
    V2 = res[1:2,21:23]
    W2 = res[1:2,24:26]

    n31 = norm(V1 - ∂g∂vel(joint1, 1, bot)[1:2,1:3])
    n32 = norm(W1 - ∂g∂vel(joint1, 1, bot)[1:2,4:6])
    n41 = norm(V2 - ∂g∂vel(joint1, 2, bot)[1:2,1:3])
    n42 = norm(W2 - ∂g∂vel(joint1, 2, bot)[1:2,4:6])

    display((n11, n12, n21, n22))
    display((n31, n32, n41, n42))
    return
end

function rottest1()
    xa = rand(3)
    qa = Quaternion(rand(RotMatrix{3}))
    xb = rand(3)
    qb = Quaternion(rand(RotMatrix{3}))

    va = rand(3)
    wa = rand(3)
    vb = rand(3)
    wb = rand(3)

    offset = Quaternion(rand(RotMatrix{3}))

    v = rand(3)

    origin = Origin{Float64}()
    link1 = Body(Box(1., 1., 1., 1.))
    link2 = Body(Box(1., 1., 1., 1.))

    oc1 = EqualityConstraint(OriginConnection(origin, link1))
    oc2 = EqualityConstraint(OriginConnection(origin, link2))
    joint1 = EqualityConstraint(ConstrainedDynamics.Rotational1(link1, link2, v, offset = offset))
    V3 = joint1.constraints[1].V3'

    bot = Mechanism(origin, [link1;link2], [oc1;oc2;joint1])

    setPosition!(bot, link1, x = xa, q = qa)
    setPosition!(bot, link2, x = xb, q = qb)
    link1.s1 = SVector([va;wa]...)
    link2.s1 = SVector([vb;wb]...)

    res = ForwardDiff.gradient(rotfunc1pos, [xa;qa;xb;qb;offset;V3])
    X1 = res[1:3]'
    Q1 = res[4:7]' * Lmat(Quaternion(qa)) * Vᵀmat()
    X2 = res[8:10]'
    Q2 = res[11:14]' * Lmat(Quaternion(qb)) * Vᵀmat()

    n11 = norm(X1 - ∂g∂pos(joint1, 1, bot)[1:3]')
    n12 = norm(Q1 - ∂g∂pos(joint1, 1, bot)[4:6]')
    n21 = norm(X2 - ∂g∂pos(joint1, 2, bot)[1:3]')
    n22 = norm(Q2 - ∂g∂pos(joint1, 2, bot)[4:6]')

    res = ForwardDiff.gradient(rotfunc1vel, [xa;qa;xb;qb;va;wa;vb;wb;offset;V3])
    V1 = res[15:17]'
    W1 = res[18:20]'
    V2 = res[21:23]'
    W2 = res[24:26]'

    n31 = norm(V1 - ∂g∂vel(joint1, 1, bot)[1:3]')
    n32 = norm(W1 - ∂g∂vel(joint1, 1, bot)[4:6]')
    n41 = norm(V2 - ∂g∂vel(joint1, 2, bot)[1:3]')
    n42 = norm(W2 - ∂g∂vel(joint1, 2, bot)[4:6]')

    display((n11, n12, n21, n22))
    display((n31, n32, n41, n42))
    return
end


text:
module Roadnet_MDP

using POMDPs, MetaGraphs, LightGraphs

#  export roadnet_with_pursuer, roadnet_pursuer_state, samenode, rewards, exit_nodes
export roadnet_with_pursuer, roadnet_pursuer_state, rewards, exit_nodes

######### state definition
struct roadnet_pursuer_state
    node::Int64 # current UGV node
    pnode::Int64 # current puruser node
end

function exit_nodes(g::MetaGraph)
    @assert :POMDPgraph in keys(g.gprops)

    return g.gprops[:exit_nodes]
end
function rewards(g::MetaGraph)
    # make sure this is a POMDP MetaGraph -- i.e. one we created with the special metadata
    @assert :POMDPgraph in keys(g.gprops)

    num_vert = nv(g)
    states = Array{roadnet_pursuer_state}(0)
    vals = Array{Float64}(0)
    rwd = g.gprops[:reward_dict]
    # may need to add current believed state in here to account for the reward of hitting the intruder
    for i in vertices(g.graph)
        if g.vprops[i][:state_property] in keys(rwd) # check to see if the state has a reward assigned
            push!(states, roadnet_pursuer_state(i,4))
            push!(vals, rwd[g.vprops[i][:state_property]])
        end
    end
    return states, vals
end

######### MDP definition
struct roadnet_with_pursuer <: MDP{roadnet_pursuer_state,Symbol}
    tprob::Float64 #probability of transitioning to desired state
    discount::Float64 #discount factor
    reward_vals::Vector{Float64}
    reward_states::Vector{roadnet_pursuer_state}
    exit_nodes::Vector{Int64} #list of states that are exits
    road_net::MetaGraph
end

function roadnet_with_pursuer(rn::MetaGraph;tp::Float64=1.0,
                               d::Float64=0.9,
                               rv::Vector{Float64}=rewards(rn)[2],
                               rs::Vector{roadnet_pursuer_state}=rewards(rn)[1],
                               es::Vector{Int64}=exit_nodes(rn))
    return roadnet_with_pursuer(tp,d,rv,rs,es,rn)
end

end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC