text:
module Guide

using Compat
using Color
using Compose
using DataStructures
using Gadfly
using Iterators
using JSON

import Gadfly: render, escape_id, default_statistic, jsdata, jsplotdata,
               svg_color_class_from_label

# Where the guide should be placed in relation to the plot.
abstract GuidePosition
immutable TopGuidePosition    <: GuidePosition end
immutable RightGuidePosition  <: GuidePosition end
immutable BottomGuidePosition <: GuidePosition end
immutable LeftGuidePosition   <: GuidePosition end
immutable UnderGuidePosition  <: GuidePosition end
immutable OverGuidePosition   <: GuidePosition end

const top_guide_position    = TopGuidePosition()
const right_guide_position  = RightGuidePosition()
const bottom_guide_position = BottomGuidePosition()
const left_guide_position   = LeftGuidePosition()
const under_guide_position  = UnderGuidePosition()
const over_guide_position   = OverGuidePosition()

# A guide graphic is a position associated with one or more contexts.
# Multiple contexts represent multiple layout possibilites that will be
# optimized over.
immutable PositionedGuide
    ctxs::Vector{Context}
    order::Int
    position::GuidePosition
end

immutable PanelBackground <: Gadfly.GuideElement
end

const background = PanelBackground

function render(guide::PanelBackground, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    back = compose!(context(order=-1),
                    rectangle(),
                    svgclass("guide background"),
                    stroke(theme.panel_stroke),
                    fill(theme.panel_fill),
                    fillopacity(theme.panel_opacity),
                    svgattribute("pointer-events", "visible"))

    return [PositionedGuide([back], 0, under_guide_position)]
end

immutable ZoomSlider <: Gadfly.GuideElement
end

const zoomslider = ZoomSlider

function render(guide::ZoomSlider, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    edge_pad = 3mm
    slide_pad = 0.5mm
    button_size = 4mm
    slider_size = 20mm
    background_color = "#eaeaea"
    foreground_color = "#6a6a6a"
    highlight_color = "#cd5c5c";

    minus_button = compose!(
        context(1w - edge_pad - 2*button_size - slider_size,
                edge_pad, button_size, button_size),
        rectangle(),
        stroke(foreground_color),
        strokeopacity(0.0),
        linewidth(0.3mm),
        {context(),
         polygon([(0.2, 0.4), (0.8, 0.4),
                  (0.8, 0.6), (0.2, 0.6)]),
         fill(foreground_color),
         svgclass("button_logo")},
        fill(background_color),
        jscall(
            """
            click(Gadfly.zoomslider_zoomout_click)
            .mouseenter(Gadfly.zoomslider_button_mouseover)
            .mouseleave(Gadfly.zoomslider_button_mouseout)
            """),
        jsdata("mouseout_color", "\"$(foreground_color)\""),
        jsdata("mouseover_color", "\"$(highlight_color)\""))

    slider_width = 2mm
    slider_xpos = 1w - edge_pad - button_size - slider_size + slide_pad

    slider_min_pos = slider_xpos + slider_width/2
    slider_max_pos = slider_xpos + slider_size - 2*slide_pad - slider_width/2

    slider = compose!(
        context(slider_xpos, edge_pad, slider_size - 2 * slide_pad, button_size),
        {context(),
         rectangle(),
         fill(background_color),
         jscall("click(Gadfly.zoomslider_track_click)"),
         jsdata("min_pos", "%x", Measure[slider_min_pos]),
         jsdata("max_pos", "%x", Measure[slider_max_pos])},
        {context(order=1),
         rectangle(0.5cx - slider_width/2, 0.0, slider_width, 1h),
         fill(foreground_color),
         svgclass("zoomslider_thumb"),
         jscall(
            """
            drag(Gadfly.zoomslider_thumb_dragmove,
                 Gadfly.zoomslider_thumb_dragstart,
                 Gadfly.zoomslider_thumb_dragend)
            .mousedown(Gadfly.zoomslider_thumb_mousedown)
            .mouseup(Gadfly.zoomslider_thumb_mouseup)
            """),
         jsdata("mouseout_color", "\"$(foreground_color)\""),
         jsdata("mouseover_color", "\"$(highlight_color)\""),
         jsdata("min_pos", "%x", Measure[slider_min_pos]),
         jsdata("max_pos", "%x", Measure[slider_max_pos])})

    plus_button = compose!(
        context(1w - edge_pad - button_size, edge_pad,
                button_size, button_size),
        rectangle(),
        stroke(foreground_color),
        strokeopacity(0.0),
        linewidth(0.3mm),
        {context(),
         polygon([(0.2, 0.4), (0.4, 0.4), (0.4, 0.2),
                  (0.6, 0.2), (0.6, 0.4), (0.8, 0.4),
                  (0.8, 0.6), (0.6, 0.6), (0.6, 0.8),
                  (0.4, 0.8), (0.4, 0.6), (0.2, 0.6)]),
         fill(foreground_color),
         svgclass("button_logo")},
        fill(background_color),
        jscall(
            """
            click(Gadfly.zoomslider_zoomin_click)
            .mouseenter(Gadfly.zoomslider_button_mouseover)
            .mouseleave(Gadfly.zoomslider_button_mouseout)
            """),
        jsdata("mouseout_color", "\"$(foreground_color)\""),
        jsdata("mouseover_color", "\"$(highlight_color)\""))

    root = compose!(
        context(withjs=true, units=UnitBox()),
        minus_button,
        slider,
        plus_button,
        stroke(nothing),
        #stroke(foreground_color),
        svgclass("guide zoomslider"),
        fillopacity(0.0))

    return [PositionedGuide([root], 0, over_guide_position)]
end

immutable ColorKey <: Gadfly.GuideElement
    title::Union(String, Nothing)

    function ColorKey(title=nothing)
        new(title)
    end
end

const colorkey = ColorKey

# A helper for render(::ColorKey) for rendering guides for discrete color
# scales.
function render_discrete_color_key(colors::Vector{ColorValue},
                                   labels::OrderedDict{ColorValue, String},
                                   aes_color_label,
                                   title_ctx::Context,
                                   title_width::Measure,
                                   theme::Gadfly.Theme)

    # only consider layouts with a reasonable number of columns
    maxcols = 4

    n = length(colors)

    extents = text_extents(theme.key_label_font,
                           theme.key_label_font_size,
                           values(labels)...)

    ypad = 1.0mm
    title_height = title_ctx.box.height
    entry_height = maximum([height for (width, height) in extents]) + ypad
    swatch_size = entry_height / 2

    # return a context with a lyout of numcols columns
    function make_layout(numcols)
        colrows = Array(Int, numcols)
        m = n
        for i in 1:numcols
            colrows[i] = min(m, ceil(Integer, (n / numcols)))
            m -= colrows[i]
        end

        xpad = 1mm
        colwidths = Array(Measure, numcols)
        m = 0
        for (i, nrows) in enumerate(colrows)
            if m == n
                colwidths[i] = 0mm
            else
                colwidth = maximum([width for (width, height) in extents[m+1:m+nrows]])
                colwidth += swatch_size + 2xpad
                colwidths[i] = colwidth
                m += nrows
            end
        end

        ctxwidth = sum(colwidths)
        ctxheight = entry_height * colrows[1] + title_height

        ctxp = ctxpromise() do draw_context
            yoff = 0.5h - ctxheight/2
            outerctx = context()

            compose!(outerctx, {context(xpad, yoff), title_ctx})

            ctx = context(0, yoff + title_height,
                          ctxwidth, ctxheight - title_height,
                          units=UnitBox(0, 0, 1, colrows[1]))

            m = 0
            xpos = 0w
            for (i, nrows) in enumerate(colrows)
                colwidth = colwidths[i]

                if theme.colorkey_swatch_shape == :square
                    swatches_shapes = rectangle(
                        [xpad], [y*cy - swatch_size/2 for y in 1:nrows],
                        [swatch_size], [swatch_size])
                elseif theme.colorkey_swatch_shape == :circle
                    swatches_shapes = circle([0.5cy], 1:nrows .- 0.5, [swatch_size/2])
                end
                cs = colors[m+1:m+nrows]
                swatches = compose!(
                    context(),
                    swatches_shapes,
                    stroke(nothing),
                    fill(cs))

                swatch_labels = compose!(
                    context(),
                    text([2xpad + swatch_size], [y*cy for y in 1:nrows],
                         collect(values(labels))[m+1:m+nrows], [hleft], [vcenter]),
                    font(theme.key_label_font),
                    fontsize(theme.key_label_font_size),
                    fill(theme.key_label_color))

                col = compose!(context(xpos, yoff), swatches, swatch_labels)
                if aes_color_label != nothing
                    classes = [svg_color_class_from_label(aes_color_label([c])[1]) for c in cs]
                    #class_jscalls = ["data(\"color_class\", \"$(c)\")"
                                     #for c in classes]
                    compose!(col,
                        svgclass(classes),
                        jscall(["""
                            data(\"color_class\", \"$(c)\")
                            .click(Gadfly.colorkey_swatch_click)
                            """ for c in classes]))
                end
                compose!(ctx, col)

                m += nrows
                xpos += colwidths[i]
            end

            return compose!(outerctx, ctx,
                            # defeat webkit's asinine default drag behavior
                            jscall("drag(function() {}, function() {}, function() {})"),
                            svgclass("guide colorkey"))
        end

        return compose!(
            context(minwidth=max(title_width, ctxwidth),
                    minheight=ctxheight,
                    units=UnitBox()),
            ctxp)
    end

    return map(make_layout, 1:maxcols)
end

# A helper for render(::ColorKey) for rendering guides for continuous color
# scales.
function render_continuous_color_key(colors::Dict,
                                     labels::OrderedDict{ColorValue, String},
                                     color_function::Function,
                                     title_context::Context,
                                     title_width::Measure,
                                     theme::Gadfly.Theme)

    entry_width, entry_height = max_text_extents(theme.key_label_font,
                                                 theme.key_label_font_size,
                                                 values(labels)...)

    numlabels = length(labels)
    title_height = title_context.box.height
    total_height = 1.5 * numlabels * entry_height + title_height
    swatch_width = entry_height / 2
    xoff = 2mm
    padding = 1mm
    entry_width += 2padding + swatch_width + xoff

    ctx = context(minwidth=max(title_width, entry_width),
                  minheight=total_height, units=UnitBox())

    yoff = 0.5h - total_height/2

    compose!(ctx, {context(xoff, yoff), title_context})

    # color bar
    compose!(ctx,
        {context(xoff, yoff + title_height,
                 1w, total_height, units=UnitBox()),
         rectangle(
             [0],
             [1*cy - i*total_height / theme.key_color_gradations
              for i in 1:theme.key_color_gradations],
             [swatch_width],
             [total_height / theme.key_color_gradations]),

         #grid lines
         {context(),
          line([{(0, 1 - y), (swatch_width, 1 - y)} for y in values(colors)]),
          linewidth(theme.grid_line_width),
          stroke(color("white"))},

         fill([color_function((i-1) / (theme.key_color_gradations - 1))
               for i in 1:theme.key_color_gradations]),
         stroke(nothing),
         svgattribute("shape-rendering", "crispEdges")})

    compose!(ctx,
        {context(xoff + swatch_width + padding, yoff + title_height,
                 1w, total_height, units=UnitBox()),
         text([0],
              [1 - y for y in values(colors)],
              [labels[c] for c in keys(colors)],
              [hleft], [vcenter]),
         fill(theme.key_label_color),
         font(theme.key_label_font),
         fontsize(theme.key_label_font_size)},
         svgclass("guide colorkey"))

    return [ctx]
end

function render_colorkey_title(title::String, theme::Gadfly.Theme)
    title_width, title_height = max_text_extents(theme.key_title_font,
                                                 theme.key_title_font_size,
                                                 title)

    if theme.guide_title_position == :left
        title_form = text(0.0w, title_height, title, hleft, vbottom)
    elseif theme.guide_title_position == :center
        title_form = text(0.5w, title_height, title, hcenter, vbottom)
    elseif theme.guide_title_position == :right
        title_form = text(1.0w, title_height, title, hright, vbottom)
    else
        error("$(theme.guide_title_position) is not a valid guide title position")
    end

    title_padding = 4mm
    title_context = compose!(
        context(0w, 0h, 1w, title_height + title_padding),
        title_form,
        stroke(nothing),
        font(theme.key_title_font),
        fontsize(theme.key_title_font_size),
        fill(theme.key_title_color))

    return title_context, title_width
end

function render(guide::ColorKey, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    if theme.key_position == :none
        return PositionedGuide[]
    end

    used_colors = Set{ColorValue}()
    colors = Array(ColorValue, 0) # to preserve ordering
    labels = OrderedDict{ColorValue, Set{String}}()

    continuous_guide = false
    guide_title = guide.title

    if guide_title === nothing && !is(aes.color_key_title, nothing)
        guide_title = aes.color_key_title
    end

    if aes.color_key_colors != nothing &&
       aes.color_key_continuous != nothing &&
       aes.color_key_continuous
        continuous_guide = true
    end

    color_key_labels = aes.color_label(keys(aes.color_key_colors))
    for (color, label) in zip(keys(aes.color_key_colors), color_key_labels)
        if !in(color, used_colors)
            push!(used_colors, color)
            push!(colors, color)
            labels[color] = Set{String}()
            push!(labels[color], label)
        else
            push!(labels[color], label)
        end
    end

    if guide_title === nothing
        guide_title = "Color"
    end

    pretty_labels = OrderedDict{ColorValue, String}()
    for (color, label) in labels
        pretty_labels[color] = join(labels[color], ", ")
    end

    title_context, title_width = render_colorkey_title(guide_title, theme)

    if theme.colorkey_swatch_shape != :circle &&
    theme.colorkey_swatch_shape != :square
        error("$(theme.colorkey_swatch_shape) is not a valid color key swatch shape")
    end

    if continuous_guide
        ctxs = render_continuous_color_key(aes.color_key_colors,
                                           pretty_labels,
                                           aes.color_function,
                                           title_context,
                                           title_width, theme)
    else
        ctxs = render_discrete_color_key(colors, pretty_labels,
                                         aes.color_label,
                                         title_context,
                                         title_width, theme)
    end

    position = right_guide_position
    if theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :right
        position = right_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [PositionedGuide(ctxs, 0, position)]
end

immutable ManualColorKey <: Gadfly.GuideElement
    title::Union(String, Nothing)
    labels::Vector{String}
    colors::Vector
end

const manual_color_key = ManualColorKey

function render(guide::ManualColorKey, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    if theme.key_position == :none
        return PositionedGuide[]
    end

    guide_title = guide.title

    if guide_title === nothing && !is(aes.color_key_title, nothing)
        guide_title = aes.color_key_title
    end

    if guide_title === nothing
        guide_title = "Color"
    end

    title_context, title_width = render_colorkey_title(guide_title, theme)

    colors = ColorValue[color(c) for c in guide.colors]
    labels = (ColorValue => String)[c => l for (c, l) in zip(colors, guide.labels)]

    ctxs = render_discrete_color_key(colors, labels, title_context, title_width, theme)

    position = right_guide_position
    if theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :right
        position = right_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [PositionedGuide(ctxs, 0, position)]
end

immutable XTicks <: Gadfly.GuideElement
    label::Bool
    ticks::Union(Nothing, AbstractArray)
    orientation::Symbol

    function XTicks(; label::Bool=true,
                      ticks::Union(Nothing, AbstractArray)=nothing,
                      orientation::Symbol=:auto)
        return new(label, ticks, orientation)
    end
end

const xticks = XTicks

function default_statistic(guide::XTicks)
    return Stat.xticks(guide.ticks)
end

function render(guide::XTicks, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    if Gadfly.issomething(aes.xtick)
        ticks = aes.xtick
        tickvisibility = aes.xtickvisible
        scale = aes.xtickscale

        T = eltype(aes.xtick)
        labels = String[]
        for scale_ticks in groupby(zip(scale, ticks), x -> x[1])
            append!(labels, aes.xtick_label(T[t for (s, t) in scale_ticks]))
        end
    else
        labels = String[]
        ticks = {}
        tickvisibility = {}
        scale = {}
    end

    if Gadfly.issomething(aes.xgrid)
        grids = aes.xgrid
        if length(grids) < length(ticks)
            gridvisibility = tickvisibility[2:end]
        else
            gridvisibility = tickvisibility
        end
    else
        grids = {}
        gridvisibility = {}
    end

    # grid lines
    static_grid_lines = compose!(
        context(withoutjs=true),
        line([{(t, 0h), (t, 1h)} for t in grids[gridvisibility]]),
        stroke(theme.grid_color),
        linewidth(theme.grid_line_width),
        strokedash(theme.grid_strokedash),
        svgclass("guide xgridlines yfixed"))

    dynamic_grid_lines = compose!(
        context(withjs=true),
        line([{(t, 0h), (t, 1h)} for t in grids]),
        visible(gridvisibility),
        stroke(theme.grid_color),
        linewidth(theme.grid_line_width),
        strokedash(theme.grid_strokedash),
        svgclass("guide xgridlines yfixed"),
        svgattribute("gadfly:scale", scale),
        jsplotdata("focused_xgrid_color",
                   "\"#$(hex(theme.grid_color_focused))\""),
        jsplotdata("unfocused_xgrid_color",
                   "\"#$(hex(theme.grid_color))\""))

    grid_lines = compose!(context(), static_grid_lines, dynamic_grid_lines)

    if !guide.label
        return [PositionedGuide([grid_lines], 0, under_guide_position)]
    end

    label_sizes = text_extents(theme.minor_label_font,
                               theme.minor_label_font_size,
                               labels...)
    label_widths = [width for (width, height) in label_sizes]
    label_heights = [height for (width, height) in label_sizes]

    padding = 1mm

    hlayout = ctxpromise() do draw_context
        static_labels = compose!(
            context(withoutjs=true),
            text(ticks[tickvisibility], [1h + padding], labels[tickvisibility],
                 [hcenter], [vbottom]),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgclass("guide xlabels"))

        dynamic_labels = compose!(
            context(withjs=true),
            text(ticks, [1h + padding], labels, [hcenter], [vbottom]),
            visible(tickvisibility),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgattribute("gadfly:scale", scale),
            svgclass("guide xlabels"))

        return compose!(context(), static_labels, dynamic_labels)
    end
    hlayout_context = compose!(context(minwidth=sum(label_widths[tickvisibility]),
                                       minheight=maximum(label_heights[tickvisibility])),
                               hlayout)

    vlayout = ctxpromise() do draw_context
        static_labels = compose!(
            context(withoutjs=true),
            text(ticks[tickvisibility],
                 [padding],
                 labels[tickvisibility],
                 [hright], [vcenter],
                 [Rotation(-0.5pi, tick, padding) for tick in ticks[tickvisibility]]),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgclass("guide xlabels"))

        dynamic_labels = compose!(
            context(withjs=true),
            text(ticks, [padding], labels, [hright], [vbottom],
                 [Rotation(-0.5pi, tick, padding) for tick in ticks]),
            visible(tickvisibility),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgattribute("gadfly:scale", scale),
            svgclass("guide xlabels"))

        return compose!(context(), static_labels, dynamic_labels)

    end
    vpenalty = 3
    vlayout_context = compose!(context(minwidth=sum(label_heights[tickvisibility]),
                                       minheight=maximum(label_widths[tickvisibility]),
                                       penalty=vpenalty), vlayout)

    if guide.orientation == :horizontal
        contexts = [hlayout_context]
    elseif guide.orientation == :vertical
        contexts = [vlayout_context]
    elseif guide.orientation == :auto
        contexts = [hlayout_context, vlayout_context]
    else
        error("$(guide.layout) is not a valid orientation for Guide.yticks")
    end

    return [PositionedGuide(contexts, 10,
                            bottom_guide_position),
            PositionedGuide([grid_lines], 0, under_guide_position)]
end

immutable YTicks <: Gadfly.GuideElement
    label::Bool
    ticks::Union(Nothing, AbstractArray)
    orientation::Symbol

    function YTicks(; label::Bool=true,
                      ticks::Union(Nothing, AbstractArray)=nothing,
                      orientation::Symbol=:horizontal)
        new(label, ticks, orientation)
    end
end

const yticks = YTicks

function default_statistic(guide::YTicks)
    Stat.yticks(guide.ticks)
end

function render(guide::YTicks, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    if Gadfly.issomething(aes.ytick)
        ticks = aes.ytick
        tickvisibility = aes.ytickvisible
        scale = aes.ytickscale
        T = eltype(aes.ytick)
        labels = String[]
        for scale_ticks in groupby(zip(scale, ticks), x -> x[1])
            append!(labels, aes.ytick_label(T[t for (s, t) in scale_ticks]))
        end
    else
        labels = String[]
        ticks = {}
        tickvisibility = {}
        scale = {}
    end

    if Gadfly.issomething(aes.ygrid)
        grids = aes.ygrid
        if length(grids) < length(ticks)
            gridvisibility = tickvisibility[2:end]
        else
            gridvisibility = tickvisibility
        end
    else
        grids = {}
        gridvisibility = {}
    end

    # grid lines
    static_grid_lines = compose!(
        context(withoutjs=true),
        line([{(0w, t), (1w, t)} for t in grids[gridvisibility]]),
        stroke(theme.grid_color),
        linewidth(theme.grid_line_width),
        strokedash(theme.grid_strokedash),
        svgclass("guide ygridlines xfixed"))

    dynamic_grid_lines = compose!(
        context(withjs=true),
        line([{(0w, t), (1w, t)} for t in grids]),
        visible(gridvisibility),
        stroke(theme.grid_color),
        linewidth(theme.grid_line_width),
        strokedash(theme.grid_strokedash),
        svgclass("guide ygridlines xfixed"),
        svgattribute("gadfly:scale", scale),
        jsplotdata("focused_ygrid_color",
                   "\"#$(hex(theme.grid_color_focused))\""),
        jsplotdata("unfocused_ygrid_color",
                   "\"#$(hex(theme.grid_color))\""))

    grid_lines = compose!(context(), static_grid_lines, dynamic_grid_lines)

    if !guide.label
        return [PositionedGuide([grid_lines], 0, under_guide_position)]
    end

    label_sizes = text_extents(theme.minor_label_font,
                               theme.minor_label_font_size,
                               labels...)
    label_widths = [width for (width, height) in label_sizes]
    label_heights = [height for (width, height) in label_sizes]
    padding = 1mm

    hlayout = ctxpromise() do draw_context
        static_labels = compose!(
            context(withoutjs=true),
            text([1.0w - padding], ticks[tickvisibility], labels[tickvisibility],
                 [hright], [vcenter]),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgclass("guide ylabels"))

        dynamic_labels = compose!(
            context(withjs=true),
            text([1.0w - padding], ticks, labels,
                 [hright], [vcenter]),
            visible(tickvisibility),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgattribute("gadfly:scale", scale),
            svgclass("guide ylabels"))

        return compose!(context(), static_labels, dynamic_labels)
    end
    hpenalty = 3
    hlayout_context =
        compose!(context(minwidth=maximum(label_widths[tickvisibility]) + 2padding,
                         minheight=sum(label_heights[tickvisibility]),
                         penalty=hpenalty), hlayout)

    vlayout = ctxpromise() do draw_context
        static_grid_lines = compose!(
            context(),
            text([1.0w - padding], ticks[tickvisibility], labels[tickvisibility],
                 [hcenter], [vbottom],
                 [Rotation(-0.5pi, (1.0w - padding, tick))
                  for tick in ticks[tickvisibility]]),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgclass("guide ylabels"))

        dynamic_grid_lines = compose!(
            context(),
            text([1.0w - padding], ticks, labels,
                 [hcenter], [vbottom],
                 [Rotation(-0.5pi, (1.0w - padding, tick))
                  for tick in ticks[tickvisibility]]),
            visible(tickvisibility),
            fill(theme.minor_label_color),
            font(theme.minor_label_font),
            fontsize(theme.minor_label_font_size),
            svgattribute("gadfly:scale", scale),
            svgclass("guide ylabels"))

        return compose!(context(), static_grid_lines, dynamic_grid_lines)
    end
    vlayout_context =
        compose!(context(minwidth=maximum(label_heights[tickvisibility]),
                         minheight=sum(label_widths[tickvisibility])),
                 vlayout)

    if guide.orientation == :horizontal
        contexts = [hlayout_context]
    elseif guide.orientation == :vertical
        contexts = [vlayout_context]
    elseif guide.orientation == :auto
        contexts = [hlayout_context, vlayout_context]
    else
        error("$(guide.layout) is not a valid orientation for Guide.yticks")
    end

    return [PositionedGuide(contexts, 10,
                            left_guide_position),
            PositionedGuide([grid_lines], 0, under_guide_position)]
end

# X-axis label Guide
immutable XLabel <: Gadfly.GuideElement
    label::Union(Nothing, String)
    orientation::Symbol

    function XLabel(label; orientation::Symbol=:auto)
        return new(label, orientation)
    end
end

const xlabel = XLabel

function render(guide::XLabel, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    if guide.label === nothing || isempty(guide.label)
        return nothing
    end

    text_width, text_height = max_text_extents(theme.major_label_font,
                                               theme.major_label_font_size,
                                               guide.label)

    padding = 3mm
    hlayout = ctxpromise() do draw_context
        return compose!(context(),
                        text(0.5w, 1h - padding, guide.label, hcenter, vbottom),
                        stroke(nothing),
                        fill(theme.major_label_color),
                        font(theme.major_label_font),
                        fontsize(theme.major_label_font_size))
    end
    hlayout_context = compose!(context(minwidth=text_width + 2padding,
                                       minheight=text_height + 2padding),
                               hlayout)

    vlayout = ctxpromise() do draw_context
        return compose!(context(),
                        text(0.5w, padding, guide.label, hright, vcenter,
                             Rotation(-0.5pi, 0.5w, padding)),
                        stroke(nothing),
                        fill(theme.major_label_color),
                        font(theme.major_label_font),
                        fontsize(theme.major_label_font_size))
    end
    vpenalty = 3
    vlayout_context = compose!(context(minwidth=text_height + 2padding,
                                       minheight=text_width + 2padding,
                                       penalty=3),
                               vlayout)

    if guide.orientation == :horizontal
        contexts = [hlayout_context]
    elseif guide.orientation == :vertical
        contexts = [vlayout_context]
    elseif guide.orientation == :auto
        contexts = [hlayout_context, vlayout_context]
    else
        error("$(guide.layout) is not a valid orientation for Guide.xlabel")
    end

    return [PositionedGuide(contexts, 0, bottom_guide_position)]
end

# Y-axis label Guide
immutable YLabel <: Gadfly.GuideElement
    label::Union(Nothing, String)
    orientation::Symbol

    function YLabel(label; orientation::Symbol=:auto)
        return new(label, orientation)
    end
end

const ylabel = YLabel

function render(guide::YLabel, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    if guide.label === nothing || isempty(guide.label)
        return nothing
    end

    text_width, text_height = max_text_extents(theme.major_label_font,
                                               theme.major_label_font_size,
                                               guide.label)

    padding = 2mm
    hlayout = ctxpromise() do draw_context
        return compose!(context(),
                        text(1.0w - padding, 0.5h, guide.label, hright, vcenter),
                        stroke(nothing),
                        fill(theme.major_label_color),
                        font(theme.major_label_font),
                        fontsize(theme.major_label_font_size))
    end
    hlayout_context = compose!(context(minwidth=text_width + 2padding,
                                       minheight=text_height + 2padding), hlayout)

    vlayout = ctxpromise() do draw_context
        return compose!(context(),
                        text(0.5w, 0.5h - padding, guide.label, hcenter, vcenter, Rotation(-0.5pi)),
                        stroke(nothing),
                        fill(theme.major_label_color),
                        font(theme.major_label_font),
                        fontsize(theme.major_label_font_size))
    end
    vlayout_context = compose!(context(minwidth=text_height + 2padding,
                                       minheight=text_width + 2padding), vlayout)

    if guide.orientation == :horizontal
        contexts = [hlayout_context]
    elseif guide.orientation == :vertical
        contexts = [vlayout_context]
    elseif guide.orientation == :auto
        contexts = [hlayout_context, vlayout_context]
    else
        error("$(guide.layout) is not a valid orientation for Guide.ylabel")
    end

    return [PositionedGuide(contexts, 0, left_guide_position)]
end

# Title Guide
immutable Title <: Gadfly.GuideElement
    label::Union(Nothing, String)
end

const title = Title

function render(guide::Title, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    if guide.label === nothing || isempty(guide.label)
        return nothing
    end

    (text_width, text_height) = max_text_extents(theme.major_label_font,
                                                 theme.major_label_font_size,
                                                 guide.label)

    padding = 2mm
    ctx = compose!(
        context(minwidth=text_width, minheight=text_height + 2padding),
        text(0.5w, 1h - padding, guide.label, hcenter, vbottom),
        stroke(nothing),
        fill(theme.major_label_color),
        font(theme.major_label_font),
        fontsize(theme.major_label_font_size))

    return [PositionedGuide([ctx], 0, top_guide_position)]
end

# Arrange a plot with its guides
#
# Args:
#   plot_canvas: A canvas containing the plot graphics.
#   guides: Tuples of guide canvases each with a GuidePosition giving
#           where the guide should be placed relative to the plot.
#
# Returns:
#   A new canvas containing the plot with guides layed out in the specified
#   manner.
function layout_guides(plot_context::Context,
                       coord::Gadfly.CoordinateElement,
                       theme::Gadfly.Theme,
                       positioned_guides::PositionedGuide...)
    # Organize guides by position
    guides = DefaultDict(() -> (Vector{Context}, Int)[])
    for positioned_guide in positioned_guides
        push!(guides[positioned_guide.position],
              (positioned_guide.ctxs, positioned_guide.order))
    end

    for (position, ordered_guides) in guides
        if position == left_guide_position || position == top_guide_position
            sort!(ordered_guides, by=x -> x[2])
        else
            sort!(ordered_guides, by=x -> -x[2])
        end
    end

    m = 1 + length(guides[top_guide_position]) +
            length(guides[bottom_guide_position])
    n = 1 + length(guides[left_guide_position]) +
            length(guides[right_guide_position])

    focus_y = 1 + length(guides[top_guide_position])
    focus_x = 1 + length(guides[left_guide_position])

    plot_units = plot_context.units

    # Populate the table

    aspect_ratio = nothing
    if isa(coord, Gadfly.Coord.cartesian)
        if coord.fixed
            aspect_ratio = abs(plot_context.units.width / plot_context.units.height)
        elseif coord.aspect_ratio != nothing
            aspect_ratio = coord.aspect_ratio
        end
    end
    tbl = table(m, n, focus_y:focus_y, focus_x:focus_x, units=plot_units,
                aspect_ratio=aspect_ratio)

    i = 1
    for (ctxs, order) in guides[top_guide_position]
        for ctx in ctxs
            if ctx.units == Compose.nil_unit_box
                ctx.units = UnitBox(plot_units, toppad=0mm, bottompad=0mm)
            end
        end

        tbl[i, focus_x] = ctxs
        i += 1
    end
    i += 1
    for (ctxs, order) in guides[bottom_guide_position]
        for ctx in ctxs
            if ctx.units == Compose.nil_unit_box
                ctx.units = UnitBox(plot_units, toppad=0mm, bottompad=0mm)
            end
        end

        tbl[i, focus_x] = ctxs
        i += 1
    end

    j = 1
    for (ctxs, order) in guides[left_guide_position]
        for ctx in ctxs
            if ctx.units == Compose.nil_unit_box
                ctx.units = UnitBox(plot_units, leftpad=0mm, rightpad=0mm)
            end
        end

        tbl[focus_y, j] = ctxs
        j += 1
    end
    j += 1
    for (ctxs, order) in guides[right_guide_position]
        for ctx in ctxs
            if ctx.units == Compose.nil_unit_box
                ctx.units = UnitBox(plot_units, leftpad=0mm, rightpad=0mm)
            end
        end

        tbl[focus_y, j] = ctxs
        j += 1
    end

    tbl[focus_y, focus_x] =
        [compose!(context(minwidth=minwidth(plot_context),
                          minheight=minheight(plot_context),
                          units=plot_units,
                          clip=true),
                  {context(order=-1),
                     [c for (c, o) in guides[under_guide_position]]...},
                  {context(order=1000),
                     [c for (c, o) in guides[over_guide_position]]...},
                  {context(order=0),
                     plot_context},
                  jscall(
                    """
                    mouseenter(Gadfly.plot_mouseover)
                    .mouseleave(Gadfly.plot_mouseout)
                    .mousewheel(Gadfly.guide_background_scroll)
                    .drag(Gadfly.guide_background_drag_onmove,
                          Gadfly.guide_background_drag_onstart,
                          Gadfly.guide_background_drag_onend)
                    """))]

    return tbl
end

end # module Guide


text:
using ProgressMeter
using Distances

abstract type ParticleFilter end

function Base.filter!(ps::ParticleFilter, ys::AbstractVector{T} where T, progress=true;
                      cb=(ps,y)->nothing)
    @showprogress (progress ? 1 : Inf) "Fitting particles..." for y in ys
        fit!(ps, y)
        cb(ps, y)
    end
    ps
end

particles(p::ParticleFilter) = p.particles

Statistics.mean(f::Function, p::ParticleFilter) =
    mean(f.(particles(p)), Weights(weight.(particles(p))))

state_entropy(p::ParticleFilter) = mean(state_entropy, p)

normalize!(x) = (x ./= sum(x))
posterior_predictive(p::ParticleFilter) =
    MixtureModel(posterior_predictive.(particles(p)), normalize!(weight.(particles(p))))

function assignments(p::ParticleFilter)
    ps = particles(p)
    a = assignments(first(ps))
    asgn = Matrix{eltype(a)}(undef, length(a), length(ps))
    for (i,p) in enumerate(ps)
        asgn[:,i] .= assignments(p)
    end
    asgn
end

ncomponents_dist(p::ParticleFilter) =
    fit(Categorical, ncomponents.(p.particles), weight.(p.particles))

function assignment_similarity(ps::ParticleFilter)
    as = assignments(ps)
    1 .- pairwise(Hamming(), as, dims=1) ./ size(as, 2)
end

"""
    randindex(ps::ParticleFilter, truth::Vector{<:Integer})

Compute the Rand Index of the clusters found by the particles in `ps`, given
ground truth cluster indices.  The Rand Index is the number of pairs of
observations on which the two clusterings agree, divided by the total number of
pairs.

"""
function randindex(ps::ParticleFilter, truth::T, type=:adjusted) where {T<:AbstractVector{<:Integer}}
    ps_sim = assignment_similarity(ps)
    truth_sim = truth .== truth'

    size(truth_sim) == size(ps_sim) ||
        throw(DimensionMismatch("Ground truth and particles have different number of obs!"))

    # number of pairs
    N = prod(size(ps_sim))

    # adjustment needs the sum of squares of the clusters sizes, which is just
    # the area of each assignment similarity matrix that's 1
    nis = sum(ps_sim)
    njs = sum(truth_sim)

    both_same = sum(ps_sim .* truth_sim)
    # sum((1 .- ps_sim) .* (1 .- truth_sim)) = sum(1 - truth_sim - ps_sim + truth_sim*ps_sim)
    both_diff = N - nis - njs + both_same

    if type == :adjusted
        n = length(truth)
        # from Clustering.jl:
        nc = (n*(n^2+1)-(n+1)*nis-(n+1)*njs+2*(nis*njs)/n)/(2*(n-1))
        adj_rand_index = (both_same + both_diff - nc) / (N - nc)
    else
        rand_index = (both_same + both_diff) / N
    end
end

ext: jl


text:
using LightGraphs
using MatrixDepot
using Base.Profile

"""Find the largest connected component of graph and return it as a vector of indices."""
function symmetrize(A)
    println("Symmetrizing ")
    tic()
    if !issymmetric(A)
        println(STDERR, "the matrix is not symmetric using A+A'")
        A = A + A'
    end
    if isa(A, Base.LinAlg.Symmetric)
        A.data.nzval = abs(A.data.nzval)
    else
        A.nzval = abs(A.nzval)
    end
    #= spA = abs(sparse(A)) =#
    toc()
    return A
end

function loadmat(matname)
    println("Reading MTX of $matname")
    tic()
    A = matrixdepot(matname, :read)
    A = symmetrize(A)
    tic()
    g = Graph(A)
    return g
end

names = ["Newman/football" , "Newman/cond-mat-2003", "SNAP/amazon0302","SNAP/roadNet-CA"]
sucesses = []
failures = []
for matname in names
    println("Working on $matname")
    g = loadmat(matname)
    println("Finding Components")
    visitor = LightGraphs.TreeBFSVisitorVector(zeros(Int, nv(g)))
    label = zeros(Int, nv(g))
    @time label = LightGraphs.connected_components!(label, g)
    @time components = LightGraphs.connected_components!(visitor, g)
    fill!(visitor.tree, 0)
    @time LightGraphs.bfs_tree!(visitor, g, 1)
    @show length(components)
    #=
    @time components_slow = LightGraphs.connected_components(g)
    @assert length(components) == length(components_slow)
    =#
end


text:
@testset "Multinetwork dimensions" begin

    # Some tests here intentionally throw errors, so we temporarily raise the FlexPlan
    # logger level to prevent them from being displayed.
    previous_FlexPlan_logger_level = Memento.getlevel(Memento.getlogger(FlexPlan))
    Memento.setlevel!(Memento.getlogger(FlexPlan), "alert")

    #=
    julia> DataFrames.DataFrame(
               nw = 1:24,
               hour = repeat(1:4; outer=6),
               scenario = repeat(1:3; inner=4, outer=2),
               sub_nw = repeat(1:2; inner=12)
           )
    24×4 DataFrame
     Row │ nw     hour   scenario  sub_nw
         │ Int64  Int64  Int64     Int64
    ─────┼────────────────────────────────
       1 │     1      1         1       1
       2 │     2      2         1       1
       3 │     3      3         1       1
       4 │     4      4         1       1
       5 │     5      1         2       1
       6 │     6      2         2       1
       7 │     7      3         2       1
       8 │     8      4         2       1
       9 │     9      1         3       1
      10 │    10      2         3       1
      11 │    11      3         3       1
      12 │    12      4         3       1
      13 │    13      1         1       2
      14 │    14      2         1       2
      15 │    15      3         1       2
      16 │    16      4         1       2
      17 │    17      1         2       2
      18 │    18      2         2       2
      19 │    19      3         2       2
      20 │    20      4         2       2
      21 │    21      1         3       2
      22 │    22      2         3       2
      23 │    23      3         3       2
      24 │    24      4         3       2
    =#

    sn_data = Dict{String,Any}(c=>Dict{String,Any}() for c in ("bus","branch","dcline","gen","load","shunt","switch","storage")) # Fake a single-network data structure
    _FP.add_dimension!(sn_data, :hour, 4)
    _FP.add_dimension!(sn_data, :scenario, Dict(s => Dict{String,Any}("probability"=>s/6) for s in 1:3))
    _FP.add_dimension!(sn_data, :sub_nw, 2; metadata = Dict{String,Any}("description"=>"sub_nws model different physical networks"))
    dt = Dict{String,Any}("dim"=>_FP.dim(sn_data), "multinetwork"=>true, "nw"=>Dict{String,Any}("1"=>sn_data)) # Fake a multinetwork data structure
    pm = _PM.instantiate_model(dt, _PM.ACPPowerModel, pm->nothing)
    dim = _FP.dim(pm)

    dim_shift = deepcopy(dim)
    _FP.shift_ids!(dim_shift, 24)
    dt_shift = Dict{String,Any}("dim"=>dim_shift, "multinetwork"=>true, "nw"=>Dict{String,Any}("1"=>sn_data))
    pm_shift = _PM.instantiate_model(dt_shift, _PM.ACPPowerModel, pm->nothing)

    @testset "add_dimension!" begin
        @test_throws ErrorException _FP.add_dimension!(sn_data, :hour, 4) # Trying to add a dimension having the same name of an existing one
        @test_throws ErrorException _FP.add_dimension!(sn_data, :newdim, Dict(s => Dict{String,Any}("prop"=>"val") for s in [1,2,4])) # Trying to add a dimension having a property Dict whose keys are not consecutive Ints starting at 1
    end

    @testset "shift_ids!" begin
        @test _FP.nw_ids(dim_shift) == collect(25:48)
        @test _FP.shift_ids!(deepcopy(sn_data), 24) == collect(25:48)
        @test_throws ErrorException _FP.shift_ids!(dt, 1) # Trying to shift ids of a multinetwork
    end

    @testset "merge_dim!" begin
        dt1 = deepcopy(dt)
        dt2 = deepcopy(dt)
        delete!(dt2, "dim")
        _FP.add_dimension!(dt2, :hour, 4)
        _FP.add_dimension!(dt2, :sub_nw, 2; metadata = Dict{String,Any}("description"=>"sub_nws model different physical networks"))
        _FP.add_dimension!(dt2, :scenario, Dict(s => Dict{String,Any}("probability"=>s/6) for s in 1:3))
        @test_throws ErrorException _FP.merge_dim!(dt1["dim"], dt2["dim"], :sub_nw) # Dimensions are not sorted in the same way
        dt1 = deepcopy(dt)
        dt2 = deepcopy(dt)
        _FP.dim_prop(dt2, :scenario, 1)["probability"] = 1/2
        @test_throws ErrorException _FP.merge_dim!(dt1["dim"], dt2["dim"], :sub_nw) # Different property along a dimension that is not being merged
        dt1 = deepcopy(dt)
        dt2 = deepcopy(dt)
        _FP.dim_meta(dt2, :sub_nw)["description"] = ""
        @test_throws ErrorException _FP.merge_dim!(dt1["dim"], dt2["dim"], :sub_nw) # Different metadata
        dt1 = deepcopy(dt)
        sn_data_shift = deepcopy(sn_data)
        _FP.shift_ids!(sn_data_shift, 23)
        dt2 = Dict{String,Any}("dim"=>sn_data_shift["dim"], "multinetwork"=>true, "nw"=>Dict{String,Any}("1"=>sn_data))
        @test_throws ErrorException _FP.merge_dim!(dt1["dim"], dt2["dim"], :sub_nw) # Ids are not contiguous
        dt1 = deepcopy(dt)
        sn_data_shift = deepcopy(sn_data)
        _FP.shift_ids!(sn_data_shift, 25)
        dt2 = Dict{String,Any}("dim"=>sn_data_shift["dim"], "multinetwork"=>true, "nw"=>Dict{String,Any}("1"=>sn_data))
        @test_throws ErrorException _FP.merge_dim!(dt1["dim"], dt2["dim"], :sub_nw) # Ids are not contiguous
        dt1 = deepcopy(dt)
        dt2 = deepcopy(dt)
        delete!(dt2, "dim")
        _FP.add_dimension!(dt2, :hour, 4)
        _FP.add_dimension!(dt2, :scenario, Dict(s => Dict{String,Any}("probability"=>s/6) for s in 1:3))
        _FP.add_dimension!(dt2, :sub_nw, 4; metadata = Dict{String,Any}("description"=>"sub_nws model different physical networks"))
        @test _FP.merge_dim!(dt1["dim"], dt_shift["dim"], :sub_nw) == dt2["dim"]
    end

    @testset "slice_dim" begin
        slice, ids = _FP.slice_dim(dim, hour=2)
        @test _FP.dim_length(slice) == 6
        @test _FP.dim_length(slice, :hour) == 1
        @test _FP.dim_length(slice, :scenario) == 3
        @test _FP.dim_length(slice, :sub_nw) == 2
        @test _FP.dim_meta(slice, :hour, "orig_id") == 2
        @test ids == [2,6,10,14,18,22]
        slice, ids = _FP.slice_dim(dim, hour=2, scenario=3)
        @test _FP.dim_length(slice) == 2
        @test _FP.dim_length(slice, :hour) == 1
        @test _FP.dim_length(slice, :scenario) == 1
        @test _FP.dim_length(slice, :sub_nw) == 2
        @test _FP.dim_meta(slice, :hour, "orig_id") == 2
        @test _FP.dim_meta(slice, :scenario, "orig_id") == 3
        @test ids == [10,22]
    end

    @testset "require_dim" begin
        @test_throws ErrorException _FP.require_dim(Dict{String,Any}()) # Missing `dim` dict
        @test_throws ErrorException _FP.require_dim(dt, :newdim) # Missing `newdim` dimension
    end

    @testset "nw_ids" begin
        @test _FP.nw_ids(dim)                         == collect(1:24)
        @test _FP.nw_ids(dim, hour=4)                 == [4,8,12,16,20,24]
        @test _FP.nw_ids(dim, scenario=2)             == [5,6,7,8,17,18,19,20]
        @test _FP.nw_ids(dim, sub_nw=1)               == [1,2,3,4,5,6,7,8,9,10,11,12]
        @test _FP.nw_ids(dim, hour=4, scenario=2)     == [8,20]
        @test _FP.nw_ids(dim, hour=2:4)               == [2,3,4,6,7,8,10,11,12,14,15,16,18,19,20,22,23,24]
        @test _FP.nw_ids(dim, hour=2:4, scenario=2)   == [6,7,8,18,19,20]
        @test _FP.nw_ids(dim, hour=[2,4])             == [2,4,6,8,10,12,14,16,18,20,22,24]
        @test _FP.nw_ids(dim, hour=[2,4], scenario=2) == [6,8,18,20]
        @test _FP.nw_ids(dim_shift)                   == collect(25:48)
        @test _FP.nw_ids(dt)                          == _FP.nw_ids(dim)
        @test _FP.nw_ids(pm)                          == _FP.nw_ids(dim)
    end

    @testset "similar_ids" begin
        @test _FP.similar_ids(dim, 7)                           == [7]
        @test _FP.similar_ids(dim, 7, hour=4)                   == [8]
        @test _FP.similar_ids(dim, 7, scenario=1)               == [3]
        @test _FP.similar_ids(dim, 7, hour=4, scenario=1)       == [4]
        @test _FP.similar_ids(dim, 7, hour=2:4)                 == [6,7,8]
        @test _FP.similar_ids(dim, 7, hour=[2,4])               == [6,8]
        @test _FP.similar_ids(dim, 7, scenario=1:3)             == [3,7,11]
        @test _FP.similar_ids(dim, 7, scenario=[1,3])           == [3,11]
        @test _FP.similar_ids(dim, 7, hour=[2,4], scenario=1:3) == [2,4,6,8,10,12]
        @test _FP.similar_ids(dim_shift, 31)                    == [31]
        @test _FP.similar_ids(dt, 7) == _FP.similar_ids(dim, 7)
        @test _FP.similar_ids(pm, 7) == _FP.similar_ids(dim, 7)
    end

    @testset "similar_id" begin
        @test _FP.similar_id(dim, 7)                     == 7
        @test _FP.similar_id(dim, 7, hour=4)             == 8
        @test _FP.similar_id(dim, 7, scenario=1)         == 3
        @test _FP.similar_id(dim, 7, hour=4, scenario=1) == 4
        @test _FP.similar_id(dim_shift, 31)              == 31
        @test _FP.similar_id(dt, 7) == _FP.similar_id(dim, 7)
        @test _FP.similar_id(pm, 7) == _FP.similar_id(dim, 7)
    end

    @testset "first_id" begin
        @test _FP.first_id(dim, 17, :hour) == 17
        @test _FP.first_id(dim, 18, :hour) == 17
        @test _FP.first_id(dim, 19, :hour) == 17
        @test _FP.first_id(dim, 20, :hour) == 17
        @test _FP.first_id(dim, 16, :scenario) == 16
        @test _FP.first_id(dim, 19, :scenario) == 15
        @test _FP.first_id(dim, 22, :scenario) == 14
        @test _FP.first_id(dim, 13, :hour, :scenario) == 13
        @test _FP.first_id(dim, 16, :hour, :scenario) == 13
        @test _FP.first_id(dim, 21, :hour, :scenario) == 13
        @test _FP.first_id(dim, 24, :hour, :scenario) == 13
        @test _FP.first_id(dim_shift, 41, :hour) == 41
        @test _FP.first_id(dt, 17, :hour) == _FP.first_id(dim, 17, :hour)
        @test _FP.first_id(pm, 17, :hour) == _FP.first_id(dim, 17, :hour)
    end

    @testset "last_id" begin
        @test _FP.last_id(dim,  8, :hour) ==  8
        @test _FP.last_id(dim,  7, :hour) ==  8
        @test _FP.last_id(dim,  6, :hour) ==  8
        @test _FP.last_id(dim,  5, :hour) ==  8
        @test _FP.last_id(dim,  9, :scenario) ==  9
        @test _FP.last_id(dim,  6, :scenario) == 10
        @test _FP.last_id(dim,  3, :scenario) == 11
        @test _FP.last_id(dim, 12, :hour, :scenario) == 12
        @test _FP.last_id(dim,  9, :hour, :scenario) == 12
        @test _FP.last_id(dim,  4, :hour, :scenario) == 12
        @test _FP.last_id(dim,  1, :hour, :scenario) == 12
        @test _FP.last_id(dim_shift, 32, :hour) == 32
        @test _FP.last_id(dt, 8, :hour) == _FP.last_id(dim, 8, :hour)
        @test _FP.last_id(pm, 8, :hour) == _FP.last_id(dim, 8, :hour)
    end

    @testset "prev_id" begin
        @test_throws BoundsError _FP.prev_id(dim, 17, :hour)
        @test _FP.prev_id(dim, 18, :hour) == 17
        @test _FP.prev_id(dim, 19, :hour) == 18
        @test _FP.prev_id(dim, 20, :hour) == 19
        @test_throws BoundsError _FP.prev_id(dim, 16, :scenario)
        @test _FP.prev_id(dim, 19, :scenario) == 15
        @test _FP.prev_id(dim, 22, :scenario) == 18
        @test _FP.prev_id(dim_shift, 42, :hour) == 41
        @test _FP.prev_id(dt, 18, :hour) == _FP.prev_id(dim, 18, :hour)
        @test _FP.prev_id(pm, 18, :hour) == _FP.prev_id(dim, 18, :hour)
    end

    @testset "prev_ids" begin
        @test _FP.prev_ids(dim, 17, :hour) == []
        @test _FP.prev_ids(dim, 18, :hour) == [17]
        @test _FP.prev_ids(dim, 20, :hour) == [17,18,19]
        @test _FP.prev_ids(dim, 16, :scenario) == []
        @test _FP.prev_ids(dim, 19, :scenario) == [15]
        @test _FP.prev_ids(dim, 22, :scenario) == [14,18]
        @test _FP.prev_ids(dim_shift, 42, :hour) == [41]
        @test _FP.prev_ids(dt, 17, :hour) == _FP.prev_ids(dim, 17, :hour)
        @test _FP.prev_ids(pm, 17, :hour) == _FP.prev_ids(dim, 17, :hour)
    end

    @testset "next_id" begin
        @test _FP.next_id(dim, 5, :hour) ==  6
        @test _FP.next_id(dim, 6, :hour) ==  7
        @test _FP.next_id(dim, 7, :hour) ==  8
        @test_throws BoundsError _FP.next_id(dim, 8, :hour)
        @test_throws BoundsError _FP.next_id(dim, 9, :scenario)
        @test _FP.next_id(dim, 6, :scenario) == 10
        @test _FP.next_id(dim, 3, :scenario) ==  7
        @test _FP.next_id(dim_shift, 29, :hour) == 30
        @test _FP.next_id(dt, 5, :hour) == _FP.next_id(dim, 5, :hour)
        @test _FP.next_id(pm, 5, :hour) == _FP.next_id(dim, 5, :hour)
    end

    @testset "next_ids" begin
        @test _FP.next_ids(dim, 5, :hour) == [6,7,8]
        @test _FP.next_ids(dim, 7, :hour) == [8]
        @test _FP.next_ids(dim, 8, :hour) == []
        @test _FP.next_ids(dim, 9, :scenario) == []
        @test _FP.next_ids(dim, 6, :scenario) == [10]
        @test _FP.next_ids(dim, 3, :scenario) == [7,11]
        @test _FP.next_ids(dim_shift, 29, :hour) == [30,31,32]
        @test _FP.next_ids(dt, 5, :hour) == _FP.next_ids(dim, 5, :hour)
        @test _FP.next_ids(pm, 5, :hour) == _FP.next_ids(dim, 5, :hour)
    end

    @testset "coord" begin
        @test _FP.coord(dim, 7, :hour) == 3
        @test _FP.coord(dim, 7, :scenario) == 2
        @test _FP.coord(dim_shift, 31, :hour) == 3
        @test _FP.coord(dim_shift, 31, :scenario) == 2
        @test _FP.coord(dt, 7, :hour) == _FP.coord(dim, 7, :hour)
        @test _FP.coord(pm, 7, :hour) == _FP.coord(dim, 7, :hour)
    end

    @testset "is_first_id" begin
        @test _FP.is_first_id(dim, 14, :hour) == false
        @test _FP.is_first_id(dim, 14, :scenario) == true
        @test _FP.is_first_id(dim, 17, :hour) == true
        @test _FP.is_first_id(dim, 17, :scenario) == false
        @test _FP.is_first_id(dim_shift, 38, :hour) == false
        @test _FP.is_first_id(dim_shift, 38, :scenario) == true
        @test _FP.is_first_id(dt, 14, :hour) == _FP.is_first_id(dim, 14, :hour)
        @test _FP.is_first_id(pm, 14, :hour) == _FP.is_first_id(dim, 14, :hour)
    end

    @testset "is_last_id" begin
        @test _FP.is_last_id(dim, 20, :hour) == true
        @test _FP.is_last_id(dim, 20, :scenario) == false
        @test _FP.is_last_id(dim, 21, :hour) == false
        @test _FP.is_last_id(dim, 21, :scenario) == true
        @test _FP.is_last_id(dim_shift, 44, :hour) == true
        @test _FP.is_last_id(dim_shift, 44, :scenario) == false
        @test _FP.is_last_id(dt, 20, :hour) == _FP.is_last_id(dim, 20, :hour)
        @test _FP.is_last_id(pm, 20, :hour) == _FP.is_last_id(dim, 20, :hour)
    end

    @testset "dim_prop" begin
        @test Set(keys(_FP.dim_prop(dim))) == Set((:hour, :scenario, :sub_nw))
        @test _FP.dim_prop(dim, :hour) == Dict(h => Dict{String,Any}() for h in 1:4)
        @test _FP.dim_prop(dim, :scenario) == Dict(s => Dict{String,Any}("probability"=>s/6) for s in 1:3)
        @test _FP.dim_prop(dim, :scenario, 1) == Dict{String,Any}("probability"=>1/6)
        @test _FP.dim_prop(dim, :scenario, 1, "probability") == 1/6
        @test _FP.dim_prop(dim, 13, :scenario) == Dict{String,Any}("probability"=>1/6)
        @test _FP.dim_prop(dim, 13, :scenario, "probability") == 1/6
        @test _FP.dim_prop(dt) == _FP.dim_prop(dim)
        @test _FP.dim_prop(pm) == _FP.dim_prop(dim)
    end

    @testset "dim_meta" begin
        @test Set(keys(_FP.dim_meta(dim))) == Set((:hour, :scenario, :sub_nw))
        @test _FP.dim_meta(dim, :hour) == Dict{String,Any}()
        @test _FP.dim_meta(dim, :sub_nw) == Dict{String,Any}("description" => "sub_nws model different physical networks")
        @test _FP.dim_meta(dim, :sub_nw, "description") == "sub_nws model different physical networks"
        @test _FP.dim_meta(dt) == _FP.dim_meta(dim)
        @test _FP.dim_meta(pm) == _FP.dim_meta(dim)
    end

    @testset "dim_length" begin
        @test _FP.dim_length(dim) == 24
        @test _FP.dim_length(dim, :hour) == 4
        @test _FP.dim_length(dim, :scenario) == 3
        @test _FP.dim_length(dim_shift) == 24
        @test _FP.dim_length(dt) == _FP.dim_length(dim)
        @test _FP.dim_length(pm) == _FP.dim_length(dim)
    end

    Memento.setlevel!(Memento.getlogger(FlexPlan), previous_FlexPlan_logger_level)

end;


text:
#=

  Euler #27 in Julia.

  Problem 27.
  """
  Euler published the remarkable quadratic formula:

  n^2 + n + 41

  It turns out that the formula will produce 40 primes for the consecutive values
  n = 0 to 39. However, when n = 40, 402 + 40 + 41 = 40(40 + 1) + 41 is divisible by
  41, and certainly when n = 41, 41^2 + 41 + 41 is clearly divisible by 41.

  Using computers, the incredible formula  n^2 − 79n + 1601 was discovered, which
  produces 80 primes for the consecutive values n = 0 to 79. The product of the
  coefficients, −79 and 1601, is −126479.

  Considering quadratics of the form:

      n^2 + an + b, where |a| < 1000 and |b| < 1000

      where |n| is the modulus/absolute value of n
      e.g. |11| = 11 and |−4| = 4

  Find the product of the coefficients, a and b, for the quadratic
  expression that produces the maximum number of primes for consecutive
  values of n, starting with n = 0.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.04097529s
# with isPrimeCached: 0.07010393s (slower)
function euler27a()
    t = 999
    bestLen = 0
    bestA = 0
    bestB = 0
    for a in -t:t
        for b in -t:t
            len = 0;
            pp = len^2 + a*len + b
            # while pp > 1 && isPrimeCached(pp) # slower
            while pp > 1 && isPrime(pp)
                len += 1
                pp = len^2 + a*len + b;
            end
            if len > bestLen
                bestLen = len
                bestA = a
                bestB = b
            end
        end
    end

    return(bestA * bestB)
end

run_euler(euler27a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module TestBootsampleNonParametric

using Bootstrap
using Test

using Statistics
using Random

using DataFrames
using RDatasets
using StatsBase

## 'city' dataset
const city = dataset("boot", "city")
const citya = convert(Matrix, city)

@testset "Non-parametric bootstraps" begin

    function test_bootsample(bs, ref, raw_data, n)

        show(IOBuffer(), bs)
        @test typeof(bs) <: NonParametricBootstrapSample
        t0 = original(bs)
        @test length(t0) == length(ref)
        [@test t ≈ r for (t, r) in zip(t0, ref)]

        t1 = straps(bs)
        @test length(t1) == length(t0)
        [@test length(t) == n for t in t1]
        [@test (minimum(t) <= tr && maximum(t) >= tr) for (t, tr) in zip(t1, t0)]
        [@test eltype(t) == eltype(tr) for (t, tr) in zip(t1, t0)]

        @test Bootstrap.data(bs) == raw_data

        @test nrun(sampling(bs)) == n
        @test nrun(bs) == n

        @test length(bias(bs)) == length(ref)
        [@test (b[1] > -Inf && b[1] < Inf) for b in bias(bs)]
        @test length(stderror(bs)) == length(ref)
        [@test s >= 0 for s in stderror(bs)]

        [@test original(bs, i) == original(bs)[i]  for i in 1:nvar(bs)]
        [@test straps(bs, i) == straps(bs)[i]  for i in 1:nvar(bs)]
        [@test bias(bs, i) == bias(bs)[i]  for i in 1:nvar(bs)]
        [@test stderror(bs, i) == stderror(bs)[i]  for i in 1:nvar(bs)]

        @test_throws MethodError model(bs)

        return Nothing
    end

    function test_confint(bs)

        cim_all = (BasicConfInt(), PercentileConfInt(), NormalConfInt(), BCaConfInt())
        for cim in cim_all
            c = confint(bs, cim)
            [@test (x[1] >= x[2]  && x[1] <= x[3]) for x in c]
            [@test x[1] ≈ t0 for (x, t0) in zip(c, original(bs))]
            @test level(cim) == 0.95
        end

        return Nothing
    end

    function test_confint_width0(bs)

        cim_all = (BasicConfInt(), PercentileConfInt(), NormalConfInt(), BCaConfInt())
        for cim in cim_all
            c = confint(bs, cim)
            # We expect a confidence interval of width 0
            [@test x[2] == x[1] == x[3] for x in c]
            [@test x[1] ≈ t0 for (x, t0) in zip(c, original(bs))]
            @test level(cim) == 0.95
        end

        return Nothing
    end

    n = 250

    city_ratio(x::AbstractArray) = mean(x[:,2]) ./ mean(x[:,1])
    city_ratio(x::AbstractDataFrame) = mean(x[!,:X]) ./ mean(x[!,:U])

    city_cor(x::AbstractArray) = cor(x[:,1], x[:,2])
    city_cor(x::AbstractDataFrame) = cor(x[!,:X], x[!,:U])

    @testset "Basic resampling" begin

        @testset "city_ratio with DataFrame input" begin
            ref = city_ratio(city)
            @test ref  ≈ 1.5203125
            bs = bootstrap(city_ratio, city, BasicSampling(n))
            test_bootsample(bs, ref, city, n)
            test_confint(bs)
        end

        @testset "city_cor with DataFrame input" begin
            ref = city_cor(city)
            bs = bootstrap(city_cor, city, BasicSampling(n))
            test_bootsample(bs, ref, city, n)
            test_confint(bs)
        end

        @testset "city_cor with DataArray input" begin
            ref = city_cor(citya)
            bs = bootstrap(city_cor, citya, BasicSampling(n))
            test_bootsample(bs, ref, citya, n)
            test_confint(bs)
        end

        @testset "mean_and_sd: Vector input, 2 output variables" begin
            r = randn(25)
            ref = mean_and_std(r)
            bs = bootstrap(mean_and_std, r, BasicSampling(n))
            test_bootsample(bs, ref, r, n)
            test_confint(bs)
        end

        @testset "mean_and_sd: Student CI" begin
            r = randn(25)
            bs = bootstrap(mean_and_std, r, BasicSampling(n))
            ## Student confint
            cim = StudentConfInt()
            c = confint(bs, straps(bs, 2), cim, 1)
            @test c[1] >= c[2]  && c[1] <= c[3]
            @test c[1] ≈ original(bs, 1)
            @test level(cim) == 0.95
        end
    end

    @testset "Antithetic resampling" begin

        @testset "mean_and_sd: Vector input, 2 output variables" begin
            r = randn(50)
            ref = mean_and_std(r)
            bs = bootstrap(mean_and_std, r, AntitheticSampling(n))
            test_bootsample(bs, ref, r, n)
            test_confint(bs)
        end

    end

    @testset "Balanced resampling" begin

        @testset "city_ratio with DataFrame input" begin
            ref = city_ratio(city)
            @test ref ≈ 1.5203125
            bs = bootstrap(city_ratio, city, BalancedSampling(n))
            test_bootsample(bs, ref, city, n)
            test_confint(bs)
        end

        @testset "city_cor with DataFrame input" begin
            ref = city_cor(city)
            bs = bootstrap(city_cor, city, BalancedSampling(n))
            test_bootsample(bs, ref, city, n)
            test_confint(bs)
        end

        @testset "city_cor with DataArray input" begin
            ref = city_cor(citya)
            bs = bootstrap(city_cor, citya, BalancedSampling(n))
            test_bootsample(bs, ref, citya, n)
            test_confint(bs)
        end

        @testset "mean_and_sd: Vector input, 2 output variables" begin
            r = randn(50)
            ref = mean_and_std(r)
            bs = bootstrap(mean_and_std, r, BalancedSampling(n))
            test_bootsample(bs, ref, r, n)
            test_confint(bs)
            ## mean should be unbiased
            @test isapprox(bias(bs)[1], 0.0, atol = 1e-8)
        end

    end

    @testset "Exact resampling" begin

        nc = Bootstrap.nrun_exact(nrow(city))

        @testset "city_ratio with DataFrame input" begin
            ref = city_ratio(city)
            @test ref ≈ 1.5203125
            bs = bootstrap(city_ratio, city, ExactSampling())
            test_bootsample(bs, ref, city, nc)
            test_confint(bs)
        end

        @testset "mean: Vector input, 1 output variables" begin
            r = randn(10)
            ref = mean(r)
            bs = bootstrap(mean, r, ExactSampling())
            test_bootsample(bs, ref, r, nc)
            test_confint(bs)
        end

    end

    @testset "Maximum Entropy resampling" begin

        # Simulated AR(1) data copied from
        # https://github.com/colintbowers/DependentBootstrap.jl/blob/master/test/runtests.jl#L8

        nobs = 100

        function test_obs(n, seed = 1234)
            Random.seed!(seed)
            e = randn(n)
            x = zeros(Float64, n)
            for i = 2:n
                x[i] = 0.8 * x[i - 1] + e[i]
            end
            return sin.(x)
        end

        r = test_obs(nobs)
        ref = mean(r)
        s = MaximumEntropySampling(n)
        bs = bootstrap(mean, r, s)
        test_bootsample(bs, ref, r, n)
        test_confint(bs)

        # Collect the samples
        samples = zeros(eltype(r), (nobs, n))
        for i in 1:n
            # For some reason the samples are only filled in if we have
            # an explicit assignment back into the matrix.
            samples[:, i] = draw!(s.cache, r, samples[:, i])
        end

        # Add some checks to ensure that our within sample variation is greater than our
        # across sample variation at any given "timestep".
        @test all(std(samples, dims = 2) .< std(r))
        @test mean(std(samples, dims = 2)) < 0.1  # NOTE: This is about 0.09 in julia and 0.08 in the R package
        @test std(r) > 0.5

    end

    @testset "ConfInt" begin
        @testset "All equal" begin
            r = rand(10)
            # In this case our count will always be the same because all values are positive
            bs = bootstrap(x->count(x .> 0.0), r, BasicSampling(n))
            test_confint_width0(bs)
        end
    end
end

end


text:
amap = Dict(stay=>[ss_state1=>ss_state1])
fsmsingle = FiniteStateMachine(amap,ss_state1,ss_state1)
@test fire!(fsmsingle,stay)
@test finished(fsmsingle)
@test can(fsmsingle,stay)

amap = Dict(
startup => [ts_state1=>ts_state1],
st1_1 => [ts_state1=>ts_state1],
st1_2 => [ts_state1=>ts_state2],
st2_2 => [ts_state2=>ts_state2]
)
fsmtwostate = FiniteStateMachine(amap,ts_state1,ts_state2)

@test !finished(fsmtwostate)

@test can(fsmtwostate,startup)
@test fire!(fsmtwostate,startup)
@test fsmtwostate.current == ts_state1
@test !finished(fsmtwostate)

@test can(fsmtwostate,st1_1)
@test fire!(fsmtwostate,st1_1)
@test fsmtwostate.current == ts_state1
@test !finished(fsmtwostate)

@test can(fsmtwostate,st1_2)
@test fire!(fsmtwostate,st1_2)
@test fsmtwostate.current == ts_state2
@test finished(fsmtwostate)

@test !can(fsmtwostate,st2_1)
@test !fire!(fsmtwostate,st2_1)
@test fsmtwostate.current == ts_state2
@test finished(fsmtwostate)

@test !can(fsmtwostate,st1_2)
@test !can(fsmtwostate,st1_1)
@test !can(fsmtwostate,startup)

@test can(fsmtwostate,st2_2)
@test fire!(fsmtwostate,st2_2)
@test fsmtwostate.current == ts_state2
@test finished(fsmtwostate)


text:
using MOCNeutronTransport
@testset "Triangle6_2D" begin
    for F in [Float32, Float64]
        @testset "Constructors" begin
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 1)
            p₃ = Point_2D(F, 1, 1)
            p₄ = Point_2D(F, 1//2)
            p₅ = Point_2D(F, 1, 1//2)
            p₆ = Point_2D(F, 1//2, 1//2)
            tri6 = Triangle6_2D(SVector(p₁, p₂, p₃, p₄, p₅, p₆))
            @test tri6.points == SVector(p₁, p₂, p₃, p₄, p₅, p₆)

            # single constructor
            tri6 = Triangle6_2D(p₁, p₂, p₃, p₄, p₅, p₆)
            @test tri6.points == SVector(p₁, p₂, p₃, p₄, p₅, p₆)
        end

        @testset "Methods" begin
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 1)
            p₃ = Point_2D(F, 0, 1)
            p₄ = Point_2D(F, 1//2)
            p₅ = Point_2D(F, 1//2, 1//2)
            p₆ = Point_2D(F, 0, 1//2)
            tri6 = Triangle6_2D(p₁, p₂, p₃, p₄, p₅, p₆)

            # interpolation
            @test tri6(0, 0) ≈ p₁
            @test tri6(1, 0) ≈ p₂
            @test tri6(0, 1) ≈ p₃
            @test tri6(1//2, 0) ≈ p₄
            @test tri6(1//2, 1//2) ≈ p₅
            @test tri6(0, 1//2) ≈ p₆

            # derivative
            r, s = derivative(tri6, 0, 0)
            @test r ≈ Point_2D(F, 1, 0)
            @test s ≈ Point_2D(F, 0, 1)
            r, s = derivative(tri6, 1, 0)
            @test r ≈ Point_2D(F, 1, 0)
            @test s ≈ Point_2D(F, 0, 1)
            r, s = derivative(tri6, 0, 1)
            @test r ≈ Point_2D(F, 1, 0)
            @test s ≈ Point_2D(F, 0, 1)
            r, s = derivative(tri6, 1//2, 1//2)
            @test r ≈ Point_2D(F, 1, 0)
            @test s ≈ Point_2D(F, 0, 1)

            # area
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 2)
            p₃ = Point_2D(F, 2, 2)
            p₄ = Point_2D(F, 3//2, 1//4)
            p₅ = Point_2D(F, 3, 1)
            p₆ = Point_2D(F, 1, 1)
            tri6 = Triangle6_2D(p₁, p₂, p₃, p₄, p₅, p₆)
            # 2D default
            @test isapprox(area(tri6), 3, atol=1.0e-6)

            # real_to_parametric
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 2)
            p₃ = Point_2D(F, 2, 2)
            p₄ = Point_2D(F, 3//2, 1//4)
            p₅ = Point_2D(F, 3, 1)
            p₆ = Point_2D(F, 1, 1)
            tri6 = Triangle6_2D(p₁, p₂, p₃, p₄, p₅, p₆)
            @test tri6(real_to_parametric(p₁, tri6)) ≈ p₁
            @test tri6(real_to_parametric(p₂, tri6)) ≈ p₂
            @test tri6(real_to_parametric(p₃, tri6)) ≈ p₃
            @test tri6(real_to_parametric(p₄, tri6)) ≈ p₄
            @test tri6(real_to_parametric(p₅, tri6)) ≈ p₅
            @test tri6(real_to_parametric(p₆, tri6)) ≈ p₆

            # in
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 2)
            p₃ = Point_2D(F, 2, 2)
            p₄ = Point_2D(F, 3//2, 1//4)
            p₅ = Point_2D(F, 3, 1)
            p₆ = Point_2D(F, 1, 1)
            tri6 = Triangle6_2D(p₁, p₂, p₃, p₄, p₅, p₆)
            @test Point_2D(F, 1, 1//2) ∈  tri6
            @test Point_2D(F, 1, 0) ∉  tri6

            # intersect
            # 0 intersection
            l = LineSegment_2D(Point_2D(F, 0, -1), Point_2D(F, 4, -1))
            n, points = l ∩ tri6
            @test n == 0

            # 2 intersection
            l = LineSegment_2D(Point_2D(F, 0, 1), Point_2D(F, 4, 1))
            n, points = l ∩ tri6
            @test n == 2
            @test points[1] ≈ Point_2D(F, 3, 1)
            @test points[2] ≈ Point_2D(F, 1, 1)

            # 4 intersection
            l = LineSegment_2D(Point_2D(F, 0, 1//10), Point_2D(F, 4, 1//10))
            n, points = l ∩ tri6
            @test n == 4
            @test points[1] ≈ Point_2D(F, 0.4254033307585166, 1//10)
            @test points[2] ≈ Point_2D(F, 1.9745966692414834, 1//10)
            @test points[3] ≈ Point_2D(F, 2.1900000000000000, 1//10)
            @test points[4] ≈ Point_2D(F, 1//10,              1//10)

            # 6 intersection
            p₁ = Point_2D(F,  1, 0)
            p₂ = Point_2D(F,  0, 0)
            p₃ = Point_2D(F, -1, 0)
            p₄ = Point_2D(F,  1//2, -1//2)
            p₅ = Point_2D(F, -1//2, -1//2)
            p₆ = Point_2D(F,   0,    -2)
            tri6 = Triangle6_2D(p₁, p₂, p₃, p₄, p₅, p₆)
            l = LineSegment_2D(Point_2D(F, -2, -1//4), Point_2D(F, 2, -1//4))
            n, points = l ∩ tri6
            @test n == 6
            @test points[1] ≈ Point_2D(F,  0.14644659, -1//4)
            @test points[2] ≈ Point_2D(F,  0.8535534,  -1//4)
            @test points[3] ≈ Point_2D(F, -0.8535534,  -1//4)
            @test points[4] ≈ Point_2D(F, -0.14644665, -1//4)
            @test points[5] ≈ Point_2D(F,  0.9354143,  -1//4)
            @test points[6] ≈ Point_2D(F, -0.9354143,  -1//4)
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"Action consisting of longitudinal acceleration and lateral speed"
struct AccelSpeed <: FieldVector{2, Float64}
    long_accel::Float64
    lat_speed::Float64
end

"Action consisting of a longitudinal acceleration"
struct LongAccel
    a::Float64
end

# state: [x, y, xdot]; y = 0 is middle of car in middle of rightmost lane

struct HighwayVehicle end

function propagate(m::HighwayVehicle, s::AbstractVector{Float64}, a::AccelSpeed, context, dt)
    xddot = a[1]
    ydot = a[2]
    x = s[1]
    xdot = s[3]
    xp = x + xdot*dt + 0.5*xddot*dt^2
    xdotp = xdot + xddot*dt
    y = s[2]
    yp = y + ydot*dt
    return SVector(xp, yp, xdotp)
end

brake_limit(v::HighwayVehicle) = 8.0
Base.length(v::HighwayVehicle) = 4.8
width(v::HighwayVehicle) = 1.8

pos(v::AbstractVector{Float64}) = v[1:2]
# XXX :/
vel(v::AbstractVector{Float64}) = v[3]