text:
#############
# Repo Type #
#############

@ghdef mutable struct Repo
    name::Union{String, Nothing}
    full_name::Union{String, Nothing}
    description::Union{String, Nothing}
    language::Union{String, Nothing}
    default_branch::Union{String, Nothing}
    owner::Union{Owner, Nothing}
    parent::Union{Repo, Nothing}
    source::Union{Repo, Nothing}
    id::Union{Int, Nothing}
    size::Union{Int, Nothing}
    subscribers_count::Union{Int, Nothing}
    forks_count::Union{Int, Nothing}
    stargazers_count::Union{Int, Nothing}
    watchers_count::Union{Int, Nothing}
    open_issues_count::Union{Int, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    clone_url::Union{HTTP.URI, Nothing}
    ssh_url::Union{HTTP.URI, Nothing}
    homepage::Union{HTTP.URI, Nothing}
    pushed_at::Union{Dates.DateTime, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    has_issues::Union{Bool, Nothing}
    has_wiki::Union{Bool, Nothing}
    has_downloads::Union{Bool, Nothing}
    has_pages::Union{Bool, Nothing}
    private::Union{Bool, Nothing}
    fork::Union{Bool, Nothing}
    permissions::Union{Dict, Nothing}
end

Repo(full_name::AbstractString) = Repo(Dict("full_name" => full_name))

namefield(repo::Repo) = repo.full_name

###############
# API Methods #
###############

# repos #
#-------#

@api_default function repo(api::GitHubAPI, repo_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo_obj))"; options...)
    return Repo(result)
end

@api_default function create_repo(api::GitHubAPI, owner, repo_name::String, params=Dict{String,Any}(); options...)
    params["name"] = repo_name
    if isorg(owner)
        result = gh_post_json(api, "/orgs/$(name(owner))/repos"; params=params, options...)
    else
        result = gh_post_json(api, "/user/repos"; params=params, options...)
    end
    return Repo(result)
end

@api_default function delete_repo(api::GitHubAPI, repo; options...)
    return gh_delete(api, "/repos/$(name(repo))"; options...)
end

# forks #
#-------#

@api_default function forks(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/forks"; options...)
    return map(Repo, results), page_data
end

@api_default function create_fork(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/forks"; options...)
    return Repo(result)
end

# contributors/collaborators #
#----------------------------#

@api_default function contributors(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/contributors"; options...)
    results = [Dict("contributor" => Owner(i), "contributions" => i["contributions"]) for i in results]
    return results, page_data
end

@api_default function collaborators(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_json(api, "/repos/$(name(repo))/collaborators"; options...)
    return map(Owner, results), page_data
end

@api_default function iscollaborator(api::GitHubAPI, repo, user; options...)
    path = "/repos/$(name(repo))/collaborators/$(name(user))"
    r = gh_get(api, path; handle_error = false, options...)
    r.status == 204 && return true
    r.status == 404 && return false
    handle_response_error(r)  # 404 is not an error in this case
    return false
end

@api_default function add_collaborator(api::GitHubAPI, repo, user; options...)
    path = "/repos/$(name(repo))/collaborators/$(name(user))"
    return gh_put(api, path; options...)
end

@api_default function remove_collaborator(api::GitHubAPI, repo, user; options...)
    path = "/repos/$(name(repo))/collaborators/$(name(user))"
    return gh_delete(api, path; options...)
end

@api_default function collaborator_permission(api::GitHubAPI, repo, user; options...)
    results, page_data = gh_get_json(api, "/repos/$(name(repo))/collaborators/$(name(user))/permission"; options...)
    return results[2]
end

# stats #
#-------#

@api_default function stats(api::GitHubAPI, repo, stat, attempts = 3; options...)
    path = "/repos/$(name(repo))/stats/$(name(stat))"
    local r
    for a in 1:attempts
        r = gh_get(api, path; handle_error = false, options...)
        r.status == 200 && return r
        sleep(2.0)
    end
    return r
end


text:
module Fortune

using ..Geometry
using ..EventQueue
using ..BeachLine
using ..Diagram

@enum SIDE LEFT RIGHT

function randf(start, finish, n)
   v = rand(n)
   return map(x -> start + x*(finish-start), v)
end

function init(points::Array{Tuple{Real, Real}, 1}, W, H)
	n = size(points)[1]

    V = Diagram.DCEL(points)
	T = BeachLine.BST()
	Q = EventQueue.Heap(3n) # 3n is enough, see page 166 on BCKO (de Berg, Cheong, Kreveld, Overmars)

    Geometry.setSafe(H)

	for region in V.regions
		EventQueue.push!(Q, EventQueue.SiteEvent(region))
	end

	return V, T, Q
end

function handleEvent(V::Diagram.DCEL, T::BeachLine.BST, Q::EventQueue.Heap, event::EventQueue.SiteEvent)
   #println("SITE EVENT: ", event.region.generator)
   #readline(stdin)
   ly = event.region.generator[2] # sweep line

   arc, arcAbove, sideOnSpecialCase = BeachLine.insert(T, event.region, ly)

   if T.root == arc
       return
   end

   if arcAbove == nothing # first couple of points have the same y
       bp = Geometry.parabolaIntersection(arc.parent.leftFocus, arc.parent.rightFocus, ly)

       if sideOnSpecialCase == BeachLine.LEFT
           left = arc
           right = arc.parent.rightChild
       else
           left = arc.parent.leftChild
           right = arc
       end

       #dir = Geometry.rotateVectorCCW(Geometry.subVector(left.region.generator, right.region.generator))
       dir = Geometry.rotateVectorCCW(Geometry.subVector(left.region.generator, right.region.generator))

       he1 = Diagram.HalfEdge(Geometry.subVector(bp, dir), false, nothing, nothing, nothing, left.region.generator)
       he2 = Diagram.HalfEdge(Geometry.addVector(bp, dir), false, nothing, nothing, nothing, right.region.generator)
       Diagram.makeTwins(he1, he2)
       Base.push!(V.halfEdges, he1)
       Base.push!(V.halfEdges, he2)
       left.parent.halfEdge = he1
       left.region.borderHead = he1 # following the convention that inside half-edges go counter-clockwise
       right.region.borderHead = he2
   else
       if arcAbove.disappearsAt != nothing # circle event for arc vertically above 'arc' was a false alarm
          EventQueue.remove(Q, arcAbove.disappearsAt)
          arcAbove.disappearsAt = nothing
       end

       rightBreakpoint = arc.parent
       leftBreakpoint = rightBreakpoint.parent

       # create half edge records for each new breakpoint
       f = Geometry.parabola(arcAbove.region.generator, ly)
       x = event.region.generator[1]

       if f != nothing
           breakpoint = (x, f(x))
       else
           breakpoint = nothing
       end

       # direction of the halfedge, rotating clockwise around the arc above (always going right)
       dir = Geometry.rotateVectorCCW(Geometry.subVector(arc.region.generator, arcAbove.region.generator))

       he1 = Diagram.HalfEdge(Geometry.subVector(breakpoint, dir), false, nothing, nothing, nothing, arcAbove.region.generator)
       he2 = Diagram.HalfEdge(Geometry.addVector(breakpoint, dir), false, nothing, nothing, nothing, event.region.generator)
       Diagram.makeTwins(he1, he2)
       Base.push!(V.halfEdges, he1)
       Base.push!(V.halfEdges, he2)
       leftBreakpoint.halfEdge = he1
       rightBreakpoint.halfEdge = he2
       event.region.borderHead = he2 # following the convention that inside half-edges go counter-clockwise
       arcAbove.region.borderHead = he1
   end

   # check for new circle events where 'arc' is the rightmost arc
   c = arc
   b = arc.prev
   if b != nothing
      a = b.prev

      # there's a triple with 'arc' on the rightmost position
      if a != nothing && Geometry.arcWillConverge(a.region.generator, b.region.generator, c.region.generator)
         O, r = Geometry.circumcircle(a.region.generator, b.region.generator, c.region.generator)

         ev = EventQueue.CircleEvent(O, r, b)
         b.disappearsAt = ev
         EventQueue.push!(Q, ev)
      end
   end

   # check for new circle events where 'arc' is the leftmost arc
   a = arc
   b = arc.next
   if b != nothing
      c = b.next

      # there's a triple with 'arc' on the leftmost position
      if c != nothing && Geometry.arcWillConverge(a.region.generator, b.region.generator, c.region.generator)
         O, r = Geometry.circumcircle(a.region.generator, b.region.generator, c.region.generator)

         ev = EventQueue.CircleEvent(O, r, b)
         b.disappearsAt = ev
         EventQueue.push!(Q, ev)
      end
   end
end

function handleEvent(V::Diagram.DCEL, T::BeachLine.BST, Q::EventQueue.Heap, event::EventQueue.CircleEvent)
   #println("CIRCLE EVENT: ", event.coordinates, "; ", event.center)
   #readline(stdin)
   arc = event.disappearingArc

   # find the Voronoi edges that are going to join together
   leftBreakpoint = rightBreakpoint = arc.parent
   while leftBreakpoint.leftFocus != arc.prev.region.generator ||
         leftBreakpoint.rightFocus != arc.region.generator
      leftBreakpoint = leftBreakpoint.parent
   end
   while rightBreakpoint.leftFocus != arc.region.generator ||
         rightBreakpoint.rightFocus != arc.next.region.generator
      rightBreakpoint = rightBreakpoint.parent
   end

   # fix the extremes of the joining Voronoi edges
   leftBreakpoint.halfEdge.origin = rightBreakpoint.halfEdge.origin = event.center
   leftBreakpoint.halfEdge.isFixed = rightBreakpoint.halfEdge.isFixed = true

   # remembering breakpoint edges, because BeachLine.remove will likely change the breakpoints
   lBhe = leftBreakpoint.halfEdge
   rBhe = rightBreakpoint.halfEdge

   arc.prev.next = arc.next
   arc.next.prev = arc.prev
   newBreakpoint = BeachLine.remove(T, event.disappearingArc)

   # direction of the new halfedge, rotating clockwise around the arc to the right of the disappearing arc (always going down)
   dir = Geometry.rotateVectorCCW(Geometry.subVector(arc.prev.region.generator, arc.next.region.generator))
   #dir = Geometry.multVector(1000, dir)

   # create new half edge for the newly formed breakpoint
   he1 = Diagram.HalfEdge(Geometry.addVector(event.center, dir), false, nothing, nothing, nothing, arc.prev.region.generator)
   he2 = Diagram.HalfEdge(event.center, false, nothing, nothing, nothing, arc.next.region.generator)
   Diagram.makeTwins(he1, he2)
   Base.push!(V.halfEdges, he1)
   Base.push!(V.halfEdges, he2)
   newBreakpoint.halfEdge = he1 # he2 is left bound to the vertex
   he2.isFixed = true

   # join the adjacent edges in the region lists
   Diagram.concat(newBreakpoint.halfEdge, lBhe)
   Diagram.concat(rBhe.twin, newBreakpoint.halfEdge.twin)
   Diagram.concat(lBhe.twin, rBhe)

   # remove circle events of triples that don't exist anymore
   if arc.prev.disappearsAt != nothing
      EventQueue.remove(Q, arc.prev.disappearsAt)
      arc.prev.disappearsAt = nothing
   end
   if arc.next.disappearsAt != nothing
      EventQueue.remove(Q, arc.next.disappearsAt)
      arc.next.disappearsAt = nothing
   end

   # check the left new triple of adjacent arcs for circle events
   if arc.prev != nothing && arc.next != nothing && arc.prev.prev != nothing
      a = arc.prev.prev
      b = arc.prev
      c = arc.next

      if Geometry.arcWillConverge(a.region.generator, b.region.generator, c.region.generator)
         O, r = Geometry.circumcircle(a.region.generator, b.region.generator, c.region.generator)

         ev = EventQueue.CircleEvent(O, r, b)
         b.disappearsAt = ev
         EventQueue.push!(Q, ev)
      end
   end

   # check the right new triple of adjacent arcs for circle events
   if arc.prev != nothing && arc.next != nothing && arc.next.next != nothing
      a = arc.prev
      b = arc.next
      c = arc.next.next

      if Geometry.arcWillConverge(a.region.generator, b.region.generator, c.region.generator)
         O, r = Geometry.circumcircle(a.region.generator, b.region.generator, c.region.generator)

         ev = EventQueue.CircleEvent(O, r, b)
         b.disappearsAt = ev
         EventQueue.push!(Q, ev)
      end
   end
end

function compute(points::Array{Tuple{Real, Real}, 1}, W, H)
   V, T, Q = init(points, W, H)

   while (event = EventQueue.pop(Q)) != nothing
      Fortune.handleEvent(V, T, Q, event)
   end

   return V
end

end # module

ext: jl

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
@testset "Differential reduction" begin
    ode = @ODEmodel(
        x'(t) = a * x(t),
        y(t) = x(t)
    )
    pbr = PBRepresentation(ode, find_ioequations(ode))
    R, (y_5, y_4) = Nemo.PolynomialRing(Nemo.QQ, ["y_5", "y_4"])
    res = diffreduce(y_5, pbr)
    @test res == str_to_var("a", parent(res))^5 * str_to_var("y_0", parent(res))
    res = diffreduce(y_4, pbr)
    @test res == str_to_var("a", parent(res))^4 * str_to_var("y_0", parent(res))
    res = diffreduce(y_4^2, pbr)
    @test res == str_to_var("a", parent(res))^8 * str_to_var("y_0", parent(res))^2

    ode = @ODEmodel(
        x1'(t) = x2(t),
        x2'(t) = -x1(t),
        y(t) = x1(t) + u(t)
    )
    pbr = PBRepresentation(ode, find_ioequations(ode))
    R, (y_5, y_4, u_10) = Nemo.PolynomialRing(Nemo.QQ, ["y_5", "y_4", "u_10"])
    res = diffreduce(y_4 + u_10, pbr)
    @test res == str_to_var("u_10", parent(res)) + str_to_var("y_0", parent(res)) - str_to_var("u_0", parent(res)) + str_to_var("u_4", parent(res))

    # Next two verified with Maple
    # Mizuka's example
    R, (y_0, y_1, y_2, y_3, u_0, u_1, u_2) = Nemo.PolynomialRing(Nemo.QQ, ["y_0", "y_1", "y_2", "y_3", "u_0", "u_1", "u_2"])
    pbr = PBRepresentation(["y"], ["u"], Array{String, 1}(), Dict("y" => 2), Dict("y" => 27*y_0^9 + 27*y_0^6*y_1^3 - 54*u_0*y_0^6 + 54*y_0^7 + 54*y_0^6*y_1 - 27*y_0^5*y_1^2 - 27*y_0^4*u_1*y_1^2 + 27*y_0^4*y_1^3 + 27*y_0^4*y_1^2*y_2 + 27*u_0^2*y_0^3 - 54*u_0*y_0^4 - 54*u_0*y_0^3*y_1 + 27*y_0^5 + 54*y_0^4*y_1 + 18*y_0^3*u_1*y_1 + 9*y_0^3*y_1^2 - 18*y_0^3*y_1*y_2 + 9*y_0^2*u_1^2*y_1 - 18*y_0^2*u_1*y_1^2 - 18*y_0^2*u_1*y_1*y_2 + 9*y_0^2*y_1^3 + 18*y_0^2*y_1^2*y_2 + 9*y_0^2*y_1*y_2^2 + 4*y_0^3 - 3*y_0*u_1^2 + 6*y_0*u_1*y_1 + 6*y_0*u_1*y_2 - 3*y_0*y_1^2 - 6*y_0*y_1*y_2 - 3*y_0*y_2^2 - u_1^3 + 3*u_1^2*y_1 + 3*u_1^2*y_2 - 3*u_1*y_1^2 - 6*u_1*y_1*y_2 - 3*u_1*y_2^2 + y_1^3 + 3*y_1^2*y_2 + 3*y_1*y_2^2 + y_2^3))
    p = 100*y_0^4 + 10*y_0^3*y_1 + 410//3*y_0*y_2 + 10*y_0*y_3 - 110*y_1^2 - 10*y_1*y_2
    @time res = diffreduce(p, pbr)
    expected = parent_ring_change(-100*y_2^2*y_1^2+70*y_2^2*y_0^4+380//3*y_2^2*y_0^2-200*y_2*y_1^3-180*y_2*y_0^7-380*y_2*y_0^5-270*y_1^4*y_0^6-1080*y_1^4*y_0^4-630*y_1^4*y_0^2+810*y_1^3*y_0^9-2520*y_1^3*y_0^7-2910*y_1^3*y_0^5-440*y_1^3*y_0^3+220//3*y_1^3*y_0+90*y_1^2*y_0^8+3690*y_1^2*y_0^6-1330*y_1^2*y_0^4+380*y_1^2*y_0^2+1080*y_1*y_0^9-6540*y_1*y_0^7-7220*y_1*y_0^5+810*u_0^2*y_0^6-3420*u_0^2*y_0^4-1620*u_0*y_0^9+5220*u_0*y_0^7+6840*u_0*y_0^5-10*u_1^3*y_1-30*u_1^3*y_0^3+380//3*u_1^3*y_0-80*u_1^2*y_1^2+10*u_1^2*y_0^4+380*u_1^2*y_0^2+190*u_1*y_1^3+180*u_1*y_0^7+380*u_1*y_0^5-3300*y_0^6-1520//3*y_0^4-20*u_2*y_2*u_1*y_0+60*u_2*y_2*y_1*y_0^3+20*u_2*y_2*y_1*y_0-60*u_2*u_1*y_1*y_0^3-20*u_2*u_1*y_1*y_0-360*y_2*u_1*y_1*y_0^5+1380*y_2*u_1*y_1*y_0^3+1580//3*y_2*u_1*y_1*y_0-30*y_2^2*y_1^2*y_0^2+90*y_2^2*y_1*y_0^5-340*y_2^2*y_1*y_0^3-380//3*y_2^2*y_1*y_0-180*y_2*y_1^3*y_0^4-660*y_2*y_1^3*y_0^2+540*y_2*y_1^2*y_0^7-1860*y_2*y_1^2*y_0^5-1380*y_2*y_1^2*y_0^3-160//3*y_2*y_1^2*y_0+240*y_2*y_1*y_0^6+1400*y_2*y_1*y_0^4+1520//3*y_2*y_1*y_0^2+10*u_2*y_2^2*y_0-20*u_2*y_2*y_0^2+10*u_2*u_1^2*y_0+20*u_2*u_1*y_0^2+90*u_2*y_1^2*y_0^5+60*u_2*y_1^2*y_0^3+10*u_2*y_1^2*y_0-60*u_2*y_1*y_0^4-20*u_2*y_1*y_0^2-10*y_2^2*u_1*y_1-30*y_2^2*u_1*y_0^3+380//3*y_2^2*u_1*y_0+20*y_2*u_1^2*y_1+60*y_2*u_1^2*y_0^3-760//3*y_2*u_1^2*y_0+180*y_2*u_1*y_1^2-80*y_2*u_1*y_0^4-1520//3*y_2*u_1*y_0^2+180*y_2*u_0*y_0^4+30*u_1^2*y_1^2*y_0^2+270*u_1^2*y_1*y_0^5-1040*u_1^2*y_1*y_0^3-400*u_1^2*y_1*y_0+90*u_1*y_1^3*y_0^4+600*u_1*y_1^3*y_0^2-810*u_1*y_1^2*y_0^7+2820*u_1*y_1^2*y_0^5+2170*u_1*y_1^2*y_0^3+200*u_1*y_1^2*y_0-60*u_1*y_1*y_0^6-2100*u_1*y_1*y_0^4-760*u_1*y_1*y_0^2-180*u_1*u_0*y_0^4-1080*y_1*u_0*y_0^6+7020*y_1*u_0*y_0^4-6030*y_0^8+810*y_0^12-1800*y_0^10-100*y_1^4, parent(res))
    @test total_degree(divexact(res, expected)) == 0

    ode = @ODEmodel(
        x1'(t) = x1(t) + 2 * x1(t) * x2(t) + u(t),
        x2'(t) = x2(t) + 3 * x1(t) * x2(t),
        y(t) = x1(t)
    )
    pbr = PBRepresentation(ode, find_ioequations(ode))
    R, (y_0, y_1, y_2, y_3, y_4, u_0, u_3) = Nemo.PolynomialRing(Nemo.QQ, ["y_0", "y_1", "y_2", "y_3", "y_4", "u_0", "u_3"])
    io_switch!(pbr)
    @time res = diffreduce(u_3, pbr)
    expected = parent_ring_change(y_0*(27*u_0*y_0^4+27*y_0^5-27*y_0^4*y_1+27*u_0*y_0^3+54*u_0*y_0^2*y_1+27*y_0^4+27*y_0^3*y_1-54*y_0^2*y_1^2+9*u_0*y_0^2+27*u_0*y_0*y_1+12*u_0*y_0*y_2+9*u_0*y_1^2+9*y_0^3+18*y_0^2*y_1+12*y_0^2*y_2-18*y_0*y_1^2-12*y_0*y_1*y_2-9*y_1^3+u_0*y_0+3*u_0*y_1+3*u_0*y_2+u_0*y_3+y_0^2+2*y_0*y_1+3*y_0*y_2+y_0*y_4-3*y_1^2-3*y_1*y_2-y_1*y_3), parent(res))
    @test total_degree(divexact(res, expected)) == 0
end


text:
abstract type AbstractLayout end

export SpectralLayout, RandomLayout, PrecomputedLayout, KnnGraphLayout

"""
    SpectralLayout <: AbstractLayout

Initializes the embedding using the spectral layout method. It could be costly in very large datasets.
"""
struct SpectralLayout <: AbstractLayout end

"""
    RandomLayout <: AbstractLayout

Initializes the embedding using a random set of points. It may converge slowly
"""
struct RandomLayout <: AbstractLayout end

"""
    PrecomputedLayout <: AbstractLayout

Initializes the embedding using a previously computed layout, i.e., (maxoutdim, n_points) matrix.
"""
struct PrecomputedLayout <: AbstractLayout
    init::Matrix{Float32}
end

"""
    KnnGraphLayout <: AbstractLayout

A lattice like + clouds of points initialization that uses the computed all-knn graph.
This layout initialization is a simple proof of concept, so please use it under this assumption.
"""
struct KnnGraphLayout <: AbstractLayout
end

function initialize_embedding(::RandomLayout, graph::AbstractMatrix, knns, dists, maxoutdim::Integer)
    n = size(knns, 2)
    rand(-10f0:eps(Float32):10f0, maxoutdim, n)
end

function initialize_embedding(layout::PrecomputedLayout, graph::AbstractMatrix, knns, dists, maxoutdim::Integer)
    k, n = size(layout.init)
    n_ = size(knns, 2)
    if maxoutdim !=  k || n != n_
        error("PrecomputedLayout matrix ($k, $n) doesn't matches ($maxoutdim, $(n_))")
    end

    layout.init
end

function initialize_embedding(::SpectralLayout, graph::AbstractMatrix{T}, knns, dists, maxoutdim::Integer) where {T}
    local embed

    try
        embed = spectral_layout(graph, maxoutdim)
        r = randn(T, size(embed))
        r .*= 0.0001
        # expand
        embed .*= 10 / maximum(embed)
        embed .+= r
    catch e
        @info "$e\nError encountered in spectral_layout; defaulting to random layout"
        embed = initialize_embedding(RandomLayout(), graph, knns, dists, maxoutdim)
    end

    embed
end

"""
    spectral_layout(graph, embed_dim) -> embedding

Initialize the graph layout with spectral embedding.
"""
function spectral_layout(graph::SparseMatrixCSC{T},
                         embed_dim::Integer) where {T<:Real}
    graph_f64 = convert.(Float64, graph)
    D_ = Diagonal(dropdims(sum(graph_f64; dims=2); dims=2))
    D = inv(sqrt(D_))
    # normalized laplacian
    L = Symmetric(I - D*graph*D)

    k = embed_dim+1
    num_lanczos_vectors = max(2k+1, round(Int, sqrt(size(L, 1))))
    # get the 2nd - embed_dim+1th smallest eigenvectors
    eigenvals, eigenvecs = eigs(L; nev=k,
                                   ncv=num_lanczos_vectors,
                                   which=:SM,
                                   tol=1e-4,
                                   v0=ones(Float64, size(L, 1)),
                                   maxiter=size(L, 1)*5)
    layout = permutedims(eigenvecs[:, 2:k])::Array{Float64, 2}
    return convert.(T, layout)
end

function revknn_frequencies_(knns)
    n = size(knns, 2)
    F = zeros(Int32, n)
    for i in 1:n
        for objID in view(knns, :, i)
            F[objID] += 1  # appoximation to indegree centrality
        end
    end

    F
end

function knn_fast_components_(knns)  ## approximated connected components
    C = Dict{Int32,Vector{Int32}}()
    for i in 1:size(knns, 2)
        prev = i
        links = @view knns[:, i]

        while (curr = minimum(links); prev > curr)
            links = @view knns[:, curr]
            prev = curr
        end

        lst = get(C, prev, nothing)
        if lst === nothing
            C[prev] = Int32[i]
        else
            push!(lst, i)
        end
    end

    C
end

function initialize_embedding(::KnnGraphLayout, graph::AbstractMatrix, knns, dists, maxoutdim::Integer)
    n = size(knns, 2)
    embed = zeros(Float32, maxoutdim, n)
    E = MatrixDatabase(embed)
    C = knn_fast_components_(knns)
    F = revknn_frequencies_(knns)

    for lst in values(C)
        _, p = findmax(p -> F[p], lst)
        i = lst[p]
        rand_point_lattice_(E[i])  # most popular elements are in the center
        for p in lst
            if p != i
                rand_point_cloud_(E[p], E[i], 0.25f0)
            end
        end
    end

    embed
end

function rand_point_lattice_(V)
    @inbounds for i in eachindex(V)
        V[i] = rand(-10f0:0.5f0:10f0)
    end
end

function rand_point_cloud_(V, center, min_dist::Float32)
    @inbounds for i in eachindex(V)
        V[i] = center[i] + rand(-min_dist:eps(Float32):min_dist)
    end
end

"""
    initialize_embedding(graph::AbstractMatrix, ref_embedding::Matrix) -> embedding

Initialize an embedding of points corresponding to the columns of the `graph`, by taking weighted means of
the columns of `ref_embedding`, where weights are values from the rows of the `graph`.

The resulting embedding will have shape `(size(ref_embedding, 1), size(graph, 2))`, where `size(ref_embedding, 1)`
is the number of components (dimensions) of the `reference embedding`, and `size(graph, 2)` is the number of
samples in the resulting embedding. Its elements will have type T.
"""
function initialize_embedding(graph::AbstractMatrix, ref_embedding::Matrix{Float32})
    (ref_embedding * graph) ./ (sum(graph, dims=1) .+ eps(Float32))
end


text:
"""
    distance(q₁, q₂)
    distance2(q₁, q₂)

Measure the "distance" between two quaternions, or the squared distance with
`distance2`.

By default, this function just returns the natural measure in the *additive*
group of quaternions:
```julia
abs2(q₁ - q₂)
```
If both arguments are `Rotor`s, the function returns the natural measure in the
*rotation* group, which is roughly
```julia
abs2(log(q₁ / q₂))
```
Note that for `Rotor`s, this method is (efficiently) independent of the
scaling of `q₁` and `q₂`, including up to factors of -1, as is appropriate for
the rotation group.

# Examples
```jldoctest example
julia> distance(imx, imy)
1.4142135623730951
julia> distance(Rotor(imx), Rotor(imy))
1.5707963267948966
julia> distance(imz, -imz)
2.0
julia> distance(Rotor(imz), Rotor(-imz))
0.0
```
"""
distance(q₁::AbstractQuaternion, q₂::AbstractQuaternion) = √distance2(q₁, q₂)
distance2(q₁::AbstractQuaternion, q₂::AbstractQuaternion) = abs2(q₁ - q₂)
distance(q₁::Rotor, q₂::Rotor) = √distance2(q₁, q₂)
distance2(q₁::Rotor, q₂::Rotor) = _abs2_small_vec_log(q₁ / q₂)

@inline function _abs2_small_vec_log(q::Rotor{T}) where {T}
    # Like `min(abs2(log(q)), abs2(log(-q)))`, but assumes the norm of `q` is 1
    atan(absvec(q), abs(q.w))^2
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/jacobian_functor.jl
# Defines methods to used to compute Jacobians of Taylor coeffients.
#############################################################################

"""
JacTaylorFunctor!

A callable structure used to evaluate the Jacobian of Taylor cofficients. This
also contains some addition fields to be used as inplace storage when computing
and preconditioning paralleliped based methods to representing enclosure of the
pODEs (Lohner's QR, Hermite-Obreschkoff, etc.). The constructor given by
`JacTaylorFunctor!(g!, nx::Int, np::Int, k::Val{K}, t::T, q::Q)` may be used
were type `T` should use type `Q` for internal computations. The order of the
TaylorSeries is `k`, the right-hand side function is `g!`, `nx` is the number
of state variables, `np` is the number of parameters.

$(TYPEDFIELDS)
"""
mutable struct JacTaylorFunctor!{F <: Function, N, T <: Real, S <: Real, NY}
    "Right-hand side function for pODE which operates in place as g!(dx,x,p,t)"
    g!::F
    "Dimensionality of x"
    nx::Int
    "Dimensionality of p"
    np::Int
    "Order of TaylorSeries"
    s::Int
    "In-place temporary storage for Taylor coefficient calculation"
    out::Vector{S}
    "Variables y = (x,p)"
    y::Vector{S}
    "State variables x"
    x::Vector{Dual{Nothing,S,NY}}
    "Decision variables p"
    p::Vector{Dual{Nothing,S,NY}}
    "Storage for sum of Jacobian w.r.t x"
    Jxsto::Matrix{S}
    "Storage for sum of Jacobian w.r.t p"
    Jpsto::Matrix{S}
    "Temporary for transpose of Jacobian w.r.t y"
    tjac::Matrix{S}
    "Storage for vector of Jacobian w.r.t x"
    Jx::Vector{Matrix{S}}
    "Storage for vector of Jacobian w.r.t p"
    Jp::Vector{Matrix{S}}
    "Jacobian Result from DiffResults"
    result::MutableDiffResult{1, Vector{S}, Tuple{Matrix{S}}}
    "Jacobian Configuration for ForwardDiff"
    cfg::JacobianConfig{Nothing,S,NY,Tuple{Vector{Dual{Nothing,S,NY}},Vector{Dual{Nothing,S,NY}}}}
    "Store temporary STaylor1 vector for calculations"
    xtaylor::Vector{STaylor1{N,Dual{Nothing,S,NY}}}
    "Store temporary STaylor1 vector for calculations"
    xaux::Vector{STaylor1{N,Dual{Nothing,S,NY}}}
    "Store temporary STaylor1 vector for calculations"
    dx::Vector{STaylor1{N,Dual{Nothing,S,NY}}}
    taux::Vector{STaylor1{N,T}}
    t::Float64
    "Intermediate storage to avoid allocations in Taylor coefficient calc"
    vnxt::Vector{Int64}
    "Intermediate storage to avoid allocations in Taylor coefficient calc"
    fnxt::Vector{Float64}
end

function JacTaylorFunctor!(g!, nx::Int, np::Int, k::Val{K}, t::T, q::Q) where {K, T <: Number, Q <: Number}
    x0 = zeros(T, nx)
    xd0 = zeros(Dual{Nothing, T, nx + np}, nx)
    out = zeros(T, nx*(K + 1))
    y = zeros(T, nx + np)
    x = zeros(Dual{Nothing, T, nx + np}, nx)
    p = zeros(Dual{Nothing, T, nx + np}, np)
    Jxsto = zeros(T, nx, nx)
    Jpsto = zeros(T, nx, np)
    tjac = zeros(T, np + nx, nx*(K + 1))
    cfg = JacobianConfig(nothing, out, zeros(T, nx + np))
    result = JacobianResult(out, zeros(T, nx + np))
    Jx = Matrix{T}[]
    Jp = Matrix{T}[]

    temp = zero(Dual{Nothing, T, nx + np})
    taux = [STaylor1(zero(Q), Val(K))]
    xtaylor = STaylor1.(xd0, Val(K))
    dx = STaylor1.(xd0, Val(K))
    xaux = STaylor1.(xd0, Val(K))
    for i in 1:(K + 1)
        push!(Jx, zeros(T,nx,nx))
        push!(Jp, zeros(T,nx,np))
    end
    t = 0.0
    vnxt = zeros(Int, nx)
    fnxt = zeros(Float64, nx)
    return JacTaylorFunctor!{typeof(g!), K+1, Q, T, nx + np}(g!, nx, np, K, out, y, x, p,
                                                             Jxsto, Jpsto, tjac, Jx, Jp,
                                                             result, cfg, xtaylor, xaux,
                                                             dx, taux, t, vnxt, fnxt)
end

function (d::JacTaylorFunctor!{F,K,T,S,NY})(out::AbstractVector{Dual{Nothing,S,NY}},
                                            y::AbstractVector{Dual{Nothing,S,NY}}) where {F <: Function,
                                                                                K, T <: Real, S, NY}

    copyto!(d.x, 1, y, 1, d.nx)
    copyto!(d.p, 1, y, d.nx + 1, d.np)
    val = Val{K-1}()
    for i=1:d.nx
        d.xtaylor[i] = STaylor1(d.x[i], val)
    end
    jetcoeffs!(d.g!, d.t, d.xtaylor, d.xaux, d.dx, K - 1, d.p, d.vnxt, d.fnxt)
    for q = 1:K
        for i = 1:d.nx
            indx = d.nx*(q - 1) + i
            out[indx] = d.xtaylor[i].coeffs[q]
        end
    end

    return nothing
end

"""
jacobian_taylor_coeffs!

Computes the Jacobian of the Taylor coefficients w.r.t. y = (x,p) storing the
output inplace to `result`. A JacobianConfig object without tag checking, cfg,
is required input and is initialized from `cfg = ForwardDiff.JacobianConfig(nothing, out, y)`.
The JacTaylorFunctor! used for the evaluation is `g` and inputs are `x` and `p`.
"""
function jacobian_taylor_coeffs!(g::JacTaylorFunctor!{F,K,T,S,NY}, X::Vector{S}, P, t::T) where {F,K,T,S,NY}

    # copyto! is used to avoid allocations
    copyto!(g.y, 1, X, 1, g.nx)
    copyto!(g.y, g.nx + 1, P, 1, g.np)
    g.t = t
    # other AD schemes may be usable as well but this is a length(g.out) >> nx + np
    # situtation typically
    jacobian!(g.result, g, g.out, g.y, g.cfg)

    # reset sum of Jacobian storage
    fill!(g.Jxsto, zero(S))
    fill!(g.Jpsto, zero(S))
    nothing
end

"""
set_JxJp!

Extracts the Jacobian of the Taylor coefficients w.r.t. x, `Jx`, and the
Jacobian of the Taylor coefficients w.r.t. p, `Jp`, from `result`. The order of
the Taylor series is `s`, the dimensionality of x is `nx`, the dimensionality of
p is `np`, and `tjac` is preallocated storage for the transpose of the Jacobian
w.r.t. y = (x,p).
"""
function set_JxJp!(g::JacTaylorFunctor!{F,K,T,S,NY}, X::Vector{S}, P, t) where {F,K,T,S,NY}

    jacobian_taylor_coeffs!(g, X, P, t)
    jac = g.result.derivs[1]
    for i = 1:(g.s + 1)
        for q = 1:g.nx
            for z = 1:g.nx
                g.Jx[i][q, z] = jac[q + g.nx*(i-1), z]
            end
            for z = 1:g.np
                g.Jp[i][q, z] = jac[q + g.nx*(i-1), g.nx + z]
            end
        end
    end
    nothing
end


text:
function generate_traindataset1d(network::AbstractNetworkTyp, training_setup::Traindata_Settings)
    @unpack input_size = network
    @unpack meshes, polydegrees, length_data = training_setup

    to = TimerOutput()

    # Preallocate some Vectors & Arrays
    index_output = 1
    num_troubled = 0
    X = Array{Float64, 2}(undef, input_size, length_data)
    Y = Array{Float64, 2}(undef, 2, length_data)

    # loop over meshs
    @timeit to "loop over meshes" for mesh in meshes

        leaf_cell_ids = leaf_cells(mesh.tree)
        n_elements = length(leaf_cell_ids)
         # Construct cell -> element mapping for easier algorithm implementation
        c2e = zeros(Int, length(mesh.tree))
         for element_id in 1:n_elements
             c2e[leaf_cell_ids[element_id]] = element_id
         end

        #loop over polydegrees
        for r in polydegrees
            n_nodes = r + 1
            nodes, _ = gauss_lobatto_nodes_weights(n_nodes)
            node_coordinates = zeros(n_nodes)

            #loop over functions
            @timeit to "loop over functions" for func in 1:26
                u(x) = trainfunction1d(func,x)

                for element_id in 1:n_elements
                    cell_id = leaf_cell_ids[element_id]
                    dx = length_at_cell(mesh.tree, cell_id)

                    for i in 1:n_nodes
                        node_coordinates[i] = (mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
                    end

                    @timeit to "generate Inputdata" generate_inputdata1d!(X, index_output, network, u, node_coordinates, element_id, mesh, c2e)

                    # Create label Yi
                    if (func > 18 && func <=26) && (  0 >= node_coordinates[1] && 0 <= node_coordinates[end])
                        Y[1,index_output] = 1
                        Y[2,index_output] = 0
                        num_troubled += 1
                    else
                        Y[1,index_output] = 0
                        Y[2,index_output] = 1
                    end

                    index_output += 1
                end
            end

            @timeit to "loop troubled cells" for t in 1:1100
                ul = rand(Uniform(-1,1))
                ur = rand(Uniform(-1,1))
                x0 = rand(Uniform(-0.75, 0.75))
                u(x) = troubledcellfunctionstep1d(x, ul, ur, x0)

                for element_id in 1:n_elements
                    cell_id = leaf_cell_ids[element_id]
                    dx = length_at_cell(mesh.tree, cell_id)

                    # Calculate node coordinates
                    for i in 1:n_nodes
                        node_coordinates[i] = (mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
                    end

                    if good_cell1d(node_coordinates, dx, x0)
                        continue
                    end

                    @timeit to "generate Inputdata" generate_inputdata1d!(X, index_output, network, u, node_coordinates, element_id, mesh, c2e)

                    # Create label Yi
                    Y[1,index_output] = 1
                    Y[2,index_output] = 0
                    num_troubled += 1

                    index_output += 1
                end
            end
        end
    end

    @timeit to "save Dataset" h5open("datasets/1d/traindata1d$(typeof(network)).h5", "w") do file
        write(file, "X", @view X[:,1:index_output-1])
        write(file, "Y", @view Y[:,1:index_output-1])
    end

    print_output(training_setup,index_output-1,num_troubled)

    show(to)

end


text:
module ByteWorks

function floattobytes(value::Float64)::Array{UInt8,1}
    return floatstobytes([value])
end

function floatstobytes(values::Array{Float64,1})::Array{UInt8,1}
    buf = IOBuffer()
    write(buf, values)
    bytes = take!(buf)
    return bytes
end

function bytestofloats(bytes::Array{UInt8,1})::Array{Float64,1}
    return reinterpret(Float64, bytes)
end

function validate(bytes::Array{UInt8,1})::Bool
    fvals = bytestofloats(bytes)
    all(map(x -> !isnan(x), fvals))
end

function checkandrestore(
    bytes::Array{UInt8, 1},
    lowerbound::Array{Float64, 1},
    upperbound::Array{Float64, 1})::Array{UInt8, 1}

    if validate(bytes)
        return bytes
    else
        return map( (L, U) -> L + rand() * (U - L), lowerbound, upperbound) |> floatstobytes
    end
end

end # end of module


text:
module DynamicalSystems

# Creates the functions for the harmonic oscillator with a sinusoidal driving
# force.
#
# x'' + beta x' + omega^2 x = gamma sin(omega t)
#
# Args:
#   omega0: The frequency term.
#   beta: The damping term (default: 0.0).
#   gamma: The driving force amplitude (default: 0.0).
#   omega: The driving force frequency (default: 1.0).
#
# Returns: An array of two functions defining the harmonic oscillator system.
function harmonicOscillator(omega0::AbstractFloat;
    beta::AbstractFloat=0.0, gamma::AbstractFloat=0.0,
    omega::AbstractFloat=1.0)

    return [(t,x) -> x[2],
            (t,x) -> (gamma*sin(omega*t) -
                 beta*x[2] -
                 (omega0^2)*x[1])];

end #Close simpleHarmonicOscillator.

# Creates the functions for a Duffing oscillator.
#
# x'' + delta x' + alpha x + beta x^3 = gamma sin(omega t)
#
# Args:
#   delta: The damping term.
#   gamma: The driving amplitude.
#   alpha: The stiffness (default: 1.0).
#   beta: The restoring force non-linearity (default: -1.0).
#   omega: The driving force frequency (default: 1.0).
#
# Returns: An array of two functions defining the Duffing oscillator system.
function duffingOscillator(delta::AbstractFloat, gamma::AbstractFloat;
    alpha::AbstractFloat=-1.0, beta::AbstractFloat=1.0,
    omega::AbstractFloat=1.0)

    return [(t,x) -> x[2],
     (t,x) -> (gamma*sin(omega*t) -
              delta*x[2] -
              alpha*x[1] -
              beta*(x[1]^3))];

end #Close duffingOscillator.

# Creates the functions for a Lorenz attractor.
#
# x' = sigma (y - x)
# y' = x (rho - z) - y
# z' = x y - beta z
#
# Args:
#   sigma: The sigma parameter of the system.
#   rho: The rho parameter of the system.
#   beta: The beta parameter of the system.
#
# Returns: An array of three function defining the Lorenz attractor system.
function lorenzAttractor(sigma::AbstractFloat, rho::AbstractFloat,
    beta::AbstractFloat)

    return [(t,x) -> sigma*(x[2] - x[1]),
            (t,x) -> x[1]*(rho - x[3]) - x[2],
            (t,x) -> x[1]*x[2] - beta*x[3]];
end

export harmonicOscillator, duffingOscillator, lorenzAttractor

end # module


text:
module LibNFFT

export
    NFCT,
    NFFT,
    NFST,
    apply_adjoint!,
    apply_adjoint,
    apply_direct!,
    apply_direct

# Types to identify raw plan structures parameterized by the floating-point
# type.  These must be defined before including "../deps/deps.jl".
abstract type _nfft_plan{T<:AbstractFloat} end
abstract type _nfct_plan{T<:AbstractFloat} end
abstract type _nfst_plan{T<:AbstractFloat} end

isfile(joinpath(@__DIR__, "..", "deps", "deps.jl")) ||
    error("LibNFFT not properly installed.  Please run Pkg.build(\"LibNFFT\")")
include(joinpath("..", "deps", "deps.jl"))

import FFTW

const OPENMP = false

const DEFAULT_FFTW_FLAGS = FFTW.ESTIMATE
const MINIMAL_FFTW_FLAGS = FFTW.DESTROY_INPUT
const ALLOWED_FFTW_FLAGS = (FFTW.ESTIMATE | FFTW.MEASURE |
                            FFTW.PATIENT  | FFTW.EXHAUSTIVE)

const DEFAULT_NFFT_FLAGS = (MALLOC_X  | MALLOC_F_HAT | MALLOC_F |
                            PRE_PHI_HUT | PRE_PSI)
const MINIMAL_NFFT_FLAGS = (_FFTW_INIT | _FFT_OUT_OF_PLACE  |
                            (OPENMP ? _OMP_BLOCKWISE_ADJOINT : 0))
const ALLOWED_NFFT_FLAGS = (MALLOC_X    | MALLOC_F | MALLOC_F_HAT |
                            PRE_PHI_HUT | FG_PSI   | PRE_LIN_PSI  |
                            PRE_FG_PSI  | PRE_PSI  | PRE_FULL_PSI |
                            SORT_NODES)

const Floats = Union{Cfloat,Cdouble}
const RealOrComplex{T} = Union{T,Complex{T}}
const SinglePrecision{T} = RealOrComplex{Cfloat}
const DoublePrecision{T} = RealOrComplex{Cdouble}

"""

`LibNFFT.NullPointerError()` is the exception thrown when a pointer has an
invalid `NULL` value.

"""
struct NullPointerError <: Exception end

"""
    LibNFFT.version(T) -> (major, minor, path)

yields the version of the NFFT library for floating-point type `T`.

""" version

"""

Abstract type `LibNFFT.Plan{T,d,E}` is the super-type of all plans implementing
one of the transforms provided by the NFFT library.

This type is parameterized by `T` the floating-point type for the computations,
`d` and `E` the number of dimensions and the type of the elements of the input
arrays.  Parameter `E` is either `T` or `Complex{T}`.  Because, `E` depends on
`T`, `T` must be specified first in the signature.

"""
abstract type Plan{T,d,E} end

# `NFFT{T,d,E}` (nonequispaced fast Fourier transform).
mutable struct NFFT{T<:Floats,d,E<:Union{T,Complex{T}}} <: Plan{T,d,E}
    ptr::Ptr{_nfft_plan{T}}    # plan to compute NFFT
    f_hat::Array{Complex{T},d} # input samples
    f::Vector{Complex{T}}      # result of transform
    x::Matrix{T}               # nodes
    ready::Bool

    # The inner constructor takes the same arguments as the nfft_init_guru
    # function except that it allocates the structure.
    function NFFT{T,d,E}(inpdims::NTuple{d,Int}, # input dimensions (N)
                         M::Int,                 # number of nodes (M)
                         ovrdims::NTuple{d,Int}, # oversampled dimensions (n)
                         m::Int,                 # window size
                         flags::Cuint,           # NFFT flags
                         fftw_flags::Cuint       # FFTW flags
                         ) where {T<:Floats,d,E<:Union{T,Complex{T}}}
        # Create object and attach a finalizer to it.
        ptr, f_hat, f, x = _create(_nfft_plan{T}, inpdims, M, ovrdims, m,
                                   flags, fftw_flags)
        return finalizer(finalize, new{T,d,E}(ptr, f_hat, f, x, false))
    end
end

# `NFCT{T,d}` (nonequispaced fast cosine transform).
mutable struct NFCT{T<:Floats,d} <: Plan{T,d,T}
    ptr::Ptr{_nfct_plan{T}}   # plan to compute NFCT
    f_hat::Array{T,d}         # input samples
    f::Vector{T}              # result of transform
    x::Matrix{T}              # nodes
    ready::Bool

    # The inner constructor takes the same arguments as the nfct_init_guru
    # function except that it allocates the structure.
    function NFCT{T,d}(inpdims::NTuple{d,Int}, # input dimensions (N)
                       M::Int,                 # number of nodes (M)
                       ovrdims::NTuple{d,Int}, # oversampled dimensions (n)
                       m::Int,                 # window size
                       flags::Cuint,           # NFCT flags
                       fftw_flags::Cuint       # FFTW flags
                       ) where {T<:Floats,d,E<:Union{T,Complex{T}}}
        # Create object and attach a finalizer to it.
        ptr, f_hat, f, x = _create(_nfct_plan{T}, inpdims, M, ovrdims, m,
                                   flags, fftw_flags)
        return finalizer(finalize, new{T,d}(ptr, f_hat, f, x, false))
    end
end

# `NFST{T,d}` (nonequispaced fast sine transform).
mutable struct NFST{T<:Floats,d} <: Plan{T,d,T}
    ptr::Ptr{_nfst_plan{T}}   # plan to compute NFST
    f_hat::Array{T,d}         # input samples
    f::Vector{T}              # result of transform
    x::Matrix{T}              # nodes
    ready::Bool

    # The inner constructor takes the same arguments as the nfst_init_guru
    # function except that it allocates the structure.
    function NFST{T,d}(inpdims::NTuple{d,Int}, # input dimensions (N)
                       M::Int,                 # number of nodes (M)
                       ovrdims::NTuple{d,Int}, # oversampled dimensions (n)
                       m::Int,                 # window size
                       flags::Cuint,           # NFST flags
                       fftw_flags::Cuint       # FFTW flags
                       ) where {T<:Floats,d,E<:Union{T,Complex{T}}}
        # Create object and attach a finalizer to it.
        ptr, f_hat, f, x = _create(_nfst_plan{T}, inpdims, M, ovrdims, m,
                                   flags, fftw_flags)
        return finalizer(finalize, new{T,d}(ptr, f_hat, f, x, false))
    end
end

const BasicPlans{T,d,E} = Union{NFFT{T,d,E},NFCT{T,d},NFST{T,d}}

function _create(::Type{T},
                 inpdims::NTuple{d,Int}, # input dimensions (N)
                 M::Int,                 # number of nodes (M)
                 ovrdims::NTuple{d,Int}, # oversampled dimensions (n)
                 m::Int,                 # window size
                 flags::Cuint,           # NFFT flags
                 fftw_flags::Cuint       # FFTW flags
                 ) where {T<:Union{_nfft_plan,_nfct_plan,_nfst_plan},
                          d,E<:Union{T,Complex{T}}}
    # Notations:
    # `d` is the number of dimensions (the rank);
    # `N = [N1,N2,...,Nd]` are the dimensions of the result;
    # `M` are the number of nodes;
    # `n = [n1,n2,...,nd]` are the oversampled dimensions for the FFT;
    # `m` is the size of the window.

    # Check rank, number of nodes and window size.
    d ≥ 1 || argument_error("invalid number of dimensions")
    M ≥ 1 || argument_error("invalid number of nodes")
    m ≥ 2 || argument_error("window size must be at least 2")

    # Check dimensions, convert them to C-int integer type and reorder
    # them in the same order as expected by FFTW/NFFT.
    N = Vector{Cint}(undef, d)
    n = Vector{Cint}(undef, d)
    for i in 1:d
        inpdim::Int = inpdims[i]
        ovrdim::Int = ovrdims[i]
        inpdim < 1 && argument_error("invalid dimension(s)")
        ovrdim ≤ inpdim && argument_error("oversampled dimension(s) too small")
        isodd(inpdim) && argument_error("dimensions must be even numbers")
        N[d + 1 - i] = inpdim
        n[d + 1 - i] = ovrdim
    end

    # Create the plan and initialize it using the "guru" interface.
    ptr = _allocate(T)
    _init_guru(ptr, d, N, M, n, m, flags, fftw_flags)

    # Retrieve the work arrays.
    f_hat = _init_array(_pointer_to_f_hat(ptr), inpdims, (flags & MALLOC_F_HAT) != 0)
    f = _init_array(_pointer_to_f(ptr), (Int(M),), (flags & MALLOC_F) != 0)
    x = _init_array(_pointer_to_x(ptr), (d, Int(M)), (flags & MALLOC_X) != 0)

    return (ptr, f_hat, f, x)
end

function _init_array(ptr::Ptr{Ptr{T}},
                     dims::NTuple{N,Int},
                     wrap::Bool) :: Array{T,N} where {T,N}
    if wrap
        return unsafe_wrap(Array, unsafe_load(ptr), dims; own = false)
    else
        A = Array{T,N}(undef, dims)
        unsafe_store!(ptr, pointer(A))
        return A
    end
end

# Outer constructors.
function NFFT(M::Integer,                          # number of nodes
              E::Union{Type{T}, Type{Complex{T}}}, # input element type
              dims::NTuple{d,Integer};             # input dimensions
              m::Integer = 6,                      # window size
              sigma::Real = 2.0,                   # oversampling factor
              flags::Integer = DEFAULT_NFFT_FLAGS,
              fftw_flags::Integer = DEFAULT_FFTW_FLAGS) where {T<:Floats,d}
    # Create NFFT plan with minimal set of options enforced.
    inpdims, ovrdims = _fix_dims(NFFT, dims, sigma)
    return NFFT{T,d,E}(inpdims, Int(M), ovrdims, Int(m),
                       Cuint(flags | MINIMAL_NFFT_FLAGS),
                       Cuint(fftw_flags | MINIMAL_FFTW_FLAGS))
end

function NFFT(x::AbstractMatrix{<:AbstractFloat},  # coordinates of nodes
              E::Union{Type{T}, Type{Complex{T}}}, # input element type
              dims::NTuple{d,Integer};             # input dimensions
              kwds...) where {T<:Floats,d}
    size(x, 1) == d || dimension_mismatch("bad first dimension for node coordinates")
    M = size(x, 2)
    return set_x!(NFFT(M, E, dims; kwds...), x)
end

function NFCT(M::Integer,                          # number of nodes
              ::Type{T},                           # floating-point type
              dims::NTuple{d,Integer};             # input dimensions
              m::Integer = 6,                      # window size
              sigma::Real = 2.0,                   # oversampling factor
              flags::Integer = DEFAULT_NFFT_FLAGS,
              fftw_flags::Integer = DEFAULT_FFTW_FLAGS) where {T<:Floats,d}
    # Create NFCT plan with minimal set of options enforced.
    inpdims, ovrdims = _fix_dims(NFCT, dims, sigma)
    return NFCT{T,d}(inpdims, Int(M), ovrdims, Int(m),
                     Cuint(flags | MINIMAL_NFFT_FLAGS),
                     Cuint(fftw_flags | MINIMAL_FFTW_FLAGS))
end

function NFCT(x::AbstractMatrix{<:AbstractFloat},  # coordinates of nodes
              ::Type{T},                           # floating-point type
              dims::NTuple{d,Integer};             # input dimensions
              kwds...) where {T<:Floats,d}
    size(x, 1) == d || dimension_mismatch("bad first dimension for node coordinates")
    M = size(x, 2)
    return set_x!(NFCT(M, T, dims; kwds...), x)
end

function NFST(M::Integer,                          # number of nodes
              ::Type{T},                           # floating-point type
              dims::NTuple{d,Integer};             # input dimensions
              m::Integer = 6,                      # window size
              sigma::Real = 2.0,                   # oversampling factor
              flags::Integer = DEFAULT_NFFT_FLAGS,
              fftw_flags::Integer = DEFAULT_FFTW_FLAGS) where {T<:Floats,d}
    # Create NFST plan with minimal set of options enforced.
    inpdims, ovrdims = _fix_dims(NFST, dims, sigma)
    return NFST{T,d}(inpdims, Int(M), ovrdims, Int(m),
                     Cuint(flags | MINIMAL_NFFT_FLAGS),
                     Cuint(fftw_flags | MINIMAL_FFTW_FLAGS))
end

function NFST(x::AbstractMatrix{<:AbstractFloat},  # coordinates of nodes
              ::Type{T},                           # floating-point type
              dims::NTuple{d,Integer};             # input dimensions
              kwds...) where {T<:Floats,d}
    size(x, 1) == d || dimension_mismatch("bad first dimension for node coordinates")
    M = size(x, 2)
    return set_x!(NFST(M, T, dims; kwds...), x)
end

function finalize(ths::Union{NFFT{T,d},NFCT{T,d},NFST{T,d}}) where {T,d}
    if (ptr = ths.ptr) != C_NULL
        ths.ptr = C_NULL # to avoid finalize more than once
        _finalize(ptr)
        _free(ptr)
    end
    ths.f_hat = Array{Complex{T},d}(undef, ntuple(x->0, Val(d)))
    ths.f = Vector{Complex{T}}(undef, 0)
    ths.x = Matrix{T}(undef, d, 0)
    ths.ready = false
    nothing
end

function _fix_dims(::Type{<:BasicPlans},
                   dims::Tuple{Vararg{Integer}}, sigma::Real)
    # Get good dimensions for the NFFT suitable fot the FFT (hence powers of 2,
    # 3 or 5) and are larger of equal `sigma` times the dimensions in `dims`.
    sigma < 1 && argument_error("oversampling factor must be at least one")
    ovrdims = map(dim -> nextprod([2,3,5], ceil(Int, sigma*dim)), dims)
    return map(Int, dims), ovrdims
end

# Yields the number of nodes.
get_M(ths::BasicPlans) = size(ths.x, 2)

# Yields the number of dimensions of the input.
get_d(ths::Plan{T,d}) where {T,d} = d

get_flags(ths::BasicPlans) = get_flags(ths.ptr)

get_fftw_flags(ths::BasicPlans) = get_fftw_flags(ths.ptr)

isready(ths::BasicPlans) = ths.ready

argument_error(x) = throw(ArgumentError(x))
dimension_mismatch(x) = throw(DimensionMismatch(x))

Base.show(io::IO, ::MIME"text/plain", ths::NFFT{T,d,E}) where {T,d,E} =
    print(io, "LibNFFT.NFFT{$T,$d,$E} ",
          "input: $E $(size(ths.f_hat)), ",
          "output: Complex{$T} $(size(ths.f))")

Base.show(io::IO, ::MIME"text/plain", ths::NFCT{T,d}) where {T,d} =
    print(io, "LibNFFT.NFCT{$T,$d} ",
          "input: $T $(size(ths.f_hat)), ",
          "output: $T $(size(ths.f))")

Base.show(io::IO, ::MIME"text/plain", ths::NFST{T,d}) where {T,d} =
    print(io, "LibNFFT.NFST{$T,$d} ",
          "input: $T $(size(ths.f_hat)), ",
          "output: $T $(size(ths.f))")

"""
    LibNFFT.get_x(ths) -> x

yields the coordinates of the nodes in NFFT plan `ths`.  The result `x` is such
that `x[t,j]` is `t`-th coordinate of `j`-th node.

If the node coordinates have not yet been set in `ths`, all coordinates in the
result have `NaN` values.

The result `x` is a new array as coordinates are stored in reverse order by
NFFT (compared to Julia).  To re-use storage, the in-place version can be
called instead:

    LibNFFT.get_x!(x, ths) -> x

overwrites `x` with the coordinates of the nodes in NFFT plan `ths`.

See also [`NFFT`](@ref), [`NFCT`](@ref), [`NFST`](@ref),
[`LibNFFT.set_x!`](@ref).

"""
get_x(ths::NFFT{T}) where{T} = get_x!(Matrix{T}(undef, size(ths.x)), ths)

function get_x!(x::AbstractMatrix{<:AbstractFloat},
                ths::NFFT{T,d,E}) where {T,d,E}
    # Check dimensions.
    src = ths.x
    dims = size(src)
    dims[1] == d || error("corrupted NFFT structure")
    M = dims[2]
    Base.has_offset_axes(x) &&
        error("array of node coordinates has non-standard indexing")
    size(x) == dims ||
        dimension_mismatch("incompatible dimensions of node coordinates")

    # Copy node coordinates (dimensions are in reverse order in NFFT).
    if ths.ready
        @inbounds begin
            for j in 1:M
                @simd for i in 1:d
                    x[i,j] = src[(d+1)-i,j]
                end
            end
        end
    else
        fill!(x, NaN)
    end
    return x
end

@doc @doc(get_x) get_x!

"""
    LibNFFT.set_x!(ths, x) -> ths

sets the coordinates of the nodes in NFFT (or NFCT or NFST) plan `ths`.
Argument `x` is such that `x[t,j]` is `t`-th coordinate of `j`-th node.  Thus
`x` must be an array of dimensions `(d,M)` with `d` the number of dimensions of
the input argument of the direct transform and `M` the number of nodes.  All
node coordinates must be in the semi-open range `[-1/2,+1/2)`.

See also [`NFFT`](@ref), [`NFCT`](@ref), [`NFST`](@ref),
[`LibNFFT.get_x`](@ref).

"""
function set_x!(ths::BasicPlans{<:AbstractFloat,d},
                x::AbstractMatrix{T}) where {T<:AbstractFloat,d}
    # Check dimensions.  Node coordinates must be an array of dimensions
    # `(d,M)` with `d` the number of dimensions of the input argument of the
    # direct transform and `M` the number of nodes.
    dest = ths.x
    dims = size(dest)
    dims[1] == d || error("corrupted NFFT structure")
    M = dims[2]
    Base.has_offset_axes(x) &&
        error("array of node coordinates has non-standard indexing")
    size(x) == dims ||
        dimension_mismatch("incompatible dimensions of node coordinates")

    # Check values of node coordinates.
    cmax = one(T)/2
    cmin = -cmax
    @inbounds begin
        for j in 1:M
            for i in 1:d
                c = x[i,j]
                isnan(c) && error("invalid NaN node coordinate")
                isinf(c) && error("invalid infinite node coordinate")
                cmin ≤ c < cmax ||
                    error("node coordinate(s) out of range [-0.5,0.5)")
            end
        end
    end

    # Copy node coordinates (dimensions are in reverse order in NFFT).
    @inbounds begin
        for j in 1:M
            @simd for i in 1:d
                dest[(d+1)-i,j] = x[i,j]
            end
        end
    end

    if ths.ptr != C_NULL
        # Pre-compute PSI if needed.
        if isa(ths, NFFT)
            if (get_flags(ths) & (PRE_LIN_PSI|PRE_FG_PSI|PRE_PSI|PRE_FULL_PSI)) != 0
                _precompute_one_psi(ths.ptr)
            end
        else
            if (get_flags(ths) & PRE_PSI) != 0
                _precompute_psi(ths.ptr)
            end
        end
        ths.ready = true
    else
        ths.ready = false
    end
    return ths
end

# Implement behavior of the NFFT operator.
#
# It is not advisable to provide a single method, say
# `apply!(dst,ths,src,adj::Bool=false)` to apply either the direct (`adj is
# false) or the adjoint (`adj` is true) because this would break type
# stability.

function apply_direct!(dst, ths::BasicPlans, src)
    isready(ths) || error("the node coordinates must have been set")
    _copyto!(ths.f_hat, src)
    _apply_direct(ths.ptr)
    return _copyto!(dst, ths.f)
end

function apply_adjoint!(dst, ths::BasicPlans, src)
    isready(ths) || error("the node coordinates must have been set")
    _copyto!(ths.f, src)
    _apply_adjoint!(ths.ptr)
    return _copyto!(dst, ths.f_hat)
end

apply_direct(ths::Plan, src) =
    apply_direct!(create_f(ths), ths, src)

apply_adjoint(ths::Plan, src) =
    apply_adjoint!(create_f_hat(ths), ths, src)

create_f(ths::NFFT{T,d,E}) where {T,d,E} =
    Vector{Complex{T}}(undef, size(ths.f))

create_f(ths::Union{NFCT{T,d},NFST{T,d}}) where {T,d} =
    Vector{T}(undef, size(ths.f))

create_f_hat(ths::NFFT{T,d,E}) where {T,d,E} =
    Array{E,d}(undef, size(ths.f_hat))

create_f_hat(ths::Union{NFCT{T,d},NFST{T,d}}) where {T,d} =
    Array{T,d}(undef, size(ths.f_hat))

"""
    unsafe_load_field(ptr, T, off)

yields the value of a field of type `T` at relative offset `off` (in bytes) in
a structure whose base address is `ptr`.  A result of type `T` is returned.  If
`ptr` is NULL, a `LibNFFT.NullPointerError` exception is thrown.

"""
unsafe_load_field(ptr::Ptr, ::Type{T}, off::Integer) where {T} =
    (ptr == C_NULL ? throw(NullPointerError()) :
     unsafe_load(Ptr{T}(ptr + off)))

# Allocate dynamic memory:
#     ptr = _allocate(nbytes)
#     ptr = _allocate(T, num=1)
# then:
#     _free(ptr)
#
function _allocate(size::Integer)
    ptr = Libc.malloc(size)
    ptr != C_NULL || thrown(OutOfMemoryError())
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t),
          ptr, 0, size)
    return ptr
end

_allocate(::Type{T}, n::Integer=1) where {T} =
    Ptr{T}(_allocate(n*sizeof(T)))

_free(ptr::Ptr) =
    ptr != C_NULL && Libc.free(ptr)

"""
    indices(A, ...)

yields an iterable object for visiting each index of its arguments.  The is is
the same as `eachindex` except that arguments must be abstract arrays and that
it is checked that arguments have the same dimensions and have standard
indexing.

"""
@inline function indices(A::AbstractArray)
    Base.has_offset_axes(A) &&
        error("argument(s) have non-standard indexing")
    return eachindex(A)
end
@inline function indices(A::AbstractArray, B::AbstractArray)
    Base.has_offset_axes(A, B) &&
        error("argument(s) have non-standard indexing")
    size(A) == size(B) ||
        dimension_mismatch("arguments have not smae dimensions")
    return eachindex(A, B)
end
@inline function indices(A::AbstractArray, B::AbstractArray, C::AbstractArray)
    Base.has_offset_axes(A, B, C) &&
        error("argument(s) have non-standard indexing")
    size(A) == size(B) == size(C) ||
        dimension_mismatch("arguments have not smae dimensions")
    return eachindex(A, B, C)
end

"""
    _copyto!(dst, src) -> dst

copies the values of the source `src` into the destination `dst`.  An error is
thrown if the source and destination do not have the same dimensions oir have
non-standard indexing.  The types of their elements may be different.  The
`copyto!` method of Julia does not impose that the dimensions be the same and
may be slower than `_copyto!`.  For instance, copying from a complex array to a
real array is about 5 times slower because of the checking of the imaginary
part (here it is jut ignored).

"""
function _copyto!(dst::AbstractArray{T,N}, src::AbstractArray{T,N}) where {T,N}
    @inbounds @simd for i in indices(A, B)
        dst[i] = src[i]
    end
    return dst
end

function _copyto!(dst::Array{<:Real,N}, src::Array{<:Real,N}) where {N}
    @inbounds @simd for i in indices(A, B)
        dst[i] = src[i]
    end
    return dst
end

function _copyto!(dst::Array{<:Complex,N}, src::Array{<:Complex,N}) where {N}
    @inbounds @simd for i in indices(A, B)
        dst[i] = src[i]
    end
    return dst
end

function _copyto!(dst::Array{<:Complex,N}, src::Array{<:Real,N}) where {N}
    @inbounds @simd for i in indices(A, B)
        dst[i] = src[i]
    end
    return dst
end

function _copyto!(dst::Array{<:Real,N}, src::Array{<:Complex,N}) where {N}
    @inbounds @simd for i in indices(A, B)
        dst[i] = src[i].re
    end
    return dst
end

end # module


text:
using Test
using Printf
using FileIO

include("../../src/image.jl")

# the tested code
include("../../src/synthPoints.jl")

include("tinyImage.jl")

@testset "SortedOffsets" begin

    testImage = createTestMaskedImage()

    testSynthPoints = generateOrderedSynthPoints(testImage)

    println(testSynthPoints)
    @debug "Synth points" testSynthPoints

    for point in testSynthPoints
        println(point)
    end

end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
"""
This is a sample input file for the CloudScat2 code.  To use it run julia
and type

julia> include("scan_depths.jl")

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia scan_depths.jl

- - -
2019 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat
using StaticArrays
using Printf
using Logging
using Dates
using Formatting

const co = CloudScat.constants

function runone(lmbd, radius_um, h_km)
    params = Params(
        # NUmber of simulated photons
        N = 5e6,

        # The wavelength of the scattered light
        λ = lmbd * co.nano,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = @SVector([0, 0, h_km * co.kilo]),
        source_b = @SVector([0, 0, h_km * co.kilo]),

        zmin = 13 * co.kilo)

    # World geometry
    # Describe the cloud geometry here.  It may consist in combinations of
    # geometrical figures (cylinders, spheres, cones, semi-planes)
    cloud = Cylinder(7 * co.kilo, 15 * co.kilo, 0, 0, 10 * co.kilo)

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = Cylinder(7 * co.kilo, 60 * co.kilo, 0, 0, 100 * co.kilo)

    # Set a homogeneous cloud droplet density and radius.
    # See variable_droplet_radius.jl for how to deal with inhomogeneous clouds.
    # Set a homogeneous cloud droplet density and radius.
    radius = 10e-6

    # Density of the scattering centers.
    nscat = 100 * co.centi^-3

    composition = Homogeneous(params.λ, nscat, radius)

    # Define the full simulation world
    world = World(cloud, domain, composition)

    # Observers record a time curve and an image of received light.  You can
    # set up as many observers as you wish but they have a significant impact
    # on performance.
    observers = [Observer(
        # Location of the observer
        position = [0, 0, 400 * co.kilo],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)]

    basename = splitext(splitdir(@__FILE__)[2])[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    CloudScat.main(params, world, observers,
                    saveto=@sprintf("%s_%d_%02dum_%02dkm.h5",
                                    basename, lmbd, radius_um, h_km))
end

function run()
    for lmbd in [337, 777]
        for radius_um in [10, 20]
            for h_km in 8:14
                runone(lmbd, radius_um, h_km)
            end
        end
    end
end

isinteractive() || run()


text:
"""
Based on
https://github.com/bjkomer/hyperopt-tutorial/blob/master/Function-Fitting-Example.ipynb
"""
module TestFunctionFitting

using Compat
import Distributions
using Test

using TreeParzen

# Define possible functions to fit

function sinusoid(domain, amplitude, frequency, offset, phase)
    return amplitude .* sin.(frequency .* domain .+ phase) .+ offset
end

function polynomial(domain, coefficients)
    result = zeros(length(domain))
    for (p, c) in enumerate(coefficients)
        result .+= (c * domain .^ p)
    end
    return result
end

# Define the objective function

function objective_ff(params)
    args = params[:args]

    # Generate set of data points from a sinusoid
    num_points = 20
    x = @compat range(0, 10; length = num_points)
    data = 1.5sin.(x .+ 1) .+ 2

    # Add noise to the data
    data += rand(Distributions.Normal(0, 0.2), num_points)
    estimate = if args[:type] == :sinusoid
        sinusoid(x, args[:amplitude], args[:frequency], args[:offset], args[:phase])
    elseif args[:type] == :polynomial
        polynomial(x, args[:coefficients])
    end
    # sum of squared error
    return (data - estimate) .^ 2 |> sum
end

# Define the search space with sinusoid and polynomial
space = Dict(
    :args => HP.Choice(:function, [
        Dict(
            :type => :sinusoid,
            :amplitude => HP.Uniform(:amplitude, 0.1, 2.0),
            :frequency => HP.Uniform(:frequency, 0.1, 2.0),
            :offset => HP.Normal(:offset, 0.0, 2.0),
            :phase => HP.Normal(:phase, 0.0, 2.0),
        ),
        Dict(
            :type => :polynomial,
            :coefficients => HP.Choice(:degree, [
                [HP.Normal(:d0_c0, 0.0, 5.0)],
                [HP.Normal(:d1_c0, 0.0, 5.0), HP.Normal(:d1_c1, 0.0, 5.0)],
                [
                    HP.Normal(:d2_c0, 0.0, 5.0), HP.Normal(:d2_c1, 0.0, 5.0),
                    HP.Normal(:d2_c2, 0.0, 5.0),
                ],
                [
                    HP.Normal(:d3_c0, 0.0, 5.0), HP.Normal(:d3_c1, 0.0, 5.0),
                    HP.Normal(:d3_c2, 0.0, 5.0), HP.Normal(:d3_c3, 0.0, 5.0),
                ],
            ]),
        ),
    ]),
)

# Run the search for the specified number of evaluations
best = fmin(objective_ff, space, 1_000)

@test isapprox(best[:args][:amplitude], 1.4762631349533661; rtol=1e1)
@test isapprox(best[:args][:frequency], 0.98526529964042; rtol=1e1)
@test isapprox(best[:args][:offset], 2.010532261653186; rtol=1e1)
@test isapprox(best[:args][:phase], 0.8694539754903817; rtol=1e1)

# Define the search space with polynomial alone
space = Dict(
    :args => HP.Choice(:function, [
        Dict(
            :type => :polynomial,
            :coefficients => HP.Choice(:degree, [
                [HP.Normal(:d0_c0, 0.0, 5.0)],
                [HP.Normal(:d1_c0, 0.0, 5.0), HP.Normal(:d1_c1, 0.0, 5.0)],
                [
                    HP.Normal(:d2_c0, 0.0, 5.0), HP.Normal(:d2_c1, 0.0, 5.0),
                    HP.Normal(:d2_c2, 0.0, 5.0),
                ],
                [
                    HP.Normal(:d3_c0, 0.0, 5.0), HP.Normal(:d3_c1, 0.0, 5.0),
                    HP.Normal(:d3_c2, 0.0, 5.0), HP.Normal(:d3_c3, 0.0, 5.0),
                ],
            ]),
        ),
    ]),
)

best = fmin(objective_ff, space, 1_000)
@test best[:args][:type] == :polynomial
@test isa(best[:args][:coefficients], Vector{Float64})

space = Dict(
    :args => HP.PChoice(:function, [
        HP.Prob(0.7, Dict(
            :type => :sinusoid,
            :amplitude => HP.Uniform(:amplitude, 0.1, 2.0),
            :frequency => HP.Uniform(:frequency, 0.1, 2.0),
            :offset => HP.Normal(:offset, 0.0, 2.0),
            :phase => HP.Normal(:phase, 0.0, 2.0),
        )),
        HP.Prob(0.3, Dict(
            :type => :polynomial,
            :coefficients => HP.Choice(:degree, [
                [HP.Normal(:d0_c0, 0.0, 5.0)],
                [HP.Normal(:d1_c0, 0.0, 5.0), HP.Normal(:d1_c1, 0.0, 5.0)],
                [
                    HP.Normal(:d2_c0, 0.0, 5.0), HP.Normal(:d2_c1, 0.0, 5.0),
                    HP.Normal(:d2_c2, 0.0, 5.0),
                ],
                [
                    HP.Normal(:d3_c0, 0.0, 5.0), HP.Normal(:d3_c1, 0.0, 5.0),
                    HP.Normal(:d3_c2, 0.0, 5.0), HP.Normal(:d3_c3, 0.0, 5.0),
                ],
            ]),
        )),
    ]),
)

# Test that evaluation completes
fmin(objective_ff, space, 1_000)

end # module TestFunctionFitting
true


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "primecount"
version = v"6.5.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/kimwalisch/primecount.git", "fe433f39fa16eaeab41d7caaef90a9caefa32474"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/primecount/
atomic_patch -p1 ../patches/0001-Allow-disabling-building-static-library-for-Windows.patch
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    -DBUILD_STATIC_LIBS=OFF \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# code contains std::string values, this is needed for GCC 4/5 version boundary compat
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("primecount", :primecount),
    LibraryProduct("libprimecount", :libprimecount),
    LibraryProduct("libprimesieve", :libprimesieve),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
### Decimation for dynamic-variables objects
function decimate!(vars::DynamicsVars)
    decimate!(vars.F)
    decimate!(vars.Fˢ)
    decimate!(vars.ζ)
end

### Static auxiliar quantities
function dynamics_weights!(w, ws, K, S, d, grid)
    w .= ws .* K.^(d + 1) .* (1 .- inv.(S)).^2
    return w
end
#
function dynamics_weights!(w::Vector{U}, ws, K, S, d, grid) where {T, U <: DProjections{0, T}}

    @inbounds for i in eachindex(K)
        t = ws[i] * K[i]^4 * (1 - inv(S[i].t))^2
        w[i] = MDProjections(t)
    end

    return w
end
#
function dynamics_weights!(w::Vector{U}, ws, K, S, d, grid) where {U <: MDProjections{2}}

    @inbounds for i in eachindex(K)
        Sᵢ = S[i]
        k² = K[i]^2
        wᵣ = ws[i] * k²
        wₜ = wᵣ * k²

        # Projections components of the weights
        t  = wₜ * (1 - inv(Sᵢ.t))^2
        r₀ = 3wₜ * (1 - inv(Sᵢ.r[1]))^2
        r₁ = 6wᵣ * ((Sᵢ.r[1] - 1) * inv(Sᵢ.r[2]))^2

        w[i] = MDProjections(t, SVector(r₀, r₁))
    end

    return w
end
#
#function dynamics_weights!(w::Vector{U}, ws, K, S, d, grid) where {U <: MDProjections}
#
#    @inbounds for i in eachindex(K)
#        Sᵢ = S[i]
#        k² = K[i]^2
#        wᵣ = ws[i] * k²
#
#        # Projections components of the weights
#        t = wᵣ * k² * (1 - inv(Sᵢ.t))^2
#        for j = 1:l
#            α = wᵣ * (2j + 1)
#            v[j] = α * k² * (1 - inv(Sᵢ.r[j]))^2
#            v[j] = α * j * (j + 1) * ((Sᵢ.r[j] - 1) * inv(Sᵢ.r[j + l]))^2
#        end
#
#        w[i] = MDProjections(t, SVector(v))
#    end
#
#    return w
#end


text:
function estimate_error(obj::GraggBulirschStoerStateInterpolator, arg0::Vector{jdouble})
    return jcall(obj, "estimateError", jdouble, (Vector{jdouble},), arg0)
end


text:
# Julia script to compute eigenvector of Laplacian

# Load indexes and construct sparse matrix
iis = readcsv(open(ARGS[1]), Int32)
jjs = readcsv(open(ARGS[2]), Int32)
vvs = readcsv(open(ARGS[3]), Int32)
A = sparse(iis[:]+1, jjs[:]+1, vvs[:])

# Compute Laplacian
Lap = spdiagm(sum(A,2)[:,1], 0, size(A)[1], size(A)[2]) - A + 1e-9*speye(Float64, size(A)[1]);

# Get largest eigenvalue
(lbdamax, ~) = eigs(Lap, nev=1);

# Construct lambda_max*I - Laplacian, its second largest eigenvalue is the fiedler value we are looking for
negLap = lbdamax[1]*speye(Float64, size(Lap)[1]);
negLap = negLap - Lap;

# Get eigenvalue and eigenvector
(lbdadiff, V2) = eigs(negLap, nev=6, ncv=20, tol=1.e-9, v0=ones((size(Lap)[1],)), maxiter=10000000);
fidval = lbdamax[1] - lbdadiff[2]
fidvec = V2[:,2];

# Sort to recover permutation (spectral ordering)
permu = sortperm(fidvec)

writecsv(ARGS[4], permu')

max_stars_repo_name: antrec/spectrassembler

max_issues_repo_name: antrec/spectrassembler

max_forks_repo_name: antrec/spectrassembler


text:
module HierarchicalStateMachineModule

using UserExceptionModule: UserError

abstract type State end

abstract type Event end

function handle_event(state::State,event::Event)
  throw(UserError())
end

mutable struct HierarchicalStateMachine
  state::State
end

function handle_event(hsm::HierarchicalStateMachine, event::Event)
  hsm.state = handle_event(hsm.state,event)
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
x = 1
x += 2
x += 3
x += 4
x += 5

function f(x)
  y = x + 1
  return 2 * y
end

x += 1


text:
using SigmaRidgeRegression
using Test
using LinearAlgebra
using Random
import StatsBase

@testset "Woodbury and Cholesky" begin
    Random.seed!(1)
    σ = 5.0
    grp = GroupedFeatures([30;30;30;30;30;30])
    n = 400
    p = grp.p
    X = randn(n, p)
    αs = sqrt.(range(4.0,12.0,length=ngroups(grp)))#r#ange(2.0, 2.75, 3.5; length=3)#1.0:5.0
    β = random_betas(grp, αs)
    Y = X*β .+ σ .* randn(n)

    tmp_chol = BasicGroupRidgeWorkspace(X=X, Y=Y, groups=grp)

    tmp_woodb = BasicGroupRidgeWorkspace(X=X, Y=Y, groups=grp,
                               XtXpΛ_chol = WoodburyRidgePredictor(X))

    @test SigmaRidgeRegression.λωλας_λ(tmp_chol) ≈ SigmaRidgeRegression.λωλας_λ(tmp_woodb)
    beta_chol = tmp_chol.XtXpΛ_chol \ tmp_woodb.XtY
    beta_wdb =  tmp_woodb.XtXpΛ_chol \ tmp_woodb.XtY
    @test beta_chol ≈ beta_wdb

    tmp_chol.XtXpΛ_chol \ tmp_chol.X' ≈ tmp_woodb.XtXpΛ_chol \ tmp_woodb.X'

    ldiv_chol = ldiv!(tmp_chol.XtXpΛ_div_Xt, tmp_chol.XtXpΛ_chol, tmp_chol.X')
    ldiv_wdb = ldiv!(tmp_woodb.XtXpΛ_div_Xt, tmp_woodb.XtXpΛ_chol, tmp_woodb.X')

    @test ldiv_chol ≈ ldiv_wdb

    loo_chol = StatsBase.fit!(tmp_chol, SigmaRidgeRegression.λωλας_λ(tmp_chol))
    loo_wdb = StatsBase.fit!(tmp_woodb, SigmaRidgeRegression.λωλας_λ(tmp_chol))

    @test tmp_woodb.β_curr ≈ tmp_chol.β_curr
end

@testset "MMI" begin
    include("test_mmi.jl")
end
#include("test_grouplasso.jl")
include("test_covariance.jl")
include("test_simulations.jl")

@testset "Theoretical risks" begin
    include("test_theoretical_risks.jl")
end
# TO ADD
# Compare: - GroupLasso AGAINST R
#          - Check group lasso lambda_max
#          - test the TunedRidgeRegressor type
#          - test σ_max  (probably won't increase coverage but...)


text:
module Samples

import Historic

Historic.@define Default
Default.isrecording() = true

Historic.@define Custom
Custom.isrecording() = true
Custom.defaultdata() = (; foo = :bar)
Custom.prefix() = "custom prefix"

default_simple() = Default.@record(:simple)
default_withdata(v) = Default.@record(:withdata, k = v)
default_noyield() = Default.@record(:noyield, {yield = false})

custom_simple() = Custom.@record(:simple)
custom_withdata(v) = Custom.@record(:withdata, k = v)

function default_throw()
    try
        error("some exception")
    catch err
        Default.@record(:catch, exception = (err, catch_backtrace()))
    end
end

end  # module


text:
using Revise
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters

using Profile
using CountPlus
using Juno
using BenchmarkTools

using ShaleDrillingLikelihood: test_serial_simloglik!

num_i = 15
M = 60
nz = 101
nψ = 101

data_theta_sm, data_theta_lg = MakeTestData(; num_i=num_i, nz=nz, nψ=nψ)

stopcount!()

data, theta = data_theta_lg

theta_peturb = 0.9 .* theta
leo = LocalEstObj(data,theta_peturb)
s = SimulationDraws(M, ShaleDrillingLikelihood.data(leo))
reo = RemoteEstObj(leo, M)
ew = EstimationWrapper(leo, reo)

# run a couple of times to compile
VFTOL = 1e-9
for i = 1:2
    test_serial_simloglik!(ew, theta_peturb, true; vftol=VFTOL)
    test_serial_simloglik!(ew, theta_peturb, true; vftol=VFTOL)
end

# profile
Profile.clear()
@profile test_serial_simloglik!(ew, theta_peturb, true; vftol=VFTOL)
Juno.profiletree()
Juno.profiler()

# different VF Tols
for VFTOL in (1e-7,) # 1e-6, 1e-7, 1e-8, 1e-9, 1e-10)
    println("VFTOL = $VFTOL")
    @btime test_serial_simloglik!($ew, $theta_peturb, true; vftol=$VFTOL)
end


text:
using Test
using BranchAndBound
using JuMP, GLPK
using MathOptInterface

const BB = BranchAndBound
const MOI = MathOptInterface

@testset "Abstract Node" begin
    struct MyNode <: BB.AbstractNode end
    node = MyNode()
    @test BB.branch(node) == []
    BB.bound!(node)
    BB.heuristics!(node)
    BB.apply_changes!(node)
end

model = Model(GLPK.Optimizer)
@variable(model, 0 <= x <= 1)
@objective(model, Min, x)

tree = BB.AbstractTree()
node = BB.JuMPNode{BB.VariableBranch}(model)

@testset "Initial Tree" begin
    @test tree.node_counter == 0
    @test length(tree.nodes) == 0
    @test length(tree.processed) == 0
end

@testset "Initial Node" begin
    @test node.id == -1
    @test isnothing(node.parent)
    @test node.depth == 0
    @test isnothing(node.branch)
    @test node.solution_status == MOI.OPTIMIZE_NOT_CALLED
    @test isinf(-node.bound)
    @test !isnothing(node.model)
end

@testset "Tree with root" begin
    BB.push!(tree, node)
    @test tree.node_counter == 1
    @test length(tree.nodes) == 1
    @test length(tree.processed) == 0
    @test tree.nodes[1] == node

    tree1 = BB.initialize_tree(model)
    @test tree1.node_counter == 1
    @test length(tree1.nodes) == 1
    @test length(tree1.processed) == 0
end

current_node = BB.next_node(tree)

@testset "Root node with model" begin
    @test !isnothing(current_node.model)
    @test JuMP.lower_bound(x) == 0.
    @test JuMP.upper_bound(x) == 1.
end

@testset "Bounding root node" begin
    BB.bound!(current_node)
    @test current_node.solution_status == MOI.OPTIMAL
    @test current_node.bound == 0.0
    @test JuMP.value(x) == 0.0
end

@testset "Marked root node as processed" begin
    BB.processed!(tree, current_node)
    @test BB.termination(tree)
    @test length(tree.processed) == 1
    @test tree.processed[1] == current_node
end

@testset "Empty branching" begin
    BB.branch!(tree, current_node)
    @test tree.node_counter == 1
    @test length(tree.nodes) == 0
    @test length(tree.processed) == 1
end

@testset "Branching at root" begin
    child1 = BB.create_child_node_with_lb(current_node, x, 0.5)
    child2 = BB.create_child_node_with_ub(current_node, x, 0.5)
    BB.push!(tree, [child1,child2])
    @test tree.node_counter == 3
    @test length(tree.nodes) == 2
    @test length(tree.processed) == 1
end

next_node = BB.next_node(tree)

@testset "Child 1 with model" begin
    @test !isnothing(next_node.model)
    print(next_node.model)
    @test JuMP.lower_bound(x) == 0.
    @test JuMP.upper_bound(x) == 0.5
    @test next_node.bound == 0.0
end

@testset "Bounding Child 1 node" begin
    BB.bound!(next_node)
    @test next_node.solution_status == MOI.OPTIMAL
    @test next_node.bound == 0.0
    @test JuMP.value(x) == 0.0
end


text:
using JuMP, GLPK

function get_α(P::Matrix{Float64}, e::Vector{Float64}, f::Vector{Float64})
    # Set up model and payoff
    model = direct_model(GLPK.Optimizer())
    @variable(model, z)
    @objective(model, Max, 1.0 * z)

    @variable(model, x[1:length(e)], lower_bound = 0.0)
    @constraint(model, sum(x) == 1.0)
    @constraint(model, x' * P .>= f)
    @constraint(model, x' * e .>= z)

    optimize!(model)
    return JuMP.value(z)
end

function get_β(O::Matrix{Float64}, e::Vector{Float64}, f::Vector{Float64})
    # Set up model and payoff
    model = direct_model(GLPK.Optimizer())
    @variable(model, z)
    @objective(model, Min, 1.0 * z)

    @variable(model, x[1:length(e)], lower_bound = 0.0)
    @constraint(model, sum(x) == 1.0)
    @constraint(model, O * x .<= f)
    @constraint(model, e * x' .<= z)

    optimize!(model)
    return JuMP.objective_value(model)
end

function SMAB(state::DynamicState, static_state::StaticState, α₀::Float64,
  β₀::Float64, depth::Int64; allow_waiting::Bool = false, sim_to_end::Bool = false, ϵ = 0.001)
    A, B = get_simultaneous_decisions(state, static_state, allow_waiting = allow_waiting)
    (length(A) == 0 || depth == 0) &&
        return sim_to_end ? NashResult(sum(battle_scores(state, static_state, 100) / 100) - 0.5,
        no_strat, no_strat) : NashResult(battle_score(state, static_state) - 0.5,
        no_strat, no_strat)

    Q = [play_turn(state, static_state, (a, b)) for a in A, b in B]
    P = vcat(map(q -> get_min_score(q, static_state), Q), repeat([α₀], 1, length(B)))
    O = hcat(map(q -> get_max_score(q, static_state), Q), repeat([β₀], length(A), 1))
    non_dominated_rows = trues(length(A) + 1)
    non_dominated_columns = trues(length(B))
    non_dominated_rows[end] = false

    for i in 1:length(A), j in 1:length(B)
        if non_dominated_rows[i] && non_dominated_columns[j]
            @inbounds if (5 in B || 7 in B) && !(5 <= B[j] <= 8) && iseven(A[i]) &&
              (5 in A || 7 in A) && !(5 <= A[i] <= 8) && iseven(B[j])
                @inbounds P[i, j] = P[i - 1, j - 1]
                @inbounds O[i, j] = O[i - 1, j - 1]
            elseif @inbounds (5 in B || 7 in B) && !(5 <= B[j] <= 8) && iseven(A[i])
                @inbounds P[i, j] = P[i - 1, j]
                @inbounds O[i, j] = O[i - 1, j]
            elseif @inbounds (5 in A || 7 in A) && !(5 <= A[i] <= 8) && iseven(B[j])
                @inbounds P[i, j] = P[i, j - 1]
                @inbounds O[i, j] = O[i, j - 1]
            else
                α = get_α(P[1:end .!= i, 1:end .!= j], P[1:end .!= i, j], O[i, 1:end-1 .!= j])
                β = get_β(O[1:end .!= i, 1:end .!= j], O[i, 1:end .!= j], P[1:end-1 .!= i, j])
                if α < β
                    @inbounds v = SMAB(Q[i, j], static_state, α, β, depth - 1,
                        allow_waiting = allow_waiting, sim_to_end = sim_to_end).payoff
                    if v <= α
                        @inbounds non_dominated_rows[i] = false
                    elseif v >= β
                        @inbounds non_dominated_columns[j] = false
                    else
                        P[i, j] = O[i, j] = v
                    end
                else
                    @inbounds v = SMAB(Q[i, j], static_state, α, α + ϵ, depth - 1,
                        allow_waiting = allow_waiting, sim_to_end = sim_to_end).payoff
                    if v <= α
                        @inbounds non_dominated_rows[i] = false
                    else
                        @inbounds non_dominated_columns[j] = false
                    end
                end
            end
        end
    end
    return nash(P[non_dominated_rows, non_dominated_columns])
end


text:
# # Sensitivity Analysis of Ridge Regression

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/sensitivity-analysis-ridge.jl)

# This example illustrates the sensitivity analysis of data points in a
# [Ridge Regression](https://en.wikipedia.org/wiki/Ridge_regression) problem.
# The general form of the problem is given below:

# ```math
# \begin{split}
# \begin{array} {ll}
# \mbox{minimize} & \sum_{i=1}^{N} (y_{i} - w x_{i} - b)^2 + \alpha (w^2 + b^2) \\
# \end{array}
# \end{split}
# ```
# where
# - `w`, `b` are slope and intercept of the regressing line
# - `x`, `y` are the N data points
# - `α` is the regularization constant
#
# which is equivalent to:
# ```math
# \begin{split}
# \begin{array} {ll}
# \mbox{minimize} & e^{\top}e + \alpha (w^2) \\
# \mbox{s.t.} & e_{i} = y_{i} - w x_{i} - b \quad \quad i=1..N  \\
# \end{array}
# \end{split}
# ```

# This tutorial uses the following packages

using JuMP
import DiffOpt
import Random
import Ipopt
import Plots
using LinearAlgebra: dot

# ## Define and solve the problem

# Construct a set of noisy (guassian) data points around a line.

Random.seed!(42)

N = 150

w = 2 * abs(randn())
b = rand()
X = randn(N)
Y = w * X .+ b + 0.8 * randn(N);

# The helper method `fit_ridge` defines and solves the corresponding model.
# The ridge regression is modeled with quadratic programming
# (quadratic objective and linear constraints) and solved in generic methods
# of Ipopt. This is not the standard way of solving the ridge regression problem
# this is done here for didactic purposes.

function fit_ridge(X, Y, alpha = 0.1)
    N = length(Y)
    ## Initialize a JuMP Model with Ipopt solver
    model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
    set_silent(model)
    @variable(model, w) # angular coefficient
    @variable(model, b) # linear coefficient
    ## expression defining approximation error
    @expression(model, e[i=1:N], Y[i] - w * X[i] - b)
    ## objective minimizing squared error and ridge penalty
    @objective(
        model,
        Min,
        1 / N * dot(e, e) + alpha * (w^2),
    )
    optimize!(model)
    return model, w, b # return model & variables
end

# Plot the data points and the fitted line for different alpha values

p = Plots.scatter(X, Y, label=nothing, legend=:topleft)
mi, ma = minimum(X), maximum(X)
Plots.title!("Fitted lines and points")

for alpha in 0.5:0.5:1.5
    local model, w, b = fit_ridge(X, Y, alpha)
    ŵ = value(w)
    b̂ = value(b)
    Plots.plot!(p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], label="alpha=$alpha", width=2)
end
p

# ## Differentiate

# Now that we've solved the problem, we can compute the sensitivity of optimal
# values of the slope `w` with
# respect to perturbations in the data points (`x`,`y`).

alpha = 0.4
model, w, b = fit_ridge(X, Y, alpha)
ŵ = value(w)
b̂ = value(b)

# We first compute sensitivity of the slope with respect to a perturbation of the independent
# variable `x`.

# Recalling that the points $(x_i, y_i)$ appear in the objective function as:
# `(yi - b - w*xi)^2`, the `DiffOpt.ForwardInObjective` attribute must be set accordingly,
# with the terms multiplying the parameter in the objective.
# When considering the perturbation of a parameter θ, `DiffOpt.ForwardInObjective()` takes in the expression in the
# objective that multiplies θ.
# If θ appears with a quadratic and a linear form: `θ^2 a x + θ b y`, then the expression to pass to
# `ForwardInObjective` is `2θ a x + b y`.

# Sensitivity with respect to x and y

∇y = zero(X)
∇x = zero(X)
for i in 1:N
    MOI.set(
        model,
        DiffOpt.ForwardInObjective(),
        2w^2 * X[i] + 2b * w - 2 * w * Y[i]
    )
    DiffOpt.forward(model)
    ∇x[i] = MOI.get(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        w
    )
    MOI.set(
        model,
        DiffOpt.ForwardInObjective(),
        (2Y[i] - 2b - 2w * X[i]),
    )
    DiffOpt.forward(model)
    ∇y[i] = MOI.get(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        w
    )
end

# Visualize point sensitivities with respect to regression points.

p = Plots.scatter(
    X, Y,
    color = [dw < 0 ? :blue : :red for dw in ∇x],
    markersize = [5 * abs(dw) + 1.2 for dw in ∇x],
    label = ""
)
mi, ma = minimum(X), maximum(X)
Plots.plot!(p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], color = :blue, label = "")
Plots.title!("Regression slope sensitivity with respect to x")

#

p = Plots.scatter(
    X, Y,
    color = [dw < 0 ? :blue : :red for dw in ∇y],
    markersize = [5 * abs(dw) + 1.2 for dw in ∇y],
    label = ""
)
mi, ma = minimum(X), maximum(X)
Plots.plot!(p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], color = :blue, label = "")
Plots.title!("Regression slope sensitivity with respect to y")

# Note the points with less central `x` values induce a greater y sensitivity of the slope.


text:
module PlotlyJSTest
using Test

using PlotlyJS
const M = PlotlyJS

# using Blink
# !Blink.AtomShell.isinstalled() && Blink.AtomShell.install()

# include("blink.jl")
include("kaleido.jl")

end


text:
unpackconstruct(Construct::Ordered_Construct) = map((x->getfield(x, :m)|>String |> (x -> parse(Int64, x))),Construct)

quad(x,Q,λ) = transpose(x)*Q*x - λ*norm(x,1)
# ref https://github.com/baptistar/BOCS/blob/master/test_problems/Quadratic/quad_mat.m
function compute_decay(rng,n_vars,Lc)
	K = (s,t) -> exp(-1*(s-t)^2/Lc)
	decay = Array{Float64,2}(undef,n_vars,n_vars)

	for i=1:n_vars
		for j=1:n_vars
			decay[i,j] = K(i,j);
		end
	end

	#Generate random quadratic model
	#and apply exponential decay to Q
	Q = randn(rng,n_vars, n_vars);
	Qa = Q.*decay;
end


text:
# Create example data
numbers = [10, 15, 35, 20]
# Plot bars with labels on top
barplot(numbers)
annotations(1:4, numbers .+1 , string.(numbers), halign="center")


text:
MOIU.@model(
    DualizableModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
        MOI.PositiveSemidefiniteConeTriangle,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,)
)

"""
    PrimalDualMap{T}

Maps information from all structures of the primal to the dual model.

* `constrained_var_idx::Dict{VI,Tuple{CI,Int}}`: maps original primal
constrained variables to their primal original constraints (the special ones
that makes them constrained variables) and their internal index (if vector
constraints, VectorOfVariables-in-Set), 1 otherwise (VariableIndex-in-Set).

* `constrained_var_dual::Dict{CI,CI}`: maps the original primal constraint index
of constrained variables to the dual model's constraint index of the
associated dual constraint.

* `constrained_var_zero::Dict{CI,Unions{SAF,VAF}}`: caches scalar affine
functions or vector affine functions associated with constrained variables
of type `VectorOfVariables`-in-`Zeros` or
`VariableIndex`-in-`EqualTo(zero(T))` as
their duals would be `func`-in-`Reals`, which are "irrelevant" to the model.
This information is cached for completeness of the `DualOptimizer` for
`get`ting `ConstraintDuals`.

* `primal_var_dual_con::Dict{VI,CI}`: maps "free" primal variables to their
associated dual constraints. Free variables as opposed to constrained
variables. Note that Dualization will select automatically which variables
are free and which are constrained.

* `primal_con_dual_var::Dict{CI,Vector{VI}}`: maps primal constraint indices to
vectors of dual variable indices. For scalar constraints those vectors will be
single element vectors.

* `primal_con_dual_con::Dict{CI,CI}`: maps primal constraints to
dual variable constraints (if there is such constraint the dual
dual variable is said to be constrained). If the primal constraint's set
is EqualTo or Zeros, no constraint is added in the dual variable (the
dual variable is said to be free).

* `primal_con_constants::Dict{CI,Vector{T}}`: maps primal constraints to
their respective constants, which might be inside the set.
This map is used in `MOI.get(::DualOptimizer,::MOI.ConstraintPrimal,ci)`
that requires extra information in the case that the scalar set constrains
a constant (`EqualtTo`, `GreaterThan`, `LessThan`).

* `primal_parameter::Dict{VI,VI}`: maps parameters in the primal to parameters
in the dual model.

* `primal_var_dual_quad_slack::Dict{VI,VI}`: maps primal variables
(that appear in quadratic objective terms) to dual "slack" variables.

"""
mutable struct PrimalDualMap{T}
    constrained_var_idx::Dict{VI,Tuple{CI,Int}}
    constrained_var_dual::Dict{CI,CI}
    constrained_var_zero::Dict{
        CI,
        Union{MOI.VectorAffineFunction{T},MOI.ScalarAffineFunction{T}},
    }
    primal_var_dual_con::Dict{VI,CI}
    primal_con_dual_var::Dict{CI,Vector{VI}}
    primal_con_dual_con::Dict{CI,CI}
    primal_con_constants::Dict{CI,Vector{T}}

    primal_parameter::Dict{VI,VI}
    primal_var_dual_quad_slack::Dict{VI,VI}

    function PrimalDualMap{T}() where {T}
        return new(
            Dict{VI,Tuple{CI,Int}}(),
            Dict{CI,CI}(),
            Dict{
                CI,
                Union{MOI.VectorAffineFunction{T},MOI.ScalarAffineFunction{T}},
            }(),
            Dict{VI,CI}(),
            Dict{CI,Vector{VI}}(),
            Dict{CI,CI}(),
            Dict{CI,Vector{T}}(),
            Dict{VI,VI}(),
            Dict{VI,VI}(),
        )
    end
end

function is_empty(primal_dual_map::PrimalDualMap{T}) where {T}
    return isempty(primal_dual_map.constrained_var_idx) &&
           isempty(primal_dual_map.constrained_var_dual) &&
           isempty(primal_dual_map.constrained_var_zero) &&
           isempty(primal_dual_map.primal_var_dual_con) &&
           isempty(primal_dual_map.primal_con_dual_var) &&
           isempty(primal_dual_map.primal_con_dual_con) &&
           isempty(primal_dual_map.primal_con_constants) &&
           isempty(primal_dual_map.primal_parameter) &&
           isempty(primal_dual_map.primal_var_dual_quad_slack)
end

function empty!(primal_dual_map::PrimalDualMap)
    Base.empty!(primal_dual_map.constrained_var_idx)
    Base.empty!(primal_dual_map.constrained_var_dual)
    Base.empty!(primal_dual_map.constrained_var_zero)
    Base.empty!(primal_dual_map.primal_var_dual_con)
    Base.empty!(primal_dual_map.primal_con_dual_var)
    Base.empty!(primal_dual_map.primal_con_dual_con)
    Base.empty!(primal_dual_map.primal_con_constants)
    Base.empty!(primal_dual_map.primal_parameter)
    Base.empty!(primal_dual_map.primal_var_dual_quad_slack)
    return
end

"""
    DualProblem{T,OT<:MOI.ModelLike}

Result of the `dualize` function. Contains the fields:

* `dual_model::OT`: contaninf a optimizer or data structure with the
`MathOptInterface` definition of the resulting dual model.

* `primal_dual_map::PrimalDualMap{T}`: a data structure to hold information of
how the primal and dual model are related in terms of indices (`VariableIndex`
and `ConstraintIndex`) and other data.
"""
struct DualProblem{T,OT<:MOI.ModelLike}
    dual_model::OT #It can be a model from an optimizer or a DualizableModel{T}
    primal_dual_map::PrimalDualMap{T}

    function DualProblem{T}(
        dual_optimizer::OT,
        pdmap::PrimalDualMap{T},
    ) where {T,OT<:MOI.ModelLike}
        return new{T,OT}(dual_optimizer, pdmap)
    end
end

function DualProblem{T}(dual_optimizer::OT) where {T,OT<:MOI.ModelLike}
    return DualProblem{T}(dual_optimizer, PrimalDualMap{T}())
end

function DualProblem(dual_optimizer::OT) where {OT<:MOI.ModelLike}
    return DualProblem{Float64}(dual_optimizer)
end

# Empty DualProblem cosntructor
function DualProblem{T}() where {T}
    return DualProblem{T}(DualizableModel{T}(), PrimalDualMap{T}())
end