text:
mutable struct PGFPropertyFrame
    # Vector properties in this frame.
    vector_properties::Dict{Type, Property}

    # True if this property frame has scalar properties. Scalar properties are
    # emitted as a group {scope} that must be closed when the frame is popped.
    has_scalar_properties::Bool
end

PGFPropertyFrame() = PGFPropertyFrame(Dict{Type, Property}(), false)

mutable struct PGF <: Backend
    # Image size in millimeters.
    width::AbsoluteLength
    height::AbsoluteLength

    # Output stream.
    out::IO

    # Fill properties cannot be "cleanly" applied to
    # multiple form primitives.  It must be applied
    # each time an object is drawn
    fill::Union{Color, Nothing}
    fill_opacity::Float64

    stroke::Union{Color, Nothing}
    stroke_opacity::Float64

    fontfamily::Union{AbstractString, Nothing}
    fontsize::Float64

    # Current level of indentation.
    indentation::Int

    # Output buffer.  We want the ability to add to the beginning of
    buf::IOBuffer

    # Skip drawing if visible is false
    visible::Bool

    # Have not found an easy way to define color as
    # a draw parameter.  Whenever we encounter a color, we add it to the
    # color_set set.  That way, we can write out all the color
    # definitions at the same time.
    color_set::Set{Color}

    # Stack of property frames (groups of properties) currently in effect.
    property_stack::Vector{PGFPropertyFrame}

    # SVG forbids defining the same property twice, so we have to keep track
    # of which vector property of which type is in effect. If two properties of
    # the same type are in effect, the one higher on the stack takes precedence.
    vector_properties::Dict{Type, Union{Property, Nothing}}

    # Clip-paths that need to be defined at the end of the document.
    # Not quite sure how to deal with clip paths yet
    clippath::Union{ClipPrimitive, Nothing}
    # clippaths::Dict{ClipPrimitive, String}

    # True when finish has been called and no more drawing should occur
    finished::Bool

    # Backend is responsible for opening/closing the file
    ownedfile::Bool

    # Filename when ownedfile is true
    filename::Union{AbstractString, Nothing}

    # Emit the graphic on finish when writing to a buffer.
    emit_on_finish::Bool

    # Emit only the tikzpicture environment
    only_tikz::Bool

    # Use default TeX fonts instead of fonts specified by the theme.
    texfonts::Bool
end

function PGF(out::IO,
             width::AbsoluteLength,
             height::AbsoluteLength,
             emit_on_finish::Bool=true,
             only_tikz = false;

             texfonts = false,
             fill = default_fill_color,
             fill_opacity = 1.0,
             stroke = default_stroke_color,
             stroke_opacity = 1.0,
             fontfamily = nothing,
             fontsize = 12.0,
             indentation = 0,
             buf = IOBuffer(),
             visible = true,
             color_set = Set{Color}([colorant"black"]),
             property_stack = Array{PGFPropertyFrame}(undef, 0),
             vector_properties = Dict{Type, Union{Property, Nothing}}(),
             clippath = nothing,
             finished = false,
             ownedfile = false,
             filename = nothing)
    PGF(width, height, out,
        fill, fill_opacity,
        stroke, stroke_opacity,
        fontfamily, fontsize,
        indentation,
        buf,
        visible,
        color_set,
        property_stack,
        vector_properties,
        clippath,
        finished,
        ownedfile,
        filename,
        emit_on_finish,
        only_tikz,
        texfonts)
end

# Write to a file.
PGF(filename::AbstractString, width=default_graphic_width, height=default_graphic_height,
            only_tikz=false; texfonts=false) =
        PGF(open(filename, "w"), width, height, true, only_tikz;
            texfonts = texfonts, ownedfile=true, filename=filename)

# Write to buffer.
PGF(width::MeasureOrNumber=default_graphic_width, height::MeasureOrNumber=default_graphic_height,
            emit_on_finish::Bool=true, only_tikz=false; texfonts=false) =
        PGF(IOBuffer(), width, height, emit_on_finish, only_tikz, texfonts=texfonts)

function (img::PGF)(x)
    draw(img, x)
end

function finish(img::PGF)
    img.finished && return

    while !isempty(img.property_stack)
        pop_property_frame(img)
    end

    # if length(img.clippaths) > 0
    #     for (clippath, id) in img.clippaths
    #         write(img.out, "\\clip")
    #         print_svg_path(img.out, clippath.points)
    #         write(img.out, ";\n")
    #     end
    # end

    writeheader(img)
    writecolors(img)
    write(img.out, take!(img.buf))
    write(img.out,
        """
        \\end{tikzpicture}
        """
        )
    !img.only_tikz &&  write(img.out,
        """
        \\end{document}
        """
        )

    hasmethod(flush, (typeof(img.out),)) && flush(img.out)
    close(img.buf)

    img.ownedfile && close(img.out)

    img.finished = true

    # If we are writing to a buffer. Collect the string and emit it.
    img.emit_on_finish && typeof(img.out) == IOBuffer && display(img)
end

isfinished(img::PGF) = img.finished
root_box(img::PGF) = BoundingBox(0mm, 0mm, img.width, img.height)

function writeheader(img::PGF)
    !img.only_tikz && write(img.out,
        """
        \\documentclass{minimal}
        \\usepackage{pgfplots}
        $(img.texfonts ? "" : "\\usepackage{fontspec}")
        \\usepackage{amsmath}
        \\usepackage[active,tightpage]{preview}
        \\PreviewEnvironment{tikzpicture}
        \\begin{document}
        """)
    write(img.out, """
        \\begin{tikzpicture}[x=1mm,y=-1mm]
        """)
    return img
end

function writecolors(img::PGF)
    for color in img.color_set
        @printf(img.out, "\\definecolor{mycolor%s}{rgb}{%s,%s,%s}\n",
            hex(color),
            svg_fmt_float(color.r),
            svg_fmt_float(color.g),
            svg_fmt_float(color.b)
            )
    end
end

function print_pgf_path(out::IO, points::Vector{AbsoluteVec2},
                        bridge_gaps::Bool=false)
    isfirst = true
    for point in points
        x, y = point[1].value, point[2].value
        if !(isfinite(x) && isfinite(y))
            isfirst = true
            continue
        end

        if isfirst
            isfirst = false
            @printf(out, " (%s,%s)",
                    svg_fmt_float(x),
                    svg_fmt_float(y))
        else
            @printf(out, " -- (%s,%s)",
                    svg_fmt_float(x),
                    svg_fmt_float(y))
        end
    end
end

function get_vector_properties(img::PGF, idx::Int)
    props_str = String[]
    modifiers = String[]
    for (propertytype, property) in img.vector_properties
        (property === nothing) && continue
        primitives = property.primitives
        idx > length(primitives) &&
                error("Vector form and vector property differ in length. Can't distribute.")
        push_property!(props_str, img, primitives[idx])
    end

    if img.fill !== nothing
        push!(props_str, string("fill=mycolor",hex(img.fill)))
        img.fill_opacity < 1.0 &&
                push!(props_str, string("fill opacity=",svg_fmt_float(img.fill_opacity)))
    end

    if img.stroke !== nothing
        push!(props_str, string("draw=mycolor",hex(img.stroke)))
        img.stroke_opacity < 1.0 &&
                push!(props_str, string("draw opacity=",svg_fmt_float(img.stroke_opacity)))
    end

    return modifiers, props_str
end

push_property!(props_str, img::PGF, property::StrokeDashPrimitive) =
        isempty(property.value) ||
            push!(props_str, string("dash pattern=", join(map( v -> join(v, " "),
                zip(Iterators.cycle(["on", "off"]),
                map(v -> string(svg_fmt_float(v.value), "mm"), property.value)) )," ")))

function push_property!(props_str, img::PGF, property::StrokePrimitive)
    if isa(property.color, TransparentColor)
        img.stroke = color(property.color)
        img.stroke_opacity = property.color.alpha
    else
        img.stroke = property.color
    end

    (img.stroke === nothing) || push!(img.color_set, convert(RGB, property.color))
end

function push_property!(props_str, img::PGF, property::FillPrimitive)
    if isa(property.color, TransparentColor)
        img.fill = color(property.color)
        img.fill_opacity = property.color.alpha
    else
        img.fill = property.color
    end

    (img.fill === nothing) || push!(img.color_set, convert(RGB, property.color))
end

push_property!(props_str, img::PGF, property::VisiblePrimitive) =
        img.visible = property.value
push_property!(props_str, img::PGF, property::LineWidthPrimitive) =
        push!(props_str, string("line width=", svg_fmt_float(property.value.value), "mm"))

pgf_fmt_linecap(::LineCapButt) = "butt"
pgf_fmt_linecap(::LineCapSquare) = "rect"
pgf_fmt_linecap(::LineCapRound) = "round"

push_property!(props_str, img::PGF, property::StrokeLineCapPrimitive) =
        push!(props_str, string("line cap=", pgf_fmt_linecap(property.value)))
push_property!(props_str, img::PGF, property::StrokeLineJoinPrimitive) =
        push!(props_str, string("line join=", svg_fmt_linejoin(property.value)))
push_property!(props_str, img::PGF, property::FillOpacityPrimitive) =
        img.fill_opacity = property.value
push_property!(props_str, img::PGF, property::StrokeOpacityPrimitive) =
        img.stroke_opacity = property.value
# Can only only work with one font family for now
push_property!(props_str, img::PGF, property::FontPrimitive) =
        img.fontfamily = strip(split(escape_string(property.family),',')[1],'\'')
push_property!(props_str, img::PGF, property::FontSizePrimitive) =
        img.fontsize = property.value.value
# Not quite sure how to handle clipping yet, stub for now
push_property!(props_str, img::PGF, property::ClipPrimitive) =
        img.clippath = property

# Stubs for SVG and JS specific properties
push_property!(props_str, img::PGF, property::JSIncludePrimitive) = nothing
push_property!(props_str, img::PGF, property::JSCallPrimitive) = nothing
push_property!(props_str, img::PGF, property::SVGClassPrimitive) = nothing
push_property!(props_str, img::PGF, property::SVGAttributePrimitive) = nothing
iswithjs(img::PGF) = false
iswithousjs(img::PGF) = true

push_property!(props_str, img::PGF, property::ArrowPrimitive) =
        push!(props_str, "arrows=->")

# Form Drawing
# ------------

function draw(img::PGF, form::Form)
    for (idx, primitive) in enumerate(form.primitives)
        draw(img, primitive, idx)
    end
end

function draw(img::PGF, prim::LinePrimitive, idx::Int)
    n = length(prim.points)
    n <= 1 && return

    modifiers, props = get_vector_properties(img, idx)
    img.visible || return

    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] ", join(props, ","));
    print_pgf_path(img.buf, prim.points, true)
    write(img.buf, ";\n")
end

function draw(img::PGF, prim::RectanglePrimitive, idx::Int)
    width = max(prim.width.value, 0.01)
    height = max(prim.height.value, 0.01)

    modifiers, props = get_vector_properties(img, idx)
    img.visible || return

    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] ", join(props, ","));
    @printf(img.buf, "(%s,%s) rectangle +(%s,%s);\n",
            svg_fmt_float(prim.corner[1].value),
            svg_fmt_float(prim.corner[2].value),
            svg_fmt_float(width),
            svg_fmt_float(height))
end

function draw(img::PGF, prim::PolygonPrimitive, idx::Int)
    n = length(prim.points)
    n <= 1 || return

    modifiers, props = get_vector_properties(img, idx)
    img.visible || return

    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] ", join(props, ","))
    print_pgf_path(img.buf, prim.points, true)
    write(img.buf, " -- cycle;\n")
end

function draw(img::PGF, prim::CirclePrimitive, idx::Int)
    modifiers, props = get_vector_properties(img, idx)
    img.visible || return
    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] ", join(props, ","))
    @printf(img.buf, "(%s,%s) circle [radius=%s];\n",
        svg_fmt_float(prim.center[1].value),
        svg_fmt_float(prim.center[2].value),
        svg_fmt_float(prim.radius.value))
end

function draw(img::PGF, prim::EllipsePrimitive, idx::Int)
    modifiers, props = get_vector_properties(img, idx)
    img.visible || return

    cx = prim.center[1].value
    cy = prim.center[2].value
    rx = sqrt((prim.x_point[1].value - cx)^2 +
              (prim.x_point[2].value - cy)^2)
    ry = sqrt((prim.y_point[1].value - cx)^2 +
              (prim.y_point[2].value - cy)^2)
    theta = rad2deg(atan(prim.x_point[2].value - cy,
                         prim.x_point[1].value - cx))

    all(isfinite,[cx, cy, rx, ry, theta]) || return

    modifiers, props = get_vector_properties(img, idx)
    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] ", join(props, ","))
    @printf(img.buf, "(%s,%s) circle [x radius=%s, y radius=%s",
        svg_fmt_float(cx),
        svg_fmt_float(cy),
        svg_fmt_float(rx),
        svg_fmt_float(ry))
    abs(theta) > 1e-4 && @printf(img.buf, " rotate=%s", svg_fmt_float(theta))
    write(img.buf, "];\n")
end

function draw(img::PGF, prim::CurvePrimitive, idx::Int)
    modifiers, props = get_vector_properties(img, idx)
    img.visible || return
    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] ", join(props, ","))
    @printf(img.buf, "(%s,%s) .. controls (%s,%s) and (%s,%s) .. (%s,%s);\n",
        svg_fmt_float(prim.anchor0[1].value),
        svg_fmt_float(prim.anchor0[2].value),
        svg_fmt_float(prim.ctrl0[1].value),
        svg_fmt_float(prim.ctrl0[2].value),
        svg_fmt_float(prim.ctrl1[1].value),
        svg_fmt_float(prim.ctrl1[2].value),
        svg_fmt_float(prim.anchor1[1].value),
        svg_fmt_float(prim.anchor1[2].value))
end

function draw(img::PGF, prim::TextPrimitive, idx::Int)
    # Rotation direction is reversed!
    modifiers, props = get_vector_properties(img, idx)
    img.visible || return
    push!(props, string(
        "rotate around={",
            svg_fmt_float(-rad2deg(prim.rot.theta)),
            ": (",
            svg_fmt_float(prim.rot.offset[1].value - prim.position[1].value),
            ",",
            svg_fmt_float(prim.rot.offset[2].value - prim.position[2].value),
            ")}"))
    if prim.halign === hcenter
        push!(props, "inner sep=0.0")
    elseif prim.halign === hright
        push!(props, "left,inner sep=0.0")
    else
        push!(props, "right,inner sep=0.0")
    end
    write(img.buf, join(modifiers))
    @printf(img.buf, "\\draw (%s,%s) node [%s]{\\fontsize{%smm}{%smm}\\selectfont \$%s\$};\n",
        svg_fmt_float(prim.position[1].value),
        svg_fmt_float(prim.position[2].value),
        replace(join(props, ","), "fill"=>"text"),
        svg_fmt_float(img.fontsize),
        svg_fmt_float(1.2*img.fontsize),
        pango_to_pgf(prim.value)
    )
end

function indent(img::PGF)
    for i in 1:img.indentation
        write(img.buf, "  ")
    end
end

function push_property_frame(img::PGF, properties::Vector{Property})
    isempty(properties) && return

    frame = PGFPropertyFrame()
    applied_properties = Set{Type}()
    scalar_properties = Array{Property}(undef, 0)
    for property in properties
        if !isrepeatable(property) && (typeof(property) in applied_properties)
            continue
        elseif isscalar(property)
            push!(scalar_properties, property)
            push!(applied_properties, typeof(property))
            frame.has_scalar_properties = true
            img.vector_properties[typeof(property)] = nothing
        else
            frame.vector_properties[typeof(property)] = property
            img.vector_properties[typeof(property)] = property
        end
    end
    push!(img.property_stack, frame)
    isempty(scalar_properties) && return

    write(img.buf, "\\begin{scope}\n")
    prop_str = AbstractString[]
    for property in scalar_properties
        push_property!(prop_str, img, property.primitives[1])
    end
    if length(prop_str) > 0
        @printf(img.buf, "[%s]\n", join(prop_str, ","))
    end
    if img.clippath !== nothing
        write(img.buf, "\\clip ")
        print_pgf_path(img.buf, img.clippath.points)
        write(img.buf, ";\n")
    end
    if (!img.texfonts && (img.fontfamily !== nothing))
        @printf(img.buf, "\\fontspec{%s}\n", img.fontfamily)
    end
end

function pop_property_frame(img::PGF)
    @assert !isempty(img.property_stack)
    frame = pop!(img.property_stack)

    if frame.has_scalar_properties
        write(img.buf, "\\end{scope}\n")
        # There should be a better way to to this:
        # Maybe put the applied properties on a stack then
        # pop them out here?
        # FIX ME!
        img.fill = default_fill_color
        img.stroke = default_stroke_color
        img.fill_opacity = 1.0
        img.stroke_opacity = 1.0
        img.fontfamily = nothing
        img.clippath = nothing
        img.visible = true
    end

    for (propertytype, property) in frame.vector_properties
        img.vector_properties[propertytype] = nothing
        for i in length(img.property_stack):-1:1
            if haskey(img.property_stack[i].vector_properties, propertytype)
                img.vector_properties[propertytype] =
                        img.property_stack[i].vector_properties[propertytype]
            end
        end
    end
end

# Horrible abuse of Latex inline math mode just to
# get something working first.
# FIX ME!
function pango_to_pgf(text::AbstractString)
    pat = r"<(/?)\s*([^>]*)\s*>"
    input = codeunits(escape_tex_chars(text))
    output = IOBuffer()
    lastpos = 1
    for mat in eachmatch(pat, text)
        write(output, "\\text{")
        write(output, input[lastpos:mat.offset-1])
        write(output, "}")

        if mat.captures[2] == "sup"
            if mat.captures[1] == "/"
                write(output, "}")
            else
                write(output, "^{")
            end
        elseif mat.captures[2] == "sub"
            if mat.captures[1] == "/"
                write(output, "}")
            else
                write(output, "_{")
            end
        elseif mat.captures[2] == "i"
            if mat.captures[1] == "/"
                write(output, "}")
            else
                write(output, "\\textit{")
            end
        elseif mat.captures[2] == "b"
            if mat.captures[1] == "/"
                write(output, "}")
            else
                write(output, "\\textbf{")
            end
        end
        lastpos = mat.offset + length(mat.match)
    end
    write(output, "\\text{")
    write(output, input[lastpos:end])
    write(output, "}")

    return String(take!(output))
end

function escape_tex_chars(text::AbstractString)
	# see http://www.cespedes.org/blog/85/how-to-escape-latex-special-characters
	escaped_str = text
	escaped_str = replace(escaped_str, "\\"=>"\\textbackslash{}")
	escaped_str = replace(escaped_str, "#"=>"\\#")
	escaped_str = replace(escaped_str, "\$"=>"\\\$")
	escaped_str = replace(escaped_str, "%"=>"\\%")
	escaped_str = replace(escaped_str, "&"=>"\\&")
	escaped_str = replace(escaped_str, "_"=>"\\_")
	escaped_str = replace(escaped_str, "{"=>"\\{")
	escaped_str = replace(escaped_str, "}"=>"\\}")
	escaped_str = replace(escaped_str, "^"=>"\\textasciicircum{}")
	escaped_str = replace(escaped_str, "~"=>"\\textasciitilde{}")
	escaped_str = replace(escaped_str, "\\textbackslash\\{\\}"=>"\\textbackslash{}")
end

function draw(img::PGF, prim::ArcPrimitive, idx::Int)

    angle2 = prim.angle2 + (prim.angle2<prim.angle1)*2π
    modifiers, props = get_vector_properties(img, idx)
    img.visible || return
    write(img.buf, join(modifiers))
    @printf(img.buf, "\\path [%s] (%s,%s) ",
        join(props, ","),
        svg_fmt_float(prim.center[1].value),
        svg_fmt_float(prim.center[2].value) )
    prim.sector && write(img.buf, "-- ")
    @printf(img.buf, "+(%s:%s) ",
        svg_fmt_float(rad2deg(prim.angle1)),
        svg_fmt_float(prim.radius.value) )
    @printf(img.buf,  "arc [radius=%s, start angle=%s, end angle=%s]",
        svg_fmt_float(prim.radius.value),
        svg_fmt_float(rad2deg(prim.angle1)),
        svg_fmt_float(rad2deg(angle2)) )
    prim.sector && write(img.buf, " -- cycle")
    write(img.buf, ";\n")
end


text:
# 2PL
"""
    Parameters2PL <: AbstractParametersBinary

Contains information of a set of item parameters (values, bounds, priors, posteiors, chains, expected Fisher information, calibrated) under the 2-parameter logistic model.
"""
mutable struct Parameters2PL <: AbstractParametersBinary
    a::Float64
    bounds_a::Vector{Float64}
    b::Float64
    bounds_b::Vector{Float64}
    prior::Distributions.MultivariateDistribution
    posterior::Distributions.MultivariateDistribution
    chain::Vector{Vector{Float64}}
    expected_information::Matrix{Float64}
    calibrated::Bool

    Parameters2PL(
        a::Float64,
        bounds_a::Vector{Float64},
        b::Float64,
        bounds_b::Vector{Float64},
        prior::Distributions.MultivariateDistribution,
        posterior::Distributions.MultivariateDistribution,
        chain::Vector{Vector{Float64}},
        expected_information::Matrix{Float64},
        calibrated::Bool,
    ) = new(
        a,
        bounds_a,
        b,
        bounds_b,
        prior,
        posterior,
        chain,
        expected_information,
        calibrated,
    )

    # Random Initializers
    function Parameters2PL(
        bivariate_dist::Distributions.MultivariateDistribution,
        bounds_a::Vector{Float64},
        bounds_b::Vector{Float64},
    )
        pars = truncate_rand(bivariate_dist, [bounds_a, bounds_b])
        return Parameters2PL(
            pars[1][1],
            bounds_a,
            pars[1][2],
            bounds_b,
            bivariate_dist,
            bivariate_dist,
            [pars[1]],
            [1.0 0.0; 0.0 1.0],
            false,
        )
    end

    function Parameters2PL()
        return Parameters2PL(
            Distributions.Product([
                Distributions.LogNormal(0, 0.25),
                Distributions.Normal(0, 1),
            ]),
            [1e-5, 5.0],
            [-6.0, 6.0],
        )
    end

    function Parameters2PL(bounds_a::Vector{Float64}, bounds_b::Vector{Float64})
        par = Parameters2PL()
        par.bounds_a = bounds_a
        par.bounds_b = bounds_b
        return par::Parameters2PL
    end
end

"""
    _empty_chain!(parameters::Parameters2PL)
"""
function _empty_chain!(parameters::Parameters2PL)
    parameters.chain = Vector{Vector{Float64}}(undef, 0)
end

"""
    _set_val!(parameters::Parameters2PL, vals::Vector{Float64})
"""
function _set_val!(parameters::Parameters2PL, vals::Vector{Float64})
    parameters.a = vals[1]
    parameters.b = vals[2]
    return nothing
end

"""
    _set_val_from_chain!(parameters::Parameters2PL)
"""
function _set_val_from_chain!(parameters::Parameters2PL)
    parameters.a = parameters.chain[end][1]
    parameters.b = parameters.chain[end][2]
    return nothing
end

"""
    _update_estimate!(parameters::Parameters2PL; sampling = true)
"""
function _update_estimate!(parameters::Parameters2PL; sampling = true)
    chain_size = size(parameters.chain, 1)
    if sampling
        chain_matrix = hcat(parameters.chain[(chain_size-min(999, chain_size - 1)):end]...)
        vals = [sum(i) / min(1000, chain_size) for i in eachrow(chain_matrix)]
    else
        chain_matrix = hcat(parameters.chain...)
        vals = [sum(i) / chain_size for i in eachrow(chain_matrix)]
    end
    parameters.a = clamp(vals[1], parameters.bounds_a[1], parameters.bounds_a[2])
    parameters.b = clamp(vals[2], parameters.bounds_b[1], parameters.bounds_b[2])
    return nothing
end

"""
    _get_parameters_vals(parameters::Parameters2PL)
"""
function _get_parameters_vals(parameters::Parameters2PL)
    return [parameters.a, parameters.b]
end

ext: jl


text:
export SimpleRandomWalkEnv, CartPoleEnv

include("simple_random_walk_env.jl")
include("cart_pole_env.jl")


text:
"""
$(TYPEDEF)

Hypograph of real symmetric or complex Hermitian root-determinant cone of
dimension `dim` in svec format.

    $(FUNCTIONNAME){T, R}(dim::Int, use_dual::Bool = false)
"""
mutable struct HypoRootdetTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d::Int
    is_complex::Bool
    rt2::T

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    di::T
    ϕ::T
    ζ::T
    ϕζidi::T
    mat::Matrix{R}
    fact_W::Cholesky{R}
    Wi::Matrix{R}
    Wi_vec::Vector{T}
    tempw::Vector{T}
    mat2::Matrix{R}
    mat3::Matrix{R}
    mat4::Matrix{R}

    function HypoRootdetTri{T, R}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert dim >= 2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.rt2 = sqrt(T(2))
        cone.is_complex = (R <: Complex)
        cone.d = svec_side(R, dim - 1)
        return cone
    end
end

function setup_extra_data!(
    cone::HypoRootdetTri{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    dim = cone.dim
    d = cone.d
    cone.di = inv(T(d))
    cone.mat = zeros(R, d, d)
    cone.Wi = zeros(R, d, d)
    cone.Wi_vec = zeros(T, dim - 1)
    cone.tempw = zeros(T, dim - 1)
    cone.mat2 = zeros(R, d, d)
    cone.mat3 = zeros(R, d, d)
    cone.mat4 = zeros(R, d, d)
    return cone
end

get_nu(cone::HypoRootdetTri) = 1 + cone.d

function set_initial_point!(
    arr::AbstractVector{T},
    cone::HypoRootdetTri{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    arr .= 0
    d = cone.d
    c1 = sqrt(T(5d^2 + 2d + 1))
    c2 = arr[1] = -sqrt((3d + 1 - c1) / T(2d + 2))
    c3 = -c2 * (d + 1 + c1) / 2d
    incr = (cone.is_complex ? 2 : 1)
    k = 2
    @inbounds for i in 1:d
        arr[k] = c3
        k += incr * i + 1
    end
    return arr
end

function update_feas(cone::HypoRootdetTri{T}) where {T <: Real}
    @assert !cone.feas_updated

    @views svec_to_smat!(cone.mat, cone.point[2:end], cone.rt2)
    fact = cone.fact_W = cholesky!(Hermitian(cone.mat, :U), check = false)
    if isposdef(fact)
        cone.ϕ = exp(logdet(fact) / cone.d)
        cone.ζ = cone.ϕ - cone.point[1]
        cone.is_feas = (cone.ζ > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::HypoRootdetTri{T}) where {T <: Real}
    u = cone.dual_point[1]

    if u < -eps(T)
        @views svec_to_smat!(cone.mat2, cone.dual_point[2:end], cone.rt2)
        fact = cholesky!(Hermitian(cone.mat2, :U), check = false)
        if isposdef(fact)
            return (logdet(fact) - cone.d * log(-u / cone.d) > eps(T))
        end
    end

    return false
end

function update_grad(cone::HypoRootdetTri)
    @assert cone.is_feas
    g = cone.grad
    ζ = cone.ζ
    cone.ϕζidi = cone.ϕ / ζ * cone.di
    ϕζidi1 = -cone.ϕζidi - 1

    g[1] = inv(ζ)
    inv_fact!(cone.Wi, cone.fact_W)
    smat_to_svec!(cone.Wi_vec, cone.Wi, cone.rt2)
    @. g[2:end] = ϕζidi1 * cone.Wi_vec

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::HypoRootdetTri)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    Wi_vec = cone.Wi_vec
    ζ = cone.ζ
    ϕζidi = cone.ϕζidi
    c1 = -ϕζidi / ζ
    ϕζidi1 = ϕζidi + 1
    c2 = ϕζidi * (ϕζidi - cone.di)

    H[1, 1] = ζ^-2
    @. H[1, 2:end] = c1 * Wi_vec

    copytri!(cone.Wi, 'U', true)
    @views symm_kron!(H[2:end, 2:end], cone.Wi, cone.rt2)

    @inbounds for j in eachindex(Wi_vec)
        j1 = 1 + j
        Wi_vecj = c2 * Wi_vec[j]
        for i in 1:j
            H[1 + i, j1] = ϕζidi1 * H[1 + i, j1] + Wi_vec[i] * Wi_vecj
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::HypoRootdetTri{T},
    ) where {T <: Real}
    @assert cone.grad_updated
    di = cone.di
    ζ = cone.ζ
    ϕζidi = cone.ϕζidi
    w_aux = cone.mat2
    FU = cone.fact_W.U
    ϕζidi1 = ϕζidi + 1

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        @views svec_to_smat!(w_aux, arr[2:end, j], cone.rt2)
        copytri!(w_aux, 'U', true)
        rdiv!(w_aux, FU)
        ldiv!(FU', w_aux)

        c0 = ϕζidi * tr(Hermitian(w_aux, :U))
        c1 = c0 - p / ζ
        c2 = ϕζidi * c1 - di * c0

        lmul!(ϕζidi1, w_aux)
        for i in diagind(w_aux)
            w_aux[i] += c2
        end
        rdiv!(w_aux, FU')
        ldiv!(FU, w_aux)

        prod[1, j] = c1 / -ζ
        @views smat_to_svec!(prod[2:end, j], w_aux, cone.rt2)
    end

    return prod
end

function update_inv_hess(cone::HypoRootdetTri)
    @assert cone.grad_updated
    @assert !cone.inv_hess_updated
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    @views w = cone.point[2:end]
    svec_to_smat!(cone.mat2, w, cone.rt2)
    W = Hermitian(cone.mat2, :U)
    Hi = cone.inv_hess.data
    ζ = cone.ζ
    ϕ = cone.ϕ
    di = cone.di
    ϕdi = ϕ * di
    c2 = inv(cone.ϕζidi + 1)
    c3 = ϕdi * c2 / ζ * di

    Hi[1, 1] = abs2(ζ) + ϕdi * ϕ
    @. Hi[1, 2:end] = ϕdi * w

    @views symm_kron!(Hi[2:end, 2:end], W, cone.rt2)

    @inbounds for j in eachindex(w)
        j1 = 1 + j
        scwj = c3 * w[j]
        for i in 1:j
            Hi[1 + i, j1] = c2 * Hi[1 + i, j1] + w[i] * scwj
        end
    end

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::HypoRootdetTri{T},
    ) where {T <: Real}
    @assert cone.grad_updated
    @views w = cone.point[2:end]
    svec_to_smat!(cone.mat4, w, cone.rt2)
    W = Hermitian(cone.mat4, :U)
    ζ = cone.ζ
    ϕ = cone.ϕ
    di = cone.di
    ϕdi = ϕ * di
    c2 = inv(cone.ϕζidi + 1)
    c3 = c2 / ζ * di
    c4 = abs2(ζ) + ϕdi * ϕ
    w_aux = cone.mat2
    w_aux2 = cone.mat3

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        @views r = arr[2:end, j]
        svec_to_smat!(w_aux, r, cone.rt2)
        copytri!(w_aux, 'U', true)

        c5 = dot(w, r)
        c6 = ϕdi * (c3 * c5 + p)
        prod[1, j] = ϕdi * c5 + c4 * p

        mul!(w_aux2, w_aux, W)
        mul!(w_aux, W, w_aux2)
        @views prod_w = prod[2:end, j]
        smat_to_svec!(prod_w, w_aux, cone.rt2)
        axpby!(c6, w, c2, prod_w)
    end

    return prod
end

function dder3(cone::HypoRootdetTri{T}, dir::AbstractVector{T}) where {T <: Real}
    @assert cone.grad_updated
    dder3 = cone.dder3
    p = dir[1]
    @views r = dir[2:end]
    ζ = cone.ζ
    ϕ = cone.ϕ
    di = cone.di
    ϕζidi = cone.ϕζidi
    FU = cone.fact_W.U
    rwi = cone.mat2
    w_aux = cone.mat3
    w_aux2 = cone.mat4

    svec_to_smat!(rwi, r, cone.rt2)
    copytri!(rwi, 'U', true)
    rdiv!(rwi, FU)
    ldiv!(FU', rwi)
    c0 = tr(Hermitian(rwi, :U)) * di
    c6 = sum(abs2, rwi) * di
    ζiχ = (p - ϕ * c0) / ζ
    c1 = ζiχ^2 + ϕ / ζ * (c6 - abs2(c0)) / 2
    c7 = ϕζidi * (c1 - c6 / 2 + c0 * (ζiχ + c0 / 2))
    c8 = -ϕζidi * (ζiχ + c0)
    c9 = ϕζidi + 1

    dder3[1] = c1 / -ζ

    copyto!(w_aux2, I)
    axpby!(c9, rwi, c8, w_aux2)
    mul!(w_aux, Hermitian(rwi, :U), w_aux2)
    @inbounds for i in diagind(w_aux)
        w_aux[i] += c7
    end
    rdiv!(w_aux, FU')
    ldiv!(FU, w_aux)
    @views smat_to_svec!(dder3[2:end], w_aux, cone.rt2)

    return dder3
end


text:
# -------------------------------------------------------
#            Check input arguments
# -------------------------------------------------------
length(ARGS) >= 1 || error("input arg: prefix.meas.xml [or prefix.in.xml]")

# support for symbol link to app/measure.jl in root of repo
__FILE__ = joinpath(dirname(@__FILE__), splitpath(dirname(@__FILE__))[end] == "app" ? "" : "app")

# -------------------------------------------------------
#            Start + Check git branch
# -------------------------------------------------------
using Dates, Printf, Pkg
const start_time = now()
println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))
println("Hostname: ", gethostname())

function is_dqmc_env_activated()
    project_file = Base.active_project()
    project = Pkg.Types.read_project(project_file)
    return !isnothing(project.name) && lowercase(project.name) == "qmc"
end

if !is_dqmc_env_activated()
  println("Activating DQMC environment.")
  haskey(ENV, "JULIA_DQMC") || error("DQMC environment not loaded and JULIA_DQMC env variable not set!")
  Pkg.activate(ENV["JULIA_DQMC"])
  is_dqmc_env_activated() || error("Invalid JULIA_DQMC env variable value.")
end

using Git
branch = Git.branch(dir=dirname(@__FILE__)).string[1:end-1]
if branch != "master"
  println("!!!Not on branch master but \"$(branch)\"!!!")
  flush(stdout)
end

# -------------------------------------------------------
#                   Load Framework
# -------------------------------------------------------
include(joinpath(__FILE__, "../src/measure_framework.jl"))

if HDF5.libversion.minor |> Int < 10
  @warn "You are using a HDF5 version < 1.10 in which space isn't freed automatically.
          Files might become very large, in particular for small values of SAVE_AFTER!"
end

# #######################################################
#                    Start script
# #######################################################
# -------------------------------------------------------
#           Parse ARGS and maybe .meas.xml
# -------------------------------------------------------

if length(ARGS) == 1
  arg = ARGS[1]

  isfile(arg) || error("Input file $arg not found.")
  if endswith(arg, ".meas.xml")
    mp = measxml2MeasParams(arg)
  elseif endswith(arg, ".in.xml")
    measxml = arg[1:end-7]*".meas.xml"
    isfile(measxml) || error("File $(measxml) not found. Input arg was: $(arg).")
    mp = measxml2MeasParams(measxml)
  else
    error("No .meas.xml or .in.xml file found. Input arg was: $(arg)")
  end
elseif length(ARGS) == 2
  # Backward compatibility
  # ARGS = ["sdwO3_L_4_B_2_dt_0.1_2", 1]
  prefix = convert(String, ARGS[1])
  idx = parse(Int, ARGS[2])           # SLURM_ARRAY_TASK_ID
  println("Prefix is ", prefix, " and idx is ", idx)
  measxml = prefix * ".task" * string(idx) * ".meas.xml"
  println("Input xml constructed as: $(measxml)")
  isfile(measxml) || error("File $(measxml) not found. Input args were: $(ARGS).")
  mp = measxml2MeasParams(measxml)
  println("Input xml loaded successfully!!!!!!")
else
  error("Incompatible input arguments. Input args were: $(ARGS)")
end

# set num threads
try
  BLAS.set_num_threads(mp.num_threads)
  FFTW.set_num_threads(mp.num_threads)
  ENV["OMP_NUM_THREADS"] = mp.num_threads
  ENV["MKL_NUM_THREADS"] = mp.num_threads
  ENV["JULIA_NUM_THREADS"] = 1
catch
end

# -------------------------------------------------------
#                     Let's go
# -------------------------------------------------------
measure(mp)

# -------------------------------------------------------
#                       Exit
# -------------------------------------------------------
end_time = now()
println("\nEnded: ", Dates.format(end_time, "d.u yyyy HH:MM"))
@printf("Duration: %.2f minutes", (end_time - start_time).value / 1000 / 60)

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
include("qPlots.jl")
include("../../conf.jl")

function renderWidget(widget::QWidget, today::DateTime, saveTo::String) #::Array{Tuple{String, Any, String, String},1}

    logger(widget, "Unknwon widget type $(typeof(widget))", true)

    throw("Unknwon widget type $(typeof(widget))")
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
#! /usr/bin/env julia

function check_spdx(dir, spdx_header="// SPDX-License-Identifier: ISC")
  for file in readdir(dir, join=true)
    if isfile(file) && endswith(file, ".v")
      content = read(file, String)
      if !contains(content, spdx_header)
        write(file, spdx_header*"\n"*content)
        println("added SPDX header to $file")
      end
    elseif isdir(file)
      check_spdx(file)
    end
  end
end

check_spdx("src")

max_stars_repo_name: opme/rapcores

max_issues_repo_name: opme/rapcores

max_forks_repo_name: opme/rapcores


text:
@testset "Rank computations using descent" begin
  E1 = EllipticCurve([1, 1, 1, -352, -2689])
  E2 = EllipticCurve([0, 1, 0, -41, -116])
  E3 = EllipticCurve([1, 0, 1, -231, -442])
  E4 = EllipticCurve([0, -1, 0, -289, 289])
  @testset "2-torsion descent" begin

    lower, upper, III_phi,_III_phi = @inferred rank_2_torsion(E1)
    @test (lower <= 1 <= upper)
    @test (III_phi == 4 &&  _III_phi == 1)
    lower, upper, III_phi,_III_phi = @inferred rank_2_torsion(E2)
    @test (lower <= 0 <= upper)
    @test (III_phi == 1 && _III_phi == 1)
    lower, upper, III_phi,_III_phi = @inferred rank_2_torsion(E3)
    @test (lower <= 3 <= upper)
    @test (III_phi == 2 && _III_phi == 1)
    lower, upper, III_phi,_III_phi = @inferred rank_2_torsion(E4)
    @test (lower <= 2 <= upper)
    @test III_phi == 2 && _III_phi == 1

  end
  @testset "local solubility" begin

    @test quartic_local_solubility(3, 2, 3, 4, 5) == false
    @test quartic_local_solubility(3437843, 24789579, 3122111, -424, 57) == true
    @test quartic_local_solubility(-888, 5767, 2110, -424, 5128673) == false

  end

end


text:
# ---
# title: 66. Plus One
# id: problem66
# author: zhwang
# date: 2022-03-17
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/plus-one/description/>
# hidden: true
# ---
#
# Given a **non-empty** array of decimal digits representing a non-negative
# integer, increment one to the integer.
#
# The digits are stored such that the most significant digit is at the head of
# the list, and each element in the array contains a single digit.
#
# You may assume the integer does not contain any leading zero, except the
# number 0 itself.
#
#
#
# **Example 1:**
#
#
#
#     Input: digits = [1,2,3]
#     Output: [1,2,4]
#     Explanation: The array represents the integer 123.
#
#
# **Example 2:**
#
#
#
#     Input: digits = [4,3,2,1]
#     Output: [4,3,2,2]
#     Explanation: The array represents the integer 4321.
#
#
# **Example 3:**
#
#
#
#     Input: digits = [0]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= digits.length <= 100`
#   * `0 <= digits[i] <= 9`
#
#
## @lc code=start
using LeetCode

## use build-in function
plus_one(nums::Vector{Int}) = reverse!(digits(foldl((i, j) -> 10 * i + j, nums) + 1))
## method 2
function digits_plus_one(nums::Vector{Int})::Vector{Int}
    pos = findlast(!=(9), nums)
    isnothing(pos) && return append!([1], zeros(Int, length(nums)))
    return cat(
        @view(nums[1:(pos - 1)]), nums[pos] + 1, zeros(Int, length(nums) - pos); dims=1
    )
end

## @lc code=end


text:
export
    validated_setdiff!

validated_setdiff!(set, subset) = begin
    length(subset) == length(unique(subset)) || return false
    subset ⊆ set || return false

    setdiff!(set, subset)
    true
end


text:
module Cat

using MacroTools

export Arrow, Identity, @category, compose, lift, composition_obj, @inverse, parent_category, @morphism, @stateful
include("arrow.jl")
include("representable.jl")

end # end module


text:
# define the generative model for the urban pomdp environment

### REWARD MODEL ##################################################################################

function POMDPs.reward(pomdp::UrbanPOMDP, s::UrbanState, a::UrbanAction, sp::UrbanState)
    r = 0.
    ego = sp[findfirst(EGO_ID, sp)]
    if is_crash(sp)
        r += pomdp.collision_cost
    end
    if ego.state.posF.s >= get_end(pomdp.env.roadway[pomdp.ego_goal]) &&
       get_lane(pomdp.env.roadway, ego).tag == pomdp.ego_goal
        r += pomdp.goal_reward
    elseif a.acc > 0.
        r += pomdp.action_cost
    else
        r += pomdp.action_cost
    end
    return r
end

### TERMINAL STATES ###############################################################################

function POMDPs.isterminal(pomdp::UrbanPOMDP, s::UrbanState)
    ego = s[findfirst(EGO_ID, s)]
    return is_crash(s) || (ego.state.posF.s >= get_end(pomdp.env.roadway[pomdp.ego_goal]) &&
                           get_lane(pomdp.env.roadway, ego).tag == pomdp.ego_goal)
end

### TRANSITION MODEL ##############################################################################

function POMDPs.gen(::DDNNode{:sp}, pomdp::UrbanPOMDP, s::UrbanState, a::UrbanAction, rng::AbstractRNG)
    # clean dict
    ids = [veh.id for veh in s]
    key_ = deepcopy(keys(pomdp.models))
    for k in key_
        if !(k ∈ ids) && k!=1
            delete!(pomdp.models, k)
        end
    end
    actions = Array{Any}(undef, length(s))
    pomdp.models[1].a = a
    ego_ind = findfirst(EGO_ID, s)
    is_ego_here =  ego_ind == nothing ? 0 : 1
    sp = deepcopy(s) #XXX bad
    if rand(rng) < pomdp.car_birth && n_cars(s) < pomdp.max_cars + is_ego_here
        new_car = initial_car(pomdp, sp, rng)
        if can_push(pomdp.env, sp, new_car)
            lane = get_lane(pomdp.env.roadway, new_car)
            route = find_route(pomdp.env, random_route(rng, pomdp.env.roadway, lane))
            pomdp.models[new_car.id] = pomdp.car_models[SVector(route[1], route[end])]
            reset_hidden_state!(pomdp.models[new_car.id])
            push!(sp, new_car)
        end
    end
    if rand(rng) < pomdp.ped_birth && n_pedestrians(s) < pomdp.max_peds
        # println("Spawning new pedestrians")
        new_ped = initial_pedestrian(pomdp, sp, rng)
        # pomdp.models[new_ped.id] = ConstantPedestrian(dt = pomdp.ΔT)#TODO parameterized
        new_ped_cw = get_lane(pomdp.env.roadway, new_ped)
        new_ped_conflict_lanes = get_conflict_lanes(new_ped_cw, pomdp.env.roadway)
        pomdp.models[new_ped.id] = IntelligentPedestrian(dt = pomdp.ΔT, crosswalk=new_ped_cw, conflict_lanes=new_ped_conflict_lanes)
        reset_hidden_state!(pomdp.models[new_ped.id])
        push!(sp, new_ped)
    end
    actions = Array{Any}(undef, length(sp))
    get_actions!(actions, sp, pomdp.env.roadway, pomdp.models)
    tick!(sp, pomdp.env.roadway, actions, pomdp.ΔT)
    clean_scene!(pomdp.env, sp)
    return sp
end

function POMDPModelTools.gen(::DDNOut{(:sp, :o, :r, :info)}, p::UrbanPOMDP, s, a, rng::AbstractRNG)
    return gen(DDNOut(:sp, :o, :r), p, s, a, rng)..., deepcopy(p.models)
end

#TODO move to helpers

function AutomotiveDrivingModels.propagate(veh::Vehicle, a::UrbanAction, roadway::Roadway, Δt::Float64)
    acc = LonAccelDirection(a.acc, 4) #TODO Parameterize
    return propagate(veh, acc, roadway, Δt)
end

### INITIAL STATES ################################################################################

function POMDPs.initialstate(pomdp::UrbanPOMDP, rng::AbstractRNG, burn_in::Int64=1)
    scene = initial_scene(pomdp, rng, true)
    for t = 1:burn_in
        scene = gen(DDNOut(:sp), pomdp, scene, UrbanAction(0.), rng)
    end
    clean_scene!(pomdp.env, scene)
    # push! ego after traffic is steady
    ego = initial_ego(pomdp, rng)
    push!(scene, ego)
    clean_scene!(pomdp.env, scene)
    return scene
end

function initial_scene(pomdp::UrbanPOMDP, rng::AbstractRNG, no_ego::Bool=false)
    empty_obstacles!(pomdp.env)
    if pomdp.max_obstacles > 0
        for i=1:pomdp.max_obstacles
            sample_obstacle!(pomdp.env, pomdp.obs_dist, rng)
        end
    end
    scene = Scene()
    if !no_ego
        ego = initial_ego(pomdp, rng)
        push!(scene, ego)
    end

    # add cars
    for i=1:pomdp.max_cars
        if rand(rng) < pomdp.car_birth && n_cars(scene) < pomdp.max_cars + 1-Int(no_ego)
            new_car = initial_car(pomdp, scene, rng, true)
            if can_push(pomdp.env, scene, new_car)
                push!(scene, new_car)
                lane = get_lane(pomdp.env.roadway, new_car)
                route = find_route(pomdp.env, random_route(rng, pomdp.env.roadway, lane))
                pomdp.models[new_car.id] = pomdp.car_models[SVector(route[1], route[end])]
                reset_hidden_state!(pomdp.models[new_car.id])
            end
            clean_scene!(pomdp.env, scene)
        end
    end

    # add pedestrians
    for i=1:pomdp.max_peds
        if rand(rng) < pomdp.ped_birth && n_pedestrians(scene) < pomdp.max_peds # pedestrian appear
            new_ped = initial_pedestrian(pomdp, scene, rng, true)
            # pomdp.models[new_ped.id] = ConstantPedestrian(dt = pomdp.ΔT)#TODO parameterized
            new_ped_cw = get_lane(pomdp.env.roadway, new_ped)
            new_ped_conflict_lanes = get_conflict_lanes(new_ped_cw, pomdp.env.roadway)
            pomdp.models[new_ped.id] = IntelligentPedestrian(dt = pomdp.ΔT, crosswalk=new_ped_cw, conflict_lanes=new_ped_conflict_lanes)
            reset_hidden_state!(pomdp.models[new_ped.id])
            push!(scene, new_ped)
        end
    end
    clean_scene!(pomdp.env, scene)
    return scene
end

function initial_car(pomdp::UrbanPOMDP, scene::Scene, rng::AbstractRNG, first_scene::Bool = false)
    env = pomdp.env

    #velocity
    v0 = rand(rng, 4.0:1.0:7.0) #XXX parameterize!
    s0 = 0.
    if first_scene
        lanes = get_lanes(pomdp.env.roadway)
        lanes = delete!(Set(lanes),pomdp.env.roadway[LaneTag(6,1)])
        lane = rand(rng, lanes) # could be precomputed
        s0 = rand(rng, 0.:0.5:get_end(lane))
    else
        lanes = get_start_lanes(pomdp.env.roadway)
        # lanes = delete!(Set(lanes),pomdp.env.roadway[LaneTag(6,1)])
        lane = rand(rng, lanes) # could be precomputed
    end
    initial_posF = Frenet(lane, s0)
    initialstate = VehicleState(initial_posF, env.roadway, v0)

    id = next_car_id(pomdp, scene)

    return Vehicle(initialstate, pomdp.car_type, id)
end

function initial_ego(pomdp::UrbanPOMDP, rng::AbstractRNG)
    posF = Frenet(get_ego_route(pomdp)[1], pomdp.ego_start)
    state = VehicleState(posF, pomdp.env.roadway, 0.)
    return Vehicle(state, pomdp.ego_type, EGO_ID)
end

"""
    initial_pedestrian(pomdp::UrbanPOMDP, scene::Scene, rng::AbstractRNG)
Create a new pedestrian entity at its initial state
"""
function initial_pedestrian(pomdp::UrbanPOMDP, scene::Scene, rng::AbstractRNG, first_scene::Bool = false)
    env = pomdp.env
    cw_ind = rand(rng, 1:length(env.params.crosswalk_pos))
    crosswalk_pos = env.params.crosswalk_pos[cw_ind]

    # position along the crosswalk
    t0 = rand(rng, Distributions.Uniform(-env.params.crosswalk_width[cw_ind]/2 + 1., env.params.crosswalk_width[cw_ind]/2 - 1.))
    s0 = rand(rng, [0., get_end(env.crosswalks[cw_ind])])
    ϕ0 = float(π)
    if s0 == 0.
        ϕ0 = 0.
    end
    if first_scene
        s0 = rand(rng, Distributions.Uniform(0., get_end(env.crosswalks[cw_ind])))
    end

    #velocity
    v0 = rand(rng, Distributions.Uniform(0., env.params.ped_max_speed))
    posF = Frenet(env.crosswalks[cw_ind], s0, t0, ϕ0)

    ped_initialstate = VehicleState(posF, env.ped_roadway, v0);

    id = next_ped_id(pomdp, scene)

    return Vehicle(ped_initialstate, pomdp.ped_type, id)
end

function POMDPs.initialstate_distribution(pomdp::UrbanPOMDP)
    return GenerativeDist(pomdp)
end

### Observations ##################################################################################

function POMDPs.gen(::DDNNode{:o}, pomdp::UrbanPOMDP, s::UrbanState, a::UrbanAction, sp::UrbanState, rng::AbstractRNG)
    if pomdp.lidar
        return measure_lidar(pomdp, s, a, sp, rng)
    else
        return measure_gaussian(pomdp, s, a, sp, rng)
    end
end

function measure_gaussian(pomdp::UrbanPOMDP, s::UrbanState, a::UrbanAction, sp::UrbanState, rng::AbstractRNG)
    ego = get_ego(sp)
    obs = measure(pomdp.sensor, ego, sp, pomdp.env.roadway, pomdp.env.obstacles)
    o = obs_to_array(pomdp, ego, obs)
    return o
end

function obs_to_array(pomdp::UrbanPOMDP, ego_veh::Vehicle, obs::Vector{Vehicle})
    sorted_vehicles = sort!(obs, by=x->x.id)
    n_features = 4
    n_obstacles = pomdp.max_obstacles
    o = zeros(n_features*(pomdp.max_cars + pomdp.max_peds + n_obstacles + 1))
    ego = ego_veh.state
    o[1] = ego.posG.x
    o[2] = ego.posG.y
    o[3] = ego.posG.θ
    o[4] = ego.v
    pos_off = get_off_the_grid(pomdp)
    for i=2:pomdp.max_cars + pomdp.max_peds + n_obstacles + 1
        o[n_features*i - 3] = pos_off.posG.x - ego.posG.x
        o[n_features*i - 2] = pos_off.posG.y - ego.posG.y
        o[n_features*i - 1] = pos_off.posG.θ
        o[n_features*i] = 0.
    end
    for veh in sorted_vehicles
        if veh.id == EGO_ID
            continue
        end
        if veh.def.class == AgentClass.CAR
            @assert veh.id <= pomdp.max_cars+1 "$(veh.id)"
            o[n_features*veh.id - 3] = veh.state.posG.x - ego.posG.x
            o[n_features*veh.id - 2] = veh.state.posG.y - ego.posG.y
            o[n_features*veh.id - 1] = veh.state.posG.θ
            o[n_features*veh.id] = veh.state.v
        end
        if veh.def.class == AgentClass.PEDESTRIAN
            o[n_features*(veh.id - 100 + pomdp.max_cars + 1) - 3] = veh.state.posG.x - ego.posG.x
            o[n_features*(veh.id - 100 + pomdp.max_cars + 1) - 2] = veh.state.posG.y - ego.posG.y
            o[n_features*(veh.id - 100 + pomdp.max_cars + 1) - 1] = veh.state.posG.θ
            o[n_features*(veh.id - 100 + pomdp.max_cars + 1)] = veh.state.v
        end
    end
    for (j,obs) in enumerate(pomdp.env.obstacles)
        o[n_features*(j + pomdp.max_cars + pomdp.max_peds + 1) - 3] = get_width(obs)
        o[n_features*(j + pomdp.max_cars + pomdp.max_peds + 1) - 2] = get_height(obs)
        o[n_features*(j + pomdp.max_cars + pomdp.max_peds + 1) - 1] = get_center(obs).x - ego.posG.x
        o[n_features*(j + pomdp.max_cars + pomdp.max_peds + 1)] = get_center(obs).y - ego.posG.y
    end
    return rescale!(o, pomdp)
end

function POMDPs.convert_s(::Type{Vector{Float64}}, s::UrbanState, pomdp::UrbanPOMDP)
    ego = get_ego(s)
    obs = [veh for veh in s if veh.id != ego]
    svec = obs_to_array(pomdp, ego, obs)
    return svec
end

function n_dims(pomdp::UrbanPOMDP)
    n_features = 4
    n_obstacles = pomdp.max_obstacles
    return n_features*(pomdp.max_cars + pomdp.max_peds + n_obstacles + 1)
end

function POMDPs.initialobs(pomdp::UrbanPOMDP, s::UrbanState, rng::AbstractRNG)
    return gen(DDNNode(:o), pomdp, s, UrbanAction(0.), s, rng::AbstractRNG)
end

function rescale!(o::UrbanObs, pomdp::UrbanPOMDP)
    # rescale
    n_features = 4
    n_obstacles = pomdp.max_obstacles
    max_ego_dist = get_end(pomdp.env.roadway[pomdp.ego_goal])
    o[1] /= max_ego_dist
    o[2] /= max_ego_dist
    o[3] /= pi/2
    o[4] /= pomdp.env.params.speed_limit
    for i=2:pomdp.max_cars + pomdp.max_peds + 1
        o[n_features*i - 3] /= max_ego_dist
        o[n_features*i - 2] /= max_ego_dist
        o[n_features*i - 1] /= pi
        o[n_features*i] /= pomdp.env.params.speed_limit # XXX parameterized
    end
    for i=pomdp.max_cars + pomdp.max_peds + 2:pomdp.max_cars + pomdp.max_peds + 1 + n_obstacles
        o[n_features*i - 3] /= max_ego_dist
        o[n_features*i - 2] /= max_ego_dist
        o[n_features*i - 1] /= max_ego_dist
        o[n_features*i] /= max_ego_dist
    end

    return o
end

function unrescale!(o::UrbanObs, pomdp::UrbanPOMDP)
    # unrescale
    n_features = 4
    n_obstacles = pomdp.max_obstacles
    max_ego_dist = get_end(pomdp.env.roadway[pomdp.ego_goal])
    o[1] *= max_ego_dist
    o[2] *= max_ego_dist
    o[3] *= pi/2
    o[4] *= pomdp.env.params.speed_limit
    for i=2:pomdp.max_cars + pomdp.max_peds + 1
        o[n_features*i - 3] *= max_ego_dist
        o[n_features*i - 2] *= max_ego_dist
        o[n_features*i - 1] *= pi
        o[n_features*i] *= pomdp.env.params.speed_limit # XXX parameterized
    end
    for i=pomdp.max_cars + pomdp.max_peds + 2:pomdp.max_cars + pomdp.max_peds + 1 + n_obstacles
        o[n_features*i - 3] *= max_ego_dist
        o[n_features*i - 2] *= max_ego_dist
        o[n_features*i - 1] *= max_ego_dist
        o[n_features*i] *= max_ego_dist
    end

    return o
end

function to_global_coordinates!(o::UrbanObs, pomdp::UrbanPOMDP)
    n_features = 4
    n_obstacles = pomdp.max_obstacles
    ego = o[1:n_features]
    for i=2:pomdp.max_cars + pomdp.max_peds + 1
        o[n_features*i - 3] += ego[1]
        o[n_features*i - 2] += ego[2]
    end
    for i=pomdp.max_cars + pomdp.max_peds + 2:pomdp.max_cars + pomdp.max_peds + 1 + n_obstacles
        o[n_features*i - 1] += ego[1]
        o[n_features*i]+= ego[2]
    end
    return o
end

function get_normalized_absent_state(pomdp::UrbanPOMDP, ego::Vector{Float64})
    ego_x, ego_y, theta, v = ego
    pos_off =  get_off_the_grid(pomdp)
    max_ego_dist = get_end(pomdp.env.roadway[pomdp.ego_goal])
    return [pos_off.posG.x/max_ego_dist - ego_x,
                    pos_off.posG.y/max_ego_dist - ego_y,
                    pos_off.posG.θ/float(pi),
                    0. ]
end

function measure_lidar(pomdp::UrbanPOMDP, s::UrbanState, a::UrbanAction, sp::UrbanState, rng::AbstractRNG)
    observe!(pomdp.sensor, sp, pomdp.env.roadway, EGO_ID)
    return vcat(pomdp.sensor.ranges/pomdp.sensor.max_range, pomdp.sensor.range_rates/pomdp.env.params.speed_limit)
end

function POMDPs.convert_o(::Type{Vector{Float64}}, o::UrbanObs, pomdp::UrbanPOMDP)
    return o
end

function POMDPs.gen(::DDNOut{(:sp, :o, :r, :info)}, p::UrbanPOMDP, s::UrbanState, a::UrbanAction, rng::AbstractRNG)
    if p.lidar
        return gen(DDNOut(:sp,:o,:r), p, s, a, rng)..., p.sensor
    else
        return gen(DDNOut(:sp,:o,:r), p, s, a, rng)..., nothing
    end
end

# uncomment for Lidar
# uncomment for LIDAR observation

#
# function POMDPs.generate_o(pomdp::OCPOMDP, s::OCState, rng::AbstractRNG)
#     observe!(pomdp.sensor, s, pomdp.env.roadway, EGO_INDEX)
#     return pomdp.sensor
# end
#
# function POMDPs.convert_o(::Type{Vector{Float64}}, o::OCObs, pomdp::OCPOMDP)
#     return [o.ranges/pomdp.sensor.max_range, o.range_rates/pomdp.env.params.speed_limit]
# end

#### Helpers
# """
#     ego_state(pomdp:ICPOMDP, s::Float64, v::Float64)
# Helper that returns the ego car state given its position on the lane and the velocity
# """
# function ego_state(pomdp::UrbanPOMDP, s::Float64, v::Float64)
#     # ego_lane = pomdp.env.lane_map["ego_left"] #TODO parameterized
#     posF = Frenet(get_ego_route(pomdp)[1], s)
#     return VehicleState(posF, pomdp.env.roadway, v)
# end

"""
    car_state(pomdp::UrbanPOMDP, lane_id::String, s::Float64, v::Float64)
    car_state(pomdp::UrbanPOMDP, lane::Lane, s::Float64, v::Float64)
    car_state(pomdp::UrbanPOMDP, x::Float64, y::Float64, θ::Float64, v::Float64)
Helper that returns the state of a car given the lane and the position on the lane
"""
function car_state(pomdp::UrbanPOMDP, lane_id::String, s::Float64, v::Float64)
    lane = pomdp.env.lane_map[lane_id]
    posF = Frenet(lane, s)
    return VehicleState(posF, pomdp.env.roadway, v)
end
function car_state(pomdp::UrbanPOMDP, lane::Lane, s::Float64, v::Float64)
    posF = Frenet(lane, s)
    return VehicleState(posF, pomdp.env.roadway, v)
end
function car_state(pomdp::UrbanPOMDP, x::Float64, y::Float64, θ::Float64, v::Float64)
    posG = VecSE2(x, y, θ)
    return VehicleState(posG, pomdp.env.roadway, v)
end

"""
    get_off_the_grid(pomdp::UrbanPOMDP)
return the off the grid state
"""
function get_off_the_grid(pomdp::UrbanPOMDP)
    posG = pomdp.off_grid
    return VehicleState(posG, pomdp.env.roadway, 0.)
end

function off_the_grid(veh::VehicleState, pomdp::UrbanPOMDP)
    return veh.posG == pomdp.off_grid
end

"""
return the list of lanes the ego car should take
"""
function get_ego_route(pomdp::UrbanPOMDP)
    tags = [LaneTag(6,1), LaneTag(13, 1), LaneTag(2,pomdp.env.params.nlanes_main)]
    lanes = Array{Lane}(undef, length(tags))
    for (i,tag) in enumerate(tags)
        lanes[i] = pomdp.env.roadway[tag]
    end
    lanes
end

"""
create a unique ID for a new car
"""
function next_car_id(pomdp::UrbanPOMDP, scene::Scene, rng::AbstractRNG)
    max_id = pomdp.max_cars+1
    current_car_ids = []
    for veh in scene
        if veh.def.class == AgentClass.CAR
            push!(current_car_ids, veh.id)
        end
    end
    possible_ids = setdiff(2:max_id, current_car_ids)
    id = rand(rng, possible_ids)
    return id
end

function next_car_id(pomdp::UrbanPOMDP, scene::Scene)
    max_id = pomdp.max_cars+1
    current_car_ids = []
    for veh in scene
        if veh.def.class == AgentClass.CAR
            push!(current_car_ids, veh.id)
        end
    end
    possible_ids = setdiff(2:max_id, current_car_ids)
    return possible_ids[1]
end

"""
create a unique ID for a new pedestrian
"""
function next_ped_id(pomdp::UrbanPOMDP, scene::Scene, rng::AbstractRNG)
    max_id = pomdp.max_peds
    current_ped_ids = []
    for veh in scene
        if veh.def.class == AgentClass.PEDESTRIAN
            push!(current_ped_ids, veh.id)
        end
    end
    possible_ids = setdiff(101:100+max_id, current_ped_ids)
    id = rand(rng, possible_ids)
    return id
end
function next_ped_id(pomdp::UrbanPOMDP, scene::Scene)
    max_id = pomdp.max_peds
    current_ped_ids = []
    for veh in scene
        if veh.def.class == AgentClass.PEDESTRIAN
            push!(current_ped_ids, veh.id)
        end
    end
    possible_ids = setdiff(101:100+max_id, current_ped_ids)
    return possible_ids[1]
end

function n_cars(scene::Scene)
    n = 0
    for veh in scene
        if veh.def.class == AgentClass.CAR
            n += 1
        end
    end
    return n
end

function n_pedestrians(scene::Scene)
    n = 0
    for veh in scene
        if veh.def.class == AgentClass.PEDESTRIAN
            n += 1
        end
    end
    return n
end

"""
Create a scene that can be rendered from an observation
"""
function obs_to_scene(pomdp::UrbanPOMDP, obs::UrbanObs)
    o = deepcopy(obs)
    o = unrescale!(o, pomdp)
    o = to_global_coordinates!(o, pomdp)
    scene = Scene()
    # extract
    ego, car_map, ped_map, obs_map = split_o(o, pomdp, false)
    # println("cars observed : ", keys(car_map), "ped observed :", keys(ped_map))
    # project to roadway
    params = pomdp.env.params
    intersection_params = TInterParams(params.x_min, params.x_max, params.y_min, params.inter_x,
                                    params.inter_y, params.inter_width, params.inter_height,
                                    params. lane_width, params.nlanes_main, params.nlanes,
                                    params. stop_line, params.speed_limit, params.car_rate)
    car_roadway = gen_T_roadway(intersection_params)
    ego_state = VehicleState(VecSE2(ego[1], ego[2], ego[3]),  car_roadway, ego[4])
    ego = Vehicle(ego_state, pomdp.ego_type, EGO_ID)
    push!(scene, ego)
    for (str_id, vec_state) in car_map
        # make sure to not project it on a crosswalk
        car_state = VehicleState(VecSE2(vec_state[1], vec_state[2], vec_state[3]), car_roadway, vec_state[4])
        if !off_the_grid(car_state, pomdp)
            id = next_car_id(pomdp, scene)
            car = Vehicle(car_state, pomdp.car_type, id)
            push!(scene, car)
        end
    end
    for (str_id, vec_state) in ped_map
        ped_state = VehicleState(VecSE2(vec_state[1], vec_state[2], vec_state[3]), pomdp.env.ped_roadway, vec_state[4])
        # println("ped_state ", ped_state)
        if !off_the_grid(ped_state, pomdp)
            id = next_ped_id(pomdp, scene)
            ped = Vehicle(ped_state, pomdp.ped_type, id)
            push!(scene, ped)
        end
    end
    @assert issorted([veh.id for veh in scene])
    return scene
end

# given a big observation vector, split to an entity-wise representation
function split_o(obs::UrbanObs, pomdp::UrbanPOMDP, normalized=true)
    car_map, ped_map, obs_map = OrderedDict{Int64, Vector{Float64}}(), OrderedDict{Int64, Vector{Float64}}(), OrderedDict{Int64, Vector{Float64}}() #XXX Dictionary creation is sloooooow
    n_features = pomdp.n_features
    ego = obs[1:n_features]
    if normalized
        absent = normalized_off_the_grid_pos(pomdp, ego[1], ego[2])
    else
        pos_off = get_off_the_grid(pomdp)
        absent= [pos_off.posG.x, pos_off.posG.y, pos_off.posG.θ, pos_off.v]
    end
    # println("ego ", ego)
    # println("absent state :", absent)
#     println("ego idx ", 1, ":", n_features)
    for (j,i) in enumerate(1:pomdp.max_cars)
        if !isapprox(obs[i*n_features+1:(i+1)*n_features], absent, atol=1.0)
            car_map[j+1] = obs[i*n_features+1:(i+1)*n_features]
        end
#         println("car$j idx ", i*n_features+1:(i+1)*n_features)
    end
    for (j,i) in enumerate(pomdp.max_cars + 1:pomdp.max_cars + pomdp.max_peds)
        if !isapprox(obs[i*n_features+1:(i+1)*n_features], absent, atol=1.0)
            # println(norm(obs[i*n_features+1:(i+1)*n_features] - absent))
            ped_map[100+j] = obs[i*n_features+1:(i+1)*n_features]
        end
#         println("ped$j idx ", i*n_features+1:(i+1)*n_features)
    end
    for (j,i)  in enumerate(pomdp.max_cars + pomdp.max_peds + 1:pomdp.max_cars + pomdp.max_peds + pomdp.max_obstacles)
        obs_map[j] = obs[i*n_features+1:(i+1)*n_features]
#         println("obs$j idx ", i*n_features+1:(i+1)*n_features)
    end
    return ego, car_map, ped_map, obs_map
end

function normalized_off_the_grid_pos(pomdp::UrbanPOMDP,normalized_ego_x::Float64, normalized_ego_y::Float64)
    pos_off = get_off_the_grid(pomdp)
    max_ego_dist = get_end(pomdp.env.roadway[pomdp.ego_goal])
    return [pos_off.posG.x/max_ego_dist - normalized_ego_x, pos_off.posG.y/max_ego_dist - normalized_ego_y, pos_off.posG.θ, 0.]
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
module Peridynamics

include("discretization/Discretization.jl")

end # module


text:
module TestJulia07

using Pkg
using LibGit2

function pull_repo(pkg)
    path = normpath(Pkg.devdir(), pkg)
    repo = LibGit2.GitRepo(path)
    LibGit2.fetch(repo)
    ann = LibGit2.GitAnnotated(repo, "refs/remotes/origin/master")
    (ma, mp) = LibGit2.merge_analysis(repo, [ann])
    if LibGit2.isset(ma, Cint(LibGit2.Consts.MERGE_ANALYSIS_UP_TO_DATE))
        println("Already up to date.")
    else
        println()
        cmt = LibGit2.GitCommit(repo, LibGit2.GitHash(ann))
        sig = LibGit2.author(cmt)
        for (idx, line) in enumerate(split(LibGit2.message(cmt), '\n'))
            if 1 == idx
                print("  - ")
            else
                print("    ")
            end
            println(line)
        end
        println("    by ", sig.name)
        LibGit2.merge!(repo, fastforward=true)
    end
    LibGit2.close(repo)
end

# julia --color=yes -e 'using TestJulia07; TestJulia07.pull(split("Metalhead Flux"))'
function pull(pkgs::Vector)
    for pkg in pkgs
        printstyled(pkg, ' ', color=:yellow)
        pull_repo(pkg)
    end
end

end # module TestJulia07


text:
# export names
# do typealiases have to be exported?  I don't think so

global const comm_type = MPI.CComm

export PETSC_NULL, PETSC_IGNORE, PETSC_DECIDE, PETSC_DETERMINE, PETSCDEFAULT, PETSC_COMM_SELF

export PetscInt, PetscScalar, PetscBool, PetscErrorCode, PetscDataType, PetscReal

export INSERT_VALUES, ADD_VALUES, PETSC_COPY_VALUESA
export PETSC_TRUE, PETSC_FALSE
export NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY, NORM_MAX

export MAT_FLUSH_ASSEMBLY, MAT_FINAL_ASSEMBLY

export KSPRICHARDSON, KSPCHEBYSHEV, KSPCG, KSPGROPPCG, KSPIPECG, KSPCGNE, KSPNASH, KSPSTCG, KSPGLTR, KSPFCG, KSPGMRES, KSPFGMRES, KSPLGMRES, KSPDGMRES, KSPPGMRES, KSPTCQMR, SKPBCGS, KSPIBCGS, KSPFBCGS, KSPFBCGSR, KSPBCGSL, KSPCGS, KSPTFQMR, KSPCR, KSPPIPECR, KSPSLQR, SKPPREONLY, KSPQCQ, KSPBICG, KSPMINRES, KSPSYMMLQ, KSPLCD, KSPPYTHON, KSPCR

export KSP_NORM_DEFAULT, KSP_NORM_NONE, KSP_NORM_PRECONDITIONED, KSP_NORM_UNPRECONGDITIONED, KSP_NORM_NATURAL

export KSP_NORM_MAX

export PetscIntNullArray, PetscScalarNullArray

export VECSEQ
export VECMPI
export VECSTANDARD
export VECSHARED
export VECSEQCUSP
export VECMPICUSP
export VECCUSP
export VECSEQVIENNACL
export VECMPIVIENNACL
export VECVIENNACL
export VECNEST
export VECSEQPTHREAD
export VECMPIPTHREAD
export VECPTHREAD

# decide which version of Petsc to use
# if no PETSC_DIR or PETSC_ARCH defined, use the one (hopefully) built
# during package installation
if !haskey(ENV, "PETSC_DIR") && !haskey(ENV, "PETSC_ARCH")
  file_path = joinpath(dirname(@__FILE__), "..", "deps/petsc_evars")
  args = open(readdlm, file_path)
  ENV["PETSC_DIR"] = args[1]
  ENV["PETSC_ARCH"] = args[2]
end

global const PETSC_DIR = ENV["PETSC_DIR"]
global const PETSC_ARCH = ENV["PETSC_ARCH"]

function finalizeMPI()
  if MPI.Initialized()
    MPI.Finalize()
  end
end
if !MPI.Initialized()
  MPI.Init()
  atexit(finalizeMPI)
end

myrank = MPI.Comm_rank(MPI.COMM_WORLD)

#PETSC_DIR = readall(`echo $PETSC_DIR`)
#PETSC_ARCH = readall(`echo $PETSC_ARCH`)

#PETSC_DIR = "/home/jared/build/petsc-3.6.0"
#PETSC_ARCH = "arch-linux2-c-debug"

#PETSC_DIR = getenv("PETSC_DIR");
#PETSC_ARCH = getenv("PETSC_ARCH");
#=
if (length(PETSC_DIR) == 0)
  disp("Must have environmental variable PETSC_DIR set")
end
if (length(PETSC_ARCH) == 0)
  disp("Must have environmental variable PETSC_ARCH set")
end
=#

global const libpetsclocation = string(PETSC_DIR, "/", PETSC_ARCH, "/lib/", "libpetsc")
global const petsc = libpetsclocation # for compatability with auto generated wrappers
global const libpetsc = Libdl.dlopen(libpetsclocation)

# definitions of Petsc types
# a way to automatically generate these would be preferable
# define all the data types that are known a priori
global const Petsc64bitInt = Int64
#global const PetscBLASInt = Int32
#global const PetscScalar = Cint
global const PetscBool = UInt32
global const PetscDataType = Cint  # C enums are Int32

global const PetscLogDouble = Cdouble
global const PetscErrorCode = Cint
#=
global const PETSC_PRECISION_SINGLE = (UInt32)(4)
global const PETSC_PRECISION_DOUBLE = (UInt32)(8)
=#

global const PetscViewer = Ptr{Void}

global const PETSC_FALSE = (UInt32)(0)
global const PETSC_TRUE = (UInt32)(1)

global const PETSC_INT = (Int32)(0)
global const PETSC_DOUBLE = (Int32)(1)
global const PETSC_COMPLEX = (Int32)(2)
global const PETSC_LONG = (Int32)(3)
global const PETSC_SHORT = (Int32)(4)
global const PETSC_FLOAT = (Int32)(5)
global const PETSC_CHAR = (Int32)(6)
global const PETSC_BIT_LOGICAL = (Int32)(7)
global const PETSC_ENUM = (Int32)(8)
global const PETSC_BOOL = (Int32)(9)
global const PETSC___FLOAT128 = (Int32)(10)
global const PETSC_OBJECT = (Int32)(11)
global const PETSC_FUNCTION = (Int32)(12)
global const PETSC_STRING = (Int32)(12)

# these functions are used for figuring out the sizes of the datatypes
# use the ones in PETSc.jl for writing programs
function PetscDataTypeFromString_(name::AbstractString)
    ptype = Array{Cint}(1)
    found = Array{PetscBool}(1)
    ccall((:PetscDataTypeFromString,petsc),PetscErrorCode,(Cstring,Ptr{PetscDataType},Ptr{PetscBool}), name, ptype, found)

    return ptype[1], convert(Bool, found[1])
end

function PetscDataTypeGetSize_(dtype::PetscDataType)
    datasize = Array{Csize_t}(1)
    ccall((:PetscDataTypeGetSize,petsc),PetscErrorCode,(PetscDataType,Ptr{Csize_t}), dtype, datasize)

    return datasize[1]
end

# define types that depend on the options Petsc was compiled with
(real_type_enum, found_real) = PetscDataTypeFromString_("Real")
(scalar_type_enum, found_scalar) = PetscDataTypeFromString_("Scalar")
int_size = PetscDataTypeGetSize_(PETSC_INT)

# confirm a value was found
@assert(found_real)
@assert(found_scalar)

# figure out what the values mean
if real_type_enum == PETSC_DOUBLE
  precision = "double"
elseif real_type_enum == PETSC_FLOAT
  precision = "single"
else
  if myrank == 0
    println("unknown type of Real")
    println("real_type_enum = ", real_type_enum)
  end
  throw(ErrorException("unsupported PetscReal ty pe"))
end

if scalar_type_enum == real_type_enum
  scalar_type = "real"
elseif scalar_type_enum == PETSC_COMPLEX
  scalar_type = "complex"
else
  if myrank == 0
    println("unknown type of Scalar")
    println("scalar_type_enum = ", scalar_type_enum)
  end
  throw(ErrorException("unsupported PetscScalar type"))
end

# figure out types

if scalar_type == "real"
  # single or double precision real
  if precision == "single"
    scalar_dtype = Float32
    real_dtype = Float32
  else
    scalar_dtype = Float64
    real_dtype = Float64
  end
else  # scalar_type = complex
  if precision == "single"
    scalar_dtype = Complex64
    real_dtype = Float32
  else
    scalar_dtype = Complex128
    real_dtype = Float64
  end
end

if int_size == 4
   int_dtype = Int32
elseif int_size == 8
   int_dtype = Int64
else
  println("unknown Int size")
  println("int_size = ", int_size)
  throw(ErrorException("unsupported integer size"))
end

if myrank == 0
  println("PetscScalar type = ", scalar_dtype)
  println("PetscReal type = ", real_dtype)
  println("PetscInt type = ", int_dtype)
end

"""
  Element type of Petsc vectors and matrices
"""
global const PetscScalar = scalar_dtype

"""
  The closest real type to `PetscScalar`
"""
global const PetscReal = real_dtype

"""
  Petsc integer type
"""
global const PetscInt = int_dtype

# some useful type unions
PetscInt_arr_or_null = Union{AbstractArray{PetscInt}, Ptr{Void}}

global const PetscIntNullArray = unsafe_wrap(Vector{PetscInt}, Ptr{PetscInt}(C_NULL), 0)
global const PetscScalarNullArray = unsafe_wrap(Vector{PetscScalar}, Ptr{PetscScalar}(C_NULL), 0)

#=
global const PETSC_PI = pi
global const PETSC_MAX_INT = 2147483647
global const PETSC_MIN_INT = -PETSC_MAX_INT - 1
global const PETSC_MAX_REAL = PETSC_MAX_INT  # made up
global const PETSC_INFINITY = PETSC_MAX_REAL / 4.0
global const PETSC_NINFINITY = -PETSC_INFINITY
#const PassiveReal = Cint  # what is this?
=#

# some useful constants
#const PassiveScalar = PetscScalar
#const MPIU_MATSCALAR = MPIU_SCALAR
#const MPIU_2INT = Cint
global const PETSC_NULL = C_NULL
global const PETSC_IGNORE = C_NULL
global const PETSC_DECIDE = convert(Int32, -1)
global const PETSC_DETERMINE = PETSC_DECIDE
global const PETSC_DEFAULT = convert(Int32, -2)
global const PETSC_COMM_SELF = MPI.COMM_SELF

global const MPI_Comm = MPI.Comm  # use MPI package communicator type
# typealias to size of MPI communicator value
# this is not defined by the MPI C standard, so it might be
# 32 bits (MPICH), or possibly 64 bits
global const comm_type = MPI.CComm

export PetscMatStructure
export DIFFERENT_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, SAME_NONZERO_PATTERN
export MATSEQAIJ, MATMPIAIJ, MATMPIBAIJ
export MAT_LOCAL, MAT_GLOBAL_MAX, MAT_GLOBAL_SUM
export MAT_INITIAL_MATRIX, MAT_REUSE_MATRIX, MAT_IGNORE_MATRIX
export KSPType
export KSPConvergedReason

# get the auto generated definitions
include("petsc_constants_gen.jl")

# map values to string for printing
export KSPConvergedReasonDict
global const KSPConvergedReasonDict = Dict{KSPConvergedReason, String}(
 KSP_CONVERGED_RTOL_NORMAL => "Converged: RTol normal",
 KSP_CONVERGED_ATOL_NORMAL => "Converged: ATol normal",
 KSP_CONVERGED_RTOL => "Converged: RTol",
 KSP_CONVERGED_ATOL => "converged: ATol",
 KSP_CONVERGED_ITS => "Converged: ITS",
 KSP_CONVERGED_CG_NEG_CURVE => "Converged: CG Negative Curvature",
 KSP_CONVERGED_CG_CONSTRAINED => "Converged: CG Constrained",
 KSP_CONVERGED_STEP_LENGTH => "Converged: step length",
 KSP_CONVERGED_HAPPY_BREAKDOWN => "Converged: happy breakdown",
 KSP_DIVERGED_NULL => "Diverged: Null",
 KSP_DIVERGED_ITS => "Diverged: iteration max",
 KSP_DIVERGED_DTOL => "Diverged divergence tolerance",
 KSP_DIVERGED_BREAKDOWN => "Diverged: Krylov breakdown",
 KSP_DIVERGED_BREAKDOWN_BICG => "Diverged: BICG breakdown",
 KSP_DIVERGED_NONSYMMETRIC => "Diverged: non symmetric",
 KSP_DIVERGED_INDEFINITE_PC => "Diverged: indefinitate PC",
 KSP_DIVERGED_NANORINF => "Diverged nan or inf",
 KSP_DIVERGED_INDEFINITE_MAT => "Diverged: indefinate matrix",
 KSP_DIVERGED_PCSETUP_FAILED => "Diverged: PC setup failed",
 KSP_CONVERGED_ITERATING => "Still running: please be patient",
 )


text:
#                        .-'''-.                               _..._
#                       '   _    \         _______          .-'_..._''.
#  /|                 /   /` '.   \        \  ___ `'.     .' .'      '.\
#  ||                .   |     \  '         ' |--.\  \   / .'
#  ||        .-,.--. |   '      |  '        | |    \  ' . '                                     .|
#  ||  __    |  .-. |\    \     / /  __     | |     |  '| |                 __                .' |_
#  ||/'__ '. | |  | | `.   ` ..' /.:--.'.   | |     |  || |              .:--.'.         _  .'     |
#  |:/`  '. '| |  | |    '-...-'`/ |   \ |  | |     ' .'. '             / |   \ |      .' |'--.  .-'
#  ||     | || |  '-             `" __ | |  | |___.' /'  \ '.          .`" __ | |     .   | / |  |
#  ||\    / '| |                  .'.''| | /_______.'/    '. `._____.-'/ .'.''| |   .'.'| |// |  |
#  |/\'..' / | |                 / /   | |_\_______|/       `-.______ / / /   | |_.'.'.-'  /  |  '.'
#  '  `'-'`  |_|                 \ \._,\ '/                          `  \ \._,\ '/.'   \_.'   |   /
#                                 `--'  `"                               `--'  `"             `'-'

using Base.Broadcast
using Base.Broadcast: Broadcasted, AbstractArrayStyle, broadcasted, materialize
using NNlib

@nograd Broadcast.combine_styles, Broadcast.result_style

# There's a saying that debugging code is about twice as hard as writing it in
# the first place. So if you're as clever as you can be when writing code, how
# will you ever debug it?

# AD faces a similar dilemma: if you write code that's as clever as the compiler
# can handle, how will you ever differentiate it? Differentiating makes clever
# code that bit more complex and the compiler gives up, usually resulting in
# 100x worse performance.

# Base's broadcasting is very cleverly written, and this makes differentiating
# it... somewhat tricky.

# Utilities
# =========

accum_sum(xs; dims = :) = reduce(accum, xs, dims = dims)

# Work around reducedim_init issue
# https://github.com/JuliaLang/julia/issues/31427
accum_sum(xs::Nothing; dims = :) = nothing
accum_sum(xs::AbstractArray{Nothing}; dims = :) = nothing
accum_sum(xs::AbstractArray{<:Number}; dims = :) = sum(xs, dims = dims)
accum_sum(xs::AbstractArray{<:AbstractArray{<:Number}}; dims = :) = sum(xs, dims = dims)
accum_sum(xs::Number; dims = :) = xs

trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

unbroadcast(x::AbstractArray, x̄) =
  size(x) == size(x̄) ? x̄ :
  length(x) == length(x̄) ? trim(x, x̄) :
    trim(x, accum_sum(x̄, dims = ntuple(i -> size(x, i) == 1 ? i : ndims(x̄)+1, Val(ndims(x̄)))))

unbroadcast(x::Number, x̄) = accum_sum(x̄)
unbroadcast(x::Tuple{<:Any}, x̄) = (accum_sum(x̄),)
unbroadcast(x::Base.RefValue, x̄) = (x=accum_sum(x̄),)

unbroadcast(x::AbstractArray, x̄::Nothing) = nothing

# Split Reverse Mode
# ==================

# TODO: use DiffRules here. It's complicated a little by the fact that we need
# to do CSE, then broadcast-ify the expression so that the closure captures the
# right arrays.

Numeric{T<:Number} = Union{T,AbstractArray{<:T}}

@adjoint broadcasted(::typeof(+), xs::Numeric...) =
  broadcast(+, xs...), ȳ -> (nothing, map(x -> unbroadcast(x, ȳ), xs)...)

@adjoint broadcasted(::typeof(*), x::Numeric, y::Numeric) = x.*y,
  z̄ -> (nothing, unbroadcast(x, z̄ .* conj.(y)), unbroadcast(y, z̄ .* conj.(x)))

@adjoint function broadcasted(::typeof(/), x::Numeric, y::Numeric)
  res = x ./ y
  res, Δ -> (nothing, unbroadcast(x, Δ ./ y), unbroadcast(y, -Δ .* res ./ y))
end

@adjoint function broadcasted(::typeof(Base.literal_pow), ::typeof(^), x::Numeric, exp::Val{p}) where p
  y = Base.literal_pow.(^, x, exp)
  y, ȳ -> (nothing, nothing, ȳ .* p .* conj.(x .^ (p - 1)), nothing)
end

@adjoint broadcasted(::typeof(identity), x::Numeric) = x, Δ -> (nothing, Δ)

@adjoint function broadcasted(::typeof(σ), x::Numeric)
  y = σ.(x)
  y, ȳ -> (nothing, ȳ .* conj.(y .* (1 .- y)))
end

@adjoint function broadcasted(::typeof(tanh), x::Numeric)
  y = tanh.(x)
  y, ȳ -> (nothing, ȳ .* conj.(1 .- y.^2))
end

@adjoint broadcasted(::typeof(conj), x::Numeric) =
  conj.(x), z̄ -> (nothing, conj.(z̄))

@adjoint broadcasted(::typeof(real), x::Numeric) =
  real.(x), z̄ -> (nothing, real.(z̄))

@adjoint broadcasted(::typeof(imag), x::Numeric) =
  imag.(x), z̄ -> (nothing, im .* real.(z̄))

# General Fallback
# ================

# The fused reverse mode implementation is the most general but currently has
# poor performance. It works by flattening the broadcast and mapping the call to
# `_pullback` over the input.

# However, the core call
# broadcast(_pullback, (cx,), f, args...)
# is already 10x slower than a simple broadcast (presumably due to inlining
# issues, or something similar) and the other operations needed take it to about
# 100x overhead.

@generated inclen(::NTuple{N,Any}) where N = Val(N+1)

# Avoid hitting special cases for `Adjoint` etc.
_broadcast(f::F, x...) where F = materialize(broadcasted(f, x...))

_get(x::Tuple, i) = x[i]
_get(::Nothing, i) = nothing
collapse_nothings(xs::Vector{Nothing}) = nothing
collapse_nothings(xs) = xs

@adjoint function broadcasted(::AbstractArrayStyle, f, args...)
  len = inclen(args)
  y∂b = _broadcast((x...) -> _pullback(__context__, f, x...), args...)
  y = map(x -> x[1], y∂b)
  ∂b = map(x -> x[2], y∂b)
  y, function (ȳ)
    dxs_zip = map((∂b, ȳ) -> ∂b(ȳ), ∂b, ȳ)
    dxs = collapse_nothings.(ntuple(i -> map(x -> _get(x, i), dxs_zip), len))
    (nothing, accum_sum(dxs[1]), map(unbroadcast, args, Base.tail(dxs))...)
  end
end

@adjoint function broadcasted(::AbstractArrayStyle{0}, f, args...)
  len = inclen(args)
  y, ∂b = _broadcast((x...) -> _pullback(__context__, f, x...), args...)
  y, function (ȳ)
    dxs = ∂b(ȳ)
    (nothing, dxs...)
  end
end

@adjoint! (b::typeof(broadcast))(f, args...) = _pullback(__context__, broadcasted, f, args...)

# Forward Mode (mainly necessary for CUDA)

import ForwardDiff
using ForwardDiff: Dual

dual(x, p) = x
dual(x::Real, p) = Dual(x, p)

dualtype(::Type{Dual{G,T,P}}) where {G,T,P} = T
dualtype(T) = T

function dual_function(f::F) where F
  function (args::Vararg{Any,N}) where N
    ds = map(args, ntuple(identity,Val(N))) do x, i
      dual(x, ntuple(j -> i==j, Val(N)))
    end
    return f(ds...)
  end
end

@inline function broadcast_forward(f, args::Vararg{Any,N}) where N
  T = Broadcast.combine_eltypes(f, args)
  out = dual_function(f).(args...)
  eltype(out) <: Dual || return (out, _ -> nothing)
  y = map(x -> x.value, out)
  _back(ȳ, i) = unbroadcast(args[i], ((a, b) -> a*b.partials[i]).(ȳ, out))
  back(ȳ) = ntuple(i -> _back(ȳ, i), N)
  return y, back
end

@init @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
  @adjoint function broadcasted(::Broadcast.ArrayStyle{CuArrays.CuArray}, f, args...)
    y, back = broadcast_forward(CuArrays.cufunc(f), args...)
    y, ȳ -> (nothing, nothing, back(ȳ)...)
  end

  @adjoint CuArrays.CuArray{N,T}(xs::Array) where {N,T} =
    CuArrays.CuArray{N,T}(xs), Δ -> (convert(Array, Δ), )

  @adjoint function sum(xs::CuArrays.CuArray; dims = :)
    placeholder = similar(xs)
    sum(xs, dims = dims), Δ -> (placeholder .= Δ,)
  end

  @adjoint function Base.convert(::Type{T}, xs::Array)  where {T<:CuArrays.CuArray}
    Base.convert(T, xs), Δ -> (nothing, Base.convert(Array, Δ),)
  end

end


text:
using Turing
using Turing: Sampler
using Test

alg = SGHMC(1000, 0.01, 0.1)
sampler = Sampler(alg)
@test isa(alg, SGHMC)
@test isa(sampler, Sampler{<:Turing.SGHMC})

alg = SGHMC(200, 0.01, 0.1, :m)
sampler = Sampler(alg)
@test isa(alg, SGHMC)
@test isa(sampler, Sampler{<:Turing.SGHMC})

alg = SGHMC(1000, 0.01, 0.1, :s)
sampler = Sampler(alg)
@test isa(alg, SGHMC)
@test isa(sampler, Sampler{<:Turing.SGHMC})


text:
function getopts(gt::GradientTolerance{T}, maxiter::Int, time_limit::Real = NaN) where T
    # TODO: Figure out a way to make passing options more flexible (time_limit etc.)
    opts = Optim.Options(f_tol = zero(T), x_tol = zero(T),
                         g_tol = gt.g0norm*gt.tol,
                         allow_f_increases = true,
                         iterations = maxiter,
                         time_limit = time_limit)#,
                         #show_trace=true,extended_trace=true, store_trace=true)
    return opts
end

function getopts(ft::FunctionTolerance{T}, maxiter::Int, time_limit::Real = NaN) where T
    # TODO: Figure out a way to make passing options more flexible (time_limit etc.)
    cb = FunctionStop(ft.fL + ft.tol*(ft.f0 - ft.fL))
    opts = Optim.Options(f_tol = zero(T), x_tol = zero(T), g_tol = zero(T),
                         allow_f_increases = true,
                         iterations = maxiter, callback = cb,
                         time_limit = time_limit)
    return opts
end

"Return time elapsed and result of expression"
macro my_elapsed(ex)
    quote
        local t0 = time_ns()
        local val = $(esc(ex))
        (time_ns()-t0)/1e9, val
    end
end

function runproblem(df, x0, solver::Optim.AbstractOptimizer, solvername::AbstractString,
                    problemname::AbstractString,
                    metrictol::StopTolerance;
                    numrecordtime::Int = 0, maxiter::Int = 1000,
                    time_limit::Real = NaN)
    opts = getopts(metrictol, maxiter, time_limit)

    f0, g0 = value_gradient!(df, x0)
    g0norm = norm(g0, Inf)
    try
        # TODO: Does this help for precompilation at all?
        # This can be unnecessarily costly if the line search goes bananas
        optimize(df, x0, solver, Optim.Options(iterations=1))
    catch e
        warn(e)
    end

    reset!(df) # Get correct f_calls numbering
    sname = isempty(solvername) ? summary(solver) : solvername
    orun = OptimizationRun(0, 0, 0, 0, 0.0, f0, f0, g0norm, g0norm,
                           sname, problemname, metrictol)
    try
        runtime, r = @my_elapsed optimize(df, x0, solver, opts)
        @show Optim.f_calls(r), Optim.g_calls(r)
        for k = 2:numrecordtime
            # TODO: Do this with BenchmarkTools?
            reset!(df) # Reset df to get f_calls etc. correct
            tmptime, tmpr = @my_elapsed optimize(df, x0, solver, opts)
            if tmptime < runtime
                runtime = tmptime
                r = tmpr # In case previous r hit time_limit
            end
        end

        orun = OptimizationRun(Optim.iterations(r), Optim.f_calls(r), Optim.g_calls(r),
                               Optim.h_calls(r), runtime, Optim.minimum(r), f0,
                               Optim.g_residual(r), g0norm, sname, problemname,
                               metrictol)
    catch e
        warn(e)
    end
    return orun
end

"""
Loop over solvers for OptimizationProblem `prob` and optimize until a given `stoptolerance`.

Return vector of `OptimizationRun`s for a given problem.
Optionally include runtime (runs optimization multiple times per solver).
"""
function createruns(prob::OptimizationProblem, problemname::AbstractString,
                    solvers::AbstractVector, solvernames::AbstractVector{<:AbstractString},
                    tol::T, stoptype::Symbol = :GradientTolRelative, timelog::Int = 0,
                    maxiter::Int = 1000; time_limit::Real = NaN,
                    verbose::Bool = true) where T <: Real
    # TODO: Is it problematic for specialized compilation to use stoptype::Symbol?
    verbose && println("Running $problemname ...")
    local x0, f0, g0norm # Ensure variables stay in method scope
    try
        x0 = initial_x(prob)
        df0 = optim_problem(prob)
        f0, g0 = value_gradient!(df0, x0)
        g0norm = norm(g0, Inf)
    catch e
        warn(e)
        return Vector{OptimizationRun{typeof(time_ns()/1e9),T}}(0)
    end

    if stoptype == :GradientTolRelative
        metrictol = GradientTolerance(tol,g0norm)
    elseif stoptype == :GradientTol
        metrictol = GradientTolerance(tol,one(T))
    elseif stoptype == :FunctionTolRelative
        fL = solver_optimum(prob)
        @assert isfinite(fL) # TODO: create function to find solver_optimum
        metrictol = FunctionTolerance(tol, f0, fL)
    elseif stoptype == :FunctionTol
        fL = solver_optimum(prob)
        @assert isfinite(fL)
        metrictol = FunctionTolerance(tol, one(fL) + fL, fL)
    else
        Base.error("The parameter `stoptype` must be one of :GradientTolRelative, :GradientTol, :FunctionTolRelative, or :FunctionTol.")
    end

    oruns = Vector{OptimizationRun{typeof(time_ns()/1e9),T}}(length(solvers))
    for (k, solver) in enumerate(solvers)
        df = optim_problem(prob)
        try
            oruns[k] = runproblem(df, x0, solver, solvernames[k],
                                  problemname, metrictol;
                                  numrecordtime=timelog, maxiter = maxiter,
                                  time_limit = time_limit)
        catch probsolve
            warn(probsolve)
            oruns[k] = OptimizationRun(0,0,0,0,zero(time_ns()/1e9),
                                       f0,f0,g0norm,g0norm,
                                       solvernames[k],problemname,false)
        end
    end
    verbose && println("Finished $problemname.")

    return oruns
end

"""
Loop over solvers for CUTEst problem `cutestname` and optimize until a given `stoptolerance`.

Return vector of `OptimizationRun`s for a given problem.
Optionally include runtime (runs optimization multiple times per solver).
"""
function createruns(cutestname::AbstractString,
                    solvers::AbstractVector, solvernames::AbstractVector{<:AbstractString},
                    tol::T, stoptype::Symbol = :GradientTolRelative, timelog::Int = 0,
                    maxiter::Int = 1000; time_limit::Real = NaN,
                    decode::Bool = false,
                    verbose::Bool = true) where T <: Real
    verbose && println("Loading $cutestname ...")
    local oruns
    try
        nlp = CUTEstModel(cutestname, decode=decode)
        prob = optimizationproblem(nlp)
        oruns = createruns(prob, prob.name, solvers, solvernames,
                           tol, stoptype, timelog, maxiter;
                           verbose = verbose, time_limit = time_limit)
        finalize(nlp)
    catch e
        warn(e)
        oruns = Vector{OptimizationRun{typeof(time_ns()/1e9),T}}(0)
    end
    return oruns
end
createruns(cutestname::AbstractString, ts::TestSetup) =
    createruns(cutestname, ts.solvers, ts.solvernames,
               ts.stoptol, ts.stoptype, ts.timelog, ts.maxiter;
               time_limit = ts.timelimit)

"""
Create performance measures of a list of solvers on a collection of OptimizationProblems

Returns a vector of `OptimizationRun`s.
"""
function createmeasures(problems::AbstractVector{<:OptimizationProblem}, problemnames::AbstractVector{<:AbstractString},
                        solvers::AbstractVector, solvernames::AbstractVector{<:AbstractString},
                        tol::T, stoptype::Symbol = :GradientTolRelative, timelog::Int = 0,
                        maxiter::Int = 1000;
                        time_limit::Real = NaN) where T <: Real
    oruns = reduce(append!, pmap(k->createruns(problems[k], problemnames[k],
                                               solvers, solvernames, tol,
                                               stoptype, timelog, maxiter;
                                               time_limit = time_limit),
                                 1:length(problems)))
    return oruns
end

createmeasures(problems::AbstractVector{<:OptimizationProblem}, problemnames::AbstractVector{<:AbstractString},
               ts::TestSetup) =
                   createmeasures(problems, problemnames, ts.solvers, ts.solvernames,
                                  ts.stoptol, ts.stoptype, ts.timelog, ts.maxiter;
                                  time_limit = ts.timelimit)

"""
Create performance measures of a list of solvers on a collection of CUTEst problems.

Currently, the CUTEst models are run with the default parameters.

Returns a vector of `OptimizationRun`s.
"""
function createmeasures(cutestnames::AbstractVector{<:AbstractString},
                        solvers::AbstractVector, solvernames::AbstractVector{<:AbstractString},
                        tol::T, stoptype::Symbol = :GradientTolRelative, timelog::Int = 0,
                        maxiter::Int = 1000; time_limit::Real = NaN) where T <: Real
    oruns = reduce(append!, pmap(k->createruns(cutestnames[k],
                                               solvers, solvernames, tol,
                                               stoptype, timelog, maxiter;
                                               decode = false,
                                               time_limit = time_limit),
                                 1:length(cutestnames)))
    return oruns
end

createmeasures(cutestnames::AbstractVector{<:AbstractString},
               ts::TestSetup) =
                   createmeasures(cutestnames, ts.solvers, ts.solvernames,
                                  ts.stoptol, ts.stoptype, ts.timelog, ts.maxiter;
                                  time_limit = ts.timelimit)

"""
Take a vector of `OptimizationRun`s representing performance measures and
convert it to a DataFrame of performance measures
"""
function createmeasuredataframe(oruns::Vector{OptimizationRun{T,Tf}}) where T where Tf
    n         = length(oruns)
    snames    = Vector{String}(n)
    pnames    = similar(snames)
    iters     = Vector{Int}(n)
    fcallarr  = Vector{Int}(n)
    gcallarr  = similar(fcallarr)
    hcallarr  = similar(fcallarr)
    cputimes  = zeros(T, n)
    fvalarr   = zeros(Tf, n)
    f0arr     = similar(fvalarr)
    gnormarr  = similar(fvalarr)
    g0normarr = similar(fvalarr)
    succarr   = Vector{Bool}(n)

    for (k, orun) in enumerate(oruns)
        pnames[k]    = orun.problemname
        snames[k]    = orun.solvername
        iters[k]     = orun.iterations
        fcallarr[k]  = orun.fevals
        gcallarr[k]  = orun.gevals
        hcallarr[k]  = orun.hevals
        cputimes[k]  = orun.cputime
        fvalarr[k]   = orun.fval
        f0arr[k]     = orun.f0
        gnormarr[k]  = orun.gnorm
        g0normarr[k] = orun.g0norm
        succarr[k]   = orun.success
    end

    dframe = DataFrame(Problem = pnames, Solver = snames, Iterations = iters,
                       fcalls = fcallarr, gcalls = gcallarr, hcalls = hcallarr,
                       CPUtime = cputimes, fval = fvalarr, f0 = f0arr,
                       gnorm = gnormarr, g0norm = g0normarr, Success = succarr)
end

" Return performance ratios x / minimum(x), but set elements where `z[k] = false` to `Inf`"
function _nanratio(x, z)
    retval = x ./ minimum(x) # ./ deals with some issues
    # Bool needed since introducing Missings.Missing data type
    retval[Bool.(.!z .| isnan.(retval) .| ismissing.(retval))] = Inf # Should we use NA here?
    return retval
end

"Create a DataFrame of performance ratios from a dataframe of performance measures."
function createratiodataframe(dframe::DataFrame, solvers::AbstractVector{<:AbstractString} = Vector{String}(0))
    if isempty(solvers)
        idx = Colon()
    else
        idx = broadcast(x -> x ∈ solvers, dframe[:Solver])
    end

    by(dframe[idx, :], :Problem) do dfr
        DataFrame(Solver = dfr[:Solver],
                  fcalls = _nanratio(dfr[:fcalls], dfr[:Success]),
                  gcalls = _nanratio(dfr[:gcalls], dfr[:Success]),
                  hcalls = _nanratio(dfr[:hcalls], dfr[:Success]),
                  CPUtime = _nanratio(dfr[:CPUtime], dfr[:Success]),
                  Iterations = _nanratio(dfr[:Iterations], dfr[:Success]),
                  Success = dfr[:Success])
    end
end

"Create a data frame with performance ratios from a vector of `OptimizationRun`s"
createratiodataframe(oruns::Vector{OptimizationRun},
                     solvers::AbstractVector{<:AbstractString} = Vector{String}(0)) =
                         createratiodataframe(createmeasuredataframe(oruns), solvers)

"""
Ensure all performance metrics/ratios are finite, by setting
non-finite value and failures to a "failure" value according to
the rule `maxfun` (default x->2*maximum(x), taken over finite values).
"""
function makefinite(rdf::DataFrame, maxfun = x-> 2*maximum(x))
    cprdf = deepcopy(rdf)
    makefinite!(cprdf, maxfun)
    return cprdf
end

"""
Ensure all performance metrics/ratios are finite, by setting
non-finite value and failures to a "failure" value according to
the rule `maxfun` (default x->2*maximum(x), taken over finite values).
"""
function makefinite!(rdf::DataFrame, maxfun = x-> 2*maximum(x))
    pltvals = [:fcalls, :gcalls,:CPUtime,:Iterations]
    for k in 1:length(pltvals)
        finidx = convert(Array{Bool}, isfinite.(rdf[pltvals[k]]))
        if any(finidx)
            rdf[.!finidx, pltvals[k]] = maxfun(rdf[finidx,pltvals[k]])
        else
            rdf[.!finidx, pltvals[k]] = -1.0
        end
    end
end

"Create violin plot of performance metric/ratios grouped by solver."
createviolins(rdf::DataFrame, yscale::Symbol = :log2) = createstatplots(rdf,violin; yscale = yscale, kwargs...)

"Create box plot of performance metric/ratios grouped by solver."
createboxplots(rdf::DataFrame, yscale::Symbol = :log2) = createstatplots(rdf,boxplot; yscale = yscale, kwargs...)

"Create plot of performance metric/ratios grouped by solver. Defaults to violin plots."
function createstatplots(rdf::DataFrame, fun::Function = violin;
                         yscale = :log2, maxfun = x -> 2*maximum(x),
                         kwargs...)
    pltlabels = ["f-calls", "g-calls", "CPU time", "Iterations"]
    plts = Vector{Plots.Plot}(length(pltlabels))

    cprdf = makefinite(rdf, maxfun)

    plts[1] = @df cprdf fun(:Solver, :fcalls, label=pltlabels[1])
    plts[2] = @df cprdf fun(:Solver, :gcalls, label=pltlabels[2])
    plts[3] = @df cprdf fun(:Solver, :CPUtime, label=pltlabels[3])
    plts[4] = @df cprdf fun(:Solver, :Iterations, label=pltlabels[4])
    return plot(plts...; yscale = yscale, kwargs...)
end

"""
Create performance profiles from a DataFrame of performance ratios.
Defaults to performance profiles of f-calls.
"""
function createperfprofiles(rdf::DataFrame, sym::Symbol = :fcalls;
                            t::Symbol = :steppost, xscale::Symbol = :log2,
                            ylims = (0,1),
                            maxfun = x -> 2*maximum(x), kwargs...)
    @assert sym in [:fcalls, :gcalls, :CPUtime, :Iterations]
    cprdf = makefinite(rdf, maxfun)
    plt = if sym == :fcalls
        @> cprdf begin
            @splitby _.Solver
            @x _.fcalls
            @y :cumulative
            @plot
        end
    elseif sym == :gcalls
        @> cprdf begin
            @splitby _.Solver
            @x _.gcalls
            @y :cumulative
            @plot
        end
    elseif sym == :CPUtime
        @> cprdf begin
            @splitby _.Solver
            @x _.CPUtime
            @y :cumulative
            @plot
        end
    elseif sym == :Iterations
        @> cprdf begin
            @splitby _.Solver
            @x _.Iterations
            @y :cumulative
            @plot
        end
    end
    return plot(plt; t=t, xscale=xscale, ylims=ylims, kwargs...)
end

"""
Create performance profiles from a vector of OptimizationRuns.
Defaults to performance profiles of f-calls.
"""
function createperfprofiles(oruns::Vector{OptimizationRun}, sym::Symbol = :fcalls,
                            solvers::AbstractVector{<:AbstractString} = Vector{String}(0);
                            t::Symbol = :steppost, xscale::Symbol = :log2,
                            ylims=(0,1),
                            kwargs...)
    createperfprofiles(createratiodataframe(oruns, solvers), sym;
                       t=t, xscale=xscale, ylims=ylims, kwargs...)
end


text:
# Autogenerated wrapper script for GR_jll for aarch64-apple-darwin
export gksqt, libGKS, libGR, libGR3, libGRM

using Bzip2_jll
using Cairo_jll
using FFMPEG_jll
using Fontconfig_jll
using GLFW_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using Pixman_jll
using Qt5Base_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("GR")
JLLWrappers.@declare_library_product(libGKS, "@rpath/libGKS.dylib")
JLLWrappers.@declare_library_product(libGR, "@rpath/libGR.dylib")
JLLWrappers.@declare_library_product(libGR3, "@rpath/libGR3.dylib")
JLLWrappers.@declare_library_product(libGRM, "@rpath/libGRM.dylib")
JLLWrappers.@declare_executable_product(gksqt)
function __init__()
    JLLWrappers.@generate_init_header(Bzip2_jll, Cairo_jll, FFMPEG_jll, Fontconfig_jll, GLFW_jll, JpegTurbo_jll, libpng_jll, Libtiff_jll, Pixman_jll, Qt5Base_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libGKS,
        "lib/libGKS.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGR,
        "lib/libGR.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGR3,
        "lib/libGR3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGRM,
        "lib/libGRM.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        gksqt,
        "bin/gksqt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Compat, ClassicalLaminateTheory

forcesAndMoments =
  @compat Dict(
    :Nx => 300,
    :Ny => 100,
    :Nxy => 0,
    :Mx => 5,
    :My => 0,
    :Mxy => 0)

materialProperties = @compat Dict(
  :E1 => 20.0*10^6,
  :E2 => 1.5*10^6,
  :G12 => 1.0*10^6,
  :ν12 => .34,
  :alpha1 => 0.2e-6,
  :alpha2 => 20e-6,
  :beta1 => 0.2e-4,
  :beta2 => 20e-4,
  :thickness => 0.0075)

laminateProperties = @compat Dict(
  :nplies => 4,
  :repeats => 1,
  :symmetric => true,
  :orientation => [0, 45, 90, -45],
  :deltaTemp => 0.0,
  :deltaMoisture => 0
)

createLaminate!(laminateProperties, materialProperties, forcesAndMoments)

ABDmat = [
  525698.0  155811.0       0.0    0.0    -0.0     0.0;
  155811.0  525698.0       0.0   -0.0     0.0     0.0;
       0.0       0.0  184943.0    0.0     0.0     0.0;
       0.0      -0.0       0.0  250.495  32.687  23.619;
      -0.0       0.0       0.0   32.687  93.036  23.619;
       0.0       0.0       0.0   23.619  23.619  41.427
];

@assert round(laminateProperties[:ABD]) == round(ABDmat)
println("Tests for CLT1 passed!")


text:
function spiral!(r)
    a = 1
    s = Nord
    while !ismarker(r)
        for i in 1:a
            if ismarker(r)
                break
            end
            move!(r, s)
        end
        s = rotate_right(s)
        for i in 1:a
            if ismarker(r)
                break
            end
            move!(r, s)
        end
        if ismarker(r)
            break
        end
        s = rotate_right(s)
        a += 1
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs