text:
import Base: rand,
             ∈,
             isempty

export Line,
       an_element,
       translate,
       translate!

"""
    Line{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a line in the form

```math
    \\{y ∈ \\mathbb{R}^n: y = p + λd, λ ∈ \\mathbb{R}\\}
```
where ``p`` is a point on the line and ``d`` is its direction vector (not necessarily
normalized).

### Fields

- `p` -- point on the line
- `d` -- direction

### Examples

The line passing through the point ``[-1, 2, 3]`` and parallel to the vector
``[3, 0, -1]``:

```jldoctest
julia> Line([-1, 2, 3.], [3, 0, -1.])
Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [3.0, 0.0, -1.0])
```
"""
struct Line{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
    p::VN
    d::VN

    # default constructor with length constraint
    function Line(p::VN, d::VN; check_direction::Bool=true) where {N, VN<:AbstractVector{N}}
        if check_direction
            @assert !iszero(d) "a line needs a non-zero direction vector"
        end
        return new{N, VN}(p, d)
    end
end

isoperationtype(::Type{<:Line}) = false
isconvextype(::Type{<:Line}) = true

"""
    Line(; from::AbstractVector, to::AbstractVector, normalize=true)

Constructor of a line given two points.

### Input

- `from`      -- point
- `to`        -- another point
- `normalize` -- (optional, default: `true`) if `true`, the direction of the line
                 has norm 1 (w.r.t the Euclidean norm)

### Output

The line which passes through `from` and `to`.

### Algorithm

Given two points ``p ∈ \\mathbb{R}^n`` and ``q ∈ \\mathbb{R}^n``, the line that
passes through these two points is
`L: `\\{y ∈ \\mathbb{R}^n: y = p + λ(q - p), λ ∈ \\mathbb{R}\\}``.
"""
function Line(; from::AbstractVector, to::AbstractVector, normalize=true)
    d = from - to
    if normalize && iszero(d)
        throw(ArgumentError("points `$from` and `$to` should be distinct"))
    end
    d_n = normalize ? d / dot(d, d) : d
    return Line(from, d_n)
end

"""
    Line(a::AbstractVector{N}, b::N; normalize=true) where {N}

Constructor of a two-dimensional line given `ax = b`.

### Input

- `a`         -- two-dimensional vector
- `b`         -- scalar
- `normalize` -- (optional, default: `true`) if `true`, the direction of the line
                 has norm 1 (w.r.t the Euclidean norm)

### Output

The `Line` that satisfies `a ⋅ x = b`.
"""
function Line(a::AbstractVector{N}, b::N; normalize=true) where {N}
    @assert length(a) == 2 "expected a normal vector of length two, but it is $(length(a))-dimensional"

    got_horizontal = iszero(a[1])
    got_vertical = iszero(a[2])

    if got_horizontal && got_vertical
        throw(ArgumentError("the vector $a must be non-zero"))
    end

    if got_horizontal
        α = b / a[2]
        p = [zero(N), α]
        q = [one(N), α]
    elseif got_vertical
        β = b / a[1]
        p = [β, zero(N)]
        q = [β, one(N)]
    else
        α = b / a[2]
        μ = a[1] / a[2]
        p = [zero(N), α]
        q = [one(N), α - μ]
    end
    return Line(from=p, to=q, normalize=normalize)
end

"""
    direction(L::Line)

Return the direction of the line.

### Input

- `L` -- line

### Output

The line's field corresponding to the direction to the line.

### Notes

The direction is not necessarily normalized.
See [`normalize(::Line, ::Real)`](@ref) / [`normalize!(::Line, ::Real)`](@ref)
for such operation.
"""
direction(L::Line) = L.d

"""
    normalize!(L::Line, p::Real=2.0)

Normalize the direction of a line storing the result in `L`.

### Input

- `L` -- line
- `p` -- (optional, default: `2.0`) vector `p`-norm used in the normalization

### Output

A line whose direction has unit norm w.r.t the given `p`-norm.
"""
function normalize!(L::Line, p::Real=2.0)
    normalize!(L.d, p)
    return L
end

"""
    normalize(L::Line, p::Real=2.0)

Normalize the direction of a line.

### Input

- `L` -- line
- `p` -- (optional, default: `2.0`) vector `p`-norm used in the normalization

### Output

A line whose direction has unit norm w.r.t the given `p`-norm.

### Notes

See also [`normalize!(::Line, ::Real)`](@ref) for the in-place version.
"""
function normalize(L::Line, p::Real=2.0)
    normalize!(copy(L), p)
end

# --- polyhedron interface functions ---

"""
    constraints_list(L::Line{N, VN}) where {N, VN}

Return the list of constraints of a line.

### Input

- `L` -- line

### Output

A list containing `2n-2` half-spaces whose intersection is `L`, where `n` is the
ambient dimension of `L`.
"""
function constraints_list(L::Line{N, VN}) where {N, VN}
    p = L.p
    n = length(p)
    d = reshape(L.d, 1, n)
    K = nullspace(d)
    m = size(K, 2)
    @assert m == n - 1 "expected $(n - 1) normal half-spaces, got $m"

    out = Vector{HalfSpace{N, VN}}(undef, 2m)
    idx = 1
    @inbounds for j in 1:m
        Kj = K[:, j]
        b = dot(Kj, p)
        out[idx] = HalfSpace(Kj, b)
        out[idx+1] = HalfSpace(-Kj, -b)
        idx += 2
    end
    return out
end

# --- LazySet interface functions ---

"""
    dim(L::Line)

Return the ambient dimension of a line.

### Input

- `L` -- line

### Output

The ambient dimension of the line.
"""
dim(L::Line) = length(L.p)

"""
    ρ(d::AbstractVector, L::Line)

Return the support function of a line in a given direction.

### Input

- `d` -- direction
- `L` -- line

### Output

The support function in the given direction.
"""
ρ(d::AbstractVector, L::Line) = _ρ(d, L)

function _ρ(d::AbstractVector, L::Line)
    if isapproxzero(dot(d, L.d))
        return dot(d, L.p)
    else
        return Inf
    end
end

"""
    σ(d::AbstractVector, L::Line)

Return the support vector of a line in a given direction.

### Input

- `d` -- direction
- `L` -- line

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector, L::Line)
    if isapproxzero(dot(d, L.d))
        return L.p
    else
        throw(ArgumentError("the support vector is undefined because the line is unbounded"))
    end
end

"""
    isbounded(L::Line)

Determine whether a line is bounded.

### Input

- `L` -- line

### Output

`false`.
"""
isbounded(::Line) = false

"""
    isuniversal(L::Line; [witness::Bool]=false)

Check whether a line is universal.

### Input

- `P`       -- line
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` is `false`: `true` if the ambient dimension is one, `false` otherwise

* If `witness` is `true`: `(true, [])` if the ambient dimension is one,
  `(false, v)` where ``v ∉ P`` otherwise
"""
isuniversal(L::Line; witness::Bool=false) = isuniversal(L, Val(witness))

# TODO implement case with witness
isuniversal(L::Line, ::Val{false}) = dim(L) == 1

"""
    an_element(L::Line)

Return some element of a line.

### Input

- `L` -- line

### Output

An element on the line.
"""
an_element(L::Line) = L.p

"""
    ∈(x::AbstractVector, L::Line)

Check whether a given point is contained in a line.

### Input

- `x` -- point/vector
- `L` -- line

### Output

`true` iff `x ∈ L`.

### Algorithm

The point ``x`` belongs to the line ``L : p + λd`` if and only if
``x - p`` is proportional to the direction ``d``.
"""
function ∈(x::AbstractVector, L::Line)
    @assert length(x) == dim(L) "expected the point and the line to have the same dimension, " *
                                "but they are $(length(x)) and $(dim(L)) respectively"
    _isapprox(x, L.p) && return true

    # TODO pass "minus" option to samedir
    return first(samedir(x - L.p, L.d)) || first(samedir(x - L.p, -L.d))
end

"""
    rand(::Type{Line}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random line.

### Input

- `Line` -- type for dispatch
- `N`    -- (optional, default: `Float64`) numeric type
- `dim`  -- (optional, default: 2) dimension
- `rng`  -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed` -- (optional, default: `nothing`) seed for reseeding

### Output

A random line.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
"""
function rand(::Type{Line};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)

    rng = reseed(rng, seed)
    d = randn(rng, N, dim)
    while iszero(d)
        d = randn(rng, N, dim)
    end
    p = randn(rng, N, dim)
    return Line(p, d)
end

"""
    isempty(L::Line)

Return if a line is empty or not.

### Input

- `L` -- line

### Output

`false`.
"""
isempty(::Line) = false

"""
    translate(L::Line, v::AbstractVector)

Translate (i.e., shift) a line by a given vector.

### Input

- `L` -- line
- `v` -- translation vector

### Output

A translated line.

### Notes

See also `translate!` for the in-place version.
"""
function translate(L::Line, v::AbstractVector)
    return translate!(copy(L), v)
end

"""
    translate!(L::Line, v::AbstractVector)

Translate (i.e., shift) a line by a given vector storing the result in `L`.

### Input

- `L` -- line
- `v` -- translation vector

### Output

A translated line, modifying `L` in-place.
"""
function translate!(L::Line, v::AbstractVector)
    @assert length(v) == dim(L) "cannot translate a $(dim(L))-dimensional " *
                                "set by a $(length(v))-dimensional vector"

    L.p .+= v
    return L
end

"""
    distance(x::AbstractVector, L::Line; [p]::Real=2.0)

Compute the distance between point `x` and the line with respect to the given
`p`-norm.

### Input

- `x` -- vector
- `L` -- line
- `p` -- (optional, default: `2.0`) the `p`-norm used; `p = 2.0` corresponds to
         the usual Euclidean norm

### Output

A scalar representing the distance between `x` and the line `L`.
"""
function distance(x::AbstractVector, L::Line; p::Real=2.0)
    d = L.d  # direction of the line
    t = dot(x - L.p, d) / dot(d, d)
    return distance(x, L.p + t*d; p=p)
end

distance(L::Line, x::AbstractVector; p::Real=2.0) = distance(x, L; p=p)

"""
    linear_map(M::AbstractMatrix, L::Line)

Concrete linear map of a line.

### Input

- `M` -- matrix
- `L` -- line

### Output

The line obtained by applying the linear map to the point and direction of `L`.
"""
function linear_map(M::AbstractMatrix, L::Line)
    @assert dim(L) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
                                 "applied to a set of dimension $(dim(L))"

    Mp = M * L.p
    Md = M * L.d
    return Line(Mp, Md)
end

function project(L::Line{N}, block::AbstractVector{Int}; kwargs...) where {N}
    d = L.d[block]
    if iszero(d)
        return Singleton(L.p[block])  # projected out all nontrivial dimensions
    elseif length(d) == 1
        return Universe{N}(1)  # special case: 1D line is a universe
    else
        return Line(L.p[block], d)
    end
end


text:
function LRect(parent::Scene; bbox = nothing, kwargs...)
    attrs = merge!(Attributes(kwargs), default_attributes(LRect))

    @extract attrs (color, visible, valign, halign, padding, strokewidth,
        strokevisible, strokecolor)

    sizeattrs = sizenode!(attrs.width, attrs.height)
    alignment = lift(tuple, halign, valign)

    suggestedbbox = create_suggested_bboxnode(bbox)

    autosizenode = Node{NTuple{2, Optional{Float32}}}((nothing, nothing))

    computedsize = computedsizenode!(sizeattrs, autosizenode)

    finalbbox = alignedbboxnode!(suggestedbbox, computedsize, alignment, sizeattrs, autosizenode)

    strokecolor_with_visibility = lift(strokecolor, strokevisible) do col, vis
        vis ? col : RGBAf0(0, 0, 0, 0)
    end

    r = poly!(parent, finalbbox, color = color, visible = visible, raw = true,
        strokecolor = strokecolor_with_visibility, strokewidth = strokewidth)[end]

    # no protrusions
    protrusions = Node(RectSides(0f0, 0f0, 0f0, 0f0))

    layoutnodes = LayoutNodes{LRect, GridLayout}(suggestedbbox, protrusions, computedsize, autosizenode, finalbbox, nothing)

    # trigger bbox
    suggestedbbox[] = suggestedbbox[]

    LRect(parent, layoutnodes, r, attrs)
end

function Base.delete!(lr::LRect)
    disconnect_layoutnodes!(lr.layoutnodes.gridcontent)
    remove_from_gridlayout!(lr.layoutnodes.gridcontent)
    empty!(lr.layoutnodes.suggestedbbox.listeners)
    empty!(lr.layoutnodes.computedbbox.listeners)
    empty!(lr.layoutnodes.computedsize.listeners)
    empty!(lr.layoutnodes.autosize.listeners)
    empty!(lr.layoutnodes.protrusions.listeners)

    # remove the plot object from the scene
    delete!(lr.parent, lr.rect)
end

ext: jl


text:
using ProfileSVG
using Test
using Profile

using FlameGraphs
using Base.StackTraces: StackFrame
using Colors

function stackframe(func, file, line; C=false, inlined=false)
    StackFrame(Symbol(func), Symbol(file), line, nothing, C, inlined, 0)
end

backtraces = UInt[0, 4, 3, 2, 1,
                  0, 4, 3, 2, 1,
                  0,    6, 5, 1,
                  0,       8, 7,
                  0]

lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                 2=>stackframe(:jl_f, :filec, 55; C=true),
                                 3=>stackframe(:jl_invoke, :file2, 1; C=true),
                                 4=>stackframe(:_ZL, Symbol("libLLVM-8.0.so"), 15),
                                 5=>stackframe(:f4, :file1, 20),
                                 6=>stackframe(:copy, Symbol(".\\expr.jl"), 1),
                                 7=>stackframe(:f1, :file1, 2),
                                 8=>stackframe(:typeinf, Symbol("./compiler/typeinfer.jl"), 10))

function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A; dims=2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B; dims=1)
        b = rand(100)
        C = B.*b
    end
end

function count_element(pattern::Regex, svg::AbstractString)
    n = 0
    i = firstindex(svg)
    while true
        r = findnext(pattern, svg, i)
        r === nothing && return n
        n += 1
        i = last(r)
    end
end

@test detect_ambiguities(ProfileSVG) == []

@testset "showable" begin
    fg = ProfileSVG.view(backtraces, lidict=lidict)
    @test Base.showable(MIME"image/svg+xml"(), fg)
    @test Base.showable(MIME"text/html"(), fg)

    @test !Base.showable(MIME"image/svg+xml"(), ProfileSVG.default_config)
    @test !Base.showable(MIME"text/html"(), ProfileSVG.default_config)
end

@testset "view" begin
    sfc = StackFrameCategory()

    Profile.init(delay=0.002)
    fg = ProfileSVG.view(sfc, backtraces,
                         C=true, lidict=lidict, width=123.4, unknown=nothing)
    @test FlameGraphs.depth(fg.g) == 5
    @test fg.fcolor isa StackFrameCategory
    @test fg.graph_options[:C] == true
    @test fg.graph_options[:lidict] == lidict
    @test fg.bgcolor == :fcolor
    @test fg.fontcolor == :fcolor
    @test fg.frameopacity == 1.0
    @test fg.yflip == false
    @test fg.width == 123.4
    @test fg.height == 0
    @test fg.roundradius == 2
    @test fg.font == "inherit"
    @test fg.fontsize == 12
    @test fg.notext == false
    @test fg.timeunit == :none
    @test fg.delay == 0.002

    Profile.init(delay=0.001)

    fg = ProfileSVG.view(backtraces,
                         C=true, lidict=lidict, fontsize=12.34, unknown=missing)
    @test FlameGraphs.depth(fg.g) == 5
    @test fg.fcolor isa FlameColors
    @test fg.graph_options[:C] == true
    @test fg.graph_options[:lidict] == lidict
    @test fg.yflip == false
    @test fg.width == 960
    @test fg.height == 0
    @test fg.roundradius == 2
    @test fg.font == "inherit"
    @test fg.fontsize == 12.34
    @test fg.notext == false
    @test fg.timeunit == :none
    @test fg.delay == 0.001

    g = flamegraph(backtraces, lidict=lidict)

    fg = ProfileSVG.view(sfc, g, C=true, height=123.4, roundradius=0, timeunit=:ms,
                         unknown=true)
    @test FlameGraphs.depth(fg.g) == 4 # `C` option does not affect the graph
    @test fg.fcolor isa StackFrameCategory
    @test fg.graph_options[:C] == true
    @test fg.yflip == false
    @test fg.width == 960
    @test fg.height == 123.4
    @test fg.roundradius == 0
    @test fg.font == "inherit"
    @test fg.fontsize == 12
    @test fg.notext == false
    @test fg.timeunit == :ms
    @test fg.delay == 0.001

    fg = ProfileSVG.view(g, C=true, yflip=true, font="serif", notext=true, delay=0.01,
                         unknown=false)
    @test FlameGraphs.depth(fg.g) == 4 # `C` option does not affect the graph
    @test fg.fcolor isa FlameColors
    @test fg.graph_options[:C] == true
    @test fg.yflip == true
    @test fg.width == 960
    @test fg.height == 0
    @test fg.roundradius == 2
    @test fg.font == "serif"
    @test fg.fontsize == 12
    @test fg.notext == true
    @test fg.timeunit == :none
    @test fg.delay == 0.01
end

@testset "save" begin
    io = IOBuffer()
    sfc = StackFrameCategory()

    function svg_size(str)
        m = match(Regex("""<svg version="1.1" width="([^"]+)" height="([^"]+)" """), str)
        m === nothing && error("svg size is unknown")
        (m.captures[1], m.captures[2])
    end
    function has_filled_rect(str, color)
        occursin(Regex("""<rect[^/]+fill="$color" """), str)
    end
    function has_filled_path(str, color)
        occursin(Regex("""<path[^/]+fill="$color" """), str)
    end

    ProfileSVG.save(sfc, io, backtraces,
                    C=true, lidict=lidict, width=123.4, unknown=nothing)
    str = String(take!(io))
    @test svg_size(str) == ("123.4", "147")
    @test has_filled_rect(str, "#FF0000")
    @test !has_filled_path(str, "#FF0000")
    filename = tempname()
    ProfileSVG.save(sfc, filename, backtraces,
                    C=true, lidict=lidict, width=123.4, unknown=nothing)
    str = read(filename, String)
    rm(filename)
    @test svg_size(str) == ("123.4", "147")
    @test has_filled_rect(str, "#FF0000")
    @test !has_filled_path(str, "#FF0000")

    ProfileSVG.save(io, backtraces,
                    C=true, lidict=lidict, fontsize=12.34, unknown=missing)
    str = String(take!(io))
    @test svg_size(str) == ("960", "151")
    @test !has_filled_rect(str, "#FF0000")
    @test !has_filled_path(str, "#FF0000")
    filename = tempname()
    ProfileSVG.save(filename, backtraces,
                    C=true, lidict=lidict, fontsize=12.34, unknown=missing)
    str = read(filename, String)
    rm(filename)
    @test svg_size(str) == ("960", "151")
    @test !has_filled_rect(str, "#FF0000")
    @test !has_filled_path(str, "#FF0000")

    g = flamegraph(backtraces, lidict=lidict)

    ProfileSVG.save(sfc, io, g, C=true, height=123.4, roundradius=0, unknown=true)
    str = String(take!(io))
    @test svg_size(str) == ("960", "123.4")
    @test !has_filled_rect(str, "#FF0000")
    @test has_filled_path(str, "#FF0000")
    filename = tempname()
    ProfileSVG.save(sfc, filename, g, C=true, height=123.4, roundradius=0, unknown=true)
    str = read(filename, String)
    rm(filename)
    @test svg_size(str) == ("960", "123.4")
    @test !has_filled_rect(str, "#FF0000")
    @test has_filled_path(str, "#FF0000")

    ProfileSVG.save(io, g,
                    C=true, yflip=true, font="serif", notext=true, unknown=false)
    str = String(take!(io))
    @test svg_size(str) == ("960", "136")
    @test !has_filled_rect(str, "#FF0000")
    @test !has_filled_path(str, "#FF0000")
    filename = tempname()
    ProfileSVG.save(filename, g,
                    C=true, yflip=true, font="serif", notext=true, unknown=false)
    str = read(filename, String)
    rm(filename)
    @test svg_size(str) == ("960", "136")
    @test !has_filled_rect(str, "#FF0000")
    @test !has_filled_path(str, "#FF0000")
end

@testset "color schemes" begin
    g = flamegraph(backtraces, lidict=lidict)

    io = IOBuffer()
    fcolor = FlameColors(colorbg=colorant"steelblue", colorfont=colorant"lightyellow")

    modcat(mod) = nothing
    loccat(sf) = occursin("f1", string(sf.func)) ? colorant"navy" : colorant"powderblue"
    sfc = StackFrameCategory(modcat, loccat, colorant"dimgray", colorant"red")

    function bgcolor(str)
        m = match(r"#fig-\w+-bg {([^}]+)", str)
        m === nothing && return m
        m = match(r"fill: ([^;]+);\n", m[1])
        m === nothing && return m
        color1 = m[1]

        m = match(r"#fig-\w+-viewport \+ rect {([^}]+)", str)
        m === nothing && return m
        m = match(r"fill: ([^;]+);\n", m[1])
        m === nothing && return m
        color2 = m[1]
        return color1 == color2 ? color1 : nothing
    end

    function fontcolor(str)
        m = match(r"#fig-\w+ text {([^}]+)", str)
        m === nothing && return m
        m = match(r"fill: ([^;]+);\n", m[1])
        m === nothing && return m
        color1 = m[1]

        m = match(r"#fig-\w+-viewport text {([^}]+)", str)
        m === nothing && return m
        m = match(r"stroke: ([^;]+);\n", m[1])
        m === nothing && return m
        color2 = m[1]

        m = match(r"#fig-\w+-viewport path:hover {([^}]+)", str)
        m === nothing && return m
        m = match(r"stroke: ([^;]+);\n", m[1])
        m === nothing && return m
        color3 = m[1]
        return color1 == color2 == color3 ? color1 : nothing
    end

    function captioncolor(str)
        m = match(r"text#fig-\w+-caption {([^}]+)", str)
        m === nothing && return m
        m = match(r"fill: ([^;]+);\n", m[1])
        m === nothing && return m
        return m[1]
    end

    function opacity(str)
        m = match(r"#fig-\w+-viewport path {([^}]+)", str)
        m === nothing && return m
        m = match(r"fill-opacity: ([^;]+);\n", m[1])
        m === nothing && return m
        opacity1 = m[1]

        m = match(r"#fig-\w+-viewport path:hover {([^}]+)", str)
        m === nothing && return m
        m = match(r"fill-opacity: ([^;]+);\n", m[1])
        m === nothing && return m
        opacity2 = m[1]
        return (opacity1, opacity2)
    end

    ProfileSVG.save(io, g)
    str = String(take!(io))
    @test bgcolor(str) == "#FFFFFF"
    @test fontcolor(str) == "#000000"
    @test captioncolor(str) == "#000000"
    @test opacity(str) == ("1", "0.75")

    ProfileSVG.save(fcolor, io, g, frameopacity=0.4f0)
    str = String(take!(io))
    @test bgcolor(str) == "#4682B4"
    @test fontcolor(str) == "#FFFFE0"
    @test captioncolor(str) == "#FFFFE0"
    @test opacity(str) == ("0.4", "0.65")

    ProfileSVG.save(sfc, io, g, bgcolor=:classic, fontcolor=:classic, frameopacity=0.75)
    str = String(take!(io))
    @test occursin("<linearGradient", str)
    @test bgcolor(str) != "#FFFFFF"
    @test fontcolor(str) == "black"
    @test captioncolor(str) == "black"
    @test opacity(str) == ("0.75", "1")

    ProfileSVG.save(io, g, bgcolor=:transparent, fontcolor=:currentcolor, frameopacity=0.9)
    str = String(take!(io))
    @test bgcolor(str) == "transparent"
    @test fontcolor(str) == "currentcolor"
    @test opacity(str) == ("0.9", "0.65")

    ProfileSVG.save(sfc, io, g, bgcolor=:unknown, fontcolor=:bw)
    str = String(take!(io))
    @test bgcolor(str) == "white"
    @test fontcolor(str) == "black"
    @test captioncolor(str) == "black"
    @test opacity(str) == ("1", "0.75")
    @test count_element(r"<text [^/]+class=\"w\">", str) == 2

    ProfileSVG.save(fcolor, io, g, fontcolor=:bw)
    str = String(take!(io))
    @test bgcolor(str) == "#4682B4"
    @test fontcolor(str) == "black"
    @test captioncolor(str) == "white"
    @test opacity(str) == ("1", "0.75")
    @test count_element(r"<text [^/]+class=\"w\">", str) == 5

    @test_throws TypeError ProfileSVG.save(sfc, io, g, bgcolor=colorant"black")
end

@testset "size limitation" begin
    io = IOBuffer()

    @test_logs(
        (:warn, r"The depth of this graph is 5, exceeding the `maxdepth` \(=4\)"),
        (:warn, r"The maximum number of frames \(`maxframes`=8\) is reached"),
        ProfileSVG.save(io, backtraces, C=true, lidict=lidict, maxdepth=4, maxframes=8))
    str = String(take!(io))
    @test occursin("""height="136" """, str)
    @test count_element(r"<rect x=[^/]+/>", str) == 8
end

@testset "show as html" begin
    io = IOBuffer()
    fg = ProfileSVG.view(backtraces, lidict=lidict)
    show(io, "text/html", fg)
    str = String(take!(io))
    @test occursin(r"<html>.+(?=</html>)"s, str)
end

@testset "set_default" begin
    sfc = StackFrameCategory()
    Profile.init(delay=0.001)
    ProfileSVG.set_default(sfc, C=true, height=567.8)

    fgc = ProfileSVG.view(sfc, backtraces,
                          C=false, lidict=lidict, width=123.4, unknown=nothing)
    @test FlameGraphs.depth(fgc.g) == 4
    @test fgc.fcolor isa StackFrameCategory
    @test fgc.graph_options[:C] == false
    @test fgc.graph_options[:lidict] == lidict
    @test fgc.yflip == false
    @test fgc.width == 123.4
    @test fgc.height == 567.8
    @test fgc.roundradius == 2
    @test fgc.font == "inherit"
    @test fgc.fontsize == 12
    @test fgc.notext == false
    @test fgc.timeunit == :none
    @test fgc.delay == 0.001

    ProfileSVG.set_default(fontsize=9, delay=0.01)
    fgc = ProfileSVG.view(backtraces,
                          lidict=lidict, width=123.4, unknown=nothing)
    @test FlameGraphs.depth(fgc.g) == 5
    @test fgc.fcolor isa StackFrameCategory
    @test fgc.graph_options[:C] == true
    @test fgc.graph_options[:lidict] == lidict
    @test fgc.yflip == false
    @test fgc.width == 123.4
    @test fgc.height == 567.8
    @test fgc.roundradius == 2
    @test fgc.font == "inherit"
    @test fgc.fontsize == 9
    @test fgc.notext == false
    @test fgc.timeunit == :none
    @test fgc.delay == 0.01

    ProfileSVG.init()
    fgc = ProfileSVG.view(backtraces,
                          C=true, lidict=lidict, width=123.4, unknown=nothing)
    @test FlameGraphs.depth(fgc.g) == 5
    @test fgc.fcolor isa FlameColors
    @test fgc.graph_options[:C] == true
    @test fgc.graph_options[:lidict] == lidict
    @test fgc.yflip == false
    @test fgc.width == 123.4
    @test fgc.height == 0
    @test fgc.roundradius == 2
    @test fgc.font == "inherit"
    @test fgc.fontsize == 12
    @test fgc.notext == false
    @test fgc.timeunit == :none
    @test fgc.delay == 0.001
end

# For these tests to work you need `rsvg-convert` installed.
# On Ubuntu this is `sudo apt install librsvg2-bin`.
@testset "profview" begin
    Profile.init(n=10000) # prevent stack overflow (on 32-bit systems)
    profile_test(1)   # to compile
    @profview profile_test(10)

    mktemp() do path, io
        ProfileSVG.save(io)
        flush(io)
        try
            # Validate the file by converting to PNG
            str = read(`rsvg-convert $path`, String)
            @test codeunits(str)[1:8] == UInt8[0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a]
        catch e
            e isa Base.IOError || rethrow()
        end
    end
end


text:
# Module with optimization functions for Julia TimeModelling
# Author: Mathias Louboutin, mloubout@eos.ubc.ca
# Date: January, 2017
#

module SLIM_optim

using PolynomialRoots, Printf, Statistics, LinearAlgebra, Random

import JUDI

#############################################################################
# Optimization algorithms
include("SPGSlim.jl")   # minConf_SPG
include("PQNSlim.jl")   # minConf_PQN
include("OptimizationFunctions.jl") # common functions
include("Constraints.jl") # Constriaints projection
include("fwi_linesearch.jl")

end


text:
function rootfind_bisection(f::Function, val::Float64, bracket1::Vector{Float64}, bracket2::Vector{Float64})

    a = bracket1[:]
    b = bracket2[:]
    fa = f(a) - val
    fb = f(b) - val
    c = (a + b)/2
    # c = (a*fb - b*fa) / (fb - fa)
    fc = f(c) - val

    while abs(fc) > 1e-8
        if fa*fc < 0
            b = c
        elseif fc*fb < 0
            a = c
        else
            @error "root finding algorithm has failed."
            break
        end
        fa = f(a) - val
        fb = f(b) - val
        c = (a + b)/2
        # c = (a*fb - b*fa) / (fb - fa)
        fc = f(c) - val
    end

    return c
end

function rk4_step(r::Vector{Float64}, f::Function, ds::Float64)

    k1 = f(r)
    k2 = f(r + ds*k1/2)
    k3 = f(r + ds*k2/2)
    k4 = f(r + ds*k3)

    return r + ds*(k1 + 2*k2 + 2*k3 + k4) / 6
end

function conservative_step(r::Vector{Float64}, b::Function, psi::Function, ds::Float64)

    # compute predictor step using rk4
    rp = rk4_step(r, b, ds)

    # compute correction using root finding algorithm
    en = zeros(length(rp))
    en[1:2] = [0 -1; 1 0] * b(rp)[1:2] / norm(b(rp)[1:2])
    bracket1 = rp + ds*en
    bracket2 = rp - ds*en
    rc = rootfind_bisection(psi, psi(r), bracket1, bracket2)

    # compute with correct displacement
    return rc, norm(rc - r)
end

function trace_fieldline(stop_condition::Function, r0::Vector{Float64}, b::Function, psi::Function, ds::Float64)

    r = r0[:]
    deltaS = 0.0

    points = Vector{Float64}[r0]
    while !stop_condition(r)
        rnew, dsnew = conservative_step(r, b, psi, ds)
        r[:] = rnew
        deltaS += dsnew
        append!(points, [rnew])
    end

    return hcat(points...), deltaS
end

function trace_fieldline_dist(r0::Vector{Float64}, b::Function, ds::Float64, dlim::Float64)
    r = r0[:]
    dS = 0.0

    points = Vector{Float64}[r0]

    while dS < dlim
        rnew = rk4_step(r, b, ds)
        dS += norm(rnew - r)
        r[:] = rnew
        append!(points, [rnew])
    end

    return hcat(points...), dS
end


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Getting started with JuMP

# This tutorial is aimed at providing a quick introduction to writing JuMP code.

# ## What is JuMP?

# JuMP ("Julia for Mathematical Programming") is an open-source modeling
# language that is embedded in Julia. It allows users to formulate various
# classes of optimization problems (linear, mixed-integer, quadratic, conic
# quadratic, semidefinite, and nonlinear) with easy-to-read code. These problems
# can then be solved using state-of-the-art open-source and commercial solvers.

# JuMP also makes advanced optimization techniques easily accessible from a
# high-level language.

# ## Installation

# JuMP is a package for Julia. From Julia, JuMP is installed by using the
# built-in package manager.

# ```julia
# import Pkg
# Pkg.add("JuMP")
# ```

# You also need to include a Julia package which provides an appropriate solver.
# One such solver is `GLPK.Optimizer`, which is provided by the
# [GLPK.jl package](https://github.com/jump-dev/GLPK.jl).
# ```julia
# import Pkg
# Pkg.add("GLPK")
# ```
# See [Installation Guide](@ref) for a list of other solvers you can use.

# ## An example

# Let's try to solve the following linear programming problem by using JuMP and
# GLPK. We will first look at the complete code to solve the problem and then go
# through it step by step.

# ```math
# \begin{aligned}
# & \min & 12x + 20y \\
# & \;\;\text{s.t.} & 6x + 8y \geq 100 \\
# & & 7x + 12y \geq 120 \\
# & & x \geq 0 \\
# & & y \in [0, 3] \\
# \end{aligned}
# ```

using JuMP
using GLPK
model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
@objective(model, Min, 12x + 20y)
@constraint(model, c1, 6x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
print(model)
optimize!(model)
@show termination_status(model)
@show primal_status(model)
@show dual_status(model)
@show objective_value(model)
@show value(x)
@show value(y)
@show shadow_price(c1)
@show shadow_price(c2)

# ## Step-by-step

# Once JuMP is installed, to use JuMP in your programs, we just need to write:

using JuMP

# We also need to include a Julia package which provides an appropriate solver.
# We want to use `GLPK.Optimizer` here which is provided by the `GLPK.jl`
# package.

using GLPK

# A model object is a container for variables, constraints, solver options, etc.
# Models are created with the [`Model`](@ref) function. The model can be created
# with an optimizer attached with default arguments by calling the constructor
# with the optimizer type, as follows:

model = Model(GLPK.Optimizer)

# Variables are modeled using [`@variable`](@ref):

@variable(model, x >= 0)

# They can have lower and upper bounds.

@variable(model, 0 <= y <= 30)

# The objective is set using [`@objective`](@ref):

@objective(model, Min, 12x + 20y)

# Constraints are modeled using [`@constraint`](@ref). Here `c1` and `c2` are
# the names of our constraint.

@constraint(model, c1, 6x + 8y >= 100)

#-

@constraint(model, c2, 7x + 12y >= 120)

#- Call `print` to display the model:

print(model)

# To solve the optimization problem, call the `optimize!` function.

optimize!(model)

# !!! info
#     The `!` after optimize is just part of the name. It's nothing special.
#     Julia has a convention that functions which mutate their arguments should
#     end in `!`. A common example is `push!`.

# Now let's see what information we can query about the solution.

# [`termination_status`](@ref) tells us why the solver stopped:

termination_status(model)

# In this case, the solver found an optimal solution. We should also check
# [`primal_status`](@ref) to see if the solver found a primal feasible point:

primal_status(model)

# and [`dual_status`](@ref) to see if the solver found a dual feasible point:

dual_status(model)

# Now we know that our solver found an optimal solution, and has a primal and a
# dual solution to query.

# Query the objective value using [`objective_value`](@ref):

objective_value(model)

# The primal solution using [`value`](@ref):

value(x)

#-

value(y)

# and the dual solution using [`shadow_price`](@ref):

shadow_price(c1)

#-

shadow_price(c2)

# ## Variable basics

model = Model()

# ### Variable bounds

# All of the variables we have created till now have had a bound. We can also
# create a free variable.

@variable(model, free_x)

# While creating a variable, instead of using the <= and >= syntax, we can also
# use the `lower_bound` and `upper_bound` keyword arguments.

@variable(model, keyword_x, lower_bound = 1, upper_bound = 2)

# We can query whether a variable has a bound using the `has_lower_bound` and
# `has_upper_bound` functions. The values of the bound can be obtained using the
# `lower_bound` and `upper_bound` functions.

has_upper_bound(keyword_x)

#-

upper_bound(keyword_x)

# Note querying the value of a bound that does not exist will result in an error.

lower_bound(free_x)

# JuMP also allows us to change the bounds on variable. We will learn this in
# the problem modification tutorial.

# ### [Containers](@id tutorial_variable_container)

# We have already seen how to add a single variable to a model using the
# [`@variable`](@ref) macro. Let's now look at more ways to add variables to a
# JuMP model.

# JuMP provides data structures for adding collections of variables to a model.
# These data structures are referred to as Containers and are of three types:
# `Arrays`, `DenseAxisArrays`, and `SparseAxisArrays`.

# #### Arrays

# JuMP arrays are created in a similar syntax to Julia arrays with the addition
# of specifying that the indices start with 1. If we do not tell JuMP that the
# indices start at 1, it will create a `DenseAxisArray` instead.

@variable(model, a[1:2, 1:2])

# An n-dimensional variable $x \in {R}^n$ having a bound $l \preceq x \preceq u$
# ($l, u \in {R}^n$) is added in the following manner.

n = 10
l = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10]
u = [10; 11; 12; 13; 14; 15; 16; 17; 18; 19]

@variable(model, l[i] <= x[i = 1:n] <= u[i])

# Note that while working with Containers, we can also create variable bounds
# depending upon the indices:

@variable(model, y[i = 1:2, j = 1:2] >= 2i + j)

# #### DenseAxisArrays

# `DenseAxisArrays` are used when the required indices are not one-based integer
# ranges. The syntax is similar except with an arbitrary vector as an index as
# opposed to a one-based range.

# An example where the indices are integers but do not start with one.

@variable(model, z[i = 2:3, j = 1:2:3] >= 0)

# Another example where the indices are an arbitrary vector.

@variable(model, w[1:5, ["red", "blue"]] <= 1)

# #### SparseAxisArrays

# `SparseAxisArrays` are created when the indices do not form a rectangular set.
# For example, this applies when indices have a dependence upon previous indices
# (called triangular indexing).

@variable(model, u[i = 1:3, j = i:5])

# We can also conditionally create variables by adding a comparison check that
# depends upon the named indices and is separated from the indices by a
# semi-colon (;).

@variable(model, v[i = 1:9; mod(i, 3) == 0])

# ### Variable types

# The last argument to the `@variable` macro is usually the variable type. Here
# we'll look at how to specify the variable type.

# #### Integer variables

# Integer optimization variables are constrained to the set $x \in {Z}$

@variable(model, integer_x, Int)

# or

@variable(model, integer_z, integer = true)

# #### Binary variables

# Binary optimization variables are constrained to the set $x \in \{0, 1\}$.

@variable(model, binary_x, Bin)

# or

@variable(model, binary_z, binary = true)

# ## Constraint basics

model = Model()
@variable(model, x)
@variable(model, y)
@variable(model, z[1:10]);

# ### Constraint references

# While calling the `@constraint` macro, we can also set up a constraint
# reference. Such a reference is useful for obtaining additional information
# about the constraint, such as its dual solution.

@constraint(model, con, x <= 4)

# ### [Containers](@id tutorial_constraint_container)

# Just as we had containers for variables, JuMP also provides `Arrays`,
# `DenseAxisArrays`, and `SparseAxisArrays` for storing collections of
# constraints. Examples for each container type are given below.

# #### Arrays

@constraint(model, [i = 1:3], i * x <= i + 1)

# #### DenseAxisArrays

@constraint(model, [i = 1:2, j = 2:3], i * x <= j + 1)

# #### SparseAxisArrays

@constraint(model, [i = 1:2, j = 1:2; i != j], i * x <= j + 1)

# ### Constraints in a loop

# We can add constraints using regular Julia loops

for i in 1:3
    @constraint(model, 6x + 4y >= 5i)
end

# or use for each loops inside the `@constraint` macro.

@constraint(model, [i in 1:3], 6x + 4y >= 5i)

# We can also create constraints such as $\sum _{i = 1}^{10} z_i \leq 1$

@constraint(model, sum(z[i] for i in 1:10) <= 1)

# ## Objective functions

# While the recommended way to set the objective is with the [`@objective`](@ref)
# macro, the functions [`set_objective_sense`](@ref) and [`set_objective_function`](@ref)
# provide an equivalent lower-level interface.

using GLPK

model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
set_objective_sense(model, MOI.MIN_SENSE)
set_objective_function(model, x + y)

optimize!(model)

#-

objective_value(model)

# To query the objective function from a model, we use the [`objective_sense`](@ref),
# [`objective_function`](@ref), and [`objective_function_type`](@ref) functions.

objective_sense(model)

#-

objective_function(model)

#-

objective_function_type(model)

# ## Vectorized syntax

# We can also add constraints and an objective to JuMP using vectorized linear
# algebra. We'll illustrate this by solving an LP in standard form i.e.

# ```math
# \begin{aligned}
# & \min & c^T x \\
# & \;\;\text{s.t.} & A x = b \\
# & & x \succeq 0 \\
# & & x \in \mathbb{R}^n
# \end{aligned}
# ```

vector_model = Model(GLPK.Optimizer)

A = [
    1 1 9 5
    3 5 0 8
    2 0 6 13
]

b = [7; 3; 5]

c = [1; 3; 5; 2]

@variable(vector_model, x[1:4] >= 0)
@constraint(vector_model, A * x .== b)
@objective(vector_model, Min, c' * x)

optimize!(vector_model)

#-

objective_value(vector_model)

max_stars_repo_name: odow/jump-training-materials

max_issues_repo_name: odow/jump-training-materials

max_forks_repo_name: odow/jump-training-materials


text:
"""
    set_param!(gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol, value)

Set the value of a trainable parameter of the generative function.

NOTE: Does not update the gradient accumulator value.
"""
function set_param!(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol, value)
    gf.params[name] = value
end

"""
    value = get_param(gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol)

Get the current value of a trainable parameter of the generative function.
"""
function get_param(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol)
    gf.params[name]
end

"""
    value = get_param_grad(gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol)

Get the current value of the gradient accumulator for a trainable parameter of the generative function.
"""
function get_param_grad(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol)
    gf.params_grad[name]
end

"""
    zero_param_grad!(gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol)

Reset the gradient accumlator for a trainable parameter of the generative function to all zeros.
"""
function zero_param_grad!(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol)
    gf.params_grad[name] = zero(gf.params[name])
end

"""
    set_param_grad!(gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol, grad_value)

Set the gradient accumlator for a trainable parameter of the generative function.
"""
function set_param_grad!(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol, grad_value)
    gf.params_grad[name] = grad_value
end

"""
    init_param!(gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol, value)

Initialize the the value of a named trainable parameter of a generative function.

Also generates the gradient accumulator for that parameter to `zero(value)`.

Example:
```julia
init_param!(foo, :theta, 0.6)
```
"""
function init_param!(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, name::Symbol, value)
    set_param!(gf, name, value)
    zero_param_grad!(gf, name)
end

get_params(gf::Union{DynamicDSLFunction,StaticIRGenerativeFunction}) = keys(gf.params)

export set_param!, get_param, get_param_grad, zero_param_grad!, set_param_grad!, init_param!

#########################################
# gradient descent with fixed step size #
#########################################

mutable struct FixedStepGradientDescentBuiltinDSLState
    step_size::Float64
    gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}
    param_list::Vector
end

function init_update_state(conf::FixedStepGradientDescent,
        gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, param_list::Vector)
    FixedStepGradientDescentBuiltinDSLState(conf.step_size, gen_fn, param_list)
end

function apply_update!(state::FixedStepGradientDescentBuiltinDSLState)
    for param_name in state.param_list
        value = get_param(state.gen_fn, param_name)
        grad = get_param_grad(state.gen_fn, param_name)
        set_param!(state.gen_fn, param_name, value + grad * state.step_size)
        zero_param_grad!(state.gen_fn, param_name)
    end
end

####################
# gradient descent #
####################

mutable struct GradientDescentBuiltinDSLState
    step_size_init::Float64
    step_size_beta::Float64
    gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}
    param_list::Vector
    t::Int
end

function init_update_state(conf::GradientDescent,
        gen_fn::Union{DynamicDSLFunction,StaticIRGenerativeFunction}, param_list::Vector)
    GradientDescentBuiltinDSLState(conf.step_size_init, conf.step_size_beta,
        gen_fn, param_list, 1)
end

function apply_update!(state::GradientDescentBuiltinDSLState)
    step_size = state.step_size_init * (state.step_size_beta + 1) / (state.step_size_beta + state.t)
    for param_name in state.param_list
        value = get_param(state.gen_fn, param_name)
        grad = get_param_grad(state.gen_fn, param_name)
        set_param!(state.gen_fn, param_name, value + grad * step_size)
        zero_param_grad!(state.gen_fn, param_name)
    end
    state.t += 1
end


text:
function del_flagged!(S::GphysData, dflag::BitArray{1}, reason::String)
  d = findall(dflag)
  L = length(d)
  if L > 0
    @warn(string("Deleting (", reason, ")"), S.id[d])
    deleteat!(S, d)
  end
  return nothing
end


text:
# The functions in this module are a direct translation of the original
# Matlab functions by E. Doland and J. J. Moré.
# See http://www.mcs.anl.gov/~more/cops.
#
# D. Orban, 2015, 2016.

"""Compute performance ratios used to produce a performance profile.

There is normally no need to call this function directly.
See the documentation of `performance_profile()` for more information.
"""
function performance_ratios(T :: Array{Float64,2}; logscale :: Bool=true)

  (np, ns) = size(T);       # Number of problems and number of solvers.

  T[isinf.(T)] = NaN;
  T[T .< 0] = NaN;
  minperf = mapslices(NaNMath.minimum, T, 2) # Minimal (i.e., best) performance per solver

  # Compute ratios and divide by smallest element in each row.
  r = zeros(np, ns);
  for p = 1 : np
    r[p, :] = T[p, :] / minperf[p];
  end

  logscale && (r = log2.(r));
  max_ratio = NaNMath.maximum(r)

  # Replace failures with twice the max_ratio and sort each column of r.
  failures = isnan.(r);
  r[failures] = 2 * max_ratio;
  r = sort(r, 1);
  return (r, max_ratio)
end

performance_ratios{Td <: Number}(T :: Array{Td,2}; logscale :: Bool=true) = performance_ratios(convert(Array{Float64,2}, T), logscale=logscale)

"""Produce a performance profile.

Each column of the matrix `T` defines the performance data for a solver
(smaller is better). Failures on a given problem are represented by a
negative value, an infinite value, or `NaN`.
The optional argument `logscale` is used to produce a logarithmic (base 2)
performance plot.
"""
function performance_profile(T :: Array{Float64,2}, labels :: Vector{AbstractString};
                             logscale :: Bool=true,
                             title :: AbstractString="",
                             kwargs...)

  (ratios, max_ratio) = performance_ratios(T, logscale=logscale)
  (np, ns) = size(ratios)

  ratios = [ratios; 2.0 * max_ratio * ones(1, ns)]
  xs = [1:np+1;] / np
  length(labels) == 0 && (labels = [@sprintf("column %d", col) for col = 1 : ns])
  profile = Plots.plot(; kwargs...)  # initial empty plot
  for s = 1 : ns
    Plots.plot!(ratios[:, s], xs, t=:steppost, label=labels[s])
  end
  Plots.xlims!(logscale ? 0.0 : 1.0, 1.1 * max_ratio)
  Plots.ylims!(0, 1.1)
  Plots.xlabel!("Within this factor of the best" * (logscale ? " (log scale)" : ""))
  Plots.ylabel!("Proportion of problems")
  Plots.title!(title)
  return profile
end

performance_profile{Tn <: Number, S <: AbstractString}(T :: Array{Tn,2}, labels :: Vector{S}; kwargs...) =
  performance_profile(convert(Array{Float64,2}, T), convert(Vector{AbstractString}, labels); kwargs...)


text:
using CeterisParibus
using Test

@testset "CeterisParibus.jl" begin
    include("f().jl")
    include("avg.jl")
    include("g().jl")
    include("Wavg().jl")
end


text:
include("load.jl")
StopWatchesTests.runtests()