text:
module Delaunay
using Random
using LightGraphs
using Plots
using LinearAlgebra
using DataStructures

function is_closer(points, p, a::Int, b::Int)
    (points[p, :] .- points[a, :]) .^ 2 < (points[p, :] .- points[b, :]) .^ 2
end

function is_aligned(points, a::Int, b::Int, c::Int)
    v1 = points[b,:] .- points[a,:]
    v2 = points[c,:] .- points[a,:]
    v1[1] * v2[2] == v1[2] * v2[1]
end

"True if c is left of a -> b where a b c refere to points in points"
function is_left_of(points, a::Int, b::Int, c::Int)
    start_to_arrival = points[b,:] - points[a, :]
    norm = [-start_to_arrival[2] start_to_arrival[1]]
    (norm * (points[c,:]-points[b,:]))[1] >= 0
end

function is_strictly_left_of(points, a::Int, b::Int, c::Int)
    start_to_arrival = points[b,:] - points[a, :]
    norm = [-start_to_arrival[2] start_to_arrival[1]]
    (norm * (points[c,:]-points[b,:]))[1] > 0
end

"True if c is right of a -> b where a b c refere to points in points"
function is_right_of(points, a::Int, b::Int, c::Int)
    is_left_of(points, b, a, c)
end

function is_strictly_right_of(points, a::Int, b::Int, c::Int)
    is_strictly_left_of(points, b, a, c)
end

"True if p is in circumcircle a b c. (a,b,c ccw)"
function is_in_circumcircle(points, a::Int, b::Int, c::Int, p::Int)
    det(
        [
            points[a,1] points[a,2] (points[a,1]^2+points[a,2]^2) 1
            points[b,1] points[b,2] (points[b,1]^2+points[b,2]^2) 1
            points[c,1] points[c,2] (points[c,1]^2+points[c,2]^2) 1
            points[p,1] points[p,2] (points[p,1]^2+points[p,2]^2) 1
        ]
    ) > 0
end

mutable struct Triangle
    v1::Int64
    v2::Int64
    v3::Int64
    t1::Triangle
    t2::Triangle
    t3::Triangle
    function Triangle(a::Int, b::Int)
        t = new(a, b, 0)
        next = new(b, a, 0)
        t.t1 = t.t2 = t.t3 = next
        next.t1 = next.t2 = next.t3 = t
        println("Creating $t")
        t
    end
    function Triangle(v1::Int, v2::Int, v3::Int, t1::Triangle, t2::Triangle, t3::Triangle)
        println("Creating Triangle($v1, $v2, $v3, $t1, $t2, $t3)")
        new(v1, v2, v3, t1, t2, t3)
    end
end

function Triangle(points, a::Int, b::Int, c::Int)
    v1, v2, v3 = sort([a b c], dims=1)
    if is_left_of(points, v1, v2, v3)
        v2, v3 = v3, v2
    end
    t1 = Triangle(v2, v3) # We want to order the point of t1 such as t1.v1 -> t1.v2 is clock wise for t
    t2 = Triangle(v3, v1)
    t3 = Triangle(v1, v2)
    t = Triangle(v1, v2, v3, t1, t2, t3)

    # since ti.v1 -> ti.v2 is always clockwise, ti.t3 is always t
    t1.t3 = t2.t3 = t3.t3 = t
    t1.t1 = t3.t2 = t2
    t1.t2 = t2.t1 = t3
    t2.t2 = t3.t1 = t1
    t
end

Base.show(io::IO, t::Triangle) = print(io, "Triangle($(t.v1), $(t.v2), $(t.v3))")

function is_hull(triangle)
    triangle.v3 == 0
end

"Returns the next triangle on the hull, clock wise"
function next(triangle)
    triangle.t1
end

"Set the next triangle on the hull, clock wise"
function next!(prev, next)
    if !is_hull(prev)
        error("$prev is not on hull.")
    end
    if !is_hull(next)
        error("$next is not on hull.")
    end
    println("Setting $prev => $next")
    prev.t1 = next
    next.t2 = prev
end

"Returns the previous triangle on the hull, clock wise"
function prev(triangle)
    triangle.t2
end

"Set the previous triangle on the hull, clock wise"
function prev!(next, prev)
    next!(prev, next)
end

function replace_comp!(triangle, neighbour, comp)
    if triangle == neighbour.t1
        println("$(neighbour).t1 = $comp")
        neighbour.t1 = comp
    elseif triangle == neighbour.t2
        println("$(neighbour).t2 = $comp")
        neighbour.t2 = comp
    else
        println("$(neighbour).t3 = $comp")
        neighbour.t3 = comp
    end
    println("$comp.t3 = $neighbour")
    comp.t3 = neighbour
end

replace_comp_t1!(triangle, comp) = replace_comp!(triangle, triangle.t1, comp)
replace_comp_t2!(triangle, comp) = replace_comp!(triangle, triangle.t2, comp)
replace_comp_t3!(triangle, comp) = replace_comp!(triangle, triangle.t3, comp)

"Delete a right set of triangle. `to_be_deleted` must be in the order of which
candidate choosing algorithm crosses the triangles.
"
function delete_right_triangles!(graph, to_be_deleted, upper_triangle)
    next_right = next(upper_triangle)
    next_next_right = next(next_right)
    println("Right deleting $to_be_deleted")
    for triangle in to_be_deleted
        if triangle == next_right # first iteration of the loop
            continue
        elseif is_hull(triangle)
            continue
        elseif triangle.t2 == next_right
            println("triangle = $triangle, next_right = $next_right, next_next_right = $next_next_right")
            println("triangle.t1 = $(triangle.t1), triangle.t2 = $(triangle.t2), triangle.t3 = $(triangle.t3)")
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            comp_t3 = Triangle(triangle.v2, triangle.v1)
            next!(upper_triangle, comp_t1)
            next!(comp_t1, comp_t3)
            next!(comp_t3, next_next_right)
            replace_comp_t1!(triangle, comp_t1)
            replace_comp_t3!(triangle, comp_t3)
            next_right = comp_t1
            next_next_right = comp_t3
            rem_edge!(graph, triangle.v1, triangle.v3)
        elseif triangle.t3 == next_right
            println("triangle = $triangle, next_right = $next_right, next_next_right = $next_next_right")
            println("triangle.t1 = $(triangle.t1), triangle.t2 = $(triangle.t2), triangle.t3 = $(triangle.t3)")
            comp_t2 = Triangle(triangle.v1, triangle.v3)
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            next!(upper_triangle, comp_t2)
            next!(comp_t2, comp_t1)
            next!(comp_t1, next_next_right)
            replace_comp_t1!(triangle, comp_t1)
            replace_comp_t2!(triangle, comp_t2)
            next_right = comp_t2
            next_next_right = comp_t1
            rem_edge!(graph, triangle.v1, triangle.v2)
        else
            println("triangle = $triangle, next_right = $next_right, next_next_right = $next_next_right")
            println("triangle.t1 = $(triangle.t1), triangle.t2 = $(triangle.t2), triangle.t3 = $(triangle.t3)")
            error("I am removing t1 on a right triangle.")
        end
    end
end

"Delete a left set of triangle. `to_be_deleted` must be in the order of which
candidate choosing algorithm crosses the triangles.
"
function delete_left_triangles!(graph, to_be_deleted, upper_triangle)
    prev_left = prev(upper_triangle)
    prev_prev_left = prev(prev_left)
    println("Left deleting $to_be_deleted")
    for triangle in to_be_deleted
        if triangle == prev_left # first iteration of the loop
            continue
        elseif is_hull(triangle)
            continue
        elseif triangle.t1 == prev_left
            comp_t3 = Triangle(triangle.v2, triangle.v1)
            comp_t2 = Triangle(triangle.v1, triangle.v3)
            replace_comp_t2!(triangle, comp_t2)
            replace_comp_t3!(triangle, comp_t3)
            prev!(upper_triangle, comp_t2)
            prev!(comp_t2, comp_t3)
            prev!(comp_t3, prev_prev_left)
            prev_left = comp_t2
            prev_prev_left = comp_t3
            rem_edge!(graph, triangle.v2, triangle.v3)
        elseif triangle.t2 == prev_left
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            comp_t3 = Triangle(triangle.v2, triangle.v1)
            replace_comp_t1!(triangle, comp_t1)
            replace_comp_t3!(triangle, comp_t3)
            prev!(upper_triangle, comp_t3)
            prev!(comp_t3, comp_t1)
            prev!(comp_t1, prev_prev_left)
            prev_left = comp_t3
            prev_prev_left = comp_t1
            rem_edge!(graph, triangle.v1, triangle.v3)
        else # triangle.t3 == next_right
            comp_t2 = Triangle(triangle.v1, triangle.v3)
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            replace_comp_t1!(triangle, comp_t1)
            replace_comp_t2!(triangle, comp_t2)
            prev!(upper_triangle, comp_t1)
            prev!(comp_t1, comp_t2)
            prev!(comp_t2, prev_prev_left)
            prev_left = comp_t1
            prev_prev_left = comp_t2
            rem_edge!(graph, triangle.v1, triangle.v2)
        end
    end
end

struct Map{T<:Real}
    delaunay::SimpleGraph
    points::Array{T,2}
end

function find_base_tangente(points, triangle_left, triangle_right)
    prev_triangle_left = prev(triangle_left)
    prev_triangle_right = prev(triangle_right)
    there_is_a_point_below_left_next = is_strictly_right_of(
        points, triangle_left.v1, triangle_right.v1, triangle_left.v2
    )
    there_is_a_point_below_left_prev = is_strictly_right_of(
        points, triangle_left.v1, triangle_right.v1, prev_triangle_left.v1
    )
    there_is_a_point_below_right_next = is_strictly_right_of(
        points, triangle_left.v1, triangle_right.v1, triangle_right.v2
    )
    there_is_a_point_below_right_prev = is_strictly_right_of(
        points, triangle_left.v1, triangle_right.v1, prev_triangle_right.v1
    )
    while there_is_a_point_below_left_next || there_is_a_point_below_left_prev || there_is_a_point_below_right_next || there_is_a_point_below_right_prev
        if there_is_a_point_below_left_next
            triangle_left, prev_triangle_left = next(triangle_left), triangle_left
        elseif there_is_a_point_below_left_prev
            triangle_left = prev_triangle_left
            prev_triangle_left = prev(triangle_left)

        end
        if there_is_a_point_below_right_next
            triangle_right, prev_triangle_right = next(triangle_right), triangle_right
        elseif there_is_a_point_below_right_prev
            triangle_right = prev_triangle_right
            prev_triangle_right = prev(triangle_right)
        end
        there_is_a_point_below_left_next = is_strictly_right_of(
            points, triangle_left.v1, triangle_right.v1, triangle_left.v2
        )
        there_is_a_point_below_left_prev = is_strictly_right_of(
            points, triangle_left.v1, triangle_right.v1, prev_triangle_left.v1
        )
        there_is_a_point_below_right_next = is_strictly_right_of(
            points, triangle_left.v1, triangle_right.v1, triangle_right.v2
        )
        there_is_a_point_below_right_prev = is_strictly_right_of(
            points, triangle_left.v1, triangle_right.v1, prev_triangle_right.v1
        )
    end

    # Corner case : aligned points
    right_next_aligned = is_aligned(points, triangle_left.v1, triangle_right.v1, triangle_right.v2)
    while right_next_aligned && is_closer(points, triangle_left.v1, triangle_right.v2, triangle_right.v1)
        prev_triangle_right, triangle_right = triangle_right, next(triangle_right)
        right_next_aligned = is_aligned(points, triangle_left.v1, triangle_right.v1, triangle_right.v2)
    end
    right_prev_aligned = is_aligned(points, triangle_left.v1, triangle_right.v1, prev_triangle_right.v1)
    while right_prev_aligned && is_closer(points, triangle_left.v1, prev_triangle_right.v1, triangle_right.v1)
        triangle_right = prev_triangle_right
        prev_triangle_right = prev(triangle_right)
        right_prev_aligned = is_aligned(points, triangle_left.v1, triangle_right.v1, prev_triangle_right.v1)
    end
    left_next_aligned = is_aligned(points, triangle_right.v1, triangle_left.v1, triangle_left.v2)
    while left_next_aligned && is_closer(points, triangle_right.v1, triangle_left.v2, triangle_left.v1)
        prev_triangle_left, triangle_left = triangle_left, next(triangle_left)
        left_next_aligned = is_aligned(points, triangle_right.v1, triangle_left.v1, triangle_left.v2)
    end
    left_prev_aligned = is_aligned(points, triangle_right.v1, triangle_left.v1, prev_triangle_left.v1)
    while left_prev_aligned && is_closer(points, triangle_right.v1, prev_triangle_left.v1, triangle_left.v1)
        triangle_left = prev_triangle_left
        prev_triangle_left = prev(triangle_left)
        left_prev_aligned = is_aligned(points, triangle_right.v1, triangle_left.v1, prev_triangle_left.v1)
    end
    triangle_left, triangle_right
end

function _in_circumcircle_left(points, triangle, base_left, base_right)
    if is_hull(triangle)
        false
    # We need to find out which vertex among v1 and v3 is base_left
    elseif triangle.v1 == base_left
        is_in_circumcircle(points, base_left, base_right, triangle.v3, triangle.v2)
    elseif triangle.v3 == base_left
        is_in_circumcircle(points, base_left, base_right, triangle.v2, triangle.v1)
    else # i.e. triangle.v2 == base_left
        is_in_circumcircle(points, base_left, base_right, triangle.v1, triangle.v3)
    end
end

function _in_circumcircle_right(points, triangle, base_left, base_right)
    if is_hull(triangle)
        false
    # We need to find out which vertex among v2 and v3 is base_right
    elseif triangle.v2 == base_right
        is_in_circumcircle(points, base_left, base_right, triangle.v3, triangle.v1)
    elseif triangle.v3 == base_right
        is_in_circumcircle(points, base_left, base_right, triangle.v1, triangle.v2)
    else # i.e triangle.v1 == base_right
        is_in_circumcircle(points, base_left, base_right, triangle.v2, triangle.v3)
    end
end

function _is_above_base_left(triangle, points, base_left, base_right)
    if is_hull(triangle) && triangle.v1 == base_left
        is_left_of(points, base_left, base_right, triangle.v2)
    elseif is_hull(triangle) || triangle.v3 == base_left
        is_left_of(points, base_left, base_right, triangle.v1)
    elseif triangle.v2 == base_left
        is_left_of(points, base_left, base_right, triangle.v3)
    else # i.e. v1 is base_left
        is_left_of(points, base_left, base_right, triangle.v2)
    end
end

function _is_above_base_right(triangle, points, base_left, base_right)
    if is_hull(triangle) && triangle.v2 == base_right
        is_left_of(points, base_left, base_right, triangle.v1)
    elseif is_hull(triangle) || triangle.v3 == base_right
        is_left_of(points, base_left, base_right, triangle.v2)
    elseif triangle.v1 == base_right
        is_left_of(points, base_left, base_right, triangle.v3)
    else # i.e. v2 is base_right
        is_left_of(points, base_left, base_right, triangle.v1)
    end
end

function _next_left(triangle, base_left)
    if is_hull(triangle) || triangle.v1 == base_left
        triangle.t3
    elseif triangle.v2 == base_left
        triangle.t1
    else # i.e. v3 is base_left
        triangle.t2
    end
end

function _next_right(triangle, base_right)
    if is_hull(triangle) || triangle.v2 == base_right
        triangle.t3
    elseif triangle.v3 == base_right
        triangle.t1
    else # i.e. v1 is base_right
        triangle.t2
    end
end

@enum Mode left=1 right=2
function find_candidate(mode::Mode, m, first_triangle, base_left, base_right)
    _in_circumcircle = if mode == left
        triangle -> _in_circumcircle_left(m.points, triangle, base_left, base_right)
    else
        triangle -> _in_circumcircle_right(m.points, triangle, base_left, base_right)
    end
    _is_above_base = if mode == left
        triangle -> _is_above_base_left(triangle, m.points, base_left, base_right)
    else
        triangle -> _is_above_base_right(triangle, m.points, base_left, base_right)
    end
    _next = if mode == left
        triangle -> _next_left(triangle, base_left)
    else
        triangle -> _next_right(triangle, base_right)
    end

    old_candidates = Deque{Triangle}()
    triangle = first_triangle
    found = false
    above_base = _is_above_base(triangle)
    while !(triangle in old_candidates) && !found && above_base
        next_triangle = _next(triangle)
        push!(old_candidates, triangle)
        if above_base && !_in_circumcircle(next_triangle)
            found=true
        else
            triangle = next_triangle
            above_base = _is_above_base(triangle)
        end
    end
    (found ? triangle : nothing), old_candidates
end

function _choose_candidate(points, base_left, base_right, left_candidate, right_candidate)
    b, c = base_left, base_right
    a = d = 0
    if is_hull(left_candidate) && left_candidate.v1 == b
        a = left_candidate.v2
    elseif left_candidate.v2 == b # also matches the cas of triangle on hull
        a = left_candidate.v1
    elseif left_candidate.v1 == b
        a = left_candidate.v3
    else # left_candidate.v3 == b
        a = left_candidate.v2
    end
    if is_hull(right_candidate) && right_candidate.v2 == c
        d = right_candidate.v1
    elseif right_candidate.v1 == c # also matches the cas of triangle on hull
        d = right_candidate.v2
    elseif right_candidate.v2 == c
        d = right_candidate.v3
    else # right_candidate.v3 == c
        d = right_candidate.v1
    end
    if is_in_circumcircle(points, a, b, c, d) # if right candidate in circumcircle(left_candidate, base_left, base_right)
        right_candidate
    else
        left_candidate
    end
end

function triangulate!(m, first=1, last=-1)
    if last == -1
        last = div(length(m.points), 2)
    end
    if first == last # never happens
        error("Trying to triangulate a single point $first.")
    elseif last - first == 1
        println("2 edges : $first and $last")
        add_edge!(m.delaunay, first, last)
        t = Triangle(first, last)
        t.t3, t
    elseif last - first == 2
        println("3 edges : $first, $(first+1) and $last.")
        if !is_aligned(m.points, first, first+1, last)
            add_edge!(m.delaunay, first, first+1)
            add_edge!(m.delaunay, first+1, last)
            add_edge!(m.delaunay, last, first)
            t = Triangle(m.points, first, first+1, last)
            if t.v2 < t.v3
                t.t3, t.t2
            else
                t.t3, t.t1
            end
        else
            add_edge!(m.delaunay, first, first+1)
            add_edge!(m.delaunay, first+1, last)
            t_first = Triangle(first, first+1)
            t_last = Triangle(first+1, last)
            next!(t_first, t_last)
            next!(t_last, t_last.t3)
            next!(t_last.t3, t_first.t3)
            next!(t_first.t3, t_first)
            if m.points[first, 1] < m.points[first+1, 1] # vertical line
                t_first, t_last.t3
            else # horizontal line
                t_first.t3, t_last
            end
        end

    else
        mid = div(first + last, 2)
        triangle_left_left, triangle_left_right = triangulate!(m, first, mid)
        triangle_right_left, triangle_right_right = triangulate!(m, mid+1, last)

        println("=================================================")
        println("triangle_left_left, triangle_left_right = $triangle_left_left, $triangle_left_right")
        println("triangle_right_left, triangle_right_right = $triangle_right_left, $triangle_right_right")
        println(collect(edges(m.delaunay)))
        base_triangle_left, base_triangle_right = find_base_tangente(m.points, triangle_left_right, triangle_right_left)
        println("$(last + 1 - first) edges : $first to $last.")
        base_left = base_triangle_left.v1
        base_right = base_triangle_right.v1
        println("Base edge is $base_left - $base_right.")

        triangle_left = prev(base_triangle_left)
        triangle_right = base_triangle_right

        # Reminder : for 2edges triangles, we force v1 -> v2 to be clockwise
        base_triangle = Triangle(base_right, base_left)
        next!(base_triangle, base_triangle_left)
        prev!(base_triangle, prev(base_triangle_right))

        upper_triangle = base_triangle.t3
        next!(upper_triangle, triangle_right)
        prev!(upper_triangle, triangle_left)

        right_candidate = next(upper_triangle)
        left_candidate = prev(upper_triangle)

        base_triangle_edge = :t3
        println("Base edge is $base_left - $base_right.")

        add_edge!(m.delaunay, base_left, base_right)

        while !isnothing(right_candidate) || !isnothing(left_candidate)
            println("Time to try to add a new edge with base triangle $base_triangle and base $base_left - $base_right")
            println("triangle_left = $triangle_left")
            println("triangle_right = $triangle_right")
            right_candidate, delete_triangle_right = find_candidate(right, m, triangle_right, base_left, base_right)
            left_candidate, delete_triangle_left = find_candidate(left, m, triangle_left, base_left, base_right)

            chosen = :none

            to_be_deleted = isnothing(right_candidate) ? delete_triangle_left : delete_triangle_right
            triangle = isnothing(right_candidate) ? left_candidate : right_candidate

            if !isnothing(left_candidate) && !isnothing(right_candidate)
                triangle = _choose_candidate(m.points, base_left, base_right, left_candidate, right_candidate)
                if triangle == left_candidate
                    to_be_deleted = delete_triangle_left
                end
            end
            if isnothing(triangle)
                chosen = :none
            elseif triangle == left_candidate
                chosen = :left
            elseif triangle == right_candidate
                chosen = :right
            end

            if chosen != :none

                println("Chose $chosen $triangle")

                new_base_left, new_base_right = base_left, base_right
                new_base_triangle_edge = :t1
                new_vertex = 0

                if chosen == :left
                    delete_left_triangles!(m.delaunay, to_be_deleted, upper_triangle)
                    triangle = prev(upper_triangle)
                    # Now we are sure that `triangle` is on hull.
                    new_vertex = triangle.v1
                else
                    delete_right_triangles!(m.delaunay, to_be_deleted, upper_triangle)
                    triangle = next(upper_triangle)
                    # Now we are sure that `triangle` is on hull.
                    new_vertex = triangle.v2
                end

                # Let's order the vertices of this new triangle
                triangle.v1, triangle.v2, triangle.v3 = sort([new_vertex, base_left, base_right])
                if is_left_of(m.points, triangle.v1, triangle.v2, triangle.v3)
                    triangle.v2, triangle.v3 = triangle.v3, triangle.v2
                end
                println("new_vertex is $new_vertex, base is $base_left, $base_right. New triangle is : $triangle.")

                # time to ensure the neighbors triangle are the right ones.

                upper_triangle.t3 = triangle
                if chosen == :left
                    # `prev(triangle)` is still valid since we have only
                    # modified its vertices for the moment.
                    new_left = prev(triangle)
                    if base_left == triangle.v1
                        upper_triangle.v1, upper_triangle.v2 = triangle.v2, triangle.v3
                        triangle.t1, triangle.t2 = upper_triangle, base_triangle
                        new_base_triangle_edge = :t1
                        new_base_left = triangle.v2
                    elseif base_left == triangle.v3
                        println("base_left == triangle.v3")
                        println("$triangle")
                        upper_triangle.v1, upper_triangle.v2 = triangle.v1, triangle.v2
                        triangle.t1, triangle.t2, triangle.t3 = base_triangle, triangle.t3, upper_triangle
                        new_base_triangle_edge = :t3
                        new_base_left = triangle.v1
                    end
                    prev!(upper_triangle, new_left)
                else
                    upper_triangle.v1, upper_triangle.v2 = triangle.v1, triangle.v2
                    next!(upper_triangle, next(triangle))
                    triangle.t1, triangle.t2, triangle.t3 = triangle.t3, base_triangle, upper_triangle
                    new_base_triangle_edge = :t3
                    new_base_right = triangle.v2
                end

                # Link base_triangle
                if base_triangle_edge == :t1
                    base_triangle.t1 = triangle
                elseif base_triangle_edge == :t2
                    base_triangle.t2 = triangle
                else
                    base_triangle.t3 = triangle
                end

                # We can finally add the new edge :D
                add_edge!(m.delaunay, new_base_left, new_base_right)

                # set up the new value for loop variables
                triangle_left = prev(upper_triangle)
                triangle_right = next(upper_triangle)
                base_triangle = triangle
                base_left, base_right = new_base_left, new_base_right
                base_triangle_edge = new_base_triangle_edge

                println("triangle = $triangle, $triangle.t1 = $(triangle.t1), $triangle.t2 = $(triangle.t2), $triangle.t3 = $(triangle.t3)")
            end
        end

        while !is_hull(triangle_left_left)
            if is_hull(triangle_left_left.t1)
                triangle_left_left = triangle_left_left.t1
            elseif is_hull(triangle_left_left.t2)
                    triangle_left_left = triangle_left_left.t2
            elseif is_hull(triangle_left_left.t3)
                triangle_left_left = triangle_left_left.t3
            else
                triangle_left_left = triangle_left_left.t1
            end
        end
        while !is_hull(triangle_right_right)
            if is_hull(triangle_right_right.t1)
                triangle_right_right = triangle_right_right.t1
            elseif is_hull(triangle_right_right.t2)
                triangle_right_right = triangle_right_right.t2
            elseif is_hull(triangle_right_right.t3)
                triangle_right_right = triangle_right_right.t3
            else
                triangle_right_right = triangle_right_right.t2
            end
        end
        triangle_left_left, triangle_right_right
    end

end

function Map(points::Array{T,2}) where T <: Real
    m = Map(SimpleGraph(div(length(points),2)), sortslices(points, dims=1))
    left,right = triangulate!(m)
    println("Got $left, $right")
    m
end

function Map(height::Number, width::Number, nb_points::Int; res::Float64=0.5)
    points = sortslices([rand(nb_points) .* width rand(nb_points) .* height], dims=1)
    println("I'm going to work on $points.")
    delaunay = SimpleGraph(nb_points)
    m = Map(delaunay, points)
    triangulate!(m)
    m
end

function plot_points(points)
    plot(
        points[:,1], points[:,2],
        marker=(20, 0.2, :orange), line=:scatter,
        series_annotations=collect(map(string, 1:length(points[:,1])))
    )
end

function plot_points!(points)
    plot!(
        points[:,1], points[:,2],
        marker=(20, 0.2, :orange), line=:scatter,
        series_annotations=collect(map(string, 1:length(points[:,1])))
    )
end

function plot_map(m)
    plot(
        hcat(map(x->[m.points[x.src, 1]; m.points[x.dst, 1]],edges(m.delaunay))...),
        hcat(map(x->[m.points[x.src, 2]; m.points[x.dst, 2]],edges(m.delaunay))...),
        color="red",
        legend=false
    )
    plot_points!(m.points)
end
function plot_map!(m)
    plot!(
        hcat(map(x->[m.points[x.src, 1]; m.points[x.dst, 1]],edges(m.delaunay))...),
        hcat(map(x->[m.points[x.src, 2]; m.points[x.dst, 2]],edges(m.delaunay))...),
        color="red",
        legend=false
    )
    plot_points!(m.points)
end

function to_geogebra(points)
    (*)(
        "Point[{",
        join(["{$x, $y}" for (x,y) in zip(points[:,1], points[:,2])], ", "),
        "}]"
    )
end

function to_geogebra(map::Map)
    convert_to_ascii(n) = Char(n + 64)
    join(
        [
            to_geogebra(map.points),
            (*)(
                "{",
                join(["Segment[$(convert_to_ascii(e.src)), $(convert_to_ascii(e.dst))]" for e in edges(map.delaunay)], ", "),
                "}"
            )
        ],
        "\n"
    )
end

end

max_stars_repo_name: Klafyvel/Land-of-Julia

max_issues_repo_name: Klafyvel/Land-of-Julia

max_forks_repo_name: Klafyvel/Land-of-Julia


text:
using DataFrames
using BenchmarkTools

function loadmgf(fname)
    FIELDS = ["TITLE=", "RTINSECONDS=", "PEPMASS=", "CHARGE=", "SCANS="]

    function format_precursor!(spectrum)
        if occursin( " ", spectrum["PEPMASS"] )
            spectrum["PEPMASS"] = map( x -> tryparse(Float64, x), split( spectrum["PEPMASS"] , " " ) )
        else
            spectrum["PEPMASS"] = [ tryparse(Float64, spectrum["PEPMASS"] ) ]
        end

        local polarity_multiplier = 1
        if spectrum["CHARGE"][end] == '-'
            polarity_multiplier = -1
        end

        if isnumeric( spectrum["CHARGE"][end] ) == false
            spectrum["CHARGE"] = polarity_multiplier * tryparse(Int32, spectrum["CHARGE"][1:end-1] )
        else
            spectrum["CHARGE"] = tryparse(Int32, spectrum["CHARGE"])
        end
        return true
    end

    function msdata_to_df!(spectrum)
        spectrum["ms_data"] = DataFrame( spectrum["ms_data"] )
        return true
    end

    local spectrum
    spectra_array = []
    open(fname) do fh
        state = false
        for line::String in eachline(fh)

            if length(line) > 0 && isnumeric( line[1] ) && state == true
                num_line = map( x -> tryparse(Float64, x), split(line, " ") )
                push!(
                    spectrum["ms_data"]["m/z"],
                    num_line[1]
                )
                push!(
                    spectrum["ms_data"]["Intensity"],
                    num_line[2]
                )
            elseif occursin("BEGIN IONS", line)
                spectrum = Dict{String,Any}( [ "ms_data" => Dict(
                            [
                                "m/z"=> Array{Float64}(undef, 0),
                                "Intensity"=> Array{Float64}(undef, 0)
                            ]
                        ) ] )
                state = true

            elseif occursin("END IONS", line)
                #println(spectrum)
                if length( spectrum["ms_data"] ) > 0
                    format_precursor!(spectrum)
                    msdata_to_df!(spectrum)
                    push!(spectra_array, spectrum)
                end
                state = false
            else
                for field_name in FIELDS
                    if occursin(field_name, line) && state == true
                        spectrum[ field_name[1:end-1] ] = split(line, field_name)[2]
                    end
                end
            end
        end
    end
    spectra_array
end

fname = "Yeast_1000spectra.mgf"

res = loadmgf(fname)
length(res)

res[501]

res[501]["ms_data"]

@benchmark loadmgf(fname)

AA_DELTAS = Dict(
    "G"=> 57.02147, "A"=> 71.03712, "S"=> 87.03203, "P"=> 97.05277, "V"=> 99.06842,
    "T"=> 101.04768, "Ccam"=> 160.03065, "Cmes"=> 148.996912, "I/L"=> 113.08407,
    "N"=> 114.04293, "D"=> 115.02695, "Q"=> 128.05858, "K"=> 128.09497, "E"=> 129.0426,
    "M"=> 131.04049, "Mox"=> 147.0354, "H"=> 137.05891, "F"=> 147.06842, "R"=> 156.10112,
    "Y"=> 163.06333, "W"=> 186.07932
)

#Flatten the values from the dictionary
single_res_Δ = collect( values(AA_DELTAS) )
println( typeof(single_res_Δ) )
#Add doubly-charged and triply-charged mass Deltas (simply divide by 2 and 3)
single_res_Δ = vcat( single_res_Δ / 3, single_res_Δ / 2, single_res_Δ )
println( length( single_res_Δ ) )
single_res_Δ[1:5]

function matches(spectra, masses_to_match, rel_tolerance)
    res_dict = Dict([
        ( "Spectrum_idx", Vector{Int64}() ),
        ( "Exp_idx", Vector{Int64}() ),
        ( "Library_idx", Vector{Int64}() ),
        ( "Rel_error", Vector{Float64}() )
    ])
    min_theo_val = minimum(masses_to_match) * (1 - rel_tolerance)
   for (idx, s) in enumerate(spectra)
        exp_Δs = filter(
            x -> x > min_theo_val,
            s["ms_data"][!, "m/z"] .- s["ms_data"][!, "m/z"]'
        )
        rel_deltas_matrix = map(
            abs,
            (masses_to_match .- exp_Δs')
        ) * 2 ./ (masses_to_match .+ exp_Δs')
        matching_inds = findall(x -> isless(x, rel_tolerance), rel_deltas_matrix)
        num_matches = size(matching_inds)[1]
        if num_matches > 0
            append!( res_dict["Spectrum_idx"], fill(idx, num_matches) )
            append!( res_dict["Library_idx"], map(x -> x[1], matching_inds) )
            append!( res_dict["Exp_idx"], map(x -> x[2], matching_inds) )
            append!( res_dict["Rel_error"], rel_deltas_matrix[matching_inds] )
        end
    end
    DataFrame(res_dict)
end

m = matches(res, single_res_Δ, 1e-5)
size(m)

@benchmark matches(res, single_res_Δ, 1e-5) samples=5

function matchesfuse(spectra, masses_to_match, rel_tolerance)
    res_dict = Dict([
        ( "Spectrum_idx", Vector{Int64}() ),
        ( "Exp_idx", Vector{Int64}() ),
        ( "Library_idx", Vector{Int64}() ),
        ( "Rel_error", Vector{Float64}() )
    ])
    min_theo_val = minimum(masses_to_match) * (1 - rel_tolerance)
   for (idx, s) in enumerate(spectra)
        exp_Δs = filter(
            x -> x > min_theo_val,
            s["ms_data"][!, "m/z"] .- s["ms_data"][!, "m/z"]'
        )
        #Fuse the vectorized operations using the dot macro
        rel_deltas_matrix = @. map(abs, (masses_to_match - exp_Δs') ) * 2 / (masses_to_match + exp_Δs')
        matching_inds = findall(x -> isless(x, rel_tolerance), rel_deltas_matrix)
        num_matches = size(matching_inds)[1]
        if num_matches > 0
            append!( res_dict["Spectrum_idx"], fill(idx, num_matches) )
            append!( res_dict["Library_idx"], map(x -> x[1], matching_inds) )
            append!( res_dict["Exp_idx"], map(x -> x[2], matching_inds) )
            append!( res_dict["Rel_error"], rel_deltas_matrix[matching_inds] )
        end
    end
    DataFrame(res_dict)
end

m = matchesfuse(res, single_res_Δ, 1e-5)
size(m)

m[ m.Spectrum_idx .== 2 , :]

@benchmark matchesfuse(res, single_res_Δ, 1e-5) samples=5

ext: jl

max_stars_repo_name: dev-ev/julia-vs-python-ms-scripting

max_issues_repo_name: dev-ev/julia-vs-python-ms-scripting

max_forks_repo_name: dev-ev/julia-vs-python-ms-scripting


text:
# # Shor's 9 qubit code

# The well-known Shor's 9 qubit code can correct any single Pauli error, but it can also correct some other
# kinds of errors. Its circuit is shown as below
# ![shor-9-circuit](assets/shor9circuit.png)
# which can be constructed by the following code

using Yao

shor(E) = chain(9,
    ## encode circuit
    cnot(1, 4), cnot(1, 7),
    put(1=>H), put(4=>H), put(7=>H),
    cnot(1,2), cnot(1,3), cnot(4,5), cnot(4,6), cnot(7,8), cnot(7,9),
    E, # the error
    ## decode circuit
    cnot(1,2), cnot(1,3), cnot((2, 3), 1),
    cnot(4,5), cnot(4,6), cnot((5, 6), 4),
    cnot(7,8), cnot(7,9), cnot((8, 9), 7),
    put(1=>H), put(4=>H), put(7=>H), cnot(1, 4), cnot(1, 7), cnot((4, 7), 1)
)

# Now we can check whether it can correct a given error by doing symbolic computation
# on an arbitrary 1-qubit pure quantum state ``α|0⟩ + β|1⟩`` and a specific weight-9 error.

@vars α β
s = α * ket"0" + β * ket"1" |> addbits!(8)
E = kron(1=>X, 2=>Z, 3=>Z, 4=>X, 5=>Z, 6=>Z, 7=>X, 8=>Z, 9=>Z);
s |> shor(E) |> partial_tr(2:9) |> expand


text:
function best_psis(psis::Vector)
  k = length(psis)
  idx = sortperm(psis, rev=true)    # sort, get the indices

  max_id = 1
  max_val = psis[idx[1]]
  for j=2:k
    v = ( sum(psis[idx[1:j]]) + j - 1 ) / j
    if v >= max_val
      max_val = v
      max_id = j
    else
      break
    end
  end

  psi_id = idx[1:max_id]

  return psi_id::Vector{Int}, max_val::Float64
end

function calc_const(psi_list::Vector, psi_id::Vector)
  j = length(psi_id)
  ret = ( sum(psi_list[psi_id]) + j - 1 ) / j
  return ret::Float64
end

function calc_cconst(psi_id::Vector)
  n_psi = length(psi_id)
  ret = (n_psi - 1 ) / n_psi
  return ret::Float64
end

## get featrure indeces for class i
function idi(m::Integer, i::Integer)
  return ((i-1)*m+1 : i*m)::UnitRange{Int64}
end

function psi_list(w::Vector, X::Matrix, y::Vector, i::Integer, c::Integer, idmi::Vector)
  psis = zeros(c)
  yi = y[i]
  for j=1:c
    if j != yi
      v1 = dot(w[idmi[j]], view(X, :, i))
      v2 = dot(w[idmi[yi]], -view(X, :, i))
      psis[j] = v1 + v2
    end
  end

  return psis::Vector{Float64}
end

## no y. for prediction
function psi_list(w::Vector, X::Matrix, i::Integer, c::Integer, idmi::Vector)
  psis = zeros(c)
  for j=1:c
    val = dot(w[idmi[j]], view(X, :, i))
    psis[j] = val
  end

  return psis::Vector{Float64}
end

# in terms of dual variable (useful for kernel methods)
function psi_list_dual(alpha::Vector, LPsi::Matrix, i::Integer, c::Integer)
  LPsi_i = view(LPsi, (i-1)*c+1 : i*c, :)
  psis = -(LPsi_i * alpha)
  return psis::Vector{Float64}
end

# in terms of dual variable (useful for kernel methods) dec
function psi_list_dual(sLa::Vector, i::Integer, c::Integer)
  psis = -(sLa[(i-1)*c+1 : i*c])
  return psis::Vector{Float64}
end

function eta_list(w::Vector, X::Matrix, i::Integer, c::Integer, idmi::Vector)
  etas = zeros(c)
  for j=1:c
    etas[j] = dot(w[idmi[j]], view(X, :, i))
  end

  return etas::Vector{Float64}
end

function calc_dot(key::Tuple{Integer,Vector,Integer,Vector}, K_ij::Float64, y::Vector)

  i = key[1]
  j = key[3]
  psi_i = key[2]
  psi_j = key[4]
  li = length(psi_i)
  lj = length(psi_j)
  yi = y[i]
  yj = y[j]

  mult = 0.0

  inii = yi in psi_i
  inij = yi in psi_j
  inji = yj in psi_i
  injj = yj in psi_j

  if yi == yj
    mult += (li - round(Int, inii)) * (lj - round(Int, injj))
  else
    if inij
      mult -= (li - round(Int, inii))
    end
    if inji
      mult -= (lj - round(Int, injj))
    end
  end

  ii = 1
  ij = 1
  while ii <= li && ij <= lj
    if psi_i[ii] > psi_j[ij]
      ij += 1
    elseif psi_i[ii] < psi_j[ij]
      ii += 1
    else  # equal
      if psi_i[ii] != yi && psi_j[ij] != yj && psi_i[ii] != yj && psi_j[ij] != yi
        mult += 1
      end
      ii += 1
      ij += 1
    end
  end

  d = (mult * K_ij) / (li * lj)

  return d::Float64
end

function calc_dot(key::Tuple{Integer,Vector,Integer,Vector}, K::Matrix, y::Vector)
  i = key[1]
  j = key[3]
  K_ij = K[i,j]

  d = calc_dot(key, K_ij, y)

  return d::Float64
end

function calc_dconst(key::Tuple{Integer,Vector}, X::Matrix, y::Vector, n_c::Integer, idmi::Vector)

    i = key[1]
    psi_i = key[2]
    li = length(psi_i)
    yi = y[i]

    xi = view(X, :, i)

    m = length(xi)
    dc = zeros(m * n_c)

    inii = yi in psi_i
    dc[idmi[yi]] = - ( (li - round(Int, inii)) * xi ) / li

    for ii = 1:li
      if psi_i[ii] != yi
        dc[idmi[psi_i[ii]]] = xi / li
      end
    end

    return dc::Vector{Float64}
end

# calc Lambda * Phi for kernel prediction
function calc_dotlphi(key::Tuple{Integer,Vector,Integer}, K_ij::Float64, y::Vector, n_c::Integer)
  i = key[1]
  psi_i = key[2]
  li = length(psi_i)
  yi = y[i]

  mults = zeros(n_c)
  for ii in psi_i
    if ii == yi
      mults[ii] = -(li - 1.0) / li
    else
      mults[ii] = 1.0 / li
    end
  end

  if !(yi in psi_i)
    mults[yi] = -1.0
  end

  return (mults * K_ij)::Vector{Float64}
end

function calc_dotlphi(key::Tuple{Integer,Vector,Integer}, K::Matrix, y::Vector, n_c::Integer)
  i = key[1]
  j = key[3]
  K_ij = K[i,j]
  return calc_dotlphi(key, K_ij, y, n_c)::Vector{Float64}
end

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
struct COSMOSACdspParam <: EoSParam
    Pnhb::SingleParam{Vector{Float64}}
    POH::SingleParam{Vector{Float64}}
    POT::SingleParam{Vector{Float64}}
    epsilon::SingleParam{Float64}
    V::SingleParam{Float64}
    A::SingleParam{Float64}
    water::SingleParam{Int64}
    COOH::SingleParam{Int64}
    hb_acc::SingleParam{Int64}
    hb_don::SingleParam{Int64}
end

abstract type COSMOSACdspModel <: COSMOSAC10Model end

struct COSMOSACdsp{c<:EoSModel} <: COSMOSACdspModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    params::COSMOSACdspParam
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel COSMOSACdsp
export COSMOSACdsp

function COSMOSACdsp(components; puremodel=PR,
    userlocations=String[],
     verbose=false)
    params = getparams(components, ["Activity/COSMOSAC/COSMOSAC10_like.csv","Activity/COSMOSAC/COSMOSACdsp_like.csv"]; userlocations=userlocations, verbose=verbose)
    Pnhb  = COSMO_parse_Pi(params["Pnhb"])
    POH  = COSMO_parse_Pi(params["POH"])
    POT  = COSMO_parse_Pi(params["POT"])
    A  = params["A"]
    V  = params["V"]
    epsilon  = params["epsilon"]
    water = params["water"]
    COOH = params["COOH"]
    hb_acc = params["hb_acc"]
    hb_don = params["hb_don"]
    icomponents = 1:length(components)

    init_puremodel = [puremodel([components[i]]) for i in icomponents]
    packagedparams = COSMOSACdspParam(Pnhb,POH,POT,epsilon,V,A,water,COOH,hb_acc,hb_don)
    references = String[]
    model = COSMOSACdsp(components,icomponents,packagedparams,init_puremodel,1e-12,references)
    return model
end

function activity_coefficient(model::COSMOSACdspModel,V,T,z)
    return exp.(@f(lnγ_comb) .+@f(lnγ_res).+@f(lnγ_dsp))
end

function lnγ_dsp(model::COSMOSACdspModel,V,T,z)
    x = z./sum(z)
    if model.params.water.values[1]==1 && model.params.hb_acc.values[2]==1 && model.params.hb_don.values[2]==0
        w = -0.27027
    elseif model.params.water.values[2]==1 && model.params.hb_acc.values[1]==1 && model.params.hb_don.values[1]==0
        w = -0.27027
    elseif model.params.COOH.values[1]==1 && (model.params.hb_acc.values[2]==1 | model.params.hb_don.values[2]==1)
        w = -0.27027
    elseif model.params.COOH.values[1]==1 && model.params.water.values[2]==1
        w = -0.27027
    elseif model.params.COOH.values[2]==1 && model.params.water.values[1]==1
        w = -0.27027
    else
        w = 0.27027
    end

    ϵ = model.params.epsilon.values

    A = w*(0.5*(ϵ[1]+ϵ[2])-√(ϵ[1]*ϵ[2]))
    return A*(1 .-x).^2
end


text:
using OffsetArrays
using Test
ambs = detect_ambiguities(Base, Core)  # in case these have ambiguities of their own
using PaddedViews
@testset "ambiguities" begin
    @test isempty(setdiff(detect_ambiguities(PaddedViews, Base, Core), ambs))
end

@testset "PaddedView" begin
    for n = 0:5
        a = @inferred(PaddedView(0, ones(Int,ntuple(d->1,n)), ntuple(x->x+1,n)))
        @test axes(a) == ntuple(x->1:x+1,n)
        @test @inferred(a[1]) == 1
        n > 0 && @test @inferred(a[2]) == 0
        @test @inferred(a[ntuple(x->1,n)...]) == 1
        n > 0 && @test @inferred(a[2, ntuple(x->1,n-1)...]) == 0
    end
    a0 = reshape([3])
    a = @inferred(PaddedView(-1, a0, ()))
    @test axes(a) == ()
    @test ndims(a) == 0
    @test a[] == 3

    a = reshape(1:9, 3, 3)
    A = @inferred(PaddedView(0.0, a, (Base.OneTo(4), Base.OneTo(5))))
    @test eltype(A) == Int
    @test ndims(A) == 2
    @test size(A) === (4,5)
    @test @inferred(axes(A)) === (Base.OneTo(4), Base.OneTo(5))
    @test @inferred(axes(A, 3)) === Base.OneTo(1)
    @test A == [1 4 7 0 0;
                2 5 8 0 0;
                3 6 9 0 0;
                0 0 0 0 0]

    A = @inferred(PaddedView(0.0, a, (0:4, -1:5)))
    @test eltype(A) == Int
    @test ndims(A) == 2
    @test size(A) == (5, 7)
    @test @inferred(axes(A)) === (0:4, -1:5)
    @test @inferred(axes(A, 3)) === 1:1
    @test A == OffsetArray([0 0 0 0 0 0 0;
                            0 0 1 4 7 0 0;
                            0 0 2 5 8 0 0;
                            0 0 3 6 9 0 0;
                            0 0 0 0 0 0 0], 0:4, -1:5)

    A = @inferred(PaddedView(0.0, a, (Base.OneTo(5), Base.OneTo(5)), (2:4, 2:4)))
    @test A == [0 0 0 0 0;
                0 1 4 7 0;
                0 2 5 8 0;
                0 3 6 9 0;
                0 0 0 0 0]
    @test A == @inferred(PaddedView(0.0, a, (5, 5), (2, 2)))

    @test PaddedView(0, ones(4), (3,), (1,)) == PaddedView(0, ones(4), (3,))
    let a = reshape(1:6, 2, 3) # [1 3 5; 2 4 6]
        @test @inferred(PaddedView(-1, a, (2, 2))) == [1 3
                                                       2 4]

        @test @inferred(PaddedView(-1, a, (1, 4)))  == [1 3 5 -1 ]

        @test @inferred(PaddedView(-1, a, (2, 4), (1, 1)))  == [1 3 5 -1
                                                                2 4 6 -1]

        @test @inferred(PaddedView(-1, a, (2, 4), (1, 2)))  == [-1 1 3 5
                                                                -1 2 4 6]

        @test @inferred(PaddedView(-1, a, (2, 3), (0, 2)))  == [-1  2  4
                                                                -1 -1 -1]

        @test @inferred(PaddedView(-1, a, (1, 4), (1, -1))) == [5 -1 -1 -1]

        # Create an OffsetArray with axes (0:1, 2:4)
        oa = OffsetArray(a, -1, 1)
        # Make a PaddedView of its elements (1:3, 1:2) (i.e pv[1, 1] == oa[1, 1] == "a[2, 0]" == -1)
        pv = @inferred PaddedView(-1, oa, (3, 2))
        @test pv[1:3, 1:2] == [-1 2; -1 -1; -1 -1]
        # Put oa[1, 1] in pv[2, 0] and make a PaddedView of the resulting (1:3, 1:2) elements
            # i.e. same as putting oa[0,2] in pv[1, 1]
        pv = @inferred PaddedView(-1, oa, (3, 2), (2, 0))
        @test pv[1:3, 1:2] == [1 3; 2 4; -1 -1]
    end
end

@testset "paddedviews" begin
    a1 = reshape([1,2], 2, 1)
    a2 = [1.0,2.0]'
    a1p, a2p = @inferred(paddedviews(0, a1, a2))
    @test a1p == [1 0; 2 0]
    @test a2p == [1.0 2.0; 0.0 0.0]
    @test eltype(a1p) == Int
    @test eltype(a2p) == Float64
    @test axes(a1p) === axes(a2p) === (Base.OneTo(2), Base.OneTo(2))

    a3 = OffsetArray([1.0,2.0]', (0,-1))
    a1p, a3p = @inferred(paddedviews(0, a1, a3))
    @test a1p == OffsetArray([0 1; 0 2], 1:2, 0:1)
    @test a3p == OffsetArray([1.0 2.0; 0.0 0.0], 1:2, 0:1)
    @test eltype(a1p) == Int
    @test eltype(a3p) == Float64
    @test axes(a1p) === axes(a3p) === (1:2, 0:1)

    @test @inferred(paddedviews(3)) == ()
    @test_throws ErrorException PaddedViews.outerinds()
    # But a zero-dimensional input should not trigger that error
    a = reshape([5])
    @test paddedviews(-1, a) == (a,)
end

nothing


text:
using DataDeps

register(DataDep(
  "blog-feedback",
  "https://archive.ics.uci.edu/ml/datasets/BlogFeedback",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/00304/BlogFeedback.zip",
  "1ba74e5ad920f7cd037502b2968581cc695146a226a73eff52fe8ad875ed4bcf",
  post_fetch_method=(path -> begin
    unpack(path)
    new_path = joinpath(dirname(path), "blogData_train.csv")
    mv(new_path, path)

    UCIData.process_dataset(path,
        target_index=281,
        feature_indices=1:280,
        categoric_indices=63:262,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@testset "piecewisepolynomial" begin
    v1 = rand(3)
    v2 = rand(3)
    m1 = rand(3, 4)
    m2 = rand(3, 4)
    maha = ones(3, 3)
    v = 3
    k = PiecewisePolynomialKernel{v}(maha)

    k2 = PiecewisePolynomialKernel(v=v, maha=maha)

    @test k2(v1, v2) ≈ k(v1, v2) atol=1e-5

    @test typeof(k(v1, v2)) <: Real
    @test size(kernelmatrix(k, m1, m2)) == (4, 4)
    @test size(kernelmatrix(k, m1)) == (4, 4)

    A1 = ones(4, 4)
    kernelmatrix!(A1, k, m1, m2)
    @test A1 ≈ kernelmatrix(k, m1, m2) atol=1e-5

    A2 = ones(4, 4)
    kernelmatrix!(A2, k, m1)
    @test A2 ≈ kernelmatrix(k, m1) atol=1e-5

    @test size(kerneldiagmatrix(k, m1)) == (4,)
    @test kerneldiagmatrix(k, m1) == ones(4)
    A3 = ones(4)
    kerneldiagmatrix!(A3, k, m1)
    @test A3 == kerneldiagmatrix(k, m1)

    @test_throws ErrorException PiecewisePolynomialKernel{4}(maha)

    @test repr(k) == "Piecewise Polynomial Kernel (v = $(v), size(maha) = $(size(maha)))"
    # test_ADs(maha-> PiecewisePolynomialKernel(v=2, maha = maha), maha)
    @test_broken "Nothing passes (problem with Mahalanobis distance in Distances)"
end


text:
@testset "FeedbackTrees" begin
    @testset "constructor" begin
        # 1. can construct a FeedbackTree like a normal chain
        c = FeedbackTree(Dense(10, 5), relu)
        @test length(c.layers) == 2
        # 2. can construct a FeedbackTree with forkers and mergers
        c = FeedbackTree(Merger("f1", identity, +), Dense(10, 10, relu), Splitter("f1"), Dense(10, 1))
        @test length(c.layers) == 4
    end # @testset "constructor"

    @testset "apply" begin
        h = Dict("f1" => rand(10))
        x = randn(10)
        # apply empty tree
        c = FeedbackTree()
        @test c(h, x) == (h, x)
        # apply tree without feedback
        c = FeedbackTree(x -> (2.0 .* x), x -> relu.(x), x -> x'x)
        state, y = c(Dict(), x)
        @test y ≈ relu.(2.0 .* x)'relu.(2.0 .* x)
        @test state isa Dict{String, Any}
        # apply tree with feedback
        l1 = Dense(10, 10, relu)
        l2 = Dense(10, 1)
        c = FeedbackTree(Merger("f1", identity, +), l1, Splitter("f1"), l2)
        state, y = c(h, x)
        @test y ≈ l2(h["f1"])
        @test haskey(state, "f1")
        @test state["f1"] ≈ l1(h["f1"] + x)
    end # @testset "apply"

    @testset "params" begin
        l1 = Dense(10, 10)
        l2 = Dense(7, 10)
        c = FeedbackTree(l1, l2)
        @test hcat(params(c)...) == hcat(params(l1)..., params(l2)...)
    end # @testset "params"

    @testset "statefull" begin
        h = Dict("f1" => randn(10))
        x = randn(10)
        l1 = Dense(10, 10, relu)
        l2 = Dense(10, 1)
        c = FeedbackTree(Merger("f1", identity, +), l1, Splitter("f1"), l2)
        c = StateFull(c, h)
        y = c(x)
        @test y ≈ l2(h["f1"])
        @test haskey(c.state, "f1")
        @test c.state["f1"] ≈ l1(h["f1"] + x)
        # check that truncation works
        #@test Flux.Tracker.istracked(c.state["f1"])
        #Flux.truncate!(c)
        #@test !Flux.Tracker.istracked(c.state["f1"])
        # check that reset works
        Flux.reset!(c)
        @test c.state["f1"] == c.init["f1"]
    end # @testset "statefull"

    @testset "indexing" begin
        c = FeedbackTree(Dense(1,2), Dense(2,3), Dense(3,1))
        @test c[1] == c.layers[1]
        @test c[2:3] == FeedbackTree(c.layers[2:3]...)
    end # @testset "indexing"

    @testset "iteration" begin
        c = FeedbackTree(Dense(1,2), Dense(2,3), Dense(3,1))
        layers = [layer for layer in c]
        @test layers == collect(c.layers)
    end # @testset "iteration"

    @testset "interface" begin
        # test that generic method `splitname` works
        c = FeedbackTree(
            Splitter("name1"),
            Merger("name2", nothing, +),
            Splitter("name3"),
            Dense(10,10)
        )
        @test splitnames(c) == ["name1", "name3"]
        # test that validation of naming works
        @test !namesvalid(c) # merger name does not fit splitter names
        c = FeedbackTree(Splitter("name1"), Merger("name1", nothing, +))
        @test namesvalid(c)
        c = FeedbackTree(Splitter("name1"), Splitter("name1"))
        @test !namesvalid(c)
    end # @testset "interface"
end # @testset "FeedbackTrees"


text:
module DockerSandbox

import Base: run, success
import Random

export DockerConfig
export run
export success
export with_container

include("types.jl")

include("base.jl")
include("container.jl")
include("convenience.jl")

end # module


text:
using JuliADAPT
using Test

@testset "JuliADAPT.jl" begin
    # Write your tests here.
    include("test_fci.jl")
end

max_stars_repo_name: hrgrimsl/JuliADAPT

max_issues_repo_name: hrgrimsl/JuliADAPT

max_forks_repo_name: hrgrimsl/JuliADAPT


text:
function sort_sections(chn::Chains)
    smap = keys(chn.name_map)
    section_list = Vector{Symbol}(undef, length(smap))
    indx = 1
    if :parameters in smap
        section_list[1] = :parameters
        indx += 1
    end
    if :internals in smap
        section_list[end] = :internals
    end
    for par in smap
        if !(par in [:parameters, :internals])
            section_list[indx] = par
            indx += 1
        end
    end
    return section_list
end

"""

# Array

Array constructor from a Chains object. Returns 3 dimensionsal
array or an Array of 2 dimensional Arrays. If only a single parameter is selected for
inclusion, a dimension is dropped in both cases, as is e.g. required by cde(), etc.

### Method
```julia
  Array(
    chn::Chains,
    sections::Vector{Symbol};
    append_chains::Bool,
    remove_missing_union::Bool
  )
```

### Required arguments
```julia
* `chn` : Chains object to convert to an Array
```

### Optional arguments
```julia
* `sections = Symbol[]` : Sections from the Chains object to be included
* `append_chains = true`  : Append chains into a single column
* `remove_missing_union = true`  : Convert Union{Missing, Real} to Float64
```

### Examples
```julia
* `Array(chns)` : Array with chain values are appended
* `Array(chns[:par])` : Array with selected parameter chain values are appended
* `Array(chns, [:parameters])`  : Array with only :parameter section
* `Array(chns, [:parameters, :internals])`  : Array includes multiple sections
* `Array(chns, append_chains=false)`  : Array of Arrays, each chain in its own array
* `Array(chns, remove_missing_union=false)` : No conversion to remove missing values
```

"""
function Base.Array(chain::Chains,
    sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters];
    append_chains = true,
    remove_missing_union = true,
    showall = false,
    sorted = false
)
    sections = showall ? keys(chain.name_map) : sections
    if remove_missing_union
        chn = concretize(Chains(chain, sections))
    else
        chn = Chains(chain, sections)
    end

    nparams = size(chn, 2)
    if append_chains
        if nparams == 1
            return to_vector(chn)
        else
            return to_matrix(chn)
        end
    else
        if nparams == 1
            return to_vector_of_vectors(chn)
        else
            return to_vector_of_matrices(chn)
        end
    end
end

function to_matrix(chain::Chains)
    return Matrix(reshape(permutedims(chain.value.data, (1, 3, 2)), :, size(chain, 2)))
end

function to_vector(chain::Chains)
    return Vector(vec(chain.value.data))
end

function to_vector_of_vectors(chain::Chains)
    data = chain.value.data
    return [Vector(vec(data[:, :, i])) for i in axes(data, 3)]
end

function to_vector_of_matrices(chain::Chains)
    data = chain.value.data
    return [Matrix(data[:, :, i]) for i in axes(data, 3)]
end


text:
Base.getindex(axis::AbstractAxis, ::Colon) = copy(axis)
Base.getindex(axis::AbstractAxis, ::Ellipsis) = copy(axis)
@propagate_inbounds Base.getindex(axis::AbstractAxis, arg::Integer) = to_index(axis, arg)
@propagate_inbounds function Base.getindex(axis::AbstractAxis, arg::AbstractUnitRange{I}) where {I<:Integer}
    return ArrayInterface.unsafe_getindex(axis, (to_index(axis, arg),))
end
@propagate_inbounds function Base.getindex(axis::AbstractAxis, arg::StepRange{I}) where {I<:Integer}
    return ArrayInterface.unsafe_getindex(axis, (to_index(axis, arg),))
end
@propagate_inbounds function Base.getindex(axis::AbstractAxis, arg::AbstractArray{I}) where {I<:Integer}
    return ArrayInterface.unsafe_getindex(axis, (to_index(axis, arg),))
end
@propagate_inbounds function Base.getindex(axis::AbstractAxis, arg)
    return ArrayInterface.unsafe_getindex(axis, (to_index(axis, arg),))
end

ArrayInterface.unsafe_get_element(axis::AbstractAxis, inds) = eltype(axis)(first(inds))

@inline function ArrayInterface.unsafe_get_collection(axis::AbstractAxis, inds::Tuple)
    return _unsafe_get_axis_collection(axis, first(inds))
end
_unsafe_get_axis_collection(axis, i::Integer) = Int(i)
@inline function _unsafe_get_axis_collection(axis, inds)
    if known_step(inds) === one(eltype(inds))
        return index_axis_to_axis(axis, inds)
    else
        return __unsafe_get_axis_collection(index_axis_to_array(axis, inds), inds)
    end
end
@inline function __unsafe_get_axis_collection(axis, inds::AbstractRange)
    T = eltype(axis)
    if eltype(inds) <: T
        return AxisArray{T,1,typeof(inds),Tuple{typeof(axis)}}(inds, (axis,); checks=NoChecks)
    else
        return __unsafe_get_axis_collection(axis, AbstractRange{T}(inds))
    end
end
@inline function __unsafe_get_axis_collection(axis, inds)
    T = eltype(axis)
    if eltype(inds) <: T
        return AxisArray{T,1,typeof(inds),Tuple{typeof(axis)}}(inds, (axis,); checks=NoChecks)
    else
        return __unsafe_get_axis_collection(axis, AbstractArray{T}(inds))
    end
end

#=
    index_axis_to_array(axis, inds)

An axis may have other axis types nested within it so we need to propagate the indexing,
but we can usually just use unsafe_reconstruct.
=#
@inline function index_axis_to_axis(axis::SimpleAxis, inds)
    return SimpleAxis(static_first(inds):static_last(inds))
end
@inline function index_axis_to_axis(axis, inds)
    return unsafe_reconstruct(axis, index_axis_to_axis(parent(axis), _sub_offset(axis, inds)))
end
@inline function index_axis_to_axis(axis::IdentityAxis, inds)
    return IdentityAxis(inds, index_axis_to_axis(parent(axis),  _sub_offset(axis, inds)))
end

# TODO this might be worth making part of official API if a better name and more thought out
# documentation/implementation accompanies it b/c eventually all axis types have to deal
# with non unit range indexing.
#=
    index_axis_to_array(axis, inds)

An axis cannot be preserved if the elements with any collection that doesn't have a step of 1.
=#
index_axis_to_array(axis::SimpleAxis, inds) = SimpleAxis(eachindex(inds))
function index_axis_to_array(axis::Axis, inds)
    if allunique(inds)  # propagate keys corresponds to inds
        return Axis(@inbounds(keys(axis)[inds]), index_axis_to_array(parent(axis), inds); checks=NoChecks)
    else  # b/c not all indices are unique it will result in non-unique keys so drop keys
        return index_axis_to_array(parent(axis), inds)
    end
end
function index_axis_to_array(axis, inds)
    return unsafe_reconstruct(axis, index_axis_to_array(parent(axis), _sub_offset(axis, inds)))
end


text:
n = rand(1:10)

r = matrixdepot("deriv2", n, false)
rf32 = matrixdepot("deriv2", Float32, n, false)

r2 = matrixdepot("deriv2", Float32, n, 2, false)
if mod(n, 2) == 0
    r3 = matrixdepot("deriv2", Float32, n, 3, false)
end
try
    matrixdepot("deriv2", Float64, n, 4, false)
catch ArgumentError
    println("deriv2: illegal value of example.")
end

A = matrixdepot("deriv2", n)

@test_approx_eq r.A*r.x  r.b
@test issym(r.A)

r = matrixdepot("shaw", 2*n, false)
A = matrixdepot("shaw", 2*n)

print(r)

@test issym(r.A)

rf32 = matrixdepot("shaw", Float32, 2*n, false)

r = matrixdepot("wing", n, false)
A = matrixdepot("wing", n)

@test r.A == matrixdepot("wing", n, 1/3, 2/3, false).A

r = matrixdepot("foxgood", n, false)
rf32 = matrixdepot("foxgood", Float32, n, false)
A = matrixdepot("foxgood", n)

r = matrixdepot("heat", 2*n, false)
rf32 = matrixdepot("heat", Float32, 2*n)

r = matrixdepot("baart", 2*n, false)
rf32 = matrixdepot("baart", Float32, 2*n, false)
A = matrixdepot("baart", 2*n)

n = rand(3:10)
r = matrixdepot("phillips", 4*n, false)
rf32 = matrixdepot("phillips", Float32, 4*n)

A = matrixdepot("gravity", n)

@test matrixdepot("gravity", n, 1, 0, 1, 0.25) == A

r1 = matrixdepot("gravity", Float32, n, 1, false)
r2 = matrixdepot("gravity", Float32, n, 2, false)
r3 = matrixdepot("gravity", Float32, n, 3, false)

A = matrixdepot("blur", n)

@test matrixdepot("blur", n, 3, 0.7) == A

r1 = matrixdepot("blur", Float32, n, false)

A = matrixdepot("spikes", n)

@test matrixdepot("spikes", n, 5) == A

n = rand(5:10)
r1 = matrixdepot("spikes", Float32, n, false)

A = matrixdepot("ursell", n)
@test matrixdepot("ursell", Float64, n) == A
r1 = matrixdepot("ursell", Float32, n, false)
print(r1)

A = matrixdepot("parallax", n)

m, n = size(A)
@test m == 26
@test matrixdepot("parallax", Float64, n) == A
r1 = matrixdepot("parallax", Float32, n, false)


text:
using graphtheoryufrj
using DataStructures

"""
Requirement (PT-BR): Compare o desempenho em termos de quantidade de memoria
utilizada pelas duas representacoes do grafo. Ou seja, determine a quantidade
de memoria (em MB) utilizada pelo seu programa quando voce representa o grafo
utilizando uma matriz de adjacencia e lista de adjacencia.
"""
function test_one(infile_name::ASCIIString)
	println("- Starting the Test One")
	try
		println("\n-- Time to fill as Adjacency List:\n")
		@time graph_list = simpleGraph(infile_name)

		println("\n-- Time to fill as Adjacency Matrix:\n")
		@time graph_matrix = simpleGraph(infile_name, "adjmatrix")

		println("\n# Graph Size report:")
		println("-- Adjacency List size   (In bytes) = ", sizeof(graph_list.Graph))
		println("-- Adjacency Matrix size (In bytes) = ", sizeof(graph_matrix.Graph))
	catch
		error("Error in \"test_one\" function.")
	end
	return "TEST_ONE_SUCESS"
end

"""
Requirement (PT-BR): Compare o desempenho em termos de tempo de execucao
das duas representacoes do grafo. Ou seja, determine o tempo necessario
para executar dez buscas em largura em cada um dos casos (utilize diferentes
vertices como ponto de partida da busca).
"""
function test_two(infile_name::ASCIIString, initial_vertex)
	println("- Starting the Test Two")
	try
		println("\n-- Performing 10 bfs on the Adjacency List:")
		graph_list = simpleGraph(infile_name)

		@time for i in initial_vertex
			fast_bfs(graph_list.Graph, s = i)
		end

		println("\n-- Performing 10 bfs on Adjacency Matrix:")
		graph_matrix = simpleGraph(infile_name, "adjmatrix")

		@time for i in initial_vertex
			fast_bfs(graph_matrix.Graph, s = i)
		end
	catch
		error("Error in \"test_two\" function.")
	end
	return "TEST_TWO_SUCESS"
end

"""
Requirement (PT-BR): Repita o item anterior para busca em profundidade (utilize
os mesmos 10 vertices iniciais)
"""
function test_three(infile_name::ASCIIString, initial_vertex)
	println("- Starting the Test Three")
	try
		# Calculate the time to perform 10 DFS in an Adjacency List
		println("\n-- Performing 10 dfs on the Adjacency List:")
		graph_list = simpleGraph(infile_name)

		@time for i in initial_vertex
			p, t = dfs(graph_list.Graph, s = i)
		end

		# Calculate the time to perform 10 DFS in an Adjacency Matrix
		println("\n-- Performing 10 dfs on Adjacency Matrix:")
		graph_matrix = simpleGraph(infile_name, "adjmatrix")

		@time for i in initial_vertex
			p, t = dfs(graph_matrix.Graph, s = i)
		end
	catch
		error("Error in \"test_three\" function.")
	end
	return "TEST_THREE_SUCESS"
end

"""
Requirement (PT-BR): Determine o pai dos vertices 10, 20, 30, 40, 50 na
arvore geradora induzida pela BFS e pela DFS quando iniciamos a busca nos
vertices 1, 2, 3, 4, 5.
"""
function test_four(infile_name::ASCIIString, outfile_name::ASCIIString)
	println("- Starting the Test Four")
	outfile = open(outfile_name, "w")

	try
		initial_vertex = [1 2 3 4 5]
		child = [10 20 30 40 50]

		println("\n-- Performing the BFS on the Adjacency List: ")
		graph_list = simpleGraph(infile_name)
		fathers = Array{Array{Int64},1}(length(child))

		# Initialize fathers with zeros
		@inbounds @simd for i in 1:length(child)
   	 		fathers[i] = Array{Int64}(0)
    	end

    	write(outfile, "\nBFS TIME\n")
    	# Do the bfs to find the parents and fill
		for i in initial_vertex
			l, p, t = bfs(graph_list.Graph, s = i)
			write(outfile, string("\nInitial vertex = ", i))
			for j in 1:length(child)
				write(outfile, "\nParent of child : ")
				for k in 1:length(p)
					if findfirst(p[k], child[j]) != 0
						fathers[j] = [child[j] k]
						write(outfile, string(fathers[j], "\n"))
						break
					end
				end
			end
		end

		println("\n-- Performing the BFS on the Adjacency List: ")
		fathers_dfs = Array{Array{Int64},1}(length(child))
		# Initialize fathers with zeros
		@inbounds @simd for i in 1:length(child)
   	 		fathers_dfs[i] = Array{Int64}(0)
    	end

    	write(outfile, "\nDFS TIME\n")
    	# Do the bfs to find the parents and fill
		for i in initial_vertex
			p, t = dfs(graph_list.Graph, s = i)
			write(outfile, string("\nInitial vertex = ", i))
			for j in 1:length(child)
				write(outfile, "\nParent of child : ")
				for k in 1:length(p)
					if findfirst(p[k], child[j]) != 0
						fathers_dfs[j] = [child[j] k]
						write(outfile, string(fathers_dfs[j], "\n"))
						break
					end
				end
			end
		end
	catch
		error("Error in \"test_four\" function.")
	finally
		close(outfile)
	end
	return "TEST_FOUR_SUCESS"
end

"""
Requirement (PT-BR): Obtenha as componentes conexas do grafo. Quantas
componentes conexas tem o grafo? Qual e o tamanho da maior e da menor
componente conexo?
"""
function test_five(infile_name::ASCIIString)
	println("- Starting the Test Five")
	#try
		println("\n-- Performing the connected components algorithm: ")
		graph_list = simpleGraph(infile_name)
		@time connected_comp = connected_components(graph_list.Graph)
		println("Number of connected components: ", length(connected_comp))
		println("Largest connected components size : ", maximum(connected_comp))
		println("Smallest connected components size: ", minimum(connected_comp))
	# catch
	# 	error("Error in \"test_five\" function.")
	# end
	return "TEST_FIVE_SUCESS"
end

"""
Requirement (PT-BR): Obtenha a distribuicao empirica do grau dos vertices.
Trace um grafico com seu resultado. Qual e o maior grau do grafo? E o menor? Como
isto se compara ao maior grau possivel?
"""
function test_six(infile_name::ASCIIString)
	println("\n- Starting the Test Six")
	try
		graph_list = simpleGraph(infile_name)
		println("-- Performing mean_degree and empirical_dist on Adjacency List")
		@time mean_degree, empirical_dist = graph_properties(graph_list)

		largest_degree = 0
		smallest_degree = graph_list.num_edge

		for i in 1:length(graph_list.Graph)
			x = length(graph_list.Graph[i])
			largest_degree =  (largest_degree < x) ? x : largest_degree
			smallest_degree = (smallest_degree > x) ? x : smallest_degree
		end

		println("--- Largest degree on the Graph = ", largest_degree)
		println("--- Smallest degree on the Graph = ", smallest_degree)
	catch
		error("""Error in \"test_six\" function.""")
	end
	return "TEST_SIX_SUCESS"
end

"""
Requirement (PT-BR): Determine o diametro do grafo. Lembrando que o diametro e
a maior distancia entre qualquer par de vertices do grafo (ou seja, o comprimento
do maior caminho minimo do grafo). Determine tambem o tempo de execucao necessario
para calcular o diametro.
"""
function test_seven(infile_name::ASCIIString)
	println("\n- Starting the Test Six")
	try
		graph_list = simpleGraph(infile_name)
		println("-- Performing mean_degree and empirical_dist on Adjacency List")
		@time max_diam = diameter(graph_list.Graph, graph_list.num_edge)
		println("Maximum diameter = ", max_diam)
	catch
		error("""Error in \"test_seven\" function.""")
	end
	return "TEST_SEVEN_SUCESS"
end


text:
"""
$(TYPEDEF)

One chemical reaction.

# Fields
$(FIELDS)

# Examples

```julia
using ModelingToolkit
@parameters t k[1:20]
@variables A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
       Reaction(k[2], [B], nothing),            # B -> 0
       Reaction(k[3],[A],[C]),                  # A -> C
       Reaction(k[4], [C], [A,B]),              # C -> A + B
       Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A
       Reaction(k[6], [A,B], [C]),              # A + B -> C
       Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A
       Reaction(k[8], [A,B], [A,C]),            # A + B -> A + C
       Reaction(k[9], [A,B], [C,D]),            # A + B -> C + D
       Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
       Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
       Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]),          # A+3B+4C -> 2C + 3D
       Reaction(k[13], [A,B], nothing, [3,1], nothing),           # 3A+B -> 0
       Reaction(k[14], nothing, [A], nothing, [2]),               # 0 -> 2A
       Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
       Reaction(k[16], [A], [B]; only_use_rate=true),             # A -> B with custom rate.
       Reaction(k[17]*A*exp(B), [C], [D], [2], [1]),              # 2C -> D with non constant rate.
       Reaction(k[18]*B, nothing, [B], nothing, [2]),             # 0 -> 2B with non constant rate.
       Reaction(k[19]*t, [A], [B]),                                # A -> B with non constant rate.
       Reaction(k[20]*t*A, [B,C], [D],[2,1],[2])                  # 2A +B -> 2C with non constant rate.
  ]
```

Notes:
- `nothing` can be used to indicate a reaction that has no reactants or no products.
  In this case the corresponding stoichiometry vector should also be set to `nothing`.
- The three-argument form assumes all reactant and product stoichiometric coefficients
  are one.
"""
struct Reaction{S, T <: Number}
    """The rate function (excluding mass action terms)."""
    rate
    """Reaction substrates."""
    substrates::Vector
    """Reaction products."""
    products::Vector
    """The stoichiometric coefficients of the reactants."""
    substoich::Vector{T}
    """The stoichiometric coefficients of the products."""
    prodstoich::Vector{T}
    """The net stoichiometric coefficients of all species changed by the reaction."""
    netstoich::Vector{Pair{S,T}}
    """
    `false` (default) if `rate` should be multiplied by mass action terms to give the rate law.
    `true` if `rate` represents the full reaction rate law.
    """
    only_use_rate::Bool
end

function Reaction(rate, subs, prods, substoich, prodstoich;
                  netstoich=nothing, only_use_rate=false, kwargs...)

    (isnothing(prods)&&isnothing(subs)) && error("A reaction requires a non-nothing substrate or product vector.")
    (isnothing(prodstoich)&&isnothing(substoich)) && error("Both substrate and product stochiometry inputs cannot be nothing.")
    if isnothing(subs)
        subs = Vector{Term}()
        !isnothing(substoich) && error("If substrates are nothing, substrate stiocihometries have to be so too.")
        substoich = typeof(prodstoich)()
    end
    if isnothing(prods)
        prods = Vector{Term}()
        !isnothing(prodstoich) && error("If products are nothing, product stiocihometries have to be so too.")
        prodstoich = typeof(substoich)()
    end
    subs = value.(subs)
    prods = value.(prods)
    ns = isnothing(netstoich) ? get_netstoich(subs, prods, substoich, prodstoich) : netstoich
    Reaction(value(rate), subs, prods, substoich, prodstoich, ns, only_use_rate)
end

# three argument constructor assumes stoichiometric coefs are one and integers
function Reaction(rate, subs, prods; kwargs...)

    sstoich = isnothing(subs) ? nothing : ones(Int,length(subs))
    pstoich = isnothing(prods) ? nothing : ones(Int,length(prods))
    Reaction(rate, subs, prods, sstoich, pstoich; kwargs...)
end

function namespace_equation(rx::Reaction, name, iv)
    Reaction(namespace_expr(rx.rate, name, iv),
             namespace_expr(rx.substrates, name, iv),
             namespace_expr(rx.products, name, iv),
             rx.substoich, rx.prodstoich,
             [namespace_expr(n[1],name,iv) => n[2] for n in rx.netstoich], rx.only_use_rate)
end

# calculates the net stoichiometry of a reaction as a vector of pairs (sub,substoich)
function get_netstoich(subs, prods, sstoich, pstoich)
    # stoichiometry as a Dictionary
    nsdict = Dict{Any, eltype(sstoich)}(sub => -sstoich[i] for (i,sub) in enumerate(subs))
    for (i,p) in enumerate(prods)
        coef = pstoich[i]
        @inbounds nsdict[p] = haskey(nsdict, p) ? nsdict[p] + coef : coef
    end

    # stoichiometry as a vector
    ns = [el for el in nsdict if el[2] != zero(el[2])]

    ns
end

"""
$(TYPEDEF)

A system of chemical reactions.

# Fields
$(FIELDS)

# Example
Continuing from the example in the [`Reaction`](@ref) definition:
```julia
rs = ReactionSystem(rxs, t, [A,B,C,D], k)
```
"""
struct ReactionSystem <: AbstractSystem
    """The reactions defining the system."""
    eqs::Vector{Reaction}
    """Independent variable (usually time)."""
    iv::Any
    """Dependent (state) variables representing amount of each species."""
    states::Vector
    """Parameter variables."""
    ps::Vector
    observed::Vector{Equation}
    """The name of the system"""
    name::Symbol
    """systems: The internal systems"""
    systems::Vector

    function ReactionSystem(eqs, iv, states, ps, observed, name, systems)
        new(eqs, value(iv), value.(states), value.(ps), observed, name, systems)
    end
end

function ReactionSystem(eqs, iv, species, params;
                        observed = [],
                        systems = [],
                        name = gensym(:ReactionSystem))

    #isempty(species) && error("ReactionSystems require at least one species.")
    ReactionSystem(eqs, iv, species, params, observed, name, systems)
end

function ReactionSystem(iv; kwargs...)
    ReactionSystem(Reaction[], iv, [], []; kwargs...)
end

function equations(sys::ModelingToolkit.ReactionSystem)
    eqs = get_eqs(sys)
    systems = get_systems(sys)
    if isempty(systems)
        return eqs
    else
        eqs = [eqs;
               reduce(vcat,
                      namespace_equations.(get_systems(sys));
                      init=[])]
        return eqs
    end
end

"""
    oderatelaw(rx; combinatoric_ratelaw=true)

Given a [`Reaction`](@ref), return the reaction rate law [`Operation`](@ref) used in
generated ODEs for the reaction. Note, for a reaction defined by

`k*X*Y, X+Z --> 2X + Y`

the expression that is returned will be `k*X(t)^2*Y(t)*Z(t)`. For a reaction
of the form

`k, 2X+3Y --> Z`

the `Operation` that is returned will be `k * (X(t)^2/2) * (Y(t)^3/6)`.

Notes:
- Allocates
- `combinatoric_ratelaw=true` uses factorial scaling factors in calculating the rate
    law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
    `combinatoric_ratelaw=false` then the ratelaw is `k*S^2`, i.e. the scaling factor is
    ignored.
"""
function oderatelaw(rx; combinatoric_ratelaw=true)
    @unpack rate, substrates, substoich, only_use_rate = rx
    rl = rate
    if !only_use_rate
        coef = one(eltype(substoich))
        for (i,stoich) in enumerate(substoich)
            coef *= factorial(stoich)
            rl   *= isone(stoich) ? substrates[i] : substrates[i]^stoich
        end
        combinatoric_ratelaw && (!isone(coef)) && (rl /= coef)
    end
    rl
end

function assemble_oderhs(rs; combinatoric_ratelaws=true)
    sts = get_states(rs)
    species_to_idx = Dict((x => i for (i,x) in enumerate(sts)))
    rhsvec         = Any[0 for i in eachindex(sts)]

    for rx in get_eqs(rs)
        rl = oderatelaw(rx; combinatoric_ratelaw=combinatoric_ratelaws)
        for (spec,stoich) in rx.netstoich
            i = species_to_idx[spec]
            if _iszero(rhsvec[i])
                signedrl  = (stoich > zero(stoich)) ? rl : -rl
                rhsvec[i] = isone(abs(stoich)) ? signedrl : stoich * rl
            else
                Δspec     = isone(abs(stoich)) ? rl : abs(stoich) * rl
                rhsvec[i] = (stoich > zero(stoich)) ? (rhsvec[i] + Δspec) : (rhsvec[i] - Δspec)
            end
        end
    end

    rhsvec
end

function assemble_drift(rs; combinatoric_ratelaws=true, as_odes=true)
    rhsvec = assemble_oderhs(rs; combinatoric_ratelaws=combinatoric_ratelaws)
    if as_odes
        D   = Differential(get_iv(rs))
        eqs = [Equation(D(x),rhs) for (x,rhs) in zip(get_states(rs),rhsvec) if (!_iszero(rhs))]
    else
        eqs = [Equation(0,rhs) for rhs in rhsvec if (!_iszero(rhs))]
    end
    eqs
end

function assemble_diffusion(rs, noise_scaling; combinatoric_ratelaws=true)
    sts  = get_states(rs)
    eqs  = Matrix{Any}(undef, length(sts), length(get_eqs(rs)))
    eqs .= 0
    species_to_idx = Dict((x => i for (i,x) in enumerate(sts)))

    for (j,rx) in enumerate(equations(rs))
        rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw=combinatoric_ratelaws)))
        (noise_scaling!==nothing) && (rlsqrt *= noise_scaling[j])
        for (spec,stoich) in rx.netstoich
            i            = species_to_idx[spec]
            signedrlsqrt = (stoich > zero(stoich)) ? rlsqrt : -rlsqrt
            eqs[i,j]     = isone(abs(stoich)) ? signedrlsqrt : stoich * rlsqrt
        end
    end
    eqs
end

"""
    jumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=true)

Given a [`Reaction`](@ref), return the reaction rate law [`Operation`](@ref) used in
generated stochastic chemical kinetics model SSAs for the reaction. Note,
for a reaction defined by

`k*X*Y, X+Z --> 2X + Y`

the expression that is returned will be `k*X^2*Y*Z`. For a reaction of
the form

`k, 2X+3Y --> Z`

the `Operation` that is returned will be `k * binomial(X,2) *
binomial(Y,3)`.

Notes:
- `rxvars` should give the `Variable`s, i.e. species and parameters, the rate depends on.
- Allocates
- `combinatoric_ratelaw=true` uses binomials in calculating the rate law, i.e. for `2S ->
  0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaw=false` then
  the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling
  factor.
"""
function jumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=true)
    @unpack rate, substrates, substoich, only_use_rate = rx
    rl = rate
    if !only_use_rate
        coef = one(eltype(substoich))
        for (i,stoich) in enumerate(substoich)
            s   = substrates[i]
            rl *= s
            isone(stoich) && continue
            for i in one(stoich):(stoich-one(stoich))
                rl *= (s - i)
            end
            combinatoric_ratelaw && (coef *= factorial(stoich))
        end
        !isone(coef) && (rl /= coef)
    end
    rl
end

# if haveivdep=false then time dependent rates will still be classified as mass action
"""
```julia
ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                              haveivdep = any(var -> isequal(get_iv(rs),var), rxvars),
                              stateset = Set(states(rs)))
```

True if a given reaction is of mass action form, i.e. `rx.rate` does not depend
on any chemical species that correspond to states of the system, and does not depend
explicitly on the independent variable (usually time).

# Arguments
- `rx`, the [`Reaction`](@ref).
- `rs`, a [`ReactionSystem`](@ref) containing the reaction.
- Optional: `rxvars`, `Variable`s which are not in `rxvars` are ignored as possible dependencies.
- Optional: `haveivdep`, `true` if the [`Reaction`](@ref) `rate` field explicitly depends on the independent variable.
- Optional: `stateset`, set of states which if the rxvars are within mean rx is non-mass action.
"""
function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                              haveivdep = any(var -> isequal(get_iv(rs),var), rxvars),
                              stateset = Set(get_states(rs)))
    # if no dependencies must be zero order
    (length(rxvars)==0) && return true
    haveivdep && return false
    rx.only_use_rate && return false
    @inbounds for i = 1:length(rxvars)
        (rxvars[i] in stateset) && return false
    end
    return true
end

@inline function makemajump(rx; combinatoric_ratelaw=true)
    @unpack rate, substrates, substoich, netstoich = rx
    zeroorder = (length(substoich) == 0)
    reactant_stoch = Vector{Pair{Any,eltype(substoich)}}(undef, length(substoich))
    @inbounds for i = 1:length(reactant_stoch)
        reactant_stoch[i] = substrates[i] => substoich[i]
    end
    #push!(rstoich, reactant_stoch)
    coef = (zeroorder || (!combinatoric_ratelaw)) ? one(eltype(substoich)) : prod(stoich -> factorial(stoich), substoich)
    (!isone(coef)) && (rate /= coef)
    #push!(rates, rate)
    net_stoch      = [Pair(p[1],p[2]) for p in netstoich]
    #push!(nstoich, net_stoch)
    MassActionJump(Num(rate), reactant_stoch, net_stoch, scale_rates=false, useiszero=false)
end

function assemble_jumps(rs; combinatoric_ratelaws=true)
    meqs = MassActionJump[]; ceqs = ConstantRateJump[]; veqs = VariableRateJump[]
    stateset = Set(get_states(rs))
    #rates = [];  rstoich = []; nstoich = []
    rxvars = []
    ivname = nameof(get_iv(rs))

    isempty(equations(rs)) && error("Must give at least one reaction before constructing a JumpSystem.")
    for rx in equations(rs)
        empty!(rxvars)
        (rx.rate isa Symbolic) && get_variables!(rxvars, rx.rate)
        haveivdep = false
        @inbounds for i = 1:length(rxvars)
            if isequal(rxvars[i], get_iv(rs))
                haveivdep = true
                break
            end
        end
        if ismassaction(rx, rs; rxvars=rxvars, haveivdep=haveivdep, stateset=stateset)
            push!(meqs, makemajump(rx, combinatoric_ratelaw=combinatoric_ratelaws))
        else
            rl     = jumpratelaw(rx, rxvars=rxvars, combinatoric_ratelaw=combinatoric_ratelaws)
            affect = Vector{Equation}()
            for (spec,stoich) in rx.netstoich
                push!(affect, spec ~ spec + stoich)
            end
            if haveivdep
                push!(veqs, VariableRateJump(rl,affect))
            else
                push!(ceqs, ConstantRateJump(rl,affect))
            end
        end
    end
    #eqs[1] = MassActionJump(rates, rstoich, nstoich, scale_rates=false, useiszero=false)
    ArrayPartition(meqs,ceqs,veqs)
end

"""
```julia
Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
```
Convert a [`ReactionSystem`](@ref) to an [`ODESystem`](@ref).

Notes:
- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate
law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
`combinatoric_ratelaws=false` then the ratelaw is `k*S^2`, i.e. the scaling factor is
ignored.
"""
function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem;
                      name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
    eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws)
    systems = map(sys -> (sys isa ODESystem) ? sys : convert(ODESystem, sys), get_systems(rs))
    ODESystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
end

"""
```julia
Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
```

Convert a [`ReactionSystem`](@ref) to an [`NonlinearSystem`](@ref).

Notes:
- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate
law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
`combinatoric_ratelaws=false` then the ratelaw is `k*S^2`, i.e. the scaling factor is
ignored.
"""
function Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem;
                      name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
    eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, as_odes=false)
    systems = convert.(NonlinearSystem, get_systems(rs))
    NonlinearSystem(eqs, get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
end

"""
```julia
Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
```

Convert a [`ReactionSystem`](@ref) to an [`SDESystem`](@ref).

Notes:
- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate
law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
`combinatoric_ratelaws=false` then the ratelaw is `k*S^2`, i.e. the scaling factor is
ignored.
- `noise_scaling=nothing::Union{Vector{Operation},Operation,Nothing}` allows for linear
scaling of the noise in the chemical Langevin equations. If `nothing` is given, the default
value as in Gillespie 2000 is used. Alternatively, an `Operation` can be given, this is
added as a parameter to the system (at the end of the parameter array). All noise terms
are linearly scaled with this value. The parameter may be one already declared in the `ReactionSystem`.
Finally, a `Vector{Operation}` can be provided (the length must be equal to the number of reactions).
Here the noise for each reaction is scaled by the corresponding parameter in the input vector.
This input may contain repeat parameters.
"""
function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
                      noise_scaling=nothing, name=nameof(rs), combinatoric_ratelaws=true, kwargs...)

    if noise_scaling isa Vector
        (length(noise_scaling)!=length(equations(rs))) &&
        error("The number of elements in 'noise_scaling' must be equal " *
              "to the number of reactions in the reaction system.")
        noise_scaling = value.(noise_scaling)
    elseif !isnothing(noise_scaling)
        noise_scaling = fill(value(noise_scaling),length(equations(rs)))
    end

    eqs      = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws)
    noiseeqs = assemble_diffusion(rs,noise_scaling;
                                  combinatoric_ratelaws=combinatoric_ratelaws)
    systems  = convert.(SDESystem, get_systems(rs))
    SDESystem(eqs, noiseeqs, get_iv(rs), get_states(rs),
              (noise_scaling===nothing) ? get_ps(rs) : union(get_ps(rs), toparam.(noise_scaling));
              name=name,
              systems=systems,
              kwargs...)
end

"""
```julia
Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
```

Convert a [`ReactionSystem`](@ref) to an [`JumpSystem`](@ref).

Notes:
- `combinatoric_ratelaws=true` uses binomials in calculating the rate law, i.e. for `2S ->
  0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaws=false` then
  the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling
  factor.
"""
function Base.convert(::Type{<:JumpSystem},rs::ReactionSystem;
                      name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
    eqs     = assemble_jumps(rs; combinatoric_ratelaws=combinatoric_ratelaws)
    systems = convert.(JumpSystem, get_systems(rs))
    JumpSystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
end

### Converts a reaction system to ODE or SDE problems ###

# ODEProblem from AbstractReactionNetwork
function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, p=DiffEqBase.NullParameters(), args...; kwargs...)
    return ODEProblem(convert(ODESystem,rs; kwargs...),u0,tspan,p, args...; kwargs...)
end

# NonlinearProblem from AbstractReactionNetwork
function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0, p=DiffEqBase.NullParameters(), args...; kwargs...)
    return NonlinearProblem(convert(NonlinearSystem,rs; kwargs...), u0, p, args...; kwargs...)
end

# SDEProblem from AbstractReactionNetwork
function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, p=DiffEqBase.NullParameters(), args...; noise_scaling=nothing, kwargs...)
    sde_sys  = convert(SDESystem,rs;noise_scaling=noise_scaling, kwargs...)
    p_matrix = zeros(length(get_states(rs)), length(get_eqs(rs)))
    return SDEProblem(sde_sys,u0,tspan,p,args...; noise_rate_prototype=p_matrix,kwargs...)
end

# DiscreteProblem from AbstractReactionNetwork
function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple, p=DiffEqBase.NullParameters(), args...; kwargs...)
    return DiscreteProblem(convert(JumpSystem,rs; kwargs...), u0,tspan,p, args...; kwargs...)
end

# JumpProblem from AbstractReactionNetwork
function DiffEqJump.JumpProblem(rs::ReactionSystem, prob, aggregator, args...; kwargs...)
    return JumpProblem(convert(JumpSystem,rs; kwargs...), prob, aggregator, args...; kwargs...)
end

# SteadyStateProblem from AbstractReactionNetwork
function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0, p=DiffEqBase.NullParameters(), args...; kwargs...)
    return SteadyStateProblem(ODEFunction(convert(ODESystem,rs; kwargs...)),u0,p, args...; kwargs...)
end

# determine which species a reaction depends on
function get_variables!(deps::Set, rx::Reaction, variables)
    (rx.rate isa Symbolic) && get_variables!(deps, rx.rate, variables)
    for s in rx.substrates
        push!(deps, s)
    end
    deps
end

# determine which species a reaction modifies
function modified_states!(mstates, rx::Reaction, sts::Set)
    for (species,stoich) in rx.netstoich
        (species in sts) && push!(mstates, species)
    end
end


text:
defaultOutputFormat(plt::Plot) = "png"

function png(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "png")
  io = open(fn, "w")
  show(io, MIME("image/png"), plt)
  close(io)
end
png(fn::AbstractString) = png(current(), fn)

function svg(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "svg")
  io = open(fn, "w")
  show(io, MIME("image/svg+xml"), plt)
  close(io)
end
svg(fn::AbstractString) = svg(current(), fn)

function pdf(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "pdf")
  io = open(fn, "w")
  show(io, MIME("application/pdf"), plt)
  close(io)
end
pdf(fn::AbstractString) = pdf(current(), fn)

function ps(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "ps")
  io = open(fn, "w")
  show(io, MIME("application/postscript"), plt)
  close(io)
end
ps(fn::AbstractString) = ps(current(), fn)

function eps(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "eps")
  io = open(fn, "w")
  show(io, MIME("image/eps"), plt)
  close(io)
end
eps(fn::AbstractString) = eps(current(), fn)

function tex(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "tex")
  io = open(fn, "w")
  show(io, MIME("application/x-tex"), plt)
  close(io)
end
tex(fn::AbstractString) = tex(current(), fn)

function html(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "html")
    io = open(fn, "w")
    _use_remote[] = true
    show(io, MIME("text/html"), plt)
    _use_remote[] = false
    close(io)
end
html(fn::AbstractString) = html(current(), fn)

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

const _savemap = Dict(
    "png" => png,
    "svg" => svg,
    "pdf" => pdf,
    "ps"  => ps,
    "eps" => eps,
    "tex" => tex,
    "html" => html,
  )

function getExtension(fn::AbstractString)
  pieces = split(fn, ".")
  length(pieces) > 1 || error("Can't extract file extension: ", fn)
  ext = pieces[end]
  haskey(_savemap, ext) || error("Invalid file extension: ", fn)
  ext
end

function addExtension(fn::AbstractString, ext::AbstractString)
  try
    oldext = getExtension(fn)
    if oldext == ext
      return fn
    else
      return "$fn.$ext"
    end
  catch
    return "$fn.$ext"
  end
end

"""
    savefig([plot,] filename)

Save a Plot (the current plot if `plot` is not passed) to file. The file
type is inferred from the file extension. All backends support png and pdf
file types, some also support svg, ps, eps, html and tex.
"""
function savefig(plt::Plot, fn::AbstractString)

  # get the extension
  local ext
  try
    ext = getExtension(fn)
  catch
    # if we couldn't extract the extension, add the default
    ext = defaultOutputFormat(plt)
    fn = addExtension(fn, ext)
  end

  # save it
  func = get(_savemap, ext) do
    error("Unsupported extension $ext with filename ", fn)
  end
  func(plt, fn)
end
savefig(fn::AbstractString) = savefig(current(), fn)

# ---------------------------------------------------------
"""
    gui([plot])

Display a plot using the backends' gui window
"""
gui(plt::Plot = current()) = display(PlotsDisplay(), plt)

# IJulia only... inline display
function inline(plt::Plot = current())
    isijulia() || error("inline() is IJulia-only")
    Main.IJulia.clear_output(true)
    display(Main.IJulia.InlineDisplay(), plt)
end

function Base.display(::PlotsDisplay, plt::Plot)
    prepare_output(plt)
    _display(plt)
end

# override the REPL display to open a gui window
Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", plt::Plot) = gui(plt)

_do_plot_show(plt, showval::Bool) = showval && gui(plt)
function _do_plot_show(plt, showval::Symbol)
    showval == :gui && gui(plt)
    showval in (:inline,:ijulia) && inline(plt)
end

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

const _mimeformats = Dict(
    "application/eps"         => "eps",
    "image/eps"               => "eps",
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg",
    "text/plain"              => "txt",
    "application/x-tex"       => "tex",
)

const _best_html_output_type = KW(
    :pyplot => :png,
    :unicodeplots => :txt,
    :glvisualize => :png
)

# a backup for html... passes to svg or png depending on the html_output_format arg
function Base.show(io::IO, ::MIME"text/html", plt::Plot)
    output_type = Symbol(plt.attr[:html_output_format])
    if output_type == :auto
        output_type = get(_best_html_output_type, backend_name(plt.backend), :svg)
    end
    if output_type == :png
        # info("writing png to html output")
        print(io, "")
    elseif output_type == :svg
        # info("writing svg to html output")
        show(io, MIME("image/svg+xml"), plt)
    elseif output_type == :txt
        show(io, MIME("text/plain"), plt)
    else
        error("only png or svg allowed. got: $output_type")
    end
end

function _show(io::IO, m, plt::Plot{B}) where B
    # Base.show_backtrace(STDOUT, backtrace())
    warn("_show is not defined for this backend. m=", string(m))
end
function _display(plt::Plot)
    warn("_display is not defined for this backend.")
end

# for writing to io streams... first prepare, then callback
for mime in keys(_mimeformats)
    @eval function Base.show(io::IO, m::MIME{Symbol($mime)}, plt::Plot{B}) where B
        prepare_output(plt)
        _show(io, m, plt)
    end
end

"Close all open gui windows of the current backend"
closeall() = closeall(backend())

# ---------------------------------------------------------
# A backup, if no PNG generation is defined, is to try to make a PDF and use FileIO to convert

if is_installed("FileIO")
    @eval import FileIO
    function _show(io::IO, ::MIME"image/png", plt::Plot)
        fn = tempname()

        # first save a pdf file
        pdf(plt, fn)

        # load that pdf into a FileIO Stream
        s = FileIO.load(fn * ".pdf")

        # save a png
        pngfn = fn * ".png"
        FileIO.save(pngfn, s)

        # now write from the file
        write(io, readstring(open(pngfn)))
    end
end

# function html_output_format(fmt)
#     if fmt == "png"
#         @eval function Base.show(io::IO, ::MIME"text/html", plt::Plot)
#             print(io, "")
#         end
#     elseif fmt == "svg"
#         @eval function Base.show(io::IO, ::MIME"text/html", plt::Plot)
#             show(io, MIME("image/svg+xml"), plt)
#         end
#     else
#         error("only png or svg allowed. got: $fmt")
#     end
# end
#
# html_output_format("svg")

# ---------------------------------------------------------
# IJulia
# ---------------------------------------------------------

const _ijulia_output = String["text/html"]

@require IJulia begin
    if IJulia.inited
        export set_ijulia_output

        function set_ijulia_output(mimestr::AbstractString)
            # info("Setting IJulia output format to $mimestr")
            global _ijulia_output
            _ijulia_output[1] = mimestr
        end
        function IJulia.display_dict(plt::Plot)
            global _ijulia_output
            Dict{String, String}(_ijulia_output[1] => sprint(show, _ijulia_output[1], plt))
        end

        # default text/plain passes to html... handles Interact issues
        function Base.show(io::IO, m::MIME"text/plain", plt::Plot)
            show(io, MIME("text/html"), plt)
        end

        ENV["MPLBACKEND"] = "Agg"
        set_ijulia_output("text/html")
    end
end

# ---------------------------------------------------------
# Atom PlotPane
# ---------------------------------------------------------
@require Juno begin
    import Hiccup, Media

    if Juno.isactive()
        Media.media(Plot, Media.Plot)

        _show(io::IO, m::MIME"text/plain", plt::Plot{B}) where {B} = print(io, "Plot{$B}()")

        function Juno.render(e::Juno.Editor, plt::Plot)
            Juno.render(e, nothing)
        end

        if get(ENV, "PLOTS_USE_ATOM_PLOTPANE", true) in (true, 1, "1", "true", "yes")
            function Juno.render(pane::Juno.PlotPane, plt::Plot)
                # temporarily overwrite size to be Atom.plotsize
                sz = plt[:size]
                jsize = Juno.plotsize()
                jsize[1] == 0 && (jsize[1] = 400)
                jsize[2] == 0 && (jsize[2] = 500)

                plt[:size] = jsize
                Juno.render(pane, HTML(stringmime(MIME("text/html"), plt)))
                plt[:size] = sz
            end
            # special handling for PlotlyJS
            function Juno.render(pane::Juno.PlotPane, plt::Plot{PlotlyJSBackend})
                display(Plots.PlotsDisplay(), plt)
            end
        else
            function Juno.render(pane::Juno.PlotPane, plt::Plot)
                display(Plots.PlotsDisplay(), plt)
                s = "PlotPane turned off.  Unset ENV[\"PLOTS_USE_ATOM_PLOTPANE\"] and restart Julia to enable it."
                Juno.render(pane, HTML(s))
            end
        end

        # special handling for plotly... use PlotsDisplay
        function Juno.render(pane::Juno.PlotPane, plt::Plot{PlotlyBackend})
            display(Plots.PlotsDisplay(), plt)
            s = "PlotPane turned off.  The plotly backend cannot render in the PlotPane due to javascript issues. Plotlyjs is similar to plotly and is compatible with the plot pane."
            Juno.render(pane, HTML(s))
        end
    end
end


text:
# Autogenerated wrapper script for NodeJS_16_jll for armv7l-linux-gnueabihf-cxx11
export node, npm, npx

JLLWrappers.@generate_wrapper_header("NodeJS_16")
JLLWrappers.@declare_executable_product(node)
JLLWrappers.@declare_file_product(npm)
JLLWrappers.@declare_file_product(npx)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        node,
        "bin/node",
    )

    JLLWrappers.@init_file_product(
        npm,
        "bin/npm",
    )

    JLLWrappers.@init_file_product(
        npx,
        "bin/npx",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()