text:
### indicator data -> string
########################################

vec2string(p::StaticVector{2}) = @sprintf("(%0.3f, %0.3f)", p[1], p[2])
vec2string(p::StaticVector{3}) = @sprintf("(%0.3f, %0.3f, %0.3f)", p[1], p[2], p[3])

position2string(p::StaticVector{2}) = @sprintf("x: %0.6f\ny: %0.6f", p[1], p[2])
position2string(p::StaticVector{3}) = @sprintf("x: %0.6f\ny: %0.6f\nz: %0.6f", p[1], p[2], p[3])

function bbox2string(bbox::Rect3)
    p0 = origin(bbox)
    p1 = p0 .+ widths(bbox)
    @sprintf(
        """
        Bounding Box:
         x: (%0.3f, %0.3f)
         y: (%0.3f, %0.3f)
         z: (%0.3f, %0.3f)
        """,
        p0[1], p1[1], p0[2], p1[2], p0[3], p1[3]
    )
end

function bbox2string(bbox::Rect2)
    p0 = origin(bbox)
    p1 = p0 .+ widths(bbox)
    @sprintf(
        """
        Bounding Box:
         x: (%0.3f, %0.3f)
         y: (%0.3f, %0.3f)
        """,
        p0[1], p1[1], p0[2], p1[2]
    )
end

color2text(c::AbstractFloat) = @sprintf("%0.3f", c)
color2text(c::Symbol) = string(c)
color2text(c) = color2text(to_color(c))
function color2text(c::RGBAf)
    if c.alpha == 1.0
        @sprintf("RGB(%0.2f, %0.2f, %0.2f)", c.r, c.g, c.b)
    else
        @sprintf("RGBA(%0.2f, %0.2f, %0.2f, %0.2f)", c.r, c.g, c.b, c.alpha)
    end
end

color2text(name, i::Integer, j::Integer, c) = "$name[$i, $j] = $(color2text(c))"
function color2text(name, i, j, c)
    idxs = @sprintf("%0.2f, %0.2f", i, j)
    "$name[$idxs] = $(color2text(c))"
end

### dealing with markersize and rotations
########################################

_to_scale(f::AbstractFloat, idx) = Vec3f(f)
_to_scale(v::Vec2f, idx) = Vec3f(v[1], v[2], 1)
_to_scale(v::Vec3f, idx) = v
_to_scale(v::Vector, idx) = _to_scale(v[idx], idx)

_to_rotation(x, idx) = to_rotation(x)
_to_rotation(x::Vector, idx) = to_rotation(x[idx])

### Selecting a point on a nearby line
########################################

function closest_point_on_line(p0::Point2f, p1::Point2f, r::Point2f)
    # This only works in 2D
    AP = P .- A; AB = B .- A
    A .+ AB * dot(AP, AB) / dot(AB, AB)
end

function view_ray(scene)
    inv_projview = inv(camera(scene).projectionview[])
    view_ray(inv_projview, events(scene).mouseposition[], pixelarea(scene)[])
end
function view_ray(inv_view_proj, mpos, area::Rect2)
    # This figures out the camera view direction from the projectionview matrix (?)
    # and computes a ray from a near and a far point.
    # Based on ComputeCameraRay from ImGuizmo
    mp = 2f0 .* (mpos .- minimum(area)) ./ widths(area) .- 1f0
    v = inv_view_proj * Vec4f(0, 0, -10, 1)
    reversed = v[3] < v[4]
    near = reversed ? 1f0 - 1e-6 : 0f0
    far = reversed ? 0f0 : 1f0 - 1e-6

    origin = inv_view_proj * Vec4f(mp[1], mp[2], near, 1f0)
    origin = origin[SOneTo(3)] ./ origin[4]

    p = inv_view_proj * Vec4f(mp[1], mp[2], far, 1f0)
    p = p[SOneTo(3)] ./ p[4]

    dir = normalize(p - origin)
    return origin, dir
end

# These work in 2D and 3D
function closest_point_on_line(A, B, origin, dir)
    closest_point_on_line(
        to_ndim(Point3f, A, 0),
        to_ndim(Point3f, B, 0),
        to_ndim(Point3f, origin, 0),
        to_ndim(Vec3f, dir, 0)
    )
end
function closest_point_on_line(A::Point3f, B::Point3f, origin::Point3f, dir::Vec3f)
    # See:
    # https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection
    AB_norm = norm(B .- A)
    u_AB = (B .- A) / AB_norm
    u_dir = normalize(dir)
    u_perp = normalize(cross(u_dir, u_AB))
    # e_RD, e_perp defines a plane with normal n
    n = normalize(cross(u_dir, u_perp))
    t = dot(origin .- A, n) / dot(u_AB, n)
    A .+ clamp(t, 0.0, AB_norm) * u_AB
end

function ray_triangle_intersection(A, B, C, origin, dir)
    # See: https://www.iue.tuwien.ac.at/phd/ertl/node114.html
    AO = A .- origin
    BO = B .- origin
    CO = C .- origin
    A1 = 0.5 * dot(cross(BO, CO), dir)
    A2 = 0.5 * dot(cross(CO, AO), dir)
    A3 = 0.5 * dot(cross(AO, BO), dir)

    e = 1e-3
    if (A1 > -e && A2 > -e && A3 > -e) || (A1 < e && A2 < e && A3 < e)
        Point3f((A1 * A .+ A2 * B .+ A3 * C) / (A1 + A2 + A3))
    else
        Point3f(NaN)
    end
end

### Surface positions
########################################

surface_x(xs::ClosedInterval, i, j, N) = minimum(xs) + (maximum(xs) - minimum(xs)) * (i-1) / (N-1)
surface_x(xs, i, j, N) = xs[i]
surface_x(xs::AbstractMatrix, i, j, N) = xs[i, j]

surface_y(ys::ClosedInterval, i, j, N) = minimum(ys) + (maximum(ys) - minimum(ys)) * (j-1) / (N-1)
surface_y(ys, i, j, N) = ys[j]
surface_y(ys::AbstractMatrix, i, j, N) = ys[i, j]

function surface_pos(xs, ys, zs, i, j)
    N, M = size(zs)
    Point3f(surface_x(xs, i, j, N), surface_y(ys, i, j, M), zs[i, j])
end

### Mapping mesh vertex indices to Vector{Polygon} index
########################################

function vertexindex2poly(polys, idx)
    counter = 0
    for i in eachindex(polys)
        step = ncoords(polys[i])
        if idx <= counter + step
            return i
        else
            counter += step
        end
    end
    return length(polys)
end

ncoords(x) = length(coordinates(x))
ncoords(mesh::Mesh) = length(coordinates(mesh))
function ncoords(poly::Polygon)
    N = length(poly.exterior) + 1
    for int in poly.interiors
        N += length(int) + 1
    end
    N
end

### Text bounding box
########################################

function Bbox_from_glyphcollection(text, gc)
    bbox = Rect2f(0, 0, 0, 0)
    bboxes = Rect2f[]
    broadcast_foreach(gc.extents, gc.fonts, gc.scales) do extent, font, scale
        b = FreeTypeAbstraction.height_insensitive_boundingbox(extent, font) * scale
        push!(bboxes, b)
    end
    for (c, o, bb) in zip(text, gc.origins, bboxes)
        c == '\n' && continue
        bbox2 = Rect2f(o[Vec(1,2)] .+ origin(bb), widths(bb))
        if bbox == Rect2f(0, 0, 0, 0)
            bbox = bbox2
        else
            bbox = union(bbox, bbox2)
        end
    end
    bbox
end

## Shifted projection
########################################

function shift_project(scene, plot, pos)
    project(
        camera(scene).projectionview[],
        Vec2f(widths(pixelarea(scene)[])),
        apply_transform(transform_func_obs(plot)[], pos)
    ) .+ Vec2f(origin(pixelarea(scene)[]))
end

################################################################################
### Base pixel-space plot for Indicator
################################################################################

@recipe(_Inspector, x) do scene
    # Attributes starting with _ are modified internally
    Attributes(
        # Text
        text_padding = Vec4f(5, 5, 3, 3), # LRBT
        text_align = (:left, :bottom),
        textcolor = :black,
        textsize = 18,
        font = theme(scene, :font),
        _display_text = " ",
        _text_position = Point2f(0),

        # Background
        background_color = :white,
        outline_color = :grey,
        outline_linestyle = nothing,
        outline_linewidth = 2,

        # pixel BBox/indicator
        indicator_color = :red,
        indicator_linewidth = 2,
        indicator_linestyle = nothing,
        _bbox2D = Rect2f(Vec2f(0), Vec2f(0)),
        _px_bbox_visible = true,

        # general
        tooltip_align = (:center, :top), # default tooltip position relative to cursor
        tooltip_offset = Vec2f(20), # default offset in alignment direction
        depth = 9e3,
        enabled = true,
        range = 10,

        _root_px_projection = Mat4f(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1),
        _model = Mat4f(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1),
        _visible = true,
        _tooltip_align = (:center, :top),
        _bbox3D = Rect3f(Vec3f(0), Vec3f(0)),
        _bbox_visible = true,
        _position = Point3f(0),
        _color = RGBAf(0,0,0,0)
    )
end

function plot!(plot::_Inspector)
    @extract plot (
        text_padding, text_align, textcolor, textsize, font,
        background_color, outline_color, outline_linestyle, outline_linewidth,
        indicator_linestyle, indicator_linewidth, indicator_color,
        tooltip_offset, depth,
        _display_text, _text_position, _bbox2D, _px_bbox_visible,
        _tooltip_align, _root_px_projection, _visible
    )

    # tooltip text
    _aligned_text_position = Observable(Point2f(0))
    id = Mat4f(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1)
    text_plot = text!(plot, _display_text,
        position = _aligned_text_position, visible = _visible, align = text_align,
        color = textcolor, font = font, textsize = textsize,
        inspectable = false,
        # with https://github.com/JuliaPlots/Makie.jl/tree/master/GLMakie/pull/183 this should
        # allow the tooltip to work in any scene.
        space = :data,
        projection = _root_px_projection, view = id, projectionview = _root_px_projection
    )

    # compute text boundingbox and adjust _aligned_text_position
    bbox = map(text_plot.plots[1][1], text_plot.position, text_padding) do gc, pos, pad
        rect = Bbox_from_glyphcollection(_display_text[], gc)
        l, r, b, t = pad
        Rect2f(
            origin(rect) .+ Vec2f(pos[1] - l, pos[2] - b),
            widths(rect) .+ Vec2f(l + r, b + t)
        )
    end
    onany(_text_position, _tooltip_align, tooltip_offset, bbox) do pos, align, offset, bbox
        halign, valign = align
        ox, oy = offset
        wx, wy = widths(bbox)
        if halign == :left
            dx = -wx - ox
        elseif halign == :center
            dx = -0.5wx
        elseif halign == :right
            dx = + ox
        end
        if valign == :top
            dy = oy
        elseif valign == :center
            dy = -0.5wy
        elseif valign == :bottom
            dy = -wy - oy
        end
        new_pos = pos .+ Point2f(dx, dy)
        if new_pos != _aligned_text_position[]
            _aligned_text_position[] = new_pos
        end
    end

    # tooltip background and frame
    background = mesh!(
        plot, bbox, color = background_color, shading = false, #fxaa = false,
        # TODO with fxaa here the text above becomes seethrough on a heatmap
        visible = _visible, inspectable = false,
        projection = _root_px_projection, view = id, projectionview = _root_px_projection
    )
    outline = wireframe!(
        plot, bbox,
        color = outline_color, visible = _visible, inspectable = false,
        linestyle = outline_linestyle, linewidth = outline_linewidth,
        projection = _root_px_projection, view = id, projectionview = _root_px_projection
    )

    # pixel-space marker for selected element (not always used)
    px_bbox = wireframe!(
        plot, _bbox2D,
        color = indicator_color, linewidth = indicator_linewidth,
        linestyle = indicator_linestyle, visible = _px_bbox_visible,
        inspectable = false,
        projection = _root_px_projection, view = id, projectionview = _root_px_projection
    )

    # To make sure inspector plots end up in front
    on(depth) do d
        translate!(background, Vec3f(0,0,d+1))
        translate!(outline,    Vec3f(0,0,d+2))
        translate!(text_plot,  Vec3f(0,0,d+3))
        translate!(px_bbox,    Vec3f(0,0,d))
    end
    depth[] = depth[]
    nothing
end

################################################################################
### Interactive selection via DataInspector
################################################################################

# TODO destructor?
mutable struct DataInspector
    root::Scene
    temp_plots::Vector{AbstractPlot}
    plot::_Inspector
    selection::AbstractPlot
    obsfuncs::Vector{Any}
end

function DataInspector(scene::Scene, plot::AbstractPlot)
    x = DataInspector(scene, AbstractPlot[], plot, plot, Any[])
    # finalizer(cleanup, x) # doesn't get triggered when this is dereferenced
    x
end

function cleanup(inspector::DataInspector)
    off.(inspector.obsfuncs)
    empty!(inspector.obsfuncs)
    delete!(inspector.root, inspector.plot)
    clear_temporary_plots!(inspector, inspector.slection)
    inspector
end

function Base.delete!(::Union{Scene, Figure}, inspector::DataInspector)
    cleanup(inspector)
end

enable!(inspector::DataInspector) = inspector.plot.enabled[] = true
disable!(inspector::DataInspector) = inspector.plot.enabled[] = false

"""
    DataInspector(figure; kwargs...)

Creates a data inspector which will show relevant information in a tooltip
when you hover over a plot. If you wish to exclude a plot you may set
`plot.inspectable[] = false`.

### Keyword Arguments:
- `range = 10`: Controls the snapping range for selecting an element of a plot.
- `enabled = true`: Disables inspection of plots when set to false. Can also be
    adjusted with `enable!(inspector)` and `disable!(inspector)`.
- `text_padding = Vec4f(5, 5, 3, 3)`: Padding for the box drawn around the
    tooltip text. (left, right, bottom, top)
- `text_align = (:left, :bottom)`: Alignment of text within the tooltip. This
    does not affect the alignment of the tooltip relative to the cursor.
- `textcolor = :black`: Tooltip text color.
- `textsize = 20`: Tooltip text size.
- `font = "Dejavu Sans"`: Tooltip font.
- `background_color = :white`: Background color of the tooltip.
- `outline_color = :grey`: Outline color of the tooltip.
- `outline_linestyle = nothing`: Linestyle of the tooltip outline.
- `outline_linewidth = 2`: Linewidth of the tooltip outline.
- `indicator_color = :red`: Color of the selection indicator.
- `indicator_linewidth = 2`: Linewidth of the selection indicator.
- `indicator_linestyle = nothing`: Linestyle of the selection indicator
- `tooltip_align = (:center, :top)`: Default position of the tooltip relative to
    the cursor or current selection. The real align may adjust to keep the
    tooltip in view.
- `tooltip_offset = Vec2f(20)`: Offset from the indicator to the tooltip.
- `depth = 9e3`: Depth value of the tooltip. This should be high so that the
    tooltip is always in front.
- `priority = 100`: The priority of creating a tooltip on a mouse movement or
    scrolling event.
"""
function DataInspector(fig_or_layoutable; kwargs...)
    DataInspector(fig_or_layoutable.scene; kwargs...)
end

function DataInspector(scene::Scene; priority = 100, kwargs...)
    parent = root(scene)
    @assert origin(pixelarea(parent)[]) == Vec2f(0)

    plot = _inspector!(
        parent, 1,
        _root_px_projection = camera(parent).pixel_space;
        kwargs...
    )
    inspector = DataInspector(parent, plot)

    e = events(parent)
    f1 = on(_ -> on_hover(inspector), e.mouseposition, priority = priority)
    f2 = on(_ -> on_hover(inspector), e.scroll, priority = priority)

    push!(inspector.obsfuncs, f1, f2)

    inspector
end

function on_hover(inspector)
    parent = inspector.root
    (inspector.plot.enabled[] && is_mouseinside(parent)) || return Consume(false)

    mp = mouseposition_px(parent)
    should_clear = true
    for (plt, idx) in pick_sorted(parent, mp, inspector.plot.range[])
        if to_value(get(plt.attributes, :inspectable, true))
            # show_data should return true if it created a tooltip
            if show_data_recursion(inspector, plt, idx)
                should_clear = false
                break
            end
        end
    end

    if should_clear
        inspector.plot._visible[] = false
        inspector.plot._bbox_visible[] = false
        inspector.plot._px_bbox_visible[] = false
    end

    return Consume(false)
end

function show_data_recursion(inspector, plot, idx)
    processed = show_data_recursion(inspector, plot.parent, idx, plot)
    if processed
        return true
    else
        return show_data(inspector, plot, idx)
    end
end
show_data_recursion(inspector, plot, idx, source) = false
function show_data_recursion(inspector, plot::AbstractPlot, idx, source)
    processed = show_data_recursion(inspector, plot.parent, idx, source)
    if processed
        return true
    else
        return show_data(inspector, plot, idx, source)
    end
end

# clears temporary plots (i.e. bboxes) and update selection
function clear_temporary_plots!(inspector::DataInspector, plot)
    inspector.selection = plot
    for i in length(inspector.obsfuncs):-1:3
        off(pop!(inspector.obsfuncs))
    end
    for p in inspector.temp_plots
        delete!(parent_scene(p), p)
    end
    empty!(inspector.temp_plots)
end

# update alignment direction
function update_tooltip_alignment!(inspector, proj_pos)
    a = inspector.plot.attributes
    a._text_position[] = proj_pos

    wx, wy = widths(pixelarea(inspector.root)[])
    px, py = proj_pos
    halign, valign = a.tooltip_align[]
    px < wx/4  && (halign = :right)
    px > 3wx/4 && (halign = :left)
    py < wy/4  && (valign = :top)
    py > 3wy/4 && (valign = :bottom)
    a._tooltip_align[] = (halign, valign)
end

################################################################################
### show_data for primitive plots
################################################################################

# TODO: better 3D scaling
function show_data(inspector::DataInspector, plot::Scatter, idx)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    proj_pos = shift_project(scene, plot, to_ndim(Point3f, plot[1][][idx], 0))
    update_tooltip_alignment!(inspector, proj_pos)
    ms = plot.markersize[]

    a._display_text[] = position2string(plot[1][][idx])
    a._bbox2D[] = Rect2f(proj_pos .- 0.5 .* ms .- Vec2f(5), Vec2f(ms) .+ Vec2f(10))
    a._px_bbox_visible[] = true
    a._bbox_visible[] = false
    a._visible[] = true

    return true
end

function show_data(inspector::DataInspector, plot::MeshScatter, idx)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    proj_pos = shift_project(scene, plot, to_ndim(Point3f, plot[1][][idx], 0))
    update_tooltip_alignment!(inspector, proj_pos)
    bbox = Rect{3, Float32}(convert_attribute(
        plot.marker[], Key{:marker}(), Key{Makie.plotkey(plot)}()
    ))

    a._model[] = transformationmatrix(
        plot[1][][idx],
        _to_scale(plot.markersize[], idx),
        _to_rotation(plot.rotations[], idx)
    )

    if inspector.selection != plot
        cc = cameracontrols(scene)
        if cc isa Camera3D
            eyeposition = cc.eyeposition[]
            lookat = cc.lookat[]
            upvector = cc.upvector[]
        end

        # To avoid putting a bbox outside the plots bbox
        a._bbox3D[] = boundingbox(plot)

        clear_temporary_plots!(inspector, plot)
        p = wireframe!(
            scene, a._bbox3D, model = a._model, color = a.indicator_color,
            linewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
            visible = a._bbox_visible, inspectable = false
        )
        push!(inspector.temp_plots, p)

        # Restore camera
        cc isa Camera3D && update_cam!(scene, eyeposition, lookat, upvector)
    end

    a._display_text[] = position2string(plot[1][][idx])
    a._bbox3D[] = bbox
    a._px_bbox_visible[] = false
    a._bbox_visible[] = true
    a._visible[] = true

    return true
end

function show_data(inspector::DataInspector, plot::Union{Lines, LineSegments}, idx)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    # cast ray from cursor into screen, find closest point to line
    p0, p1 = plot[1][][idx-1:idx]
    origin, dir = view_ray(scene)
    pos = closest_point_on_line(p0, p1, origin, dir)
    lw = plot.linewidth[] isa Vector ? plot.linewidth[][idx] : plot.linewidth[]

    proj_pos = shift_project(scene, plot, to_ndim(Point3f, pos, 0))
    update_tooltip_alignment!(inspector, proj_pos)

    a._display_text[] = position2string(typeof(p0)(pos))
    a._bbox2D[] = Rect2f(proj_pos .- 0.5 .* lw .- Vec2f(5), Vec2f(lw) .+ Vec2f(10))
    a._px_bbox_visible[] = true
    a._bbox_visible[] = false
    a._visible[] = true

    return true
end

function show_data(inspector::DataInspector, plot::Mesh, idx)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    bbox = boundingbox(plot)
    proj_pos = Point2f(mouseposition_px(inspector.root))
    update_tooltip_alignment!(inspector, proj_pos)

    a._bbox3D[] = bbox

    if inspector.selection != plot
        cc = cameracontrols(scene)
        if cc isa Camera3D
            eyeposition = cc.eyeposition[]
            lookat = cc.lookat[]
            upvector = cc.upvector[]
        end

        clear_temporary_plots!(inspector, plot)
        p = wireframe!(
            scene, a._bbox3D, color = a.indicator_color,
            linewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
            visible = a._bbox_visible, inspectable = false
        )
        push!(inspector.temp_plots, p)

        # Restore camera
        cc isa Camera3D && update_cam!(scene, eyeposition, lookat, upvector)
    end

    a._text_position[] = proj_pos
    a._display_text[] = bbox2string(bbox)
    a._px_bbox_visible[] = false
    a._bbox_visible[] = true
    a._visible[] = true

    return true
end

function show_data(inspector::DataInspector, plot::Surface, idx)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    proj_pos = Point2f(mouseposition_px(inspector.root))
    update_tooltip_alignment!(inspector, proj_pos)

    xs = plot[1][]
    ys = plot[2][]
    zs = plot[3][]
    w, h = size(zs)
    _i = mod1(idx, w); _j = div(idx-1, w)

    # This isn't the most accurate so we include some neighboring faces
    origin, dir = view_ray(scene)
    pos = Point3f(NaN)
    for i in _i-1:_i+1, j in _j-1:_j+1
        (1 <= i <= w) && (1 <= j < h) || continue

        if i - 1 > 0
            pos = ray_triangle_intersection(
                surface_pos(xs, ys, zs, i, j),
                surface_pos(xs, ys, zs, i-1, j),
                surface_pos(xs, ys, zs, i, j+1),
                origin, dir
            )
        end

        if i + 1 <= w && isnan(pos)
            pos = ray_triangle_intersection(
                surface_pos(xs, ys, zs, i, j),
                surface_pos(xs, ys, zs, i, j+1),
                surface_pos(xs, ys, zs, i+1, j+1),
                origin, dir
            )
        end

        isnan(pos) || break
    end

    if !isnan(pos)
        a._display_text[] = position2string(pos)
        a._text_position[] = proj_pos
        a._bbox2D[] = Rect2f(proj_pos .- Vec2f(5), Vec2f(10))
        a._bbox_visible[] = false
        a._px_bbox_visible[] = true
        a._visible[] = true
    else
        a._bbox_visible[] = false
        a._px_bbox_visible[] = false
        a._visible[] = false
    end

    return true
end

function show_data(inspector::DataInspector, plot::Heatmap, idx)
    show_imagelike(inspector, plot, "H", true)
end

function show_data(inspector::DataInspector, plot::Image, idx)
    show_imagelike(inspector, plot, "img", false)
end

function show_imagelike(inspector, plot, name, edge_based)
    a = inspector.plot.attributes
    scene = parent_scene(plot)
    mpos = mouseposition(scene)

    i, j, z = if plot.interpolate[]
        _interpolated_getindex(plot[1][], plot[2][], plot[3][], mpos)
    else
        _pixelated_getindex(plot[1][], plot[2][], plot[3][], mpos, edge_based)
    end

    a._color[] = if z isa AbstractFloat
        interpolated_getindex(
            to_colormap(plot.colormap[]), z,
            to_value(get(plot.attributes, :colorrange, (0, 1)))
        )
    else
        z
    end

    a._position[] = to_ndim(Point3f, mpos, 0)
    proj_pos = Point2f(mouseposition_px(inspector.root))
    update_tooltip_alignment!(inspector, proj_pos)

    if plot.interpolate[]
        if inspector.selection != plot || !(inspector.temp_plots[1] isa Scatter)
            clear_temporary_plots!(inspector, plot)
            p = scatter!(
                scene, map(p -> [p], a._position), color = a._color,
                visible = a._bbox_visible,
                inspectable = false,
                marker=Rect, markersize = map(r -> 2r - 4, a.range),
                strokecolor = a.indicator_color,
                strokewidth = a.indicator_linewidth #, linestyle = a.indicator_linestyle no?
            )
            translate!(p, Vec3f(0, 0, a.depth[]))
            push!(inspector.temp_plots, p)
            # Hacky?
            push!(
                inspector.obsfuncs,
                Observables.ObserverFunction(a._position.listeners[end], a._position, false)
            )
        end
        a._display_text[] = color2text(name, mpos[1], mpos[2], z)
    else
        a._bbox2D[] = _pixelated_image_bbox(plot[1][], plot[2][], plot[3][], i, j, edge_based)
        if inspector.selection != plot || !(inspector.temp_plots[1][1][] isa Rect2)
            clear_temporary_plots!(inspector, plot)
            p = wireframe!(
                scene, a._bbox2D, model = a._model, color = a.indicator_color,
                strokewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
                visible = a._bbox_visible, inspectable = false
            )
            translate!(p, Vec3f(0, 0, a.depth[]))
            push!(inspector.temp_plots, p)
        end
        a._display_text[] = color2text(name, i, j, z)
    end

    a._bbox_visible[] = true
    a._px_bbox_visible[] = false
    a._visible[] = true
    return true
end

function _interpolated_getindex(xs, ys, img, mpos)
    x0, x1 = extrema(xs)
    y0, y1 = extrema(ys)
    x, y = clamp.(mpos, (x0, y0), (x1, y1))

    i = clamp((x - x0) / (x1 - x0) * size(img, 1) + 0.5, 1, size(img, 1))
    j = clamp((y - y0) / (y1 - y0) * size(img, 2) + 0.5, 1, size(img, 2))
    l = clamp(floor(Int, i), 1, size(img, 1)-1);
    r = clamp(l+1, 2, size(img, 1))
    b = clamp(floor(Int, j), 1, size(img, 2)-1);
    t = clamp(b+1, 2, size(img, 2))
    z = ((r-i) * img[l, b] + (i-l) * img[r, b]) * (t-j) +
        ((r-i) * img[l, t] + (i-l) * img[r, t]) * (j-b)

    # float, float, value (i, j are no longer used)
    return i, j, z
end
function _pixelated_getindex(xs, ys, img, mpos, edge_based)
    x0, x1 = extrema(xs)
    y0, y1 = extrema(ys)
    x, y = clamp.(mpos, (x0, y0), (x1, y1))

    i = clamp(round(Int, (x - x0) / (x1 - x0) * size(img, 1) + 0.5), 1, size(img, 1))
    j = clamp(round(Int, (y - y0) / (y1 - y0) * size(img, 2) + 0.5), 1, size(img, 2))

    # int, int, value
    return i, j, img[i,j]
end

function _interpolated_getindex(xs::Vector, ys::Vector, img, mpos)
    # x, y = mpos
    # i, j, _ = _pixelated_getindex(xs, ys, img, mpos, false)
    # w = (xs[i+1] - xs[i]); h = (ys[j+1] - ys[j])
    # z = ((xs[i+1] - x) / w * img[i, j]   + (x - xs[i]) / w * img[i+1, j])   * (ys[j+1] - y) / h +
    #     ((xs[i+1] - x) / w * img[i, j+1] + (x - xs[i]) / w * img[i+1, j+1]) * (y - ys[j]) / h
    # return i, j, z
    _interpolated_getindex(minimum(xs)..maximum(xs), minimum(ys)..maximum(ys), img, mpos)
end
function _pixelated_getindex(xs::Vector, ys::Vector, img, mpos, edge_based)
    if edge_based
        x, y = mpos
        i = max(1, something(findfirst(v -> v >= x, xs), length(xs))-1)
        j = max(1, something(findfirst(v -> v >= y, ys), length(ys))-1)
        return i, j, img[i, j]
    else
        _pixelated_getindex(minimum(xs)..maximum(xs), minimum(ys)..maximum(ys), img, mpos, edge_based)
    end
end

function _pixelated_image_bbox(xs, ys, img, i::Integer, j::Integer, edge_based)
    x0, x1 = extrema(xs)
    y0, y1 = extrema(ys)
    nw, nh = ((x1 - x0), (y1 - y0)) ./ size(img)
    Rect2f(x0 + nw * (i-1), y0 + nh * (j-1), nw, nh)
end
function _pixelated_image_bbox(xs::Vector, ys::Vector, img, i::Integer, j::Integer, edge_based)
    if edge_based
        Rect2f(xs[i], ys[j], xs[i+1] - xs[i], ys[j+1] - ys[j])
    else
        _pixelated_image_bbox(
            minimum(xs)..maximum(xs), minimum(ys)..maximum(ys),
            img, i, j, edge_based
        )
    end
end

function show_data(inspector::DataInspector, plot, idx, source=nothing)
    return false
end

################################################################################
### show_data for Combined/recipe plots
################################################################################

function show_data(inspector::DataInspector, plot::BarPlot, idx, ::Lines)
    return show_data(inspector, plot, div(idx-1, 6)+1)
end

function show_data(inspector::DataInspector, plot::BarPlot, idx, ::Mesh)
    return show_data(inspector, plot, div(idx-1, 4)+1)
end

function show_data(inspector::DataInspector, plot::BarPlot, idx)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    pos = plot[1][][idx]
    proj_pos = shift_project(scene, plot, to_ndim(Point3f, pos, 0))
    update_tooltip_alignment!(inspector, proj_pos)
    a._model[] = plot.model[]
    a._bbox2D[] = plot.plots[1][1][][idx]

    if inspector.selection != plot
        clear_temporary_plots!(inspector, plot)
        p = wireframe!(
            scene, a._bbox2D, model = a._model, color = a.indicator_color,
            strokewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
            visible = a._bbox_visible, inspectable = false
        )
        translate!(p, Vec3f(0, 0, a.depth[]))
        push!(inspector.temp_plots, p)
    end

    a._display_text[] = position2string(pos)
    a._bbox_visible[] = true
    a._px_bbox_visible[] = false
    a._visible[] = true

    return true
end

function show_data(inspector::DataInspector, plot::Arrows, idx, ::LineSegments)
    return show_data(inspector, plot, div(idx+1, 2), nothing)
end
function show_data(inspector::DataInspector, plot::Arrows, idx, source)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    pos = plot[1][][idx]
    proj_pos = shift_project(scene, plot, to_ndim(Point3f, pos, 0))

    mpos = Point2f(mouseposition_px(inspector.root))
    update_tooltip_alignment!(inspector, mpos)

    p = vec2string(pos)
    v = vec2string(plot[2][][idx])

    a._text_position[] = mpos
    a._display_text[] = "Position:\n  $p\nDirection\n  $v"
    a._bbox2D[] = Rect2f(proj_pos .- Vec2f(5), Vec2f(10))
    a._bbox_visible[] = false
    a._px_bbox_visible[] = true
    a._visible[] = true

    return true
end

# This should work if contourf would place computed levels in colors and let the
# backend handle picking colors from a colormap
function show_data(inspector::DataInspector, plot::Contourf, idx, source::Mesh)
    a = inspector.plot.attributes
    scene = parent_scene(plot)
    idx, ext = show_poly(inspector, plot.plots[1], idx, source)
    level = plot.plots[1].color[][idx]

    a._text_position[] = Point2f(mouseposition_px(inspector.root))
    a._display_text[] = @sprintf("level = %0.3f", level)
    return true
end

# What should this display?
# function show_data(
#         inspector::DataInspector, plot::Poly{<: Tuple{<: AbstractVector}},
#         idx, source::Mesh
#     )
#     @info "PolyMesh"
#     idx, ext = show_poly(inspectable, plot, idx, source)
#     return true
# end

function show_poly(inspector, plot, idx, source)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    idx = vertexindex2poly(plot[1][], idx)
    m = GeometryBasics.mesh(plot[1][][idx])

    clear_temporary_plots!(inspector, plot)
    ext = plot[1][][idx].exterior
    p = lines!(
        scene, ext, color = a.indicator_color,
        strokewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
        visible = a._visible, inspectable = false
    )
    translate!(p, Vec3f(0,0,a.depth[]))
    push!(inspector.temp_plots, p)

    for int in plot[1][][idx].interiors
        p = lines!(
            scene, int, color = a.indicator_color,
            strokewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
            visible = a._visible, inspectable = false
        )
        translate!(p, Vec3f(0,0,a.depth[]))
        push!(inspector.temp_plots, p)
    end

    a._px_bbox_visible[] = false
    a._bbox_visible[] = true
    a._visible[] = true

    return idx, ext
end

function show_data(inspector::DataInspector, plot::VolumeSlices, idx, child::Heatmap)
    a = inspector.plot.attributes
    scene = parent_scene(plot)

    proj_pos = Point2f(mouseposition_px(inspector.root))
    update_tooltip_alignment!(inspector, proj_pos)

    qs = extrema(child[1][])
    ps = extrema(child[2][])
    data = child[3][]
    trans = child.transformation.translation[]

    child_idx = findfirst(isequal(child), plot.plots)
    if child_idx == 2
        vs = [ # clockwise
            Point3f(trans[1], qs[1], ps[1]),
            Point3f(trans[1], qs[1], ps[2]),
            Point3f(trans[1], qs[2], ps[2]),
            Point3f(trans[1], qs[2], ps[1])
        ]
    elseif child_idx == 3
        vs = [ # clockwise
            Point3f(qs[1], trans[2], ps[1]),
            Point3f(qs[1], trans[2], ps[2]),
            Point3f(qs[2], trans[2], ps[2]),
            Point3f(qs[2], trans[2], ps[1])
        ]
    else
        vs = [ # clockwise
            Point3f(qs[1], ps[1], trans[3]),
            Point3f(qs[1], ps[2], trans[3]),
            Point3f(qs[2], ps[2], trans[3]),
            Point3f(qs[2], ps[1], trans[3])
        ]
    end

    origin, dir = view_ray(scene)
    pos = Point3f(NaN)
    pos = ray_triangle_intersection(vs[1], vs[2], vs[3], origin, dir)
    if isnan(pos)
        pos = ray_triangle_intersection(vs[3], vs[4], vs[1], origin, dir)
    end

    if !isnan(pos)
        if child_idx == 2
            x = pos[2]; y = pos[3]
        elseif child_idx == 3
            x = pos[1]; y = pos[3]
        else
            x = pos[1]; y = pos[2]
        end
        i = clamp(round(Int, (x - qs[1]) / (qs[2] - qs[1]) * size(data, 1) + 0.5), 1, size(data, 1))
        j = clamp(round(Int, (y - ps[1]) / (ps[2] - ps[1]) * size(data, 2) + 0.5), 1, size(data, 2))
        val = data[i, j]

        a._display_text[] = @sprintf(
            "x: %0.6f\ny: %0.6f\nz: %0.6f\n%0.6f0",
            pos[1], pos[2], pos[3], val
        )
        a._text_position[] = proj_pos
        a._bbox2D[] = Rect2f(proj_pos .- Vec2f(5), Vec2f(10))
        a._bbox_visible[] = false
        a._px_bbox_visible[] = true
        a._visible[] = true
    else
        a._bbox_visible[] = false
        a._px_bbox_visible[] = false
        a._visible[] = false
    end

    return true
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test
using StatsBase
using DoubleDouble

@testset "onlineuvstats" begin
    n = 10
    data  = [(-1)^x*exp(x-n/2) for x in 1:n]
    w = abs.(data)

    T = Float64
    wK = ProbabilityWeights
    mdata = mean(data, wK(w))
    vdata = var(data, wK(w); corrected=true)
    maxdata = maximum(data)
    mindata = minimum(data)

    @testset "BAT.OnlineUvMean" begin
        @test typeof(@inferred BAT.OnlineUvMean()) <: BAT.OnlineUvMean{Float64}
        @test typeof(@inferred BAT.OnlineUvMean{Float32}()) <: BAT.OnlineUvMean{Float32}
        ouvm = OnlineUvMean()

        res = append!(ouvm, data, w)
        @test res[] ≈ mdata

        numMeans = 3
        means = Array{BAT.OnlineUvMean{T}}(numMeans)
        for i in indices(means, 1)
            means[i] = OnlineUvMean()
        end

        for i in indices(data, 1)
            x = (i % numMeans) + 1
            means[x] = push!(means[x], data[i], w[i]);
        end

        @test merge!(means...)[] ≈ mdata

    end

    @testset "BAT.OnlineUvVar" begin
        @test typeof(@inferred BAT.OnlineUvVar()) <: BAT.OnlineUvVar{Float64, ProbabilityWeights}
        @test typeof(@inferred BAT.OnlineUvVar{Float32, FrequencyWeights}()) <: BAT.OnlineUvVar{Float32, FrequencyWeights}

        for wKind in [ProbabilityWeights, FrequencyWeights, AnalyticWeights, Weights]
            res = @inferred BAT.OnlineUvVar{Float64, wKind}()
            res = append!(res, data, w)
            @test res[] ≈ var(data, wKind(w); corrected=(wKind != Weights))
        end

        numVars = 3
        vars = Array{BAT.OnlineUvVar{T, wK}}(numVars)
        for i in indices(vars, 1)
            vars[i] = OnlineUvVar()
        end

        for i in indices(data, 1)
            x = (i % numVars) + 1
            vars[x] = push!(vars[x], data[i], w[i]);
        end

        merge!(vars...)
        @test vars[1][] ≈ var(data, wK(w); corrected=true)
        @test push!(vars[1], data[1], zero(Float64))[] ≈ var(data, wK(w); corrected=true)
    end

    @testset "BAT.BasicUvStatistics" begin
        @test typeof(@inferred BAT.BasicUvStatistics{Float32, FrequencyWeights}()) <: BAT.BasicUvStatistics{Float32, FrequencyWeights}

        res = BAT.BasicUvStatistics{T, wK}()
        res = append!(res, data, w)

        @test res.mean[] ≈ mdata
        @test res.var[] ≈ vdata # 100
        @test res.maximum ≈ maxdata
        @test res.minimum ≈ mindata

        numStats = 3
        stats = Array{BAT.BasicUvStatistics{T, wK}}(numStats)
        for i in indices(stats, 1)
            stats[i] = BasicUvStatistics{T, wK}()
        end

        for i in indices(data, 1)
            x = (i % numStats) + 1
            stats[x] = push!(stats[x], data[i], w[i]);
        end

        res = merge(stats...)
        @test res.mean[] ≈ mdata
        @test res.var[] ≈ vdata
        @test res.maximum ≈ maxdata
        @test res.minimum ≈ mindata

        res = BAT.BasicUvStatistics{T, wK}()
        res = append!(res, data)
        @test res.mean[] ≈ mean(data)
        @test res.var[] ≈ cov(data)
    end

end

ext: jl


text:
function read_data_24(inputfilepath, inputfilename)
    nb_dem=0
    nb_prod=0
    nb_mag=0
    nb_noeuds=0
    matrice_S=Array{Float64}[]
    matrice_Q=Array{Float64}[]
    matrice_R=Array{Float64}[]

    filed = open(inputfilepath*"/"*inputfilename)
    while !eof(filed)
        data=readline(filed)
        data=rstrip(data)
        if sizeof(data) >= 1
            out=split(data, " = ")
            if(out[1]=="nbDEMANDES")
                nb_dem = parse(Int, out[2])
            elseif (out[1]=="nbPRODUITS")
                nb_prod = parse(Int, out[2])
            elseif (out[1]=="nbMAGASINS")
                nb_mag = parse(Int, out[2])
            elseif (out[1]=="nbNOEUDS")
                nb_noeuds = parse(Int, out[2])
            elseif (out[1]=="S[m,p]")
                for count in 1:nb_mag
                    data=readline(filed)
                    data=rstrip(data)
                    out=split(data, " ;")
                    out=split(out[1], " ")
                    #println(length(out))
                    #println(out)
                    if length(out) != nb_prod
                        println("erreur de lecture ou donnees, difference entre Smp et nbprod")
                    end
                    tmp=[]
                    for i in 1:nb_prod
                        push!(tmp,parse(Int, out[i]))
                    end
                    push!(matrice_S,tmp)
                end
                #println(matrice_S)
                #println(matrice_S[1][2])
                data=readline(filed)
                data=rstrip(data)
                if(data != "]")
                    println("erreur de lecture ou donnees a la fin de Smp")
                end
            elseif (out[1]=="R[i,j]")
                for count in 1:nb_noeuds
                    data=readline(filed)
                    data=rstrip(data)
                    out=split(data, " ;")
                    out=split(out[1], " ")
                    #println(length(out))
                    #println(out)
                    if length(out) != nb_noeuds
                        println("erreur de lecture ou donnees, difference entre Rij et nb_noeuds")
                    end
                    tmp=[]
                    for i in 1:nb_noeuds
                        push!(tmp,parse(Int, out[i]))
                    end
                    push!(matrice_R,tmp)
                end
                #println(matrice_R)
                #println(matrice_R[1][2])
                data=readline(filed)
                data=rstrip(data)
                if(data != "]")
                    println(data)
                    println("erreur de lecture ou donnees a la fin de Rij")
                end
            elseif (out[1]=="Q[d,p]")
                for count in 1:nb_dem
                    data=readline(filed)
                    data=rstrip(data)
                    out=split(data, " ;")
                    out=split(out[1], " ")
                    if length(out) != nb_prod
                        println("erreur de lecture ou donnees, difference entre Qdp et nbdem")
                    end
                    tmp=[]
                    for i in 1:nb_prod
                        push!(tmp,parse(Int, out[i]))
                    end
                    push!(matrice_Q,tmp)
                end
                #println(matrice_Q)
                #println(matrice_Q[1][2])
                data=readline(filed)
                data=rstrip(data)
                if(data != "]")
                    println("erreur de lecture ou donnees a la fin de Qdp")
                end
            end
        end
    end
    nb_dem, nb_prod, nb_mag, nb_noeuds, matrice_S, matrice_Q, matrice_R
end

max_stars_repo_name: MOUDDENEHamza/ENSEEIHT

max_issues_repo_name: MOUDDENEHamza/ENSEEIHT

max_forks_repo_name: MOUDDENEHamza/ENSEEIHT


text:
using SIIPExamples
using PowerSimulationsDynamics
PSID = PowerSimulationsDynamics
using PowerSystems
using Sundials
using Plots

file_dir = joinpath(
    dirname(dirname(pathof(SIIPExamples))),
    "script",
    "4_PowerSimulationsDynamics_examples",
    "Data",
)
threebus_sys = System(joinpath(file_dir, "threebus_sys.json"));

threebus_sys_dyn = deepcopy(threebus_sys);

#Make a copy of the original system
sys2 = deepcopy(threebus_sys)
#Triplicates the impedance of the line named "BUS 1-BUS 3-i_1"
fault_branches = get_components(ACBranch, sys2)
for br in fault_branches
    if get_name(br) == "BUS 1-BUS 3-i_1"
        br.r = 3 * br.r
        br.x = 3 * br.x
        b_new = (from = br.b.from / 3, to = br.b.to / 3)
        br.b = b_new
    end
end
#Obtain the new Ybus
Ybus_fault = Ybus(sys2).data
#Define Fault: Change of YBus
Ybus_change = NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault, #New YBus
);

#Time span of our simulation
tspan = (0.0, 30.0)

#Define Simulation
sim = PSID.Simulation(
    pwd(), #folder to output results
    threebus_sys, #system
    tspan, #time span
    Ybus_change, #Type of perturbation
)

#Will print the initial states. It also give the symbols used to describe those states.
print_device_states(sim)
#Will export a dictionary with the initial condition values to explore
x0_init = PSID.get_initial_conditions(sim)

#Run the simulation
PSID.execute!(
    sim, #simulation structure
    IDA(), #Sundials DAE Solver
    dtmax = 0.02, #Maximum step size
)

series2 = get_voltagemag_series(sim, 102)
zoom = [
    (series2[1][ix], series2[2][ix]) for
    (ix, s) in enumerate(series2[1]) if (s > 0.90 && s < 1.6)
];

dyn_branch = DynamicBranch(get_component(Branch, threebus_sys_dyn, "BUS 2-BUS 3-i_3"))

add_component!(threebus_sys_dyn, dyn_branch)

#Make a copy of the original system
sys3 = deepcopy(threebus_sys);
#Remove Line "BUS 2-BUS 3-i_3"
remove_component!(Line, sys3, "BUS 2-BUS 3-i_3")
#Triplicates the impedance of the line named "BUS 1-BUS 2-i_1"
fault_branches2 = get_components(Line, sys3)
for br in fault_branches2
    if get_name(br) == "BUS 1-BUS 3-i_1"
        br.r = 3 * br.r
        br.x = 3 * br.x
        b_new = (from = br.b.from / 3, to = br.b.to / 3)
        br.b = b_new
    end
end
#Obtain the new Ybus
Ybus_fault_dyn = Ybus(sys3).data
#Define Fault: Change of YBus
Ybus_change_dyn = PowerSimulationsDynamics.NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault_dyn, #New YBus
)

tspan = (0.0, 30.0)

sim_dyn = PSID.Simulation(
    pwd(), #folder to output results
    threebus_sys_dyn, #system
    tspan, #time span
    Ybus_change_dyn, #Type of perturbation
)

PSID.execute!(
    sim_dyn, #simulation structure
    IDA(), #Sundials DAE Solver
    dtmax = 0.02, #Maximum step size
)

#Will print the initial states. It also give the symbols used to describe those states.
print_device_states(sim_dyn)
#Will export a dictionary with the initial condition values to explore
x0_init_dyn = PSID.get_initial_conditions(sim_dyn)

series2_dyn = get_voltagemag_series(sim_dyn, 102)
zoom_dyn = [
    (series2_dyn[1][ix], series2_dyn[2][ix]) for
    (ix, s) in enumerate(series2_dyn[1]) if (s > 0.90 && s < 1.6)
];

Plots.plot(series2_dyn, label = "V_gen_dyn")
Plots.plot!(series2, label = "V_gen_st", xlabel = "Time [s]", ylabel = "Voltage [pu]")

Plots.plot(zoom_dyn, label = "V_gen_dyn")
Plots.plot!(zoom, label = "V_gen_st", xlabel = "Time [s]", ylabel = "Voltage [pu]")

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
@testset "fantome" begin

    n = 10
    ϵ = 0.1

    rng = MersenneTwister(42)

    p = SquareParticles(rng, n, ϵ)

    q = [[0.5       , 0.5       ],
         [0.20006862, 0.38139954],
         [0.32267135, 0.72687046],
         [0.68044674, 0.78095695],
         [0.77473372, 0.51989917],
         [0.4949439 , 0.23201753],
         [0.74204115, 0.28248482],
         [0.63101615, 0.3650337 ],
         [0.5846235 , 0.64552143],
         [0.45656816, 0.76644607]]

    v = [[ 0.,   1.] ,
         [ 1.,   0.] ,
         [-1.,   0.] ,
         [ 0.,  -1.] ,
         [ 1.,   0.] ,
         [ 0.,   1.] ,
         [ 0.,   1.] ,
         [-1.,   0.] ,
         [-1.,   0.] ,
         [ 0.,  -1.]]

    for i in 1:n
        p.q[i] = q[i]
        p.v[i] = v[i]
    end

    c = PeriodicCollisions( p )

    ref = [0 9 9 1 9 9 9 9 1 1;
           0 0 9 1 9 1 9 1 9 1;
           0 0 0 9 9 9 3 9 9 6;
           0 0 0 0 9 1 1 9 9 9;
           0 0 0 0 0 9 9 2 2 9;
           0 0 0 0 0 0 9 1 9 1;
           0 0 0 0 0 0 0 7 9 9;
           0 0 0 0 0 0 0 0 9 9;
           0 0 0 0 0 0 0 0 0 1;
           0 0 0 0 0 0 0 0 0 0]

    @test ref ≈ c.fantome

end


text:
@testset "PathTracker" begin

    @testset "Regauging" begin
        A = rand(3, 3)
        b = rand(3)

        @polyvar x[1:3]
        F = A * x - b

        tracker, start_sols = pathtracker_startsolutions(F, max_steps=10)
        s = first(start_sols)
        result = track(tracker, s)
        @test result.return_code == :success
        @test result.solution ≈ A \ b atol=1e-6

        @polyvar x y
        result = solve([x^2+y^2-1, 3x+2y-1], threading=false)
        u, v = solutions(result)
        @test norm(u) ≈ 1.0 atol=1e-7
        @test norm(v) ≈ 1.0 atol=1e-7
    end

    @testset "Affine tracking" begin
        @polyvar x a y b
        F = [x^2-a, x*y-a+b]
        p = [a, b]

        # parameter homotopy
        tracker = pathtracker(F, parameters=p, p₁=[1, 0], p₀=[2, 4], affine_tracking=true)
        @test tracker.problem isa Problem{AffineTracking}
        @test affine_tracking(tracker.core_tracker) == true
        @test HC.type_of_x(tracker) == Vector{ComplexF64}
        res = track(tracker, [1, 1])
        @test res.return_code == :success
        @test isa(solution(res), Vector{ComplexF64})
        @test length(solution(res)) == 2
        @test !isprojective(res)
        @test isaffine(res)

        x = @SVector [1.0, 1.0 + 0.0*im]
        tracker, starts = pathtracker_startsolutions(F, x;
                                parameters=p, p₁=[1, 0], p₀=[2, 4], affine_tracking=true)
        @test tracker.problem isa Problem{AffineTracking}
        @test length(starts) == 1
        res = track(tracker, starts[1])
        @test isa(solution(res), Vector{ComplexF64})
        @test length(solution(res)) == 2

        # start target homotopy
        F₁ = subs.(F, Ref([a, b] => [1, 0]))
        F₀ = subs.(F, Ref([a, b] => [2, 4]))
        tracker = pathtracker(F₁, F₀, affine_tracking=true)
        @test tracker.problem isa Problem{AffineTracking}
        res = track(tracker, [1, 1])
        @test res.return_code == :success
        @test isa(solution(res), Vector{ComplexF64})
        @test length(solution(res)) == 2
        @test !isprojective(res)
        @test isaffine(res)

        # total degree
        tracker, starts = pathtracker_startsolutions(equations(katsura(5)), affine_tracking=true)
        res = track(tracker, first(starts))
        @test res.return_code == :success
        @test isa(solution(res), Vector{ComplexF64})
        @test length(solution(res)) == 6
        @test !isprojective(res)
        @test isaffine(res)
    end

    @testset "Details Level" begin
        @polyvar x a y b
        F = [x^2-a, x*y-a+b]
        p = [a, b]

        tracker = pathtracker(F, parameters=p, p₁=[1, 0], p₀=[2, 4])

        res = track(tracker, [1, 1]; details=:minimal)
        test_show_juno(res)
        @test res.residual === nothing

        res = track(tracker, [1, 1]; details=:default)
        @test res.residual !== nothing

        res = track(tracker, [1, 1]; details=:extensive)
        @test res.valuation !== nothing

        @test isnonsingular(res)
    end

    @testset "Parameter homotopy with AbstractSystem" begin
        # affine
        @polyvar x a y b
        F = SPSystem([x^2-a, x*y-a+b]; parameters=[a, b])

        tracker, starts = pathtracker_startsolutions(F, [1.0, 1.0 + 0.0*im],
                        p₁=[1, 0], p₀=[2, 4])
        res = track(tracker, starts[1])
        @test res.return_code == :success

        # affine without explicit flag (numerical degree check fails)
        tracker, starts = pathtracker_startsolutions(F, [1.0, 1.0 + 0.0*im], p₁=[1, 0], p₀=[2, 4])
        res = track(tracker, starts[1])
        @test res.return_code == :success

        # projective
        @polyvar x y z a b
        F = SPSystem([x^2-a*z^2, x*y-a*z^2+b*z^2]; parameters=[a, b])

        tracker, starts = pathtracker_startsolutions(F, [1.0, 1.0 + 0.0*im],
                    p₁=[1, 0], p₀=[2, 4])
        res = track(tracker, starts[1])
        @test res.return_code == :success
        @test isa(res.solution, ProjectiveVectors.PVector)

        # projective with homvar
        tracker, starts = pathtracker_startsolutions(F, [1.0, 1.0 + 0.0*im],
                    p₁=[1, 0], p₀=[2, 4], homvar=z)
        @test isa(tracker.core_tracker.affine_patch, EmbeddingPatch)
        res = track(tracker, starts[1])
        @test res.return_code == :success
        @test isa(res.solution, Vector)
    end

    @testset "Parameter homotopy, change parameters" begin
        @polyvar x y z p[1:3]

        F = [
            x + 3 + 2y + 2y^2 - p[1],
            (x - 2 + 5y)*z + 4 - p[2] * z,
            (x + 2 + 4y)*z + 5 - p[3] * z
        ]
        # Generate generic parameters by sampling complex numbers from the normal distribution
        p₀ = randn(ComplexF64, 3)
        # Substitute p₀ for p
        F_p₀ = subs(F, p => p₀)
        # Compute all solutions for F_p₀
        result_p₀ = solve(F_p₀; affine_tracking=true)

        # Let's store the generic solutions
        S_p₀ = solutions(nonsingular(result_p₀))
        # Construct the PathTracker
        tracker = pathtracker(F; parameters=p, generic_parameters=p₀)

        q = randn(3)
        result = track(tracker, first(S_p₀); target_parameters=q)
        @test issuccess(result)

        result = track(tracker, first(S_p₀); start_parameters=p₀, target_parameters=q)
        @test issuccess(result)

        # Construct the PathTracker
        tracker = pathtracker(F; parameters=p, generic_parameters=p₀)
        set_parameters!(tracker; target_parameters=q)
        result2 = track(tracker, first(S_p₀))
        @test result.solution ≈ result2.solution atol=1e-7
    end

    @testset "Real Jacobian" begin
        @polyvar v w
        @test nfinite(solve([v - 2w , -v + 2w, v + w - 3])) == 1
        @test nfinite(solve([v - 2w , -v + 2w, v + w - 3], affine_tracking=false)) == 1
    end
end


text:
using Compat.Test
import Roots.fzero

## Test `fzero` interface to `find_zero`
## test `fzeros` interface for functions

orders = [0,1,2,5,8,16]

### derivative free
fn, xstar, x0, br = x -> sin(x) - x/2, 1.89549426703398094714, 2.0, [pi/2, pi]
@test fzero(fn, x0)  ≈ xstar
for o in orders
    @test fzero(fn, x0, order=o)  ≈ xstar
end
### bisection
@test fzero(fn, br)  ≈ xstar
@test fzero(x->5-x,5,10)==5

### test tolerances
fn, xstar, x0, br = x -> x^5 - x - 1, 1.1673039782614187, 1.0, [1.0, 2.0]
@test fzero(fn, x0, order=1)  ≈ xstar

@test_throws Roots.ConvergenceFailed fzero(fn, x0, order=1, maxevals=2)
#@test abs(fzero(fn, x0, order=1, ftol=1e-2) - xstar) > 1e-5
#@test abs(fzero(fn, x0, order=1, xtol=1e-2) - xstar) > 1e-10

@test abs(find_zero(fn, br, Roots.A42(), xatol=1e-1) - xstar) > 1e-5

## various tests
## issue #29, basically f(a) or f(b) so big we get NaN for guess from bracket.
f =  x -> x + exp(x)
@test fzero(f, [-1e6, 1e6])  ≈ -0.5671432904097838

f =  x -> 1/x - 1
@test fzero(f, [0, 2])  ≈ 1.0

## test infinite range
@test fzero(x -> x, [-Inf, Inf])  ≈ 0.0

##################################################
## fzeros function
rts = 1:5
@test all((abs.(fzeros(x -> prod([x-r for r in rts]),0,10)) .- collect(1:5)) .<= 1e-15)

fn = x -> cos(10*pi*x)
@test length(fzeros(fn, 0, 1)) == 10

### issue with fzeros and roots near 'b'
@test 0 <  maximum(fzeros(x -> sin(x) - 1/1000*x, 0, pi)) < pi


text:
using EcologicalNetworks
using EcologicalNetworksPlots
using Plots
using Test

using Random
Random.seed!(42);

# Prepare the output
figpath = joinpath(pwd(), "gallery")
ispath(figpath) && rm(figpath; recursive=true)
ispath(figpath) || mkdir(figpath)

Umod = web_of_life("M_PA_003")
Unes = web_of_life("M_SD_033")

@info "Bipartite -- nested"
for al in [true, false], re in [true, false]
    I = initial(BipartiteInitialLayout, Unes)
    position!(NestedBipartiteLayout(al, re, 0.4), I, Unes)
    plot(I, Unes, aspectratio=1)
    scatter!(I, Unes, bipartite=true)
    savefig(joinpath(figpath, "bip_nest_al_$(al)_re_$(re).png"))
    @test isfile(joinpath(figpath, "bip_nest_al_$(al)_re_$(re).png"))
end

@info "Bipartite -- circular"
I = initial(CircularInitialLayout, Unes)
position!(CircularLayout(), I, Unes)
plot(I, Unes, aspectratio=1, framestyle=:grid, legend=true)
savefig(joinpath(figpath, "bip_circular_plot.png"))
scatter!(I, Unes, bipartite=true)
savefig(joinpath(figpath, "bip_circular_full.png"))
@test isfile(joinpath(figpath, "bip_circular_full.png"))

@info "Circular layer"
I = initial(CircularInitialLayout, Umod)
position!(CircularLayout(), I, Umod)
plot(I, Umod, aspectratio=1, framestyle=:grid, legend=false)
savefig(joinpath(figpath, "bip_circular_plot_2.png"))
scatter!(I, Umod, bipartite=true, framestyle=:box, legend=true)
savefig(joinpath(figpath, "bip_circular_full_2.png"))
@test isfile(joinpath(figpath, "bip_circular_full_2.png"))

Unes = web_of_life("M_SD_033")
I = initial(BipartiteInitialLayout, Unes)
position!(NestedBipartiteLayout(0.4), I, Unes)
plot(I, Unes, aspectratio=1)
scatter!(I, Unes, bipartite=true, nodefill=degree(Unes))

Fweb = simplify(nz_stream_foodweb()[5])
I = initial(FoodwebInitialLayout, Fweb)
for step in 1:300
  position!(ForceDirectedLayout(2.5, 1.5; gravity=0.7), I, Fweb)
end
tl = trophic_level(Fweb)
for s in species(Fweb)
  I[s].y = tl[s]
end
plot(I, Fweb)
scatter!(I, Fweb)

@info "Heatmap"
heatmap(Unes)
savefig(joinpath(figpath, "heatmap.png"))
@test isfile(joinpath(figpath, "heatmap.png"))


text:
# Python buffer protocol: this is a NumPy-array-like facility to get raw
# pointers to contiguous buffers of data underlying other objects,
# with support for describing multiple dimensions, strides, etc.
#     (thanks to @jakebolewski for his work on this)

#############################################################################
# mirror of Py_buffer struct in Python Include/object.h

struct Py_buffer
    buf::Ptr{Cvoid}
    obj::PyPtr
    len::Cssize_t
    itemsize::Cssize_t

    readonly::Cint
    ndim::Cint
    format::Ptr{Cchar}
    shape::Ptr{Cssize_t}
    strides::Ptr{Cssize_t}
    suboffsets::Ptr{Cssize_t}

    # some opaque padding fields to account for differences between
    # Python versions (the structure changed in Python 2.7 and 3.3)
    internal0::Ptr{Cvoid}
    internal1::Ptr{Cvoid}
    internal2::Ptr{Cvoid}
end

mutable struct PyBuffer
    buf::Py_buffer
    PyBuffer() = begin
        b = new(Py_buffer(C_NULL, PyPtr_NULL, 0, 0,
                          0, 0, C_NULL, C_NULL, C_NULL, C_NULL,
                          C_NULL, C_NULL, C_NULL))
        finalizer(pydecref, b)
        return b
    end
end

"""
`pydecref(o::PyBuffer)`
Release the reference to buffer `o`
N.b. As per https://docs.python.org/3/c-api/buffer.html#c.PyBuffer_Release,
It is an error to call this function on a PyBuffer that was not obtained via
the python c-api function `PyObject_GetBuffer()`, unless o.obj is a PyPtr(C_NULL)
"""
function pydecref(o::PyBuffer)
    # note that PyBuffer_Release sets o.obj to NULL, and
    # is a no-op if o.obj is already NULL
    # TODO change to `Ref{PyBuffer}` when 0.6 is dropped.
    _finalized[] || ccall(@pysym(:PyBuffer_Release), Cvoid, (Any,), o)
    o
end

#############################################################################
# Array-like accessors for PyBuffer.

Base.ndims(b::PyBuffer) = UInt(b.buf.ndim)

# from the Python docs: If shape is NULL as a result of a PyBUF_SIMPLE
# or a PyBUF_WRITABLE request, the consumer must disregard itemsize
# and assume itemsize == 1.
Base.length(b::PyBuffer) = b.buf.shape == C_NULL ? Int(b.buf.len) : Int(div(b.buf.len, b.buf.itemsize))

Base.sizeof(b::PyBuffer) = Int(b.buf.len)
Base.pointer(b::PyBuffer) = b.buf.buf

function Base.size(b::PyBuffer)
    b.buf.ndim <= 1 && return (length(b),)
    @assert b.buf.shape != C_NULL
    return tuple(Int[unsafe_load(b.buf.shape, i) for i=1:b.buf.ndim]...)
end
# specialize size(b, d) for efficiency (avoid tuple construction)
function Base.size(b::PyBuffer, d::Integer)
    d > b.buf.ndim && return 1
    d < 0 && throw(BoundsError())
    b.buf.ndim <= 1 && return length(b)
    @assert b.buf.shape != C_NULL
    return Int(unsafe_load(b.buf.shape, d))
end

# stride in bytes for i-th dimension
function Base.stride(b::PyBuffer, d::Integer)
    d > b.buf.ndim && return length(b) # as in base
    d < 0 && throw(BoundsError())
    if b.buf.strides == C_NULL
        if b.buf.ndim == 1
            return b.buf.shape == C_NULL ? 1 : Int(b.buf.itemsize)
        else
            error("unknown buffer strides")
        end
    end
    return Int(unsafe_load(b.buf.strides, d))
end

# Strides in bytes
Base.strides(b::PyBuffer) = ((stride(b,i) for i in 1:b.buf.ndim)...,)

# TODO change to `Ref{PyBuffer}` when 0.6 is dropped.
iscontiguous(b::PyBuffer) =
    1 == ccall((@pysym :PyBuffer_IsContiguous), Cint,
               (Any, Cchar), b, 'A')

#############################################################################
# pybuffer constant values from Include/object.h
const PyBUF_SIMPLE    = convert(Cint, 0)
const PyBUF_WRITABLE  = convert(Cint, 0x0001)
const PyBUF_FORMAT    = convert(Cint, 0x0004)
const PyBUF_ND        = convert(Cint, 0x0008)
const PyBUF_STRIDES        = convert(Cint, 0x0010) | PyBUF_ND
const PyBUF_C_CONTIGUOUS   = convert(Cint, 0x0020) | PyBUF_STRIDES
const PyBUF_F_CONTIGUOUS   = convert(Cint, 0x0040) | PyBUF_STRIDES
const PyBUF_ANY_CONTIGUOUS = convert(Cint, 0x0080) | PyBUF_STRIDES
const PyBUF_INDIRECT       = convert(Cint, 0x0100) | PyBUF_STRIDES
const PyBUF_ND_STRIDED    = Cint(PyBUF_WRITABLE | PyBUF_FORMAT | PyBUF_ND |
                                 PyBUF_STRIDES)
const PyBUF_ND_CONTIGUOUS = PyBUF_ND_STRIDED | PyBUF_ANY_CONTIGUOUS

# construct a PyBuffer from a PyObject, if possible
function PyBuffer(o::Union{PyObject,PyPtr}, flags=PyBUF_SIMPLE)
    return PyBuffer!(PyBuffer(), o, flags)
end

function PyBuffer!(b::PyBuffer, o::Union{PyObject,PyPtr}, flags=PyBUF_SIMPLE)
    # TODO change to `Ref{PyBuffer}` when 0.6 is dropped.
    pydecref(b) # ensure b is properly released
    @pycheckz ccall((@pysym :PyObject_GetBuffer), Cint,
                     (PyPtr, Any, Cint), o, b, flags)
    return b
end

"""
`isbuftype(o::Union{PyObject,PyPtr})`
Returns true if the python object `o` supports the buffer protocol as a strided
array. False if not.
"""
function isbuftype(o::Union{PyObject,PyPtr})
    # PyObject_CheckBuffer is defined in a header file here: https://github.com/python/cpython/blob/ef5ce884a41c8553a7eff66ebace908c1dcc1f89/Include/abstract.h#L510
    # so we can't access it easily. It basically just checks if PyObject_GetBuffer exists
    # So we'll just try call PyObject_GetBuffer and check for success/failure
    b = PyBuffer()
    ret = ccall((@pysym :PyObject_GetBuffer), Cint,
                     (PyPtr, Any, Cint), o, b, PyBUF_ND_STRIDED)
    if ret != 0
        pyerr_clear()
    else
        # handle pointer types
        T, native_byteorder = array_format(b)
        T <: Ptr && (ret = 1)
    end
    return ret == 0
end

#############################################################################

# recursive function to write buffer dimension by dimension, starting at
# dimension d with the given pointer offset (in bytes).
function writedims(io::IO, b::PyBuffer, offset, d)
    n = 0
    s = stride(b, d)
    if d < b.buf.ndim
        for i = 1:size(b,d)
            n += writedims(io, b, offset, d+1)
            offset += s
        end
    else
        @assert d == b.buf.ndim
        p = convert(Ptr{UInt8}, pointer(b)) + offset
        for i = 1:size(b,d)
            # would be nicer not to write this one byte at a time,
            # but the alternative seems to be to create an Array
            # object on each loop iteration.
            for j = 1:b.buf.itemsize
                n += write(io, unsafe_load(p))
                p += 1
            end
            p += s
        end
    end
    return n
end

function Base.write(io::IO, b::PyBuffer)
    b.buf.obj != C_NULL || error("attempted to write NULL buffer")

    if iscontiguous(b)
        # (note that 0-dimensional buffers are always contiguous)
        return write(io, pointer_to_array(convert(Ptr{UInt8}, pointer(b)),
                                          sizeof(b)))
    else
        return writedims(io, b, 0, 1)
    end
end

# ref: https://github.com/numpy/numpy/blob/v1.14.2/numpy/core/src/multiarray/buffer.c#L966

const standard_typestrs = Dict{String,DataType}(
                           "?"=>Bool,        "P"=>Ptr{Cvoid},
                           "b"=>Int8,        "B"=>UInt8,
                           "h"=>Int16,       "H"=>UInt16,
                           "i"=>Int32,       "I"=>UInt32,
                           "l"=>Int32,       "L"=>UInt32,
                           "q"=>Int64,       "Q"=>UInt64,
                           "e"=>Float16,     "f"=>Float32,
                           "d"=>Float64,     "g"=>Nothing, # Float128?
                           # `Nothing` indicates no equiv Julia type
                           "Z8"=>ComplexF32, "Z16"=>ComplexF64,
                           "Zf"=>ComplexF32, "Zd"=>ComplexF64)

const native_typestrs = Dict{String,DataType}(
                           "?"=>Bool,        "P"=>Ptr{Cvoid},
                           "b"=>Int8,        "B"=>UInt8,
                           "h"=>Cshort,      "H"=>Cushort,
                           "i"=>Cint,        "I"=>Cuint,
                           "l"=>Clong,       "L"=>Culong,
                           "q"=>Clonglong,   "Q"=>Culonglong,
                           "e"=>Float16,     "f"=>Cfloat,
                           "d"=>Cdouble,     "g"=>Nothing, # Float128?
                           # `Nothing` indicates no equiv Julia type
                           "Z8"=>ComplexF32, "Z16"=>ComplexF64,
                           "Zf"=>ComplexF32, "Zd"=>ComplexF64)

const typestrs_native =
    Dict{DataType, String}(zip(values(native_typestrs), keys(native_typestrs)))

get_format_str(pybuf::PyBuffer) = unsafe_string(convert(Ptr{UInt8}, pybuf.buf.format))

function array_format(pybuf::PyBuffer)
    # a NULL-terminated format-string ... indicating what is in each element of memory.
    # TODO: handle more cases: https://www.python.org/dev/peps/pep-3118/#additions-to-the-struct-string-syntax
    # refs: https://github.com/numpy/numpy/blob/v1.14.2/numpy/core/src/multiarray/buffer.c#L966
    #       https://github.com/numpy/numpy/blob/v1.14.2/numpy/core/_internal.py#L490
    #       https://docs.python.org/3/library/struct.html#byte-order-size-and-alignment

    # "NULL implies standard unsigned bytes ("B")" --pep 3118
    pybuf.buf.format == C_NULL && return UInt8, true

    fmt_str = get_format_str(pybuf)
    native_byteorder = true
    use_native_sizes = true
    type_start_idx = 1
    if length(fmt_str) > 1
        type_start_idx = 2
        if fmt_str[1] == '@' || fmt_str[1] == '^'
            # defaults to native_byteorder: true, use_native_sizes: true
        elseif fmt_str[1] == '<'
            native_byteorder = ENDIAN_BOM == 0x04030201
            use_native_sizes = false
        elseif fmt_str[1] == '>' || fmt_str =='!'
            native_byteorder = ENDIAN_BOM == 0x01020304
            use_native_sizes = false
        elseif fmt_str[1] == '='
            use_native_sizes = false
        elseif fmt_str[1] == 'Z'
            type_start_idx = 1
        else
            error("Unsupported format string: \"$fmt_str\"")
        end
    end
    strs2types = use_native_sizes ? native_typestrs : standard_typestrs
    strs2types[fmt_str[type_start_idx:end]], native_byteorder
end

#############################################################################


text:
# Compute total squared displacement of a set of coords
function sd{T<:AbstractFloat}(C::AbstractArray{T})
    t = zero(T)
    @inbounds @simd for i in eachindex(C)
        t += C[i] ^ 2
    end
    t
end

# Compute total squared displacement between two sets of coords
function sd{T<:AbstractFloat}(C1::AbstractArray{T}, C2::AbstractArray{T})
    t = zero(T)
    @inbounds @simd for i in eachindex(C1)
        t += (C1[i] - C2[i])^2
    end
    t
end

# Compute total squared displacement between two sets of coords in a periodic box
function sd{T<:AbstractFloat}(bw::T, C1::AbstractArray{T}, C2::AbstractArray{T})
    t = zero(T)
    @inbounds @simd for i in size(C1, 2)
        t += distancesq(bw, C1[1, i], C1[2, i], C1[3, i], C2[1, i], C2[2, i], C2[3, i])
    end
    t
end

# Compute mean squared displacements
msd{T<:AbstractFloat}(C::AbstractArray{T}) = sd(C) / length(C)
msd{T<:AbstractFloat}(C1::AbstractArray{T}, C2::AbstractArray{T}) = sd(C1, C2) / length(C1)
msd{T<:AbstractFloat}(bw::T, C1::AbstractArray{T}, C2::AbstractArray{T}) = sd(bw, C1, C2) / length(C1)

# Compute root mean squared deviations
rmsd{T<:AbstractFloat}(C::AbstractArray{T}) = sqrt(msd(C))
rmsd{T<:AbstractFloat}(C1::AbstractArray{T}, C2::AbstractArray{T}) = sqrt(msd(C1, C2))
rmsd{T<:AbstractFloat}(bw::T, C1::AbstractArray{T}, C2::AbstractArray{T}) = sqrt(msd(bw, C1, C2))

# Compute the magnitude of a vector
mag{T<:AbstractFloat}(C::AbstractVector{T}) = sqrt(sd(C))

# Wrap a set of coords into a periodic box
function wrapcoords!{T<:AbstractFloat}(bw::T, C::AbstractMatrix{T})
    hbw = bw / 2
    @inbounds @simd for i in 1:size(C, 2)
        if C[1] < -hbw C[1] += bw elseif C[1] > hbw C[1] -= bw end
        if C[2] < -hbw C[2] += bw elseif C[2] > hbw C[2] -= bw end
        if C[3] < -hbw C[3] += bw elseif C[3] > hbw C[3] -= bw end
    end
end

# Unwrap each coord in C into its image in IC
function unwrapcoords!{T<:AbstractFloat}(bw::T, C::AbstractMatrix{T}, IC::AbstractMatrix{T})
    hbw = bw / 2
    @inbounds @simd for i in 1:size(C, 2)
        if (C[1, i] - IC[1, i]) < -hbw C[1, i] += bw elseif C[1, i] - IC[1, i] > hbw C[1, i] -= bw end
        if (C[2, i] - IC[2, i]) < -hbw C[2, i] += bw elseif C[2, i] - IC[2, i] > hbw C[2, i] -= bw end
        if (C[3, i] - IC[3, i]) < -hbw C[3, i] += bw elseif C[3, i] - IC[3, i] > hbw C[3, i] -= bw end
    end
end

# Compute the squared distance between two coords
function distancesq{T<:AbstractFloat}(bw::T, x1::T, y1::T, z1::T, x2::T, y2::T, z2::T)
    hbw = bw / 2
    xd = x1 - x2
    yd = y1 - y2
    zd = z1 - z2
    if xd < -hbw xd += bw elseif xd > hbw xd -= bw end
    if yd < -hbw yd += bw elseif yd > hbw yd -= bw end
    if zd < -hbw zd += bw elseif zd > hbw zd -= bw end

    return xd*xd + yd*yd + zd*zd
end

function distancesq{T<:AbstractFloat}(xbw::T, ybw::T, zbw::T, x1::T, y1::T, z1::T, x2::T, y2::T, z2::T)
    xhbw = xbw / 2
    yhbw = ybw / 2
    zhbw = zbw / 2
    xd = x1 - x2
    yd = y1 - y2
    zd = z1 - z2
    if xd < -xhbw xd += xbw elseif xd > xhbw xd -= xbw end
    if yd < -yhbw yd += ybw elseif yd > yhbw yd -= ybw end
    if zd < -zhbw zd += zbw elseif zd > zhbw zd -= zbw end

    return xd*xd + yd*yd + zd*zd
end

distancesq{T<:AbstractFloat}(bw::T, C::AbstractMatrix{T}, i::Int, j::Int) = distancesq(bw, C[1, i], C[2, i], C[3, i], C[1, j], C[2, j], C[3, j])
distancesq{T<:AbstractFloat}(bw::T, C::AbstractMatrix{T}, i::Int, x::T, y::T, z::T) = distancesq(bw, C[1, i], C[2, i], C[3, i], x, y, z)
distancesq{T<:AbstractFloat}(bw::T, C1::AbstractMatrix{T}, C2::AbstractMatrix{T}, i1::Int, i2::Int) = distancesq(bw, C1[1, i1], C1[2, i1], C1[3, i1], C2[1, i2], C2[2, i2], C2[3, i2])

distancesq{T<:AbstractFloat}(xbw::T, ybw::T, zbw::T, C::AbstractMatrix{T}, i::Int, j::Int) = distancesq(xbw, ybw, zbw, C[1, i], C[2, i], C[3, i], C[1, j], C[2, j], C[3, j])

# Compute the distances between a set of coords and a particular coord
function distances!{T<:AbstractFloat}(bw::T, D::AbstractVector{T}, C::AbstractMatrix{T}, x::T, y::T, z::T)
    hbw = bw / 2
    @inbounds @simd for i in 1:size(D, 2)
        xd = C[1, i] - x
        yd = C[2, i] - y
        zd = C[3, i] - z
        if xd < -hbw xd += bw elseif xd > hbw xd -= bw end
        if yd < -hbw yd += bw elseif yd > hbw yd -= bw end
        if zd < -hbw zd += bw elseif zd > hbw zd -= bw end
        D[i] = sqrt(xd^2 + yd^2 + zd^2)
    end
end

function distances!{T<:AbstractFloat}(bw::T, D::AbstractArray{T, 3}, ci::Int, C::AbstractMatrix{T}, indices::AbstractMatrix{Int}, i::Int)
    hbw = bw / 2
    @inbounds begin
        x = C[1, i]
        y = C[2, i]
        z = C[3, i]
        @simd for j in 1:size(indices, 1)
            xd = C[1, indices[j, ci]] - x
            yd = C[2, indices[j, ci]] - y
            zd = C[3, indices[j, ci]] - z
            if xd < -hbw xd += bw elseif xd > hbw xd -= bw end
            if yd < -hbw yd += bw elseif yd > hbw yd -= bw end
            if zd < -hbw zd += bw elseif zd > hbw zd -= bw end
            D[1, j, ci] = xd
            D[2, j, ci] = yd
            D[3, j, ci] = zd
        end
    end
end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
# Docstring signature syntax highlighting tests.
module HighlightSig
    using Test
    import Markdown
    import Documenter.Expanders: highlightsig!

    @testset "highlightsig!" begin
        s = """
                foo(bar::Baz)
            ---
                foo(bar::Baz)
            """
        original = Markdown.parse(s)
        md = Markdown.parse(s)
        highlightsig!(md)
        @test isempty(original.content[1].language)
        @test md.content[1].language == "julia"
        @test original.content[end].language == md.content[end].language

        s = """
            ```lang
             foo(bar::Baz)
            ```
            """
        original = Markdown.parse(s)
        md = Markdown.parse(s)
        highlightsig!(md)
        @test original == md
    end
end


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

struct SecurityRulesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Get the specified default network security rule.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: defaultSecurityRuleName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SecurityRule
"""
function _swaggerinternal_defaultSecurityRulesGet(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, defaultSecurityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SecurityRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}/defaultSecurityRules/{defaultSecurityRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "defaultSecurityRuleName", defaultSecurityRuleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function defaultSecurityRulesGet(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, defaultSecurityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_defaultSecurityRulesGet(_api, resourceGroupName, networkSecurityGroupName, defaultSecurityRuleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function defaultSecurityRulesGet(_api::SecurityRulesApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, defaultSecurityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_defaultSecurityRulesGet(_api, resourceGroupName, networkSecurityGroupName, defaultSecurityRuleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all default security rules in a network security group.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SecurityRuleListResult
"""
function _swaggerinternal_defaultSecurityRulesList(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SecurityRuleListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}/defaultSecurityRules", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function defaultSecurityRulesList(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_defaultSecurityRulesList(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function defaultSecurityRulesList(_api::SecurityRulesApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_defaultSecurityRulesList(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Creates or updates a security rule in the specified network security group.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: securityRuleName::String (required)
Param: securityRuleParameters::SecurityRule (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SecurityRule
"""
function _swaggerinternal_securityRulesCreateOrUpdate(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, securityRuleParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", SecurityRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}/securityRules/{securityRuleName}", ["azure_auth"], securityRuleParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "securityRuleName", securityRuleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function securityRulesCreateOrUpdate(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, securityRuleParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesCreateOrUpdate(_api, resourceGroupName, networkSecurityGroupName, securityRuleName, securityRuleParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function securityRulesCreateOrUpdate(_api::SecurityRulesApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, securityRuleParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesCreateOrUpdate(_api, resourceGroupName, networkSecurityGroupName, securityRuleName, securityRuleParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified network security rule.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: securityRuleName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_securityRulesDelete(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}/securityRules/{securityRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "securityRuleName", securityRuleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function securityRulesDelete(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesDelete(_api, resourceGroupName, networkSecurityGroupName, securityRuleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function securityRulesDelete(_api::SecurityRulesApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesDelete(_api, resourceGroupName, networkSecurityGroupName, securityRuleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Get the specified network security rule.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: securityRuleName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SecurityRule
"""
function _swaggerinternal_securityRulesGet(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SecurityRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}/securityRules/{securityRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "securityRuleName", securityRuleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function securityRulesGet(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesGet(_api, resourceGroupName, networkSecurityGroupName, securityRuleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function securityRulesGet(_api::SecurityRulesApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, securityRuleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesGet(_api, resourceGroupName, networkSecurityGroupName, securityRuleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all security rules in a network security group.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SecurityRuleListResult
"""
function _swaggerinternal_securityRulesList(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SecurityRuleListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}/securityRules", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function securityRulesList(_api::SecurityRulesApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesList(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function securityRulesList(_api::SecurityRulesApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_securityRulesList(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export defaultSecurityRulesGet, defaultSecurityRulesList, securityRulesCreateOrUpdate, securityRulesDelete, securityRulesGet, securityRulesList


text:
include("abstractfeatures.jl")

"""
    PairFeatureDescriptor{A}

A feature object that encodes features associated with pairs of atoms in a structure that may or may not be bonded to each other. For bond features, see `BondFeatureDescriptor`.

Type parameter represents the structure representation(s) from which this feature descriptor is able to compute features.

## Fields
- `name::String`: the name of the feature
- `compute_f`: function that takes in a structure of type `A` and returns values of this feature for every pair of atoms in the structure organized in a matrix
- `encoder_decoder::AbstractCodec`: codec that encodes/decodes values of this feature
- `categorical::Bool`: flag for whether the feature is categorical or continuous-valued
- `encodable_elements::Vector{String}`: list of elements (by symbol) that can be encoded by this feature
"""
struct PairFeatureDescriptor{A,C<:AbstractCodec} <: AbstractPairFeatureDescriptor
    name::String
    compute_f::Function
    encoder_decoder::C
    categorical::Bool
    encodable_elements::Vector{String}
end

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", fd::PairFeatureDescriptor{A}) where {A}
    st = "$(typeof(fd)) $(fd.name):\n   categorical: $(fd.categorical)\n   works on: $(A)"
    print(io, st)
end

encodable_elements(fd::PairFeatureDescriptor) = fd.encodable_elements

function get_value(pfd::PairFeatureDescriptor{A}, a::AbstractAtoms{<:A}) where {A}
    @assert all([el in encodable_elements(pfd) for el in elements(a)]) "Feature $(pfd.name) cannot encode some element(s) in this structure!"
    pfd.compute_f(a.structure)
end


text:
module Events

import ..Core
import ..Currencies
using Dates

#
# TODO: this content could be generated using macros
#

#
# :AMORT -> Amortization
# :INTEREST -> Interest payment
# :CAP -> Capitalization
# :PAYMENT -> Generic cashflow event
#

const EVN_AMORT = :AMORT
const EVN_INTEREST = :INTEREST
const EVN_CAPITALIZATION = :CAP
const EVN_PAYMENT = :PAYMENT

const EVENT_SYMS = [ EVN_AMORT, EVN_INTEREST, EVN_CAPITALIZATION, EVN_PAYMENT ]

is_valid(evn::Symbol) = evn ∈ EVENT_SYMS

function Amort(contract::C) where {C<:Core.WhenAt}
    return Core.FixedIncomeEvent(EVN_AMORT, contract)
end

function Amort(maturity::Date, value::T, currency::Currencies.Currency) where {T<:Union{Number, Core.Observable}}
    return Core.FixedIncomeEvent(EVN_AMORT, maturity, value, currency)
end

function Interest(contract::C) where {C<:Core.WhenAt}
    return Core.FixedIncomeEvent(EVN_INTEREST, contract)
end

function Interest(maturity::Date, value::T, currency::Currencies.Currency) where {T<:Union{Number, Core.Observable}}
    return Core.FixedIncomeEvent(EVN_INTEREST, maturity, value, currency)
end

function Payment(contract::C) where {C<:Core.WhenAt}
    return Core.FixedIncomeEvent(EVN_PAYMENT, contract)
end

function Payment(maturity::Date, value::T, currency::Currencies.Currency) where {T<:Union{Number, Core.Observable}}
    return Core.FixedIncomeEvent(EVN_PAYMENT, maturity, value, currency)
end

end


text:
using EzXML

"""
Read the contents of the document as a string.

  * Header text will be clustered
  * Paragraphs will be separated by new lines.
  * Footnotes will be clustered

**NOTE**: This will internally mutate `docx` but multiple calls to this
function with the same docx will return the same string.
"""
function read(docx::Document, ::Type{String})::String
    plaintext = ""

    for (i, file) in enumerate(docx._files)
        if file.extension != "xml"
            continue
        end

        str = strip(_extract_text(read!(file, String)))
        if !isempty(str)
            plaintext *= str

            if i < length(docx._files)
                plaintext *= "\n"
            end
        end
    end

    plaintext
end

"""
A type that indicates to return the contents as XML strings
"""
abstract type XML end

"""
    read(docx, XML)

Reads the contents of the document and returns them as a plain string

  * Header text will be clustered
  * Paragraphs will be separated by new lines.
  * Footnotes will be clustered
"""
function read(docx::Document, ::Type{XML})::String
    xml_nodes = ""
    for (i, file) in enumerate(docx._files)
        if file.extension == "xml"
            xml_nodes *= read!(file, String)
        end
    end

    xml_nodes
end

"""
Reads the docfile contents as a string
"""
function read!(docx_file::DocxFile, ::Type{String})::String
    if docx_file.string_contents !== nothing
        return docx_file.string_contents::String
    end

    str = Base.read(docx_file._zipFile, String)
    docx_file.string_contents = str
    return str
end

"""
Reads the docfile contents as a string
"""
function read_binary!(docx_file::DocxFile)::Vector{UInt8}
    if docx_file.binary_contents !== nothing
        return docx_file.binary_contents
    end

    bin = read(docx_file._zipFile, UInt8)
    docx_file.binary_contents
    return bin
end

"""
Extracts text elements from the given XML structure
"""
function _extract_text(xml_string::String)::String
    text = ""

    doc = parsexml(xml_string)
    root_elem = root(doc)

    elems = _flatten(root_elem)
    for (i, elem) in enumerate(elems)
        if elem.name == "p"
            text *= "\n"
        end

        if _is_docx_text_node(elem)
            text *= nodecontent(elem)
            text = replace(text, "\u00a0" => "\u0020")
        end
    end

    text
end

const word_processing_xmlns_2006 = "http://schemas.openxmlformats.org/wordprocessingml/2006/main"

const text_content_nodes = [
    "t"
]

"""
Determines whether the XML node element contains text data or not.
"""
_is_docx_text_node(elem) =
    elem.namespace == word_processing_xmlns_2006 && nodename(elem) in text_content_nodes

"""
Creates an array of all the XML nodes by pushing into an array as the root node and children are
traversed.
"""
function _flatten(node::EzXML.Node)::Array{EzXML.Node}
    elems = [node]
    for elem in eachelement(node)
        for subelem in _flatten(elem)
            push!(elems, subelem)
        end
    end

    elems
end


text:
"Log( data::SingleCellExpr )

Log the expression data using log(x+1), not log2 or log10
This will not touch the 0 or -1 values

"

function Log( data::SingleCellExpr )
  function logthis(x::SparseVector)
    ok = find( x .> 0 );
    x[ok] .= log(x[ok]+1);
    return x
  end
  rws = 1:size(data.data)[1]
  for  i =  1:size(data.data)[2]
    data.data[rws,i] = logthis( data.data[rws,i] )
  end
  return data
end


text:
using Test

include("card-generator.jl")

@testset "Speedrunning 101" begin
    @test generate_card("»Speedrunning 101« – Sasha Duda Krall (they/them), from GDQU. Start: 13:00, Q&A: 13:20, End: 13:30") == """
        - Our next speaker is Sasha Duda Krall, from GDQU
        - Their talk is called »Speedrunning 101«
        - They will answer your questions in the Q&A session at the end of the talk, starting at 13:20

        13:00 - 13:20 - 13:30
        """
end

@testset "How I learned to say Farewell" begin
    @test generate_card("»How I learned to say Farewell« – Madeline (she/her). Start: 13:40, Q&A: 14:00, End: 14:05") == """
        - Our next speaker is Madeline
        - Her talk is called »How I learned to say Farewell«
        - She will answer your questions in the Q&A session at the end of the talk, starting at 14:00

        13:40 - 14:00 - 14:05
        """
end

@testset "How To Find a Good Title For Your Conference Talk" begin
    @test generate_card("»How To Find a Good Title For Your Conference Talk« – Vítězslav Kruse (he/him). Start: 14:15, End: 14:40") == """
        - Our next speaker is Vítězslav Kruse
        - His talk is called »How To Find a Good Title For Your Conference Talk«
        - There will not be a Q&A session

        14:15 - NO Q&A - 14:40
        """
end

@testset "Why Emoji Matter❣" begin
    @test generate_card("»Why Emoji Matter❣« – Ash van Amelsvoort, from the University of 🧬🧪⚛. Start: 14:50, Q&A: 15:05, End: 15:10") == """
        - Our next speaker is Ash van Amelsvoort, from the University of 🧬🧪⚛
        - Ash's talk is called »Why Emoji Matter❣«
        - Ash will answer your questions in the Q&A session at the end of the talk, starting at 15:05

        14:50 - 15:05 - 15:10
        """
end

@testset "»Can dogs look up?" begin
    @test generate_card("»Can dogs look up?« – Kira \"k1ralli\" Sørensen. Start: 20:50, Q&A: 21:05, End: 21:10") == """
        - Our next speaker is Kira "k1ralli" Sørensen
        - Kira's talk is called »Can dogs look up?«
        - Kira will answer your questions in the Q&A session at the end of the talk, starting at 21:05

        20:50 - 21:05 - 21:10
        """
end


text:
@testset "Distances" begin
  d₁ = aniso2distance([1.,1.], [0.])
  d₂ = aniso2distance([1.,2.], [0.])
  @test evaluate(d₁, [1.,0.], [0.,0.]) == evaluate(d₁, [0.,1.], [0.,0.])
  @test evaluate(d₂, [1.,0.], [0.,0.]) != evaluate(d₂, [0.,1.], [0.,0.])

  d₃ = aniso2distance([1.,.5,.5], [π/4,0.,0.])
  @test evaluate(d₃, [1.,1.,0.], [0.,0.,0.]) ≈ √2
  @test evaluate(d₃, [-1.,1.,0.], [0.,0.,0.]) ≈ √8

  # intrinsic conventions
  gslib = aniso2distance([50.,25.,5.], [30.,-30.,30.], convention=:GSLIB)
  tait = aniso2distance([25.,50.,5.], [-π/6,-π/6,π/6], convention=:TaitBryanIntr)
  euler = aniso2distance([50.,25.,5.], [-deg2rad(78),-deg2rad(41),-deg2rad(50)], convention=:EulerIntr)
  lpf = aniso2distance([50.,25.,5.], [78.,41.,50.], convention=:Leapfrog)
  dm = aniso2distance([50.,25.,5.], [78.,41.,50.], convention=:Datamine)

  @test evaluate(gslib, [1.,0.,0.], [0.,0.,0.]) ≈ 0.1325707358356285
  @test evaluate(gslib, [0.,1.,0.], [0.,0.,0.]) ≈ 0.039051248379533283
  @test evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) ≈ 0.15132745950421558

  @test evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(tait, [0.,0.,1.], [0.,0.,0.])
  @test evaluate(euler, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(dm, [0.,0.,1.], [0.,0.,0.])
  @test evaluate(euler, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(lpf, [0.,0.,1.], [0.,0.,0.])
  @test evaluate(euler, [0.,0.,1.], [0.,0.,0.]) - evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) < 10^-3

  # extrinsic conventions
  xtait = aniso2distance([50.,25.,5.], [π,0,π/2], convention=:TaitBryanExtr)
  xeuler = aniso2distance([50.,25.,5.], [-π/2,-π/2,-π/2], convention=:EulerExtr)

  @test evaluate(xtait, [1.,0.,0.], [0.,0.,0.]) ≈ 0.20
  @test evaluate(xtait, [0.,1.,0.], [0.,0.,0.]) ≈ 0.04
  @test evaluate(xtait, [0.,0.,1.], [0.,0.,0.]) ≈ 0.02

  @test evaluate(xtait, [1.,0.,0.], [0.,0.,0.]) ≈ evaluate(xeuler, [1.,0.,0.], [0.,0.,0.])
  @test evaluate(xtait, [0.,1.,0.], [0.,0.,0.]) ≈ evaluate(xeuler, [0.,1.,0.], [0.,0.,0.])
  @test evaluate(xtait, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(xeuler, [0.,0.,1.], [0.,0.,0.])

end


text:
@testset "200.number-of-islands.jl" begin
    @test lands_num([
        ["1","1","1","1","0"],
        ["1","1","0","1","0"],
        ["1","1","0","0","0"],
        ["0","0","0","0","0"]
      ]) == 1
    @test lands_num([
        ["1","1","0","0","0"],
        ["1","1","0","0","0"],
        ["0","0","1","0","0"],
        ["0","0","0","1","1"]
      ]) == 3
end


text:
using DrWatson
@quickactivate "StatisticsWithJuliaPlutoNotebooks"
DrWatson.greet()


text:
using Stan
using Compat
using Mamba
import Mamba: describe

function describe_str(c::Chains; q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975],
           etype=:bm, args...)
    out = string(
        Mamba.header(c),
        "\nEmpirical Posterior Estimates:\n",
        sprint(io -> showall(io, Mamba.summarystats(c; etype=etype, args...), false)),
        "\nQuantiles:\n",
        sprint(io -> showall(io, Mamba.quantile(c, q=q), false))
    )

    monitor = ["theta", "lp__", "accept_stat__"]

    sim = c[:, monitor, :]
    p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);

    vbox(codemirror(out, readonly=true, linenumbers=false), p...)
end

function run_model(input)
    data = input[:data]
    model_str = input[:code]
    d = parse(data) |> eval
    stanmodel = Stanmodel(name="_model", model=model_str)
    sim = stan(stanmodel, d, "/tmp", CmdStanDir=CMDSTAN_HOME)
    describe_str(sim)
end

function main(window)

    push!(window.assets, "codemirror")
    push!(window.assets, "widgets")
    input = Input(Dict())
    s = Escher.sampler()

    input_pane = (vbox(
        hbox(
            size(50vw, 50vh, codemirror(name=:code) |> watch!(s)),
            hskip(1em),
            size(50vw, 50vh, codemirror(name=:data) |> watch!(s))),
        vskip(1em),
        hbox(button("Run") |> trigger!(s)),
        vskip(1em),
    ) |> Escher.sample(s)) >>> input

    output_pane =
        consume(run_model, input; init="Ready to take inputs.", typ=Any)

    vbox(
        title(2, "StanBox"),
        vskip(1em),
        vbox(input_pane, vskip(1em), output_pane),
    ) |> pad(1em)
end