text:
module DelaunayMaps

using Random
using LightGraphs
using MetaGraphs
using Plots
using LinearAlgebra
using DataStructures
using Logging

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

"Compute the intersection between a + λ * dir_a and b + γ * dir_b"
function line_intersection(a::Array{<:Real,1}, dir_a::Array{<:Real,1}, b::Array{<:Real,1}, dir_b::Array{<:Real,1})
    λ = [-dir_a dir_b]\(a .- b)
    a .+ λ[1] .* dir_a
end

line_intersection(points::Array{<:Real,1}, a::Int, b::Int, c::Int, d::Int) = line_intersection(points[a,:], points[b,:], points[c,:], points[d,:])

function circumcenter(points, a, b, c)
    xa = points[a, 1]
    ya = points[a, 2]
    xb = points[b, 1]
    yb = points[b, 2]
    xc = points[c, 1]
    yc = points[c, 2]
    M = [
        yb-ya  yb-yc
        xa-xb  xc-xb
    ]
    B = 1/2 * [xa-xc;ya-yc]
    λ = M\B
    1/2 * (points[a,:] + points[b,:]) + λ[1] * [ya-yb; xb-xa]
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
        t
    end
    function Triangle(v1::Int, v2::Int, v3::Int, t1::Triangle, t2::Triangle, t3::Triangle)
        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))")

circumcenter(points, triangle::Triangle) = circumcenter(points, triangle.v1, triangle.v2, triangle.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
    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
        neighbour.t1 = comp
    elseif triangle == neighbour.t2
        neighbour.t2 = comp
    else
        neighbour.t3 = comp
    end
    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)

mutable struct Map{T<:Real}
    delaunay::MetaGraph
    delaunay_points::Array{T,2}
    voronoi::MetaGraph
    voronoi_points::Array{T,2}
    width::Real
    height::Real
end

"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!(map::Map, map_triangles::Set{Triangle}, to_be_deleted::Deque{Triangle}, upper_triangle::Triangle)
    next_right = next(upper_triangle)
    next_next_right = next(next_right)
    @debug "Starting deletion routine on right."
    for triangle in to_be_deleted
        pop!(map_triangles, triangle)
        @debug "Processing $triangle"
        if triangle == next_right # first iteration of the loop
            @debug "First iteration."
            continue
        elseif is_hull(triangle)
            @debug "Is hull, should be last iteration."
            continue
        elseif triangle.t2 == next_right
            if triangle.t2 ∈ map_triangles # triangle.t2 was added during last iteration
                pop!(map_triangles, triangle.t2)
            end
            @debug "triangle.t2 is current next_right."
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            comp_t3 = Triangle(triangle.v2, triangle.v1)
            @debug "Created $comp_t1 and $comp_t3."
            push!(map_triangles, comp_t1)
            push!(map_triangles, comp_t3)
            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!(map.delaunay, triangle.v1, triangle.v3)
        elseif triangle.t3 == next_right
            if triangle.t3 ∈ map_triangles # triangle.t3 was added during last iteration
                pop!(map_triangles, triangle.t3)
            end
            @debug "triangle.t3 is current next_right."
            comp_t2 = Triangle(triangle.v1, triangle.v3)
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            @debug "Created $comp_t1 and $comp_t2."
            push!(map_triangles, comp_t1)
            push!(map_triangles, comp_t2)
            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!(map.delaunay, triangle.v1, triangle.v2)
        else
            error("I am removing t1 on a right triangle.")
        end
    end
    # We might have deleted the next right triangle, especially if it was on hull
    push!(map_triangles, next(upper_triangle))
    @debug "Right deletion routine is over."
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!(map::Map, map_triangles::Set{Triangle}, to_be_deleted::Deque{Triangle}, upper_triangle::Triangle)
    prev_left = prev(upper_triangle)
    prev_prev_left = prev(prev_left)
    @debug "Starting deletion routine on left."
    for triangle in to_be_deleted
        pop!(map_triangles, triangle)
        @debug "Processing $triangle"
        if triangle == prev_left # first iteration of the loop
            @debug "First iteration."
            continue
        elseif is_hull(triangle)
            @debug "Is hull, should be last iteration."
            continue
        elseif triangle.t1 == prev_left
            if triangle.t1 ∈ map_triangles # triangle.t1 was added during last iteration
                pop!(map_triangles, triangle.t1)
            end
            @debug "triangle.t1 is current prev_left."
            comp_t3 = Triangle(triangle.v2, triangle.v1)
            comp_t2 = Triangle(triangle.v1, triangle.v3)
            @debug "Created $comp_t2 and $comp_t3."
            push!(map_triangles, comp_t2)
            push!(map_triangles, comp_t3)
            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!(map.delaunay, triangle.v2, triangle.v3)
        elseif triangle.t2 == prev_left
            if triangle.t2 ∈ map_triangles # triangle.t2 was added during last iteration
                pop!(map_triangles, triangle.t2)
            end
            @debug "triangle.t2 is current prev_left."
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            comp_t3 = Triangle(triangle.v2, triangle.v1)
            @debug "Created $comp_t1 and $comp_t3."
            push!(map_triangles, comp_t1)
            push!(map_triangles, comp_t3)
            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!(map.delaunay, triangle.v1, triangle.v3)
        else # triangle.t3 == prev_left
            if triangle.t3 ∈ map_triangles # triangle.t3 was added during last iteration
                pop!(map_triangles, triangle.t3)
            end
            @debug "triangle.t3 is current prev_left."
            comp_t2 = Triangle(triangle.v1, triangle.v3)
            comp_t1 = Triangle(triangle.v3, triangle.v2)
            @debug "Created $comp_t1 and $comp_t2."
            push!(map_triangles, comp_t1)
            push!(map_triangles, comp_t2)
            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!(map.delaunay, triangle.v1, triangle.v2)
        end
    end
    # We might have deleted the next right triangle, especially if it was on hull
    push!(map_triangles, prev(upper_triangle))
    @debug "Left deletion routine is over."
end

function find_base_tangente(points, triangle_left, triangle_right)
    @debug "Starting find_base_tangente routine."
    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
        @debug "Current tangente" triangle_left.v1 triangle_right.v1
        @debug "Working on" prev_triangle_left triangle_left prev_triangle_right triangle_right
        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
    @debug "End of find_base_tangente routine."
    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.delaunay_points, triangle, base_left, base_right)
    else
        triangle -> _in_circumcircle_right(m.delaunay_points, triangle, base_left, base_right)
    end
    _is_above_base = if mode == left
        triangle -> _is_above_base_left(triangle, m.delaunay_points, base_left, base_right)
    else
        triangle -> _is_above_base_right(triangle, m.delaunay_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 # Never happens if the hull is correctly build
        a = left_candidate.v2
    elseif is_hull(left_candidate) && left_candidate.v2 == b
        a = left_candidate.v1
    elseif left_candidate.v2 == b
        a = left_candidate.v3
    elseif left_candidate.v1 == b
        a = left_candidate.v2
    else # left_candidate.v3 == b
        a = left_candidate.v1
    end
    if is_hull(right_candidate) && right_candidate.v2 == c  # Never happens if the hull is correctly build
        d = right_candidate.v1
    elseif is_hull(right_candidate) && right_candidate.v1 == c
        d = right_candidate.v2
    elseif right_candidate.v1 == c
        d = right_candidate.v3
    elseif right_candidate.v2 == c
        d = right_candidate.v1
    else # right_candidate.v3 == c
        d = right_candidate.v2
    end
    @debug "Have to choose between $left_candidate, $right_candidate with test is_in_circumcircle(points, $a, $b, $c, $d)"
    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, map_triangles=Set{Triangle}(), first=1, last=-1)
    if last == -1
        last = div(length(m.delaunay_points), 2)
    end
    if first == last # never happens
        error("Trying to triangulate a single point $first.")
    elseif last - first == 1
        @debug "Triangulate on 2 vertices : $first - $last."
        add_edge!(m.delaunay, first, last)
        t = Triangle(first, last)
        push!(map_triangles, t)
        push!(map_triangles, t.t3)
        @debug "Created $t."
        @debug to_geogebra(m, :novoronoi)
        t.t3, t
    elseif last - first == 2
        @debug "Triangulate on 3 vertices : $first to $last."
        if !is_aligned(m.delaunay_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.delaunay_points, first, first+1, last)
            push!(map_triangles, t)
            push!(map_triangles, t.t1)
            push!(map_triangles, t.t2)
            push!(map_triangles, t.t3)
            @debug "Created $t."
            @debug to_geogebra(m, :novoronoi)
            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)
            push!(map_triangles, t_first)
            push!(map_triangles, t_first.t3)
            push!(map_triangles, t_last)
            push!(map_triangles, t_last.t3)
            @debug "Created $t_first and $t_last."
            @debug to_geogebra(m, :novoronoi)
            if m.delaunay_points[first, 1] < m.delaunay_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, map_triangles, first, mid)
        triangle_right_left, triangle_right_right = triangulate!(m, map_triangles, mid+1, last)

        @debug "Triangulate on $(last - first + 1) vertices : $first to $last."

        base_triangle_left, base_triangle_right = find_base_tangente(m.delaunay_points, triangle_left_right, triangle_right_left)
        base_left = base_triangle_left.v1
        base_right = base_triangle_right.v1

        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)
        push!(map_triangles, base_triangle)
        next!(base_triangle, base_triangle_left)
        prev!(base_triangle, prev(base_triangle_right))
        lower_triangle = base_triangle

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

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

        base_triangle_edge = :t3

        add_edge!(m.delaunay, base_left, base_right)

        while !isnothing(right_candidate) || !isnothing(left_candidate)
            @debug "Working with" base_triangle upper_triangle
            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)

            @debug "Left candidate is " left_candidate
            @debug "Right candidate is " right_candidate

            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.delaunay_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
                @debug "Chose $chosen."

                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, map_triangles, 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, map_triangles, 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.delaunay_points, triangle.v1, triangle.v2, triangle.v3)
                    triangle.v2, triangle.v3 = triangle.v3, triangle.v2
                end

                # 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
                        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
            end
        end

        @debug to_geogebra(m, :novoronoi)
        lower_triangle, upper_triangle, map_triangles
    end

end

function project_bissection_on_nearest_border(points::Array{<:Real,2}, triangle::Triangle, width::Real, height::Real)
    center = (points[triangle.v1,:] .+ points[triangle.v2,:]) ./ 2
    direction = [0 -1;1 0] * (points[triangle.v1,:] .- points[triangle.v2,:])
    sort([
            line_intersection(center, direction, [0;0], [0;height]),
            line_intersection(center, direction, [width;0], [0;height]),
            line_intersection(center, direction, [0;0], [width;0]),
            line_intersection(center, direction, [0;height], [width;0])
        ],
        lt=(x,y)-> norm(center .- x) < norm(center .- y)
    )[1]
end

function create_voronoi_from_triangles!(map::Map, map_triangles::Set{Triangle})
    processed = Set{Triangle}()
    stack = Stack{Triangle}()
    indices = Dict()
    for (i,triangle) in enumerate(map_triangles)
        push!(indices, triangle => i)
    end
    push!(stack, iterate(indices)[1].first)
    map.voronoi_points = zeros(Float64, length(indices), 2)
    map.voronoi = MetaGraph(length(indices))
    while !isempty(stack)
        triangle = pop!(stack)
        push!(processed, triangle)
        index = get(indices, triangle, 0)
        map.voronoi_points[index,:] = if is_hull(triangle)
            project_bissection_on_nearest_border(map.delaunay_points, triangle, map.width, map.height)
        else
            max.([0; 0], min.([map.width; map.height], circumcenter(map.delaunay_points, triangle)))
        end
        for t_neighbour in [triangle.t1, triangle.t2, triangle.t3]
            if t_neighbour ∉ processed
                if !(is_hull(triangle) && is_hull(t_neighbour))
                    t_index = get(indices, t_neighbour, 0)
                    add_edge!(map.voronoi, index, t_index)
                    set_prop!(map.voronoi, index, t_index, :delaunay_edge, Edge(triangle.v2, triangle.v3))
                    set_prop!(map.delaunay, triangle.v2, triangle.v3, :voronoi_edge, Edge(index, t_index))
                end
                push!(stack, t_neighbour)
            end
        end
    end
end

function Map(points::Array{T,2}) where T <: Real
    m = Map(MetaGraph(div(length(points),2)), sortslices(points, dims=1), MetaGraph(), zeros(0,0), maximum(points[:,1]), maximum(points[:,2]))
    map_triangles = Set{Triangle}()
    triangulate!(m, map_triangles)
    create_voronoi_from_triangles!(m, map_triangles)
    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)
    @debug "Generating from" points
    delaunay = MetaGraph(nb_points)
    m = Map(delaunay, points, MetaGraph(), zeros(0,0), width, height)
    map_triangles = Set{Triangle}()
    triangulate!(m, map_triangles)
    create_voronoi_from_triangles!(m, map_triangles)
    m
end

function plot_points(points, type=:delaunay; show_nums=true)
    plot()
    plot_points!(points, type, show_nums=show_nums)
end

function plot_points!(points, type=:delaunay; show_nums=true)
    marker = if type == :delaunay
        if show_nums
            (:circle, 20, 0.2, :orange)
        else
            (:circle, 100, 1, :red)
        end
    else
        if show_nums
            (:circle, 10, 0.2, :blue)
        else
            (:circle, 100, 1, :blue)
        end
    end
    if show_nums
        plot!(
            points[:,1], points[:,2],
            marker=marker, line=:scatter,
            series_annotations=collect(map(string, 1:length(points[:,1]))),
            aspect_ratio=:equal
        )
    else
        plot!(
            points[:,1], points[:,2],
            marker=marker, line=:scatter, aspect_ratio=:equal
        )
    end
end

function plot_map(m; show_nums=true, window=:fit)
    plot()
    plot_map!(m, show_nums=show_nums, window=window)
end
function plot_map!(m; show_nums=true, window=:fit)
    if window == :full
        plot!(
            aspect_ratio=:equal
        )
    else
        plot!(
            aspect_ratio=:equal,
            xlims=(0,m.width),
            ylims=(0,m.height),
        )
    end
    plot!(
        hcat(map(x->[m.delaunay_points[x.src, 1]; m.delaunay_points[x.dst, 1]],edges(m.delaunay))...),
        hcat(map(x->[m.delaunay_points[x.src, 2]; m.delaunay_points[x.dst, 2]],edges(m.delaunay))...),
        color="red",
        legend=false,
    )
    plot!(
        hcat(map(x->[m.voronoi_points[x.src, 1]; m.voronoi_points[x.dst, 1]],edges(m.voronoi))...),
        hcat(map(x->[m.voronoi_points[x.src, 2]; m.voronoi_points[x.dst, 2]],edges(m.voronoi))...),
        color="blue",
        legend=false,
    )
    plot_points!(m.voronoi_points, :voronoi, show_nums=show_nums)
    plot_points!(m.delaunay_points, show_nums=show_nums)
end

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

function to_geogebra(map::Map, voronoi=:voronoi)
    convert_to_ascii(n) = Char(n + 64)
    n_delaunay = div(length(map.delaunay_points), 2)
    if voronoi == :voronoi
    join(
        [
            to_geogebra(map.delaunay_points),
            (*)(
                "{",
                join(["Segment[$(convert_to_ascii(e.src)), $(convert_to_ascii(e.dst))]" for e in edges(map.delaunay)], ", "),
                "}"
            ),
            (*)(
                "{",
                join(["Polygone[$(convert_to_ascii(t.v1)), $(convert_to_ascii(t.v2)), $(convert_to_ascii(t.v3))]" for t in map.triangles if !is_hull(t)], ", "),
                "}"
            ),
            to_geogebra(map.voronoi_points),
            (*)(
                "{",
                join(["Segment[$(convert_to_ascii(e.src+n_delaunay)), $(convert_to_ascii(e.dst+n_delaunay))]" for e in edges(map.voronoi)], ", "),
                "}"
            ),
        ],
        "\n"
    )
    else
    join(
        [
            to_geogebra(map.delaunay_points),
            (*)(
                "{",
                join(["Segment[$(convert_to_ascii(e.src)), $(convert_to_ascii(e.dst))]" for e in edges(map.delaunay)], ", "),
                "}"
            ),
            (*)(
                "{",
                join(["Polygone[$(convert_to_ascii(t.v1)), $(convert_to_ascii(t.v2)), $(convert_to_ascii(t.v3))]" for t in map.triangles if !is_hull(t)], ", "),
                "}"
            ),
        ],
        "\n"
    )
    end
end

"Dummy function to test if a Map is delaunay."
function is_delaunay(map)
    result = true
    for p in vertices(map.delaunay)
        n = neighbors(map.delaunay, p)
        for q in vertices(map.delaunay)
            if q in n
                continue
            end
            ok = !is_in_circumcircle(map.delaunay_points, t.v1, t.v3, t.v2, p)
            if !ok
                @debug "failed on triangle $t point $p"
                result = false
            end
        end
    end
    result
end

export Map, to_geogebra, plot_points, plot_points!, plot_map, plot_map!

end # module


text:
using VLConstraintBasedModelGenerationUtilities

# setup path to protein sequence file -
path_to_vff_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/Test.vff"
path_to_system_model_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/Test.bson"

# let's build the reaction table -
metabolic_reaction_table = build_metabolic_reaction_table(path_to_vff_file) |> check

# build the stm -
stm = build_stoichiometric_matrix(metabolic_reaction_table) |> check

# build the species bounds array -
species_table = build_species_table(metabolic_reaction_table) |> check
species_bounds_array = build_species_bounds_array(species_table) |> check

# build the flux bounds array -
flux_bounds_array = build_flux_bounds_array(metabolic_reaction_table) |> check

# write the system model file -
result = write_system_model_file(path=path_to_system_model_file, stoichiometric_matrix=stm,
    flux_bounds_array=flux_bounds_array, species_bounds_array=species_bounds_array)

ext: jl


text:
# SPDX-License-Identifier: MIT

"""
    abstract type Distribution end

Supertype for distributions member of the exponential family.
"""
abstract type Distribution end

function Base.getproperty(dist::Distribution, sym::Symbol)
    # Check first the fields to avoid computing the std-params for
    # every property access.
    if sym in fieldnames(typeof(dist))
        return getfield(dist, sym)
    end

    param = stdparam(dist)
    if sym ∈ keys(param)
        return getfield(param, sym)
    end

    # Raise an error.
    getfield(dist, sym)
end

function Base.show(io::IO, ::MIME"text/plain", dist::Distribution)
    print(io, typeof(dist), ":")

    params = stdparam(dist)
    for (prop, val) in zip(keys(params), params)
        println(io)
        print(io, "  $prop = $val")
    end
end

"""
    loglikelihood(p, x)

Returns the log-likelihood of `x` for the distribution `p`.
"""
function loglikelihood(p, x)
    Tx = stats(p, x)
    η = naturalparam(p)
    dot(η, Tx) - lognorm(p) + basemeasure(p, x)
end

"""
    basemeasure(p, x)

Returns the base measure of `x` for the distribution `p`.
"""
basemeasure

"""
    gradlognorm(p)

Returns the gradient of the log-normalizer of `p` w.r.t. its natural
parameters.
"""
gradlognorm

"""
    kldiv(q::T, p::T[, μ = gradlognorm(q)]) where T<:Distribution

Compute the KL-divergence between two distributions of the same type
(i.e. `kldiv(Normal, Normal)`, `kldiv(Dirichlet, Dirichlet)`, ...). You
can specify directly the expectation of the sufficient statistics `μ`.
"""
function kldiv(q::T, p::T; μ = gradlognorm(q)) where T<:Distribution
    q_η, p_η = naturalform(q.param), naturalform(p.param)
    lognorm(p) - lognorm(q) - dot(p_η .- q_η, μ)
end

"""
    lognorm(p)

Returns the log-normalization constant of `p`.
"""
lognorm

"""
    sample(p, n=1)

Draw `n` samples from the distribution `p`.
"""
sample

"""
    splitgrad(p, μ)

Split the gradient of the log-normalizer into its "standard" components.
For instance, for the Normal distribution, the output will be the
expected value of \$x\$ and \$xxᵀ\$.
"""
splitgrad

"""
    stats(p, x)

Returns the sufficient statistics of `x` for the distribution `p`.
"""
stats

"""
    stdparam(p, η)

Returns the standard parameters corresponding to the natural parameters
`η` for distributions with the same type of `p`.
"""
stdparam

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions


text:
my_f(x,y)=2x+3y

der_of_my_f(x,y) = ForwardDiff.derivative(x->my_f(x,y), x)


text:
module Test_conversions

using PeriodicSystems
using DescriptorSystems
using Symbolics
using Test
using LinearAlgebra

println("Test_conversions")

@testset "psmrc2d" begin

# example Longhi IEEETAC, 1994
sys = dss([0 -pi/6; pi/6 0],eye(2),eye(2),0)
@time psys = psmrc2d(sys,1; ki = [2,3], ko = [1,2])
@test psys.period == 6

@time psys = psmrc2d(sys,1);
@test psys.period == 1

@time psys = psmrc2d(sys,1, ki = [2,3]);
@test psys.period == 6

@time psys = psmrc2d(sys,1, ko = [1,2]);
@test psys.period == 2

sys = rss(2,2,2,disc=true)
@time psys = psmrc2d(sys,1; ki = [2,3], ko = [1,2])
@test psys.period == 6

end

@testset "psc2d" begin

# periodic time-varying matrices from IFAC2005 paper
@variables t
A = PeriodicSymbolicMatrix([cos(t) 1; 1 1-sin(t)],2*pi);
println("A(t) = $(A.F)")
B = HarmonicArray([0;1],[[1;0]],[[1;-1]],2*pi);
println("B(t) = $(convert(PeriodicSymbolicMatrix,B).F)")
C = PeriodicFunctionMatrix(t-> [sin(2*t)+cos(2*t) 1],pi);
println("C(t) = $(C.f(t))")
D = PeriodicFunctionMatrix(zeros(1,1),pi)
K = 100
Δ = 2*pi/K
ts = (0:K-1)*Δ

psysc = PeriodicStateSpace(convert(PeriodicFunctionMatrix,A),
                          convert(PeriodicFunctionMatrix,B),
                          convert(PeriodicFunctionMatrix,C),
                          convert(PeriodicFunctionMatrix,D));

@time psys = psc2d(psysc,Δ);
@time sys = psaverage(psysc);

psysc1 = PeriodicStateSpace(convert(PeriodicSymbolicMatrix,A),
                          convert(PeriodicSymbolicMatrix,B),
                          convert(PeriodicSymbolicMatrix,C),
                          convert(PeriodicSymbolicMatrix,D));
@time psys1 = psc2d(psysc1,Δ);
@test all(psys.A.M .≈ psys1.A.M) && all(psys.B.M .≈ psys1.B.M) &&
      all(psys.C.M .≈ psys1.C.M) && all(psys.D.M .≈ psys1.D.M)
@time sys1 = psaverage(psysc1);
@test iszero(sys-sys1,atol=1.e-7)

psysc2 = PeriodicStateSpace(convert(HarmonicArray,A),
                          convert(HarmonicArray,B),
                          convert(HarmonicArray,C),
                          convert(HarmonicArray,D));
@time psys2 = psc2d(psysc2,Δ);
@test all(psys.A.M .≈ psys2.A.M) && all(psys.B.M .≈ psys2.B.M) &&
      all(psys.C.M .≈ psys2.C.M) && all(psys.D.M .≈ psys2.D.M)
@time sys2 = psaverage(psysc2);
@test iszero(sys-sys2,atol=1.e-7)

psysc3 = PeriodicStateSpace(convert(FourierFunctionMatrix,A),
                          convert(FourierFunctionMatrix,B),
                          convert(FourierFunctionMatrix,C),
                          convert(FourierFunctionMatrix,D));
@time psys3 = psc2d(psysc3,Δ);
@test all(psys.A.M .≈ psys3.A.M) && all(psys.B.M .≈ psys3.B.M) &&
      all(psys.C.M .≈ psys3.C.M) && all(psys.D.M .≈ psys3.D.M)
@time sys3 = psaverage(psysc3);
@test iszero(sys-sys3,atol=1.e-7)

At = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,A).f.(ts),A.period);
Bt = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,B).f.(ts),B.period);
Ct = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,C).f.(ts),C.period);
Dt = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,D).f.(ts),D.period);
psysc4 = PeriodicStateSpace(At, Bt, Ct, Dt);
@time psys4 = psc2d(psysc4,Δ);
@test norm(psys.A.M .- psys4.A.M) < 1.e-7 && norm(psys.B.M .- psys4.B.M) < 1.e-7 &&
      all(psys.C.M .≈ psys4.C.M) && all(psys.D.M .≈ psys4.D.M)
@time sys4 = psaverage(psysc4);
@test iszero(sys-sys4,atol=1.e-7)

# Pitelkau's example
ω = 0.00103448
period = 2*pi/ω
a = [  0            0     5.318064566757217e-02                         0
       0            0                         0      5.318064566757217e-02
       -1.352134256362805e-03   0             0    -7.099273035392090e-02
       0    -7.557182037479544e-04     3.781555288420663e-02             0
];
b = t->[0; 0; 0.1389735*10^-6*sin(ω*t); -0.3701336*10^-7*cos(ω*t)];
c = [1 0 0 0;0 1 0 0];
d = zeros(2,1);
psysc = ps(a,PeriodicFunctionMatrix(b,period),c,d);

K = 120;
@time psys = psc2d(psysc,period/K,reltol = 1.e-10);
@test all(abs.(pseig(psys.A)) .≈ 1)
@test psys.A.M[1] ≈ exp(a*period/K) && psys.C.M[1] == c && psys.D.M[1] == d

# lti system
sys = rss(3,2,3);
psysc = ps(sys,10);
@time sysd = c2d(sys,1)[1];
@time psys = psc2d(psysc,1);
@test psys.A.M[1] ≈ sysd.A && psys.B.M[1] == sysd.B && psys.C.M[1] == sysd.C && psys.D.M[1] == sysd.D

# only dynamic gains
psysc = ps(C);
Ts = pi/10;
@time psys = psc2d(psysc,Ts);
@test all([psys.D.M[i] == C.f((i-1)*Ts) for i in 1:10])

# A and B constant matrices
a = rand(2,2); b = rand(2,1); c = rand(1,2); d = rand(1,1);
sys = dss(a,b,c,d);
psysc = ps(a,b,C,d);
Ts = pi/10;
@time sysd = c2d(sys,Ts)[1];
@time psys = psc2d(psysc,Ts);
@test psys.A.M[1] ≈ sysd.A && psys.B.M[1] == sysd.B &&
      all([psys.C.M[i] == C.f((i-1)*Ts) for i in 1:10]) && psys.D.M[1] == sysd.D

end # psc2d
end # module


text:
export Device, @setup_devices, alloc, name, arraytype, get_devices, get_device

abstract type Device end

name(device::Device) = ()

arraytype(device::Device) = ()

"Alloc to the memory type"
alloc(device::Device, tup::Tuple) = alloc(device, [i for i in tup])
alloc(device::Device, array::AbstractArray) = ()
alloc(device::Device, T::Type, indxs...) = ()

include("CPUImpl.jl")
include("CUDAImpl.jl")

const CPUDev = CPUDevice()

devices = Device[CPUDev]

get_devices() = devices

function get_device(device_name::Symbol)
    for i in 1:length(devices)
        (name(devices[i]) == device_name) && return devices[i]
    end
    return nothing
end

macro setup_devices()
    for device in devices
        deviceName = "$(device.name).jl"
        Core.eval(__module__, esc(:(include($deviceName))))
    end
end

function __init_device__()
    println("========Initializing HyperSphere Deep Learning======")
    (Threads.nthreads() < 8) && println("WARNING: THREADS LESS THEN 8. RECOMMENDED TO BE ATLEAST 8")
    println("CPU Thread Count:", Threads.nthreads())
    println("GPU Count:", CUDA.functional(false) ? length(CUDA.devices()) : 0)
    if CUDA.functional(false)
        push!(devices, CUDADevice())
    end
    println("Devices Available:$devices")
    println("==============")
end


text:
import AutomotiveVisualization.render
import AutomotiveVisualization.render!

function render_multipoint_drivers!(
    rendermodel::RenderModel,
    model::MultiPtsDriver,
    vehid::Int,
    car_colors::Dict{Int,Colorant}
    )

    add_instruction!(rendermodel, render_point_trail, (model.Pts,model.color,0.25))
    car_colors[vehid]=model.color
    rendermodel
end

function render_urban_drivers!(
    rendermodel::RenderModel,
    model::UrbanDriver,
    vehid::Int,
    car_colors::Dict{Int,Colorant}
    )

    car_colors[vehid]=model.color
    rendermodel
end
#=
function render_trafficlights!(rendermodel::RenderModel,trafficlights::Vector{TrafficLight})
    for trafficlight in trafficlights
        P = convert(VecE2,trafficlight.pos)
        radius = 1.0
        color_fill = colorant"green"
        if trafficlight.color == 1
            color_fill = colorant"red"
        elseif trafficlight.color == 2
            color_fill = colorant"yellow"
        end
        add_instruction!(rendermodel,render_circle,(P,radius,color_fill))
    end
    rendermodel
end
=#
function render(ctx::Cairo.CairoContext, scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway,models::Dict{Int, DriverModel};
    text::Vector{String}=["Nothing"],
    #trafficlights::Vector{TrafficLight}=TrafficLight[],
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{Int,Colorant}=Dict{Int,Colorant}()
    )

    canvas_width = floor(Int, Cairo.width(ctx))
    canvas_height = floor(Int, Cairo.height(ctx))

    clear_setup!(rendermodel)

    render!(rendermodel, roadway)

    for (i,veh) in enumerate(scene)
        if haskey(models,veh.id)
            model = models[veh.id]
            if model <: MultiPtsTurningDriver
                render_multipoint_drivers!(rendermodel,model,veh.id,car_colors)
            end
            if model <: UrbanDriver
                render_urban_drivers!(rendermodel,model,veh.id,car_colors)
            end
        end
    end

    render!(rendermodel, scene, car_colors=car_colors)

    textoverlay=TextOverlay(text,colorant"white",20,VecE2(10, 20),1.5,false)
    render!(rendermodel, textoverlay, scene, roadway)

    #render_trafficlights!(rendermodel,trafficlights)

    camera_set!(rendermodel, cam, scene, roadway, canvas_width, canvas_height)

    render(rendermodel, ctx, canvas_width, canvas_height)
    ctx
end

function render(scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway,models::Dict{Int, DriverModel};
    text::Vector{String}=["Nothing"],
    #trafficlights::Vector{TrafficLight}=TrafficLight[],
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{Int,Colorant}=Dict{Int,Colorant}(), # id
    )

    s, ctx = get_surface_and_context(canvas_width, canvas_height)
    render(ctx, scene, roadway,models, rendermodel=rendermodel, cam=cam, car_colors=car_colors,text=text#=,trafficlights=trafficlights=#)
    s
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
module BioLizard

include("l0.jl")
export ZeroNorm

include("ordinalregr.jl")
export OrdinalRegr, fitordinalsvm, score, predict

end # module


text:
# This file was generated, do not modify it. # hide
data = dataset("ISLR", "OJ")

feature_names = [
    :PriceCH, :PriceMM, :DiscCH, :DiscMM, :SalePriceMM, :SalePriceCH,
    :PriceDiff, :PctDiscMM, :PctDiscCH
]

X = select(data, feature_names);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#Load the Distributions package. Use `Pkg.install("Distributions")` to install first time.
using Distributions: TDist, ccdf

type regress_results
    coefs
    yhat
    res
    vcv
    tstat
    pval
end

# Keyword arguments are placed after semicolon.
# Symbols start with colon, e.g. `:symbol`.
function ols(y, X; corr=:none, lags::Int=Int(floor(size(X,1)^(1/4))))

    # β̂ = X \ y is more stable than β̂ = inv(X'*X) * X' \ y
    # see notes at bottom of case 1 notebook
    β̂ = X \ y
    ŷ = X * β̂
    μ̂ = y - ŷ

    T, K = size(X)
    σ̂² = dot(μ̂, μ̂) / (T - K)

    #use correction for variance covariance
    if corr == :none
        vcv = σ̂² * inv(X'*X)
    elseif corr == :white
        vcv = newey_west(X, μ̂, 0)
    elseif corr == :newey_west
        vcv = newey_west(X, μ̂, lags)
    else
        error("wrong argument for correction keyword")
    end

    # T statistics for H₀: βᵢ = 0
    tstat = β̂ ./ sqrt(diag(vcv))

    # absolute value and times two for double sided test
    pval  = 2 * ccdf(TDist(T-K), abs(tstat))

    regress_results(β̂, ŷ, μ̂, vcv, tstat, pval)
end

function newey_west(X, μ̂, lags::Integer)

    XtXInv = inv(X'*X)
    T, K = size(X)

    if lags==0 # White estimator
        return XtXInv * X' * diagm(μ̂.^2) * X * XtXInv
    end

    vcv = zeros(K, K)
    for t = 1:T
        vcv += μ̂[t]^2 * (X[t,:] * X[t,:]')
    end
    for lag in 1:lags
        w = 1 - lag / (lags + 1)
        for t in (lag + 1):T # Calculates the off-diagonal terms
            vcv += w * μ̂[t] * μ̂[t-lag] * (X[t-lag,:]*X[t,:]' + X[t,:]*X[t-lag,:]')
        end
    end
    vcv = XtXInv * vcv * XtXInv
end

function gls(y, X, Ω)

  P = chol(inv(Ω))
  return ols(P*y, P*X)

end

function gmm(y, X, Z; corr=:none, lags=nothing)

  T, Kx = size(X)
  T, Kz = size(Z)

  if corr==:none

    # Generalized 1-step IV estimator
    W = inv(Z'*Z)

  elseif corr==:white | corr==:newey_west

    if corr==:white
      lags=0
    end

    gmm1_res = gmm(y,X,Z;corr=:none)
    μ̂ = gmm1_res.res

    W = zeros(Kz, Kz)
    for lag in 0:lags
        w = 1 - lag / (lags + 1)
        for t in (lag + 1):T
            # Calculates the off-diagonal terms
            update = w * μ̂[t] * μ̂[t-lag] * (Z[t-lag,:]*Z[t,:]' + Z[t,:]*Z[t-lag,:]')
            W = W + update
        end
    end

  else
      error("wrong argument for correction keyword")
  end

  ZtX = Z'*X
  XtZ = X'*Z

  XtZ_W_ZtXInv = inv(XtZ*W*ZtX)
  β̂  = XtZ_W_ZtXInv*(XtZ*W*Z'*y)
  ŷ = X * β̂
  μ̂ = y - ŷ
  σ̂² = dot(μ̂, μ̂) / (T - Kz)
  vcv = σ̂² * XtZ_W_ZtXInv

  # T statistics for H₀: β₀ = 0
  tstat = β̂ ./ sqrt(diag(vcv))

  # absolute value and times two for double sided test
  pval  = 2 * ccdf(TDist(T-K), abs(tstat))

  return regress_results(β̂, ŷ, μ̂, vcv, tstat, pval)

end

max_stars_repo_name: Ken-B/AdvEconometrics

max_issues_repo_name: Ken-B/AdvEconometrics

max_forks_repo_name: Ken-B/AdvEconometrics


text:
using Serialization

function julia_func(z,frame)
    k = frame*0.002
    c = complex(-0.5, 0.66-k)
    z = z^2 + c
    return z
end

function gen_julia(julia_set,frame,im_width,im_height)

    max_it = 300
    max_z2 = 100
    min_x, max_x = -1, 1
    min_y, max_y = -1, 1

    for ix in 1:im_width
        for iy in 1:im_height
            real_part = ((max_x-min_x)/im_width)*ix + min_x
            im_part = ((max_y-min_y)/im_height)*iy + min_y
            z = complex(real_part, im_part)
            it = 0
            while(abs2(z) < max_z2 && it < max_it)
                z = julia_func(z,frame)
                it += 1
            ratio = it/max_it
            julia_set[ix,iy,frame] = ratio
            end
        end
    end
    return julia_set
end

function main()
    frames = 50
    im_width, im_height = 1000, 1000
    julia_set = zeros(im_width,im_height,frames)
    for frame in 1:frames
        gen_julia(julia_set,frame,im_width,im_height)
    end

    serialize("./data/julia_data",julia_set)
end

main()

max_stars_repo_name: ofbennett/julia-set-fractal

max_issues_repo_name: ofbennett/julia-set-fractal

max_forks_repo_name: ofbennett/julia-set-fractal


text:
include("solve_SDP.jl")

DATA_PATH =  joinpath("data_sdp")

instances = [
"case1354pegase","case1888rte", "case1951rte", "case_ACTIVSg2000", "case2383wp","case2736sp", "case2737sop", "case2746wop",
"case2746wp", "case2848rte", "case2868rte", "case2869pegase", "case3012wp", "case3120sp", "case3375wp","case6468rte",
"case6470rte", "case6495rte", "case6515rte",  "case9241pegase","case_ACTIVSg10k", "case13659pegase", "case_ACTIVSg25k"]

instances = ["case14"]

FORMULATIONS = ["NO_BLOCKS", "cholesky"]

isdir("Mosek_runs") || mkpath("Mosek_runs")
f = open("Mosek_runs\\summary_runs_all_Matpower.csv", "w")
write(f, "Instance ; ")
for FORMULATION in FORMULATIONS
    write(f, "$FORMULATION nb iter ; $FORMULATION time ;")
end
write(f, "\n")
for INSTANCE in instances
    println(INSTANCE)
    write(f, "$INSTANCE ; ")
    for FORMULATION in FORMULATIONS
        time, nb_iter = solve_sdp(INSTANCE, FORMULATION, DATA_PATH)
        write(f, "$nb_iter ; $time ; ")
    end
    write(f, "\n")
end
close(f)


text:
using Swizzles
using Swizzles.PermutedArrays

@testset "PermutedArrays" begin
    A = reshape(1:9, 3, 3)
    C = copy(A)
    B = Permute(map(reverse, axes(A))...).(C)
    @test B == reshape(9:-1:1, 3, 3)
    B = Permute(map(reverse, axes(A))...)(C)
    @test B == reshape(9:-1:1, 3, 3)
    B .= reshape(1:9, 3, 3)
    @test C == reshape(9:-1:1, 3, 3)

    A = reshape(1:9, 3, 3)
    B = Permute(reverse(1:3))(reshape(1:9, 3, 3))
    @test_throws PermutationMismatch A * B
    A = Permute(1:3, reverse(1:3))(reshape(1:9, 3, 3))
    A * B
    @test A * B == reshape(1:9, 3, 3) * reshape(1:9, 3, 3)
    @test Sum().(A) == sum(1:9)
end


text:
using PyPlot
using JLD

function test(path_acceleration::AbstractString,path_deceleration::AbstractString)

	files_acc = readdir(path_acceleration)
    files_dec = readdir(path_deceleration)

    a = zeros(length(files_acc) + length(files_dec))
    println("number of tests= ",length(a))

    ## Loop for acceleration

    v_acceleration1 = []  #zeros(length(files))
    v_acceleration2 = []

    commands_acc = []  #zeros(length(files))

    for index = 1:length(files_acc)

        file = files_acc[index]

        log_path = "$(homedir())/open_loop/acceleration/$file"

        d_rec=load(log_path)

        cmd_pwm_log = d_rec["cmd_pwm_log"]
        vel_est     = d_rec["vel_est"]

        t0 = max(cmd_pwm_log.t[1],vel_est.t[1])
        t_end = min(cmd_pwm_log.t[end],vel_est.t[end])

        t = t0+0.1:.02:t_end-0.1
        v = zeros(length(t))
        cmd = zeros(length(t))

        for i=1:length(t)
            v[i] = vel_est.z[t[i].>vel_est.t,1][end]
            #v[i] = pos_info.z[t[i].>pos_info.t,15][end]
            cmd[i] = cmd_pwm_log.z[t[i].>cmd_pwm_log.t,1][end]
        end

        v_opt   = v[1:end]
        cmd_opt = cmd[1:end]
        v_opt   = v[cmd.>90]
        cmd_opt = cmd[cmd.>90]

        figure(1)

        plot(v_opt)
        grid("on")

        # println("length= ",length(files))
        # println("v_opt= ",length(v_opt))
        # println("cmd_opt= ",length(cmd_opt))
        # println("v_opt[end]= ",v_opt[end])
        # println("a[i]= ",a[i])
        # println("v_opt[end]/4= ",v_opt[end]/4)

        v_lms1 = zeros(length(v_opt)-1)
        v_lms2 = zeros(length(v_opt)-1)

        for j = 1:length(v_opt)-1

            v_lms1[j] = v_opt[j+1] - v_opt[j]
            v_lms2[j] = v_opt[j]

        end

        v_acceleration1 = vcat(v_acceleration1,v_lms1)
        v_acceleration2 = vcat(v_acceleration2,v_lms2)
        a[index] = v_opt[end]/4
        commands_acc = vcat(commands_acc,cmd_opt[1:end-1])

        # println("a= ",a)

        figure(2)

        plot(cmd_opt[1],a[index], "*")
        grid("on")

    end

    ## Loop for deceleration

    v_deceleration1 = []  #zeros(length(files))
    v_deceleration2 = []

    commands_dec = []  #zeros(length(files))

    for index = 1:length(files_dec)

        file = files_dec[index]

        log_path = "$(homedir())/open_loop/deceleration/$file"

        d_rec=load(log_path)

        cmd_pwm_log = d_rec["cmd_pwm_log"]
        vel_est     = d_rec["vel_est"]

        t0 = max(cmd_pwm_log.t[1],vel_est.t[1])
        t_end = min(cmd_pwm_log.t[end],vel_est.t[end])

        t = t0+0.1:.02:t_end-0.1
        v = zeros(length(t))
        cmd = zeros(length(t))

        for i=1:length(t)
            v[i] = vel_est.z[t[i].>vel_est.t,1][end]
            #v[i] = pos_info.z[t[i].>pos_info.t,15][end]
            cmd[i] = cmd_pwm_log.z[t[i].>cmd_pwm_log.t,1][end]
        end

        v_opt   = v[1:end]
        cmd_opt = cmd[1:end]
        v_opt   = v[cmd.<100]
        cmd_opt = cmd[cmd.<100]

        v_equal = 0
        flag = 0

        for i = 2:length(v_opt)

            if v_opt[i]==v_opt[i-1]
                v_equal = v_equal + 1
            else
                v_equal = 0
            end

            if v_equal > 10
                flag = i
                break
            end
        end

        v_opt = v_opt[1:flag]
        cmd_opt = cmd_opt[1:flag]

        a[length(files_acc) + index] = (- v_opt[1])/4

        # plot(cmd_opt[5],a[length(files_acc) + index], "*")
        # grid("on")

        figure(3)

        plot(v_opt)
        grid("on")

        v_lms1 = zeros(length(v_opt)-1)
        v_lms2 = zeros(length(v_opt)-1)

        for j = 1:length(v_opt)-1

            v_lms1[j] = v_opt[j+1] - v_opt[j]
            v_lms2[j] = v_opt[j]

        end

        v_deceleration1 = vcat(v_deceleration1,v_lms1)
        v_deceleration2 = vcat(v_deceleration2,v_lms2)
        commands_dec = vcat(commands_dec,cmd_opt[1:end-1])

    end

    # Least Mean Square

    # println("length of v_acceleration1= ",length(v_acceleration1))
    # println("length of commands= ",length(commands))

    A = zeros(length(v_acceleration1)+length(v_deceleration1),5)
    B = zeros(length(v_acceleration1)+length(v_deceleration1),1)
    x = zeros(5,1)     #[c1_+;c2_+,cM]

    A[1:length(commands_acc),1] = commands_acc
    A[1:length(commands_acc),2] = 1
    A[1:length(v_acceleration2),3] = -v_acceleration2
    A[length(v_acceleration2)+1:end,3] = -v_deceleration2
    A[length(commands_acc)+1:end,4] = commands_dec
    A[length(commands_acc)+1:end,5] = 1

    B[1:length(v_acceleration1)] = v_acceleration1
    B[length(v_acceleration1)+1:end] = v_deceleration1

    A = 0.02*A # 0.02 is the discretization

    x = inv(A'*A)*A'*B

    println("the computed constants are: ", x)

    # x_axis = (findmin(commands)[1]):0.02:(findmax(commands)[1])
    # y_axis = x[1]*x_axis + x[2]
    # plot(x_axis,y_axis)

end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
using AzSessions, Dates, HTTP, JSON, JSONWebTokens, Test

credentials = JSON.parse(ENV["AZURE_CREDENTIALS"])
AzSessions.write_manifest(;client_id=credentials["clientId"], client_secret=credentials["clientSecret"], tenant=credentials["tenantId"])

function running_on_azure()
    try
        HTTP.request(
            "GET",
            "http://169.254.169.254/metadata/instance?api-version=2017-08-01",
            Dict("Metadata"=>"true"))
        return true
    catch
        return false
    end
end

if running_on_azure()
    # TODO - not sure why this doesn't work on CI
    @test_skip @testset "AzSessions, VM" begin
        session = AzSession(;protocal=AzVMCredentials)
        @test now(Dates.UTC) <= session.expiry
        t = token(session)
        @test isa(t,String)
        t2 = token(session)
        @test t2 == t

        session.token = "x"
        session.expiry = now(Dates.UTC) - Dates.Second(1)
        t2 = token(session)
        @test t2 != "x"
    end
end

@testset "AzSessions, Client Credentials" begin
    session = AzSession(;protocal=AzClientCredentials, client_id=credentials["clientId"], client_secret=credentials["clientSecret"])
    @test now(Dates.UTC) < session.expiry
    t = token(session)
    @test isa(t,String)
    t2 = token(session)
    @test t2 == t

    session.token = "x"
    session.expiry = now(Dates.UTC) - Dates.Second(1)
    t2 = token(session)
    @test t2 != "x"
end

# TODO: requires user interaction (can we use Mocking.jl)
@test_skip @testset "AzSessions, Device code flow credentials" begin
    session = AzSession(;protocal=AzDeviceCodeFlowCredentials)
    @test now(Dates.UTC) <= session.expiry
    t = token(session)
    @test isa(t,String)
    t2 = token(session)
    @test t2 == t

    session.token == "x"
    session.expiry = now(Dates.UTC) - Dates.Second(1)
    t2 = token(session)
    @test t2 != "x"

    session2 = AzSession(session;scope="https://storage.azure.com/user_impersonation")
    t = token(session2)
    decodedJWT = JSONWebTokens.decode(JSONWebTokens.None(), t)
    @test decodedJWT["aud"] == "https://storage.azure.com"
end

# TODO - the following testset will only work if the machine can start a web-browser
@test_skip @testset "AzSessions, Authorization code flow credentials" begin
    session = AzSession(;protocal=AzAuthCodeFlowCredentials)
    @test now(Dates.UTC) <= session.expiry
    t = token(session)
    @test isa(t,String)
    t2 = token(session)
    @test t2 == t

    session.token == "x"
    session.expiry = now(Dates.UTC) - Dates.Second(1)
    t2 = token(session)
    @test t2 != "x"

    session2 = AzSession(session;scope="https://storage.azure.com/user_impersonation")
    t = token(session2)
    decodedJWT = JSONWebTokens.decode(JSONWebTokens.None(), t)
    @test decodedJWT["aud"] == "https://storage.azure.com"
end

# TODO: requires user interaction, can we use Mocking.jl?
@test_skip @testset "AzSessions, Device code flow credentials is the default" begin
    session = AzSession()
    @test now(Dates.UTC) <= session.expiry
    t = token(session)
    @test isa(t,String)
    t2 = token(session)
    @test t2 == t

    session.token == "x"
    session.expiry = now(Dates.UTC) - Dates.Second(1)
    t2 = token(session)
    @test t2 != "x"
end

@testset "AzSessions, Client Credentials, serialize" begin
    session = AzSessions.AzClientCredentialsSession(
        "AzClientCredentials",
        "myclientid",
        "myclientsecret",
        now(),
        "myresource",
        "mytenant",
        "mytoken")

    jsonsession = json(session)
    _session = AzSession(jsonsession)

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.client_secret == _session.client_secret
    @test session.expiry == _session.expiry
    @test session.resource == _session.resource
    @test session.tenant == _session.tenant
    @test session.token == _session.token

    _session = AzSession(JSON.parse(jsonsession))
    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.client_secret == _session.client_secret
    @test session.expiry == _session.expiry
    @test session.resource == _session.resource
    @test session.tenant == _session.tenant
    @test session.token == _session.token
end

@testset "AzSessions, VM Credentials, serialize" begin
    session = AzSessions.AzVMSession(
        "AzVMCredentials",
        now(),
        "myresource",
        "mytoken")

    jsonsession = json(session)
    _session = AzSession(jsonsession)

    @test session.protocal == _session.protocal
    @test session.expiry == _session.expiry
    @test session.resource == _session.resource
    @test session.token == _session.token
end

@testset "AzSessions, Auth Code Credentials, serialize" begin
    session = AzSessions.AzAuthCodeFlowSession(
        "AzAuthCodeFlowCredentials",
        "clientid",
        now(),
        "myidtoken",
        false,
        "redirecturi",
        "refreshtoken",
        "scopeauth",
        "scopetoken",
        "tenant",
        "token")

    jsonsession = json(session)
    _session = AzSession(jsonsession)

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.expiry == _session.expiry
    @test session.id_token == _session.id_token
    @test session.lock == _session.lock
    @test session.redirect_uri == _session.redirect_uri
    @test session.refresh_token == _session.refresh_token
    @test session.scope_auth == _session.scope_auth
    @test session.scope == _session.scope
    @test session.tenant == _session.tenant
    @test session.token == _session.token

    _session = AzSession(JSON.parse(jsonsession))

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.expiry == _session.expiry
    @test session.id_token == _session.id_token
    @test session.lock == _session.lock
    @test session.redirect_uri == _session.redirect_uri
    @test session.refresh_token == _session.refresh_token
    @test session.scope_auth == _session.scope_auth
    @test session.scope == _session.scope
    @test session.tenant == _session.tenant
    @test session.token == _session.token
end

@testset "AzSessions, Device Code Credentials, serialize" begin
    session = AzSessions.AzDeviceCodeFlowSession(
        "AzDeviceCodeFlowCredentials",
        "myclientid",
        now(),
        "myidtoken",
        true,
        "myrefreshtoken",
        "myscope",
        "myscope_auth",
        "mytenant",
        "mytoken")

    jsonsession = json(session)
    _session = AzSession(jsonsession)

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.expiry == _session.expiry
    @test session.id_token == _session.id_token
    @test session.lock == _session.lock
    @test session.refresh_token == _session.refresh_token
    @test session.scope == _session.scope
    @test session.scope_auth == _session.scope_auth
    @test session.tenant == _session.tenant
    @test session.token == _session.token

    _session = AzSession(JSON.parse(jsonsession))

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.expiry == _session.expiry
    @test session.id_token == _session.id_token
    @test session.lock == _session.lock
    @test session.refresh_token == _session.refresh_token
    @test session.scope == _session.scope
    @test session.scope_auth == _session.scope_auth
    @test session.tenant == _session.tenant
    @test session.token == _session.token
end

@testset "AzSesions, Client Credentials, copy" begin
    session = AzSessions.AzClientCredentialsSession(
        "AzClientCredentials",
        "myclientid",
        "myclientsecret",
        now(),
        "myresource",
        "mytenant",
        "mytoken")

    _session = copy(session)
    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.client_secret == _session.client_secret
    @test session.expiry == _session.expiry
    @test session.resource == _session.resource
    @test session.tenant == _session.tenant
    @test session.token == _session.token
end

@testset "AzSessions, VM Credentials, copy" begin
    session = AzSessions.AzVMSession(
        "AzVMCredentials",
        now(),
        "myresource",
        "mytoken")

    _session = copy(session)

    @test session.protocal == _session.protocal
    @test session.expiry == _session.expiry
    @test session.resource == _session.resource
    @test session.token == _session.token
end

@testset "AzSessions, Auth code flow credentials, copy" begin
    session = AzSessions.AzAuthCodeFlowSession(
        "AzAuthCodeFlowCredentials",
        "clientid",
        now(),
        "myidtoken",
        false,
        "redirecturi",
        "refreshtoken",
        "scopeauth",
        "scopetoken",
        "tenant",
        "token")

    _session = copy(session)

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.expiry == _session.expiry
    @test session.id_token == _session.id_token
    @test session.lock == _session.lock
    @test session.redirect_uri == _session.redirect_uri
    @test session.refresh_token == _session.refresh_token
    @test session.scope_auth == _session.scope_auth
    @test session.scope == _session.scope
    @test session.tenant == _session.tenant
    @test session.token == _session.token
end

@testset "AzSessions, Device Code Credentials, copy" begin
    session = AzSessions.AzDeviceCodeFlowSession(
        "AzDeviceCodeFlowCredentials",
        "myclientid",
        now(),
        "myidtoken",
        true,
        "myrefreshtoken",
        "myscope",
        "myscope_auth",
        "mytenant",
        "mytoken")

    _session = copy(session)

    @test session.protocal == _session.protocal
    @test session.client_id == _session.client_id
    @test session.expiry == _session.expiry
    @test session.id_token == _session.id_token
    @test session.lock == _session.lock
    @test session.refresh_token == _session.refresh_token
    @test session.scope == _session.scope
    @test session.scope_auth == _session.scope_auth
    @test session.tenant == _session.tenant
    @test session.token == _session.token
end

@testset "AzSessions, merge scopes" begin
    session = AzSession(;
        scope="openid+offline_access+https://management.azure.com/user_impersonation",
        scope_auth = "openid+offline_access+https://storage.azure.com/user_impersonation",
        lazy = true)

    scopes = split(session.scope_auth, '+')
    @test length(scopes) == 4
    @test length(unique(scopes)) == 4
    for scope in scopes
        @test scope  ∈ ["openid","offline_access","https://management.azure.com/user_impersonation","https://storage.azure.com/user_impersonation"]
    end
end

@testset "AzSessions, VM credentials, scrub" begin
    session = AzSessions.AzClientCredentialsSession(
            "AzClientCredentials",
            "myclientid",
            "myclientsecret",
            now(),
            "myresource",
            "mytenant",
            "mytoken")

    scrub!(session)

    @test session.client_secret == ""
    @test session.token == ""
end

@testset "AzSessions, Auth code flow, scrub" begin
    session = AzSessions.AzAuthCodeFlowSession(
        "AzAuthCodeFlowCredentials",
        "clientid",
        now(),
        "myidtoken",
        false,
        "redirecturi",
        "refreshtoken",
        "scopeauth",
        "scopetoken",
        "tenant",
        "token")

    scrub!(session)

    @test session.id_token == ""
    @test session.refresh_token == ""
    @test session.token == ""
end

@testset "AzSessions, Device code flow, scrub" begin
    session = AzSessions.AzDeviceCodeFlowSession(
        "AzDeviceCodeFlowCredentials",
        "myclientid",
        now(),
        "myidtoken",
        true,
        "myrefreshtoken",
        "myscope",
        "myscope_auth",
        "mytenant",
        "mytoken")

    scrub!(session)

    @test session.id_token == ""
    @test session.refresh_token == ""
    @test session.token == ""
end

@testset "AzSessions, write_manifest" begin
    AzSessions.write_manifest(;client_id="myclientid", client_secret="myclientsecret", tenant="mytenant")

    manifest = JSON.parse(read(AzSessions.manifestfile(), String))
    @test manifest["client_id"] == "myclientid"
    @test manifest["client_secret"] == "myclientsecret"
    @test manifest["tenant"] == "mytenant"
end

AzSessions.write_manifest(;client_id=credentials["clientId"], client_secret=credentials["clientSecret"], tenant=credentials["tenantId"], protocal="AzClientCredentials")
@testset "AzSessions, Client Credentials is the default" begin
    session = AzSession()
    @test session.protocal == "AzClientCredentials"
    @test now(Dates.UTC) < session.expiry
    t = token(session)
    @test isa(t,String)
    t2 = token(session)
    @test t2 == t

    session.token = "x"
    session.expiry = now(Dates.UTC) - Dates.Second(1)
    t2 = token(session)
    @test t2 != "x"
end
AzSessions.write_manifest(;client_id=credentials["clientId"], client_secret=credentials["clientSecret"], tenant=credentials["tenantId"], protocal="")


text:
#=
For details, please visit:
https://computationalmindset.com/en/neural-networks/experiments-with-neural-odes-in-julia.html#exp2
=#

using Flux, DiffEqFlux, DifferentialEquations, Plots

math_law(u) = sin.(2. * u) + cos.(2. * u)

function true_ode(du,u,p,t)
    true_A = [-0.15 2.10; -2.10 -0.10]
    du .= ((math_law(u))'true_A)'
end

tbegin = 0.0
tend = 4.0
datasize = 51
t = range(tbegin,tend,length=datasize)
u0 = [2.5; 0.5]
tspan = (tbegin,tend)
trange = range(tbegin,tend,length=datasize)
prob = ODEProblem(true_ode, u0, tspan)
dataset_ts = Array(solve(prob, Tsit5(), saveat=trange))

dudt = Chain(u -> math_law(u),
             Dense(2, 50, tanh),
             Dense(50, 2))

reltol = 1e-7
abstol = 1e-9
n_ode = NeuralODE(dudt, tspan, Tsit5(), saveat=trange, reltol=reltol,abstol=abstol)
ps = Flux.params(n_ode.p)

function loss_n_ode()
  pred = n_ode(u0)
  loss = sum(abs2, dataset_ts[1,:] .- pred[1,:]) +
         sum(abs2, dataset_ts[2,:] .- pred[2,:])
  loss
end

n_epochs = 400
learning_rate = 0.01
data = Iterators.repeated((), n_epochs)
opt = ADAM(learning_rate)

cb = function ()
  loss = loss_n_ode()
  println("Loss: ", loss)
end

println();
cb()

Flux.train!(loss_n_ode, ps, data, opt, cb=cb)

pl = plot(
  trange,
  dataset_ts[1,:],
  linewidth=2, ls=:dash,
  title="Neural ODE for forecasting",
  xaxis="t",
  label="original timeseries x(t)",
  legend=:right)
display(pl)

pl = plot!(
  trange,
  dataset_ts[2,:],
  linewidth=2, ls=:dash,
  label="original timeseries y(t)")
display(pl)

pred = n_ode(u0)

pl = plot!(
  trange,
  pred[1,:],
  linewidth=1,
  label="predicted timeseries x(t)")
display(pl)

pl = plot!(
    trange,
    pred[2,:],
    linewidth=1,
    label="predicted timeseries y(t)")
display(pl)

tbegin_forecast = tend
tend_forecast = tbegin_forecast + 5.0
tspan_forecast = (tbegin_forecast, tend_forecast)
datasize_forecast = 351
trange_forecast = range(tspan_forecast[1], tspan_forecast[2], length=datasize_forecast)
u0_forecast = [dataset_ts[1,datasize], dataset_ts[2,datasize]]

n_ode_forecast = NeuralODE(
  n_ode.model, tspan_forecast;
  p=n_ode.p, saveat=trange_forecast, reltol=reltol, abstol=abstol)
forecast = n_ode_forecast(u0_forecast)

pl = plot!(
  trange_forecast,
  forecast[1,:],
  linewidth=1,
  label="forecast x(t)")
display(pl)

pl = plot!(
    trange_forecast,
    forecast[2,:],
    linewidth=1,
    label="forecast y(t)")
display(pl)

max_stars_repo_name: ettoremessina/differential-equations

max_issues_repo_name: ettoremessina/differential-equations

max_forks_repo_name: ettoremessina/differential-equations


text:
using StructJuMP

m = StructuredModel()

@defStochasticVar(m, xh0[ASSETS] >= 0)
@defStochasticVar(m, xb0[ASSETS] >= 0)
@addConstraint(m, (1+Ct)*sum{val[j]*xb0[j], j=ASSETS} <= BUDGET)

for n1 in NODES
    if PARENT[n1] == 0
        sub = StructuredModel(parent=m)
        @defStochasticVar(sub, xh1[ASSETS] >= 0)
        @defStochasticVar(sub, xb1[ASSETS] >= 0)
        @defStochasticVar(sub, xs1[ASSETS] >= 0)
        for j in ASSETS
            @addConstraint(sub, xh1[j] == (1+Ret[n1,j])*xh0[j]+xb1[j]-xs1[j])
        end
        @addConstraint(sub, (1-Ct)*sum{val[j]*xs1[j], j=ASSETS} == Liability1 + (1+Ct)*sum{val[j]*xb1[j]})
    end

    for n2 in NODES
        if PARENT[n2] == n1
            sub = StructuredModel(parent=m)
            @defStochasticVar(sub, xh2[ASSETS] >= 0)
            @defStochasticVar(sub, xb2[ASSETS] >= 0)
            @defStochasticVar(sub, xs2[ASSETS] >= 0)
            for j in ASSETS
                @addConstraint(sub, xh2[j] == (1+Ret2[n2,j])*xh1[j] + xb2[j] - xs2[j])
            end
            @addConstraint(sub, (1-Ct)*sum{val[j]*xs2[j], j=ASSETS} == Liability2 + (1+Ct)*sum{val[j]*xb2[j]})
            @defStochasticVar(sub, wealth)
            @addConstraint(sub, wealth == Prob[n1]*Prob[n2]*sum{val[j]*xh2[j], j=ASSETS})
            @setObjective(sub, mu - Rho*((wealth^2-mu^2))*Prob[n1]*Prob[n2])
        end
    end
end

@addConstraint(m, sum{Prob[n1]*Prob[n2]*wealth[(n1,n2)], n1=NODES, n2=NODES; PARENT[n1]==0 && PARENT[n2]==n1})

# things of note:
# Block(parent::Union(Model,Block), index) constructor
# should do name mangling so that variables get associated with their block
# add function block(parent, child1, child2, ...) that returns reference to subproblem. This should then be passable as index to subproblem variables


text:
# This file constains code copied from https://github.com/jeanpauphilet/SubsetSelection.jl
# as at commit e6b3258ba394f6152853498d81d3cecb9d67e097

const sgtol = 1e-6

mutable struct PolyakCache
  best_inds::Vector{Int}
  best_upper::Float64
  best_lower::Float64
end
PolyakCache(inds::Vector{Int}) = PolyakCache(inds, Inf, -Inf)

abstract type SteppingRule end
struct ConstantStepping <: SteppingRule
  stepsize::Float64
end
mutable struct PolyakStepping <: SteppingRule
  initial_factor::Float64        # initial factor for scaling (ub-lb)/norm^2
  wait::Int64                    # currently unused TODO
  function PolyakStepping(initial_factor, wait)
      if 0.0 <= initial_factor <= 2.0
          new(initial_factor, wait)
      else
          error("Initial stepping parameter must be between 0 and 2. Current
                    choice is $initial_factor.")
      end
  end
end
PolyakStepping() = PolyakStepping(2.0, 30)

function resetindices(::PolyakStepping)
  true
end
function resetindices(::SteppingRule)
  false
end

function getdelta!(sr::SteppingRule, ::Any, ::Array{Float64,2}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Int}, ::Int, ::Float64, ::SubsetSelection.Cache, ::PolyakCache)
  error("Need to define `getdelta` for stepping rule $(sr).")
end
function getdelta!(sr::ConstantStepping, ::Any, ::Array{Float64,2}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Int}, ::Int, ::Float64, ::SubsetSelection.Cache, ::PolyakCache)
  sr.stepsize
end
function getdelta!(sr::PolyakStepping,
      Y::YVector, X::Array{Float64,2},
      α::Vector{Float64}, ∇::Vector{Float64},
      indices::Vector{Int}, n_indices::Int, γ::Float64,
      cache::SubsetSelection.Cache, pc::PolyakCache
     )
  @timeit TimeStats "Dual Bound" lower_bound = dual_bound(SubsetSelection.OLS(), Y, X, α, indices, n_indices, γ, cache)
  @timeit TimeStats "Primal Bound" upper_bound = primal_bound(SubsetSelection.OLS(), Y, X, γ, indices, n_indices)
  if upper_bound < pc.best_upper
    pc.best_upper = upper_bound
    pc.best_inds .= indices
  end
  (lower_bound > pc.best_lower) && (pc.best_lower = lower_bound)
  if lower_bound - sgtol > upper_bound + sgtol
      error("Bounds overlap. LB = $lowerbound and UB = $upper_bound.")
  end
  sr.initial_factor * (pc.best_upper - lower_bound) / sum(abs2.(∇))
end

##############################################
##DUAL SUB-GRADIENT ALGORITHM
##############################################
""" Function to compute a sparse regressor/classifier. Solve an optimization problem of the form
        min_s max_α f(α, s)
by gradient ascent on α:        α_{t+1} ← α_t + δ ∇_α f(α_t, s_t)
and partial minimization on s:  s_{t+1} ← argmin_s f(α_t, s)
INPUTS
- ℓ           Loss function used
- Card        Model to enforce sparsity (constraint or penalty)
- Y (n×1)     Vector of outputs. For classification, use ±1 labels
- X (n×p)     Array of inputs.
- indInit     (optional) Initial subset of features s
- αInit       (optional) Initial dual variable α
- γ           (optional) ℓ2 regularization penalty
- intercept   (optional) Boolean. If true, an intercept term is computed as well
- maxIter     (optional) Total number of Iterations
- δ           (optional) Gradient stepsize
- gradUp      (optional) Number of gradient updates of dual variable α performed per update of primal variable s
- anticycling (optional) Boolean. If true, the algorithm stops as soon as the support is not unchanged from one iteration to another
- averaging   (optional) Boolean. If true, the dual solution is averaged over past iterates
OUTPUT
- SparseEstimator """
function subsetSelection_bm(ℓ::LossFunction, Card::Sparsity,
    Y::YVector,
    X::Array{Float64,2};
    indInit = SubsetSelection.ind_init(Card, size(X,2)),
    αInit = SubsetSelection.alpha_init(ℓ, Y),
    γ = 1/sqrt(size(X,1)),  intercept = false,
    maxIter = 100, sr::SteppingRule = ConstantStepping(1e-3), gradUp = 1,
    anticycling = false, averaging = true)

  @timeit TimeStats "in subgradient descent" begin

  n,p = size(X)
  cache = SubsetSelection.SubsetSelection.Cache(n, p)

  #Add sanity checks
  if size(Y,1) != n
    throw(DimensionMismatch("X and Y must have the same number of rows"))
  end
  if isa(ℓ, SubsetSelection.Classification)
    levels = sort(unique(Y))
    if length(levels) != 2
      throw(ArgumentError("subsetSelection only supports two-class classification"))
    elseif (levels[1] != -1) || (levels[2] != 1)
      throw(ArgumentError("Class labels must be ±1's"))
    end
  end

  indices = indInit #Support
  n_indices = length(indices)
  best_inds = copy(indices)

  n_indices_max = SubsetSelection.max_index_size(Card, p)
  resize!(indices, n_indices_max)

  indices_old = Vector{Int}(n_indices_max)
  α = αInit[:]  #Dual variable α
  a = αInit[:]  #Past average of α

  # We will compute the dual bound incorrectly unless indicies match initial α
  if resetindices(sr)
    n_indices = SubsetSelection.partial_min!(indices, Card, X, α, γ, cache)
  end

  pc = PolyakCache(copy(indices))

  bounds_match = false

  ##Dual Sub-gradient Algorithm
  for iter in 2:maxIter

    #Gradient ascent on α
    for _ in 1:min(gradUp, div(p, n_indices))
      ∇ = SubsetSelection.grad_dual(ℓ, Y, X, α, indices, n_indices, γ, cache)
      δ = getdelta!(sr, Y, X, α, ∇, indices, n_indices, γ, cache, pc)
      α .+= δ*∇
      α = SubsetSelection.proj_dual(ℓ, Y, α)
      α = SubsetSelection.proj_intercept(intercept, α)
    end

    if !all(isfinite.(α))
        warn("Algorithm diverges! Did you normalize your data? Otherwise, try reducing stepsize δ.")
    end
    #Update average a
    @__dot__ a = (iter - 1) / iter * a + 1 / iter * α
    # a *= (iter - 1)/iter; a .+= α/iter

    #Minimization w.r.t. s
    indices_old[1:n_indices] = indices[1:n_indices]
    n_indices = SubsetSelection.partial_min!(indices, Card, X, α, γ, cache)

    #Anticycling rule: Stop if indices_old == indices
    if anticycling && SubsetSelection.indices_same(indices, indices_old, n_indices)
      averaging = false #If the algorithm stops because of cycling, averaging is not needed
      break
    end
    # If lower and upper bounds have closed, skip to the end
    if abs(pc.best_upper - pc.best_lower) / (1 + abs(pc.best_upper)) < 1e-4
      indices .= pc.best_inds
      break
    end
  end

  ##Compute sparse estimator
  #Subset of relevant features
  n_indices = SubsetSelection.partial_min!(indices, Card, X, averaging ? a : α, γ, cache)
  #Regressor
  # w = [-γ * dot(X[:, indices[j]], a) for j in 1:n_indices]
  w = SubsetSelection.recover_primal(ℓ, Y, X[:,indices], γ)
  #Bias
  b = SubsetSelection.compute_bias(ℓ, Y, X, a, indices, n_indices, γ, intercept, cache)

  #Resize final indices vector to only have relevant entries
  resize!(indices, n_indices)

   end # timing

  return SubsetSelection.SparseEstimator(ℓ, Card, indices, w, a, b, maxIter)
end

function ax_squared(X, α::Vector{Float64}, indices::Vector{Int}, n_indices::Int)
  # TODO update ax in cache and use it rather than recomputing
  axsum = 0.0
  for j = 1:n_indices
    x = @view(X[:, indices[j]])
    axsum += dot(α, x)^2
  end
  axsum
end

function primal_bound(ℓ::SubsetSelection.OLS, Y::YVector, X::Array{Float64,2}, γ::Float64, indices::Vector{Int}, n_indices::Int)
  @timeit TimeStats "Inverting matrix" αstar = SubsetSelectionCIO.sparse_inverse(ℓ, Y, X[:, indices], γ) # TODO could do this less often. also don't return vector after function.
  axsum = ax_squared(X, αstar, indices, n_indices)
  bound = -0.5 * dot(αstar, αstar) - dot(Y, αstar) - γ * 0.5 * axsum
  # Normalize TODO normalize this and grad dual
  # bound / size(X, 1)
end

function dual_bound(ℓ::SubsetSelection.OLS, Y::YVector, X::Array{Float64,2}, α::Vector{Float64}, indices::Vector{Int}, n_indices::Int, γ, cache::SubsetSelection.Cache)
  axsum = ax_squared(X, α, indices, n_indices)
  bound = -0.5 * dot(α, α) - dot(Y, α) - γ * 0.5 * axsum
  # Normalize
  # bound / size(X, 1)
end