text:
"""
$(SIGNATURES)

"""
function engineer_caustics(source_image)
    imageBW = Float64.(Gray.(source_image))

    # Set global size parameters
    height, width = size(imageBW)
    global N_Pixel_Height = height
    global N_Pixel_Width = width
    global Meters_Per_Pixel = Caustics_Long_Side / N_Pixel_Height
    println("Image size: $(N_Pixel_Height) x $(N_Pixel_Width)")

    # mesh is the same size as the image with an extra row/column to have coordinates to
    # cover each image corner with a triangle.
    mesh = FaceMesh(N_Pixel_Height, N_Pixel_Width)

    # imageBW is normalised to have 1 unit per pixel on average.
    imageBW /= average(imageBW)

    marginal_change_top = marginal_change_bot = nothing
    max_update_top = marginal_change_top = Inf
    counter = 0
    while (abs(max_update_top+marginal_change_bot) > 1e-6 && counter < 10_000)
        counter += 1

        print(
            """
          ================================================================================================================
          STARTING ITERATION $(counter):
              starting ϕ TOP field = $(field_summary(mesh.corners.ϕ_top))
              starting ϕ BOT field = $(field_summary(mesh.corners.ϕ_bot))

            """,
        )

        ϕ_top = mesh.corners.ϕ_top
        ε_top, max_update_top = solve_velocity_potential!(mesh.corners.ϕ_top, imageBW/2., "it$(counter)_top")
        marginal_change_top = mesh.corners.ϕ_top - ϕ_top

        ϕ_bot = mesh.corners.ϕ_bot
        ε_bot, max_update_bot = solve_velocity_potential!(mesh.corners.ϕ_bot, imageBW/2., "it$(counter)_bot")
        marginal_change_bot = mesh.corners.ϕ_bot - ϕ_bot

        print(
            """

            RESULT AT ITERATION $(counter):
                TOP
                Luminosity error = $(field_summary(ε_top))
                Vertical move max update = $(max_update_top)
                Marginal change = $(field_summary(marginal_change_top))
                end ϕ field = $(field_summary(mesh.corners.ϕ_top))

                BOT
                end ϕ field = $(field_summary(mesh.corners.ϕ_bot))
                Luminosity error = $(field_summary(ε_bot))
                Vertical move max update = $(max_update_bot)
                Marginal change = $(field_summary(marginal_change_bot))
                end ϕ field = $(field_summary(mesh.corners.ϕ_bot))

            ================================================================================================================

          """,
        )
        plot_as_quiver(mesh.corners.ϕ_top, n_steps = 50, scale = height / 10, max_length = height / 20)
        average_absolute(marginal_change_top + marginal_change_bot) < 1e-4 && break
    end

    println("\nSTARTING HORIZONTAL ITERATION ---")

    println("\t Horizontal max update TOP = $(max_update_top) /  BOT = $(max_update_bot)")

    # Move the around a nil average.
    mesh.corners.ϕ_top .-= average(mesh.corners.ϕ_top)
    mesh.corners.ϕ_bot .-= average(mesh.corners.ϕ_bot)

    return mesh, imageBW
end

"""
$(SIGNATURES)

ϕ is the _velocity potential_. The velocity represents the direction towards which the mesh vertices should move.
Zones of high luminosity should be covered with more triangles and should attract more vertices.
"""
function solve_velocity_potential!(mesh, image, prefix)
    # Start with a new clean field
    height, width = size(mesh)

    # Illumination only depends on the position of the corners, not their heights. ϕ is not relevant.
    # Get the area of each individual pixel as stretch/shrunk on the lens. Area = energy.
    # _FENCES_SIZED_
    # Illumination only depends on the position of the corners, not their heights. ϕ is not relevant.
    lens_pixels_area_top, lens_pixels_area_bot = get_lens_pixels_area(mesh)

    # Positive error => the triangle needs to shrink (less light). Enforce nil average error.
    ε_top = zeros(Float64, height, width)
    ε_bot = zeros(Float64, height, width)

    # Each triangle should have the luminosity of half a pixel
    ε_top[1:end-1, 1:end-1] = Float64.(lens_pixels_area_top - image/2.)
    ε_top .-= average(ε_top)

    ε_bot[1:end-1, 1:end-1] = Float64.(lens_pixels_area_bot - image/2.)
    ε_bot .-= average(ε_bot)

    println("""
            solve_velocity_potential! before loop:
                $(field_summary(lens_pixels_area_top, "Pixel area"))
                $(field_summary(ε_top, "luminosity error"))
                $(field_summary(lens_pixels_area_bot, "Pixel area"))
                $(field_summary(ε_bot, "luminosity error"))
                """)

    # Save the loss image as a png.
    save_plot_scalar_field!(ε_top+ε_bot, "error_ε_$(prefix)")
    save_plot_scalar_field!(ε_top, "error_εtop_$(prefix)")
    save_plot_scalar_field!(ε_bot, "error_εbot_$(prefix)")

    # Save the generated image as a png.
    lens_pixels_area = lens_pixels_area_top + lens_pixels_area_bot
    save(
        "./examples/img_noluminosityratio_$(prefix).png",
        Gray.(clamp.(lens_pixels_area, 0.0, 1.0)),
    )

    luminosity_ratio = sum(image) / sum(lens_pixels_area)
    save(
        "./examples/img_$(prefix).png",
        Gray.(clamp.(lens_pixels_area * luminosity_ratio, 0.0, 1.0)),
    )

    luminosity_ratio = sum(image/2.) / sum(lens_pixels_area_top)
    save(
        "./examples/img_$(prefix)_top.png",
        Gray.(clamp.(lens_pixels_area_top * luminosity_ratio, 0.0, 1.0)),
    )

    luminosity_ratio = sum(image/2.) / sum(lens_pixels_area_bot)
    save(
        "./examples/img_$(prefix)_bot.png",
        Gray.(clamp.(lens_pixels_area * lens_pixels_area_bot, 0.0, 1.0)),
    )

    # Start with a clean, flat potential field.
    mesh.corners.ϕ_top = zeros(Float64, height, width)
    mesh.corners.ϕ_bot = zeros(Float64, height, width)

    ϕ_b4_top = copy(mesh.corners.ϕ_top)
    ϕ_b4_bot = copy(mesh.corners.ϕ_bot)

    count = 0
    min_update_top = min_update_bot = new_update_top = new_update_bot = 10_000
    old_update_top = 2 * new_update_top
    old_update_bot = 2 * new_update_bot
    while (
        1e-5 < new_update_top + new_update_bot < 1.5 * min_update_top + 1.5 * min_update_bot &&
        1e-4 < (old_update_top - new_update_top) / old_update_top &&
        1e-4 < (old_update_bot - new_update_bot) / old_update_bot &&
        count < 10_000
    )
        count += 1
        old_update_top = new_update_top
        old_update_bot = new_update_bot
        min_update_top = min(min_update_top, new_update_top)
        min_update_bot = min(min_update_bot, new_update_bot)

        new_update_top, ∇²ϕ_est_top, δ_top = propagate_poisson!(mesh.corners.ϕ_top,ε_top)
        new_update_bot, ∇²ϕ_est_bot, δ_bot = propagate_poisson!(mesh.corners.ϕ_bot,ε_bot)

        count % 1_000 == 1 && println("""
                                        TOP
                                        Iteration $(count) - max_update = $(round(new_update_top, sigdigits=4))
                                            $(field_summary(mesh.corners.ϕ_top, "ϕ"))
                                            $(field_summary(∇²ϕ_est_top, "∇²ϕ_est"))
                                            $(field_summary(δ_top, "δ"))
                                        BOTTOM
                                        Iteration $(count) - max_update = $(round(new_update_bot, sigdigits=4))
                                            $(field_summary(mesh.corners.ϕ_bot, "ϕ"))
                                            $(field_summary(∇²ϕ_est_bot, "∇²ϕ_est"))
                                            $(field_summary(δ_bot, "δ"))
                                            """)

        # Save the loss image as a png.
        if count % 1_000 == 1
            save_plot_scalar_field!(∇²ϕ_est_top, "∇²ϕ_est_$(prefix)")
            save_plot_scalar_field!(δ_top, "delta_$(prefix)")
            save_plot_scalar_field!(
                mesh.corners.ϕ_bot + mesh.corners.ϕ_top -ϕ_b4_top -ϕ_b4_bot ,
                "change_mesh.corners.ϕ_$(prefix)",
            )

            save_plot_scalar_field!(∇²ϕ_est_top, "∇²ϕ_est_$(prefix)_top")
            save_plot_scalar_field!(δ_top, "delta_$(prefix)_top")
            save_plot_scalar_field!(
                mesh.corners.ϕ_top - ϕ_b4_top,
                "change_mesh.corners.ϕ_$(prefix)_top",
            )

            save_plot_scalar_field!(∇²ϕ_est_bot, "∇²ϕ_est_$(prefix)_bot")
            save_plot_scalar_field!(δ_bot, "delta_$(prefix)_bot")
            save_plot_scalar_field!(
                mesh.corners.ϕ_bot - ϕ_b4_bot,
                "change_mesh.corners.ϕ_$(prefix)_bot",
            )

            ϕ_b4_top = copy(mesh.corners.ϕ_top)
            ϕ_b4_bot = copy(mesh.corners.ϕ_bot)
        end
    end

    # Now we need to march the mesh row,col corner locations according to this gradient.
    δ_top = march_mesh!(mesh, :top)
    δ_bot = march_mesh!(mesh, :bottom)
    return ε, new_update
end

"""
$(SIGNATURES)

ϕ is the height map.

`march_mesh` flexes the mesh
"""
function march_mesh(mesh, mesh_ϕ)

    # Calculate the gradient of the velocity potential and reverse its direction.
    vr, vc = ∇(mesh_ϕ)
    mesh.corners.vr .= -vr
    mesh.corners.vc .= -vc

    # Clean up the mesh borders
    # reset_border_values!(mesh.corners)

    # For each point in the mesh we need to figure out its velocity
    # However all the nodes located at a border will never move
    # I.e. velocity (Vx, Vy) = (0, 0) and the square of acrylate will remain of the same size.

    # Get the time, at that velocity, for the area of the triangle to be nil.
    # We are only interested by positive values to only move in the direction of the (opposite) gradient
    height, width = size(mesh)
    list_triangles = vcat(
        [triangle3D(mesh, row, col, :top) for row ∈ 1:height, col ∈ 1:width],
        [triangle3D(mesh, row, col, :bottom) for row ∈ 1:height, col ∈ 1:width],
    )

    list_maximum_t = [t for t ∈ find_maximum_t.(list_triangles) if !isnothing(t) && t > 0.0]

    δ = minimum(list_maximum_t) / 1.5
    new_r = mesh.corners.r - δ * vr
    new_c = mesh.corners.c - δ * vc

    return new_r, new_c, δ
end

function march_mesh!(mesh::FaceMesh, side::Union{:top, :bottom})

    # Calculate the gradient of the velocity potential and reverse its direction.
    if side == :top
        mesh.corners.vr, mesh.corners.vc = ∇(mesh.corners.ϕ_top)
    elseif side == :bottom
        mesh.corners.vr, mesh.corners.vc = ∇(mesh.corners.ϕ_bot)
    end
    mesh.corners.vr .*= -1.0
    mesh.corners.vc .*= -1.0

    # Clean up the mesh borders
    reset_border_values!(mesh.corners)

    # For each point in the mesh we need to figure out its velocity
    # However all the nodes located at a border will never move
    # I.e. velocity (Vx, Vy) = (0, 0) and the square of acrylate will remain of the same size.
    mesh_r = copy(mesh.corners.r)
    mesh_c = copy(mesh.corners.c)

    # Get the time, at that velocity, for the area of the triangle to be nil.
    # We are only interested by positive values to only move in the direction of the (opposite) gradient
    height, width = size(mesh)
    list_triangles = [triangle3D(mesh, row, col, side) for row ∈ 1:height, col ∈ 1:width]
    list_maximum_t = [t for t ∈ find_maximum_t.(list_triangles) if !isnothing(t) && t > 0.0]

    δ = minimum(list_maximum_t) / 2.0
    if side == :top
        mesh.corners.r[1:end-1, 1:end-1] .-= δ * mesh.corners.vr[1:end-1, 1:end-1]
        mesh.corners.c[1:end-1, 1:end-1] .-= δ * mesh.corners.vc[1:end-1, 1:end-1]
    elseif
        mesh.corners.r[2:end, 2:end] .-= δ * mesh.corners.vr[1:end-1, 1:end-1]
        mesh.corners.c[2:end, 2:end] .-= δ * mesh.corners.vc[1:end-1, 1:end-1]
    end

    println(
        """

        March mesh with correction_ratio δ = $(δ)
            $(field_summary(mesh.corners.vr, "∇u"))
            $(field_summary(mesh.corners.vc, "∇v"))
            $(field_summary(mesh_r - mesh.corners.r , "new mesh changes on row"))
            $(field_summary(mesh_c - mesh.corners.c, "new mesh changes on col"))
            $(field_summary(mesh_r - mesh.corners.rows_numbers , "total mesh changes on row"))
            $(field_summary(mesh_c - mesh.corners.cols_numbers, "total mesh changes on col"))
            """,
    )

    return δ
end

"""
$(SIGNATURES)

This function implements successive over relaxation for a matrix and its associated error matrix
There is a hardcoded assumption of Neumann boundary conditions -- that the derivative across the
boundary must be zero in all cases. See:
https://math.stackexchange.com/questions/3790299/how-to-iteratively-solve-poissons-equation-with-no-boundary-conditions

- ϕ is the potential to be solved subject to the Poisson equation and is the same size as the corners (_POSTS_SIZED_)
- target is ∇ϕ and is of the size in pixels (_FENCES_SIZED_)

In order of how to think about the flow of the program for the velocity potential,
    - if δ > 0, the Laplacian is too high at that point (∇²ϕ is too low).
    - To decrease the Laplacian ∇²ϕ, note that ∇² is a divergence (of a gradient) that represents net flow going outside.

    - Decrease the Laplacian
    - => decrease a divergence
    - => decrease flow towards outside

    - The flow is a gradient of the potential ϕ.
    - => decrease flow towards outside
    - => either decrease the values of the potential around that point, or increase the value of the potential at that point.
    - => Increase the value of ϕ at that point means a positive change when δ > 0.
    - δ > 0 => correction of the SAME sign.
"""
function propagate_poisson(ϕ::Matrix{Float64}, ε::Matrix{Float64})

    height, width = size(ϕ)

    ∇²ϕ_est = zeros(Float64, height, width)
    δ = zeros(Float64, height, width)
    new_ϕ = ϕ

    ω = 1.99

    for row = 1:height, col = 1:width
        val_up = row == 1 ? 0.0 : new_ϕ[row-1, col]
        val_down = row == height ? 0.0 : new_ϕ[row+1, col]
        val_left = col == 1 ? 0.0 : new_ϕ[row, col-1]
        val_right = col == width ? 0.0 : new_ϕ[row, col+1]

        ∇²ϕ_est[row, col] = val_up + val_down + val_left + val_right - 4 * new_ϕ[row, col]
        δ[row, col] = ω / 4.0 * (∇²ϕ_est[row, col] - ε[row, col])
        new_ϕ[row, col] += δ[row, col]
    end

    return new_ϕ, maximum(abs.(δ))
end

"""
$(SIGNATURES)
"""
function solve_horizontal_dispersion(mesh::FaceMesh, image; f = Focal_Length)
    # height indexes rows, width indexes columns
    height, width = size(mesh)

    # Coordinates on the caustics = simple rectangular values in corners.
    # Coordinates on the lens face comes from corner field.
    # d_row/d_col are _POSTS-SIZED_
    d_row = mesh.corners.rows_numbers - mesh.corners.r
    d_col = mesh.corners.cols_numbers - mesh.corners.c

    # The focal length is in meters. H is in corners.
    H = f / Meters_Per_Pixel

    # true_H is the real distance beteen the surface of the lens and the projection
    # H is the initial distance from the surface of the carved face to the projection screen.
    # H is constant and does not account for the change in heights due to carving.
    # Note: Higher location means closer to the caustics => negative sign
    # true_H is _POSTS_SIZED_
    true_H = zeros(Float64, height + 1, width + 1)
    true_H = -mesh.corners.ϕ .+ H

    # Normals.
    # N_row/N_col are _POSTS_SIZED_
    N_row = zeros(Float64, size(true_H))
    N_col = zeros(Float64, size(true_H))
    N_row = tan.(atan.(d_row ./ true_H) / (n₁ - n₂))
    N_col = tan.(atan.(d_col ./ true_H) / (n₁ - n₂))

    # We need to find the divergence of the Vector field described by Nx and Ny
    # div_row/div_col are _FENCES_SIZED_
    div_row = zeros(Float64, size(true_H))
    div_col = zeros(Float64, size(true_H))
    div_row[1:end-1, 1:end-1] = N_row[2:end, 1:end-1] .- N_row[1:end-1, 1:end-1]
    div_col[1:end-1, 1:end-1] = N_col[1:end-1, 2:end] .- N_col[1:end-1, 1:end-1]

    # Normalised divergence
    # divergence_direction is _FENCES_SIZED_
    divergence_direction = div_row + div_col

    ϕ_next = mesh.corners.ϕ
    max_update = Inf
    for i ∈ 1:1e6
        ϕ_next, max_update = propagate_poisson(ϕ_next, -divergence_direction)

        i % 250 == 1 &&
            println("Convergence horiz. Δ  at counter = $(i) max update = $(max_update)")

        if abs(max_update) <= 1e-6
            println(
                "Convergence horiz. Δ stopped at counter = $(i) with max_update of $(max_update)",
            )
            break
        end
    end

    return ϕ_next, max_update
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/samplefmt.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_sample_fmt_name,
    av_get_sample_fmt,
    av_get_sample_fmt_string,
    av_get_bits_per_sample_fmt,
    av_get_bytes_per_sample,
    av_sample_fmt_is_planar,
    av_samples_get_buffer_size,
    av_samples_fill_arrays,
    av_samples_alloc

function av_get_sample_fmt_name(sample_fmt::AVSampleFormat)
    ccall((:av_get_sample_fmt_name,libavutil),Ptr{UInt8},(AVSampleFormat,),sample_fmt)
end

function av_get_sample_fmt(name)
    ccall((:av_get_sample_fmt,libavutil),Cint,(Ptr{UInt8},),name)
end

function av_get_sample_fmt_string(buf,buf_size::Integer,sample_fmt::AVSampleFormat)
    ccall((:av_get_sample_fmt_string,libavutil),Ptr{UInt8},(Ptr{UInt8},Cint,AVSampleFormat),buf,buf_size,sample_fmt)
end

function av_get_bits_per_sample_fmt(sample_fmt::AVSampleFormat)
    ccall((:av_get_bits_per_sample_fmt,libavutil),Cint,(AVSampleFormat,),sample_fmt)
end

function av_get_bytes_per_sample(sample_fmt::AVSampleFormat)
    ccall((:av_get_bytes_per_sample,libavutil),Cint,(AVSampleFormat,),sample_fmt)
end

function av_sample_fmt_is_planar(sample_fmt::AVSampleFormat)
    ccall((:av_sample_fmt_is_planar,libavutil),Cint,(AVSampleFormat,),sample_fmt)
end

function av_samples_get_buffer_size(linesize,nb_channels::Integer,nb_samples::Integer,sample_fmt::AVSampleFormat,align::Integer)
    ccall((:av_samples_get_buffer_size,libavutil),Cint,(Ptr{Cint},Cint,Cint,AVSampleFormat,Cint),linesize,nb_channels,nb_samples,sample_fmt,align)
end

function av_samples_fill_arrays(audio_data,linesize,buf,nb_channels::Integer,nb_samples::Integer,sample_fmt::AVSampleFormat,align::Integer)
    ccall((:av_samples_fill_arrays,libavutil),Cint,(Ptr{Ptr{UInt8}},Ptr{Cint},Ptr{UInt8},Cint,Cint,AVSampleFormat,Cint),audio_data,linesize,buf,nb_channels,nb_samples,sample_fmt,align)
end

function av_samples_alloc(audio_data,linesize,nb_channels::Integer,nb_samples::Integer,sample_fmt::AVSampleFormat,align::Integer)
    ccall((:av_samples_alloc,libavutil),Cint,(Ptr{Ptr{UInt8}},Ptr{Cint},Cint,Cint,AVSampleFormat,Cint),audio_data,linesize,nb_channels,nb_samples,sample_fmt,align)
end

ext: jl

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
@test Psat(300.0) ≈ 0.353_658_941E-2
@test Psat(300.0u"K") ≈ 0.353_658_941E-2u"MPa"

@test Psat(500.0) ≈ 0.263_889_776E1
@test Psat(500.0u"K") ≈ 0.263_889_776E1u"MPa"

@test Psat(600.0) ≈ 0.123_443_146E2
@test Psat(600.0u"K") ≈ 0.123_443_146E2u"MPa"


text:
using Libdl

Libdl.dlopen("/System/Library/Frameworks/Foundation.framework/Foundation")

export YES, NO, nil, hostname

for c in :[NSObject
           NSBundle
           NSString
           NSArray
           NSHost].args
  @eval $(Expr(:export, c))
  @eval const $c = Class($(Expr(:quote, c)))
end

const YES = true
const NO  = false
const nil = C_NULL

toobject(s::String) = @objc [[NSString alloc] initWithUTF8String:s]

hostname() =
  (@objc [[[NSHost currentHost] localizedName] UTF8String]) |> unsafe_string

release(obj) = @objc [obj release]

function GC.gc(obj::Object)
  finalizer(obj, release)
  obj
end

function loadbundle(path)
  bundle = @objc [NSBundle bundleWithPath:path]
  bundle.ptr |> Int |> int2bool || error("Bundle $path not found")
  loadedStuff = @objc [bundle load]
  loadedStuff |> int2bool || error("Couldn't load bundle $path")
  return
end

framework(name) = loadbundle("/System/Library/Frameworks/$name.framework")


text:
# This file was generated, do not modify it. # hide
X, y = make_regression(1000, 20; sparse=0.75, noise=0.1, rng=123);


text:
using JuKeBOX
using Test

@testset "JuKeBOX.jl" begin
    # Write your tests here.
end


text:
# sys def root
const _GW_SYSTEM_DFLT_ROOT = homedir()

# depot dir
const _GW_DEPOT_NAME = ".gitworkers"
gw_depot_dir(sys_root::String) = joinpath(sys_root, _GW_DEPOT_NAME)
gw_depot_dir(gw::GitWorker) = gw_depot_dir(sys_root(gw))

# devland
const _GW_DEVLAND_FOLDER_NAME = ".devland"
gw_devland(sys_root::String) = joinpath(gw_depot_dir(sys_root), _GW_DEVLAND_FOLDER_NAME)

# worker_dir
_format_url(url::String) = replace(url, r"[^a-zA-Z0-9-_]"=> "_")
worker_dir(gw::GitWorker) = joinpath(gw_depot_dir(gw), _format_url(remote_url(gw)))
worker_relpath(gw::GitWorker, path) = _relbasepath(path, worker_dir(gw))
worker_abspath(gw::GitWorker, path) = joinpath(worker_dir(gw), worker_relpath(gw, path))

# gitlink folder
const _GL_FOLDER_NAME = ".gitlink"
gitlink_dir(gw::GitWorker) = joinpath(worker_dir(gw), _GL_FOLDER_NAME)

# repo_dir
import GitLinks.repo_dir
repo_dir(gw::GitWorker) = repo_dir(gitlink(gw))

# stage_dir
import GitLinks.stage_dir
stage_dir(gw::GitWorker) = stage_dir(gitlink(gw))

# procs dir
const _GL_PROCS_FOLDER_NAME = ".gw-procs"
gw_procs_dir(gw::GitWorker) = joinpath(worker_dir(gw), _GL_PROCS_FOLDER_NAME)

# log dir
const _GL_LOGS_FOLDER_NAME = ".gw-logs"
gw_logs_dir(gw::GitWorker) = joinpath(worker_dir(gw), _GL_LOGS_FOLDER_NAME)

# tasks dir
const _GL_TASKS_FOLDER_NAME = ".gw-tasks"
gw_tasks_dir(gw::GitWorker) = joinpath(worker_dir(gw), _GL_TASKS_FOLDER_NAME)
gw_task_dir(gw::GitWorker, tid::String) = joinpath(gw_tasks_dir(gw), tid)

# repo version
function gw_repo_mirpath(gw::GitWorker, path::String)
    wrpath = worker_relpath(gw, path)
    return joinpath(repo_dir(gw), wrpath)
end

# stage version
function gw_stage_mirpath(gw::GitWorker, path::String)
    wrpath = worker_relpath(gw, path)
    return joinpath(stage_dir(gw), wrpath)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
function make_power_curve(consumer::Consumer, days::Int64)

    d_d = Truncated(Normal(0.0, 0.5), -0.10, 0.10)
    demand_curve = Array{Float64,1}()
    for days in 1:days
        var = rand(d_d, 24).+1.0
        daily_p = (consumer.load_curve.*var)*consumer.peak_power[3]
        demand_curve = vcat(demand_curve, daily_p)
    end

    return demand_curve
end

function make_pv_curve(pvsys::PVSystem, days::Int64)

    d_s = Truncated(Normal(-0.2, 0.5), -0.10, 0.15)
    pv_curve = Array{Float64,1}()
    for days in 1:days
        solar_day_base = pvsys.time_series[(15+days)*24:(((16+days)*24) -1)]
        solar_day = [(s+rand(d_s, 1)[1])*(s > 0.1) for s in solar_day_base]*pvsys.capacity
        pv_curve = vcat(pv_curve,solar_day)
    end

    return pv_curve
end

function simulate(consumer::Consumer, pvsys::PVSystem, days::Int64)

    demand_curve = make_power_curve(consumer, days)
    pv_curve = make_pv_curve(pvsys, days)

    time_series = Dict{String,Array{Float64}}()

    available_energy_h = Array{Float64,1}()
    grid_energy_h = Array{Float64,1}()
    withdrawn_energy_h = Array{Float64,1}()
    withdrawal_h = Array{Float64,1}()
    consumer_energy_h = Array{Float64,1}()
    PV_energy_h = Array{Float64,1}()
    injection_grid_h = Array{Float64,1}()
    stored_energy_h = Array{Float64,1}()
    carry_over_h = Array{Float64,1}()

    #memory
    available_energy = 0.0
    grid_energy = 0.0
    withdrawn_energy = 0.0
    consumer_energy = 0.0
    PV_energy = 0.0
    injection_grid = 0.0
    carry_over = 0.0

    for day in 1:days
        daily_p = demand_curve[1+(day-1)*24:24+(day-1)*24]
        solar_day = pv_curve[1+(day-1)*24:24+(day-1)*24]

        #Loop for every hour of the day
        for t in 1:24

           consumer_energy += daily_p[t]
           push!(consumer_energy_h,consumer_energy)

           PV_energy +=  solar_day[t]
           push!(PV_energy_h,PV_energy)

           balance = daily_p[t] - solar_day[t]
           push!(injection_grid_h, balance)

           injection_grid += max(0.0, -1*balance)
           push!(stored_energy_h, injection_grid)

           available_energy =  carry_over + max(0.0, -1*balance)
           push!(available_energy_h, available_energy)

           withdrawal = 0.0

           if balance >= 0.0

                withdrawal = available_energy - balance

                if withdrawal > 0.0 #There is enough stored to meet energy demand in t

                    #Book Keeping

                    withdrawn_energy += balance
                    push!(withdrawal_h, balance)
                    push!(withdrawn_energy_h, withdrawn_energy)

                    grid_energy += 0.0
                    push!(grid_energy_h,grid_energy)

                    carry_over = withdrawal
                    push!(carry_over_h, withdrawal)

                elseif withdrawal <= 0.0 #There is not enough stored to meet energy demand in t

                    withdrawn_energy += available_energy
                    push!(withdrawal_h, available_energy)
                    push!(withdrawn_energy_h, withdrawn_energy)

                    carry_over = 0.0
                    push!(carry_over_h,0.0)

                    grid_energy -= withdrawal
                    push!(grid_energy_h,grid_energy)

                end

           elseif balance < 0.0

                carry_over = available_energy
                push!(carry_over_h,carry_over)
                push!(grid_energy_h,grid_energy)
                push!(withdrawal_h, 0.0)
                push!(withdrawn_energy_h, withdrawn_energy)

           end

        end #end of hourly loop

    end

    time_series["available_energy_h"] = available_energy_h
    time_series["grid_energy_h"] = grid_energy_h
    time_series["withdrawn_energy_h"] = withdrawn_energy_h
    time_series["withdrawal_h"] = withdrawal_h
    time_series["consumer_energy_h"] = consumer_energy_h
    time_series["PV_energy_h"] = PV_energy_h
    time_series["injection_grid_h"] = injection_grid_h
    time_series["stored_energy_h"] = stored_energy_h
    time_series["carry_over_h"] = carry_over_h
    time_series["demand_curve"] = demand_curve
    time_series["pv_curve"] = pv_curve

    return time_series

end

#=
plot(consumer_energy_h)
plot!(PV_energy_h+grid_energy_h-stored_energy_h+withdrawn_energy_h, fill=(0, :blue))
plot!(PV_energy_h+grid_energy_h-stored_energy_h, fill=(0, :red))
plot!(grid_energy_h, fill=(0, :green))
=#

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
"""
    sigma_matrix(λ::T, domain, u, N; num_boundary, num_interior)

Return the matrix which is used for computing `σ(λ)` in the MPS.

The matrix is of size `(num_boundary + num_interior, N)`. The element
type of the matrix is `float(T)`.

The first `num_boundary` rows consists of the first `N` basis
functions of `u` evaluated on boundary points of the domain. The
remaining `num_interior` rows corresponds to the basis functions
evaluated on interior points.
"""
function sigma_matrix(
    λ::T,
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    N::Integer;
    num_boundary::Integer = 2N,
    num_interior::Integer = 2N,
) where {T}
    @timeit_debug "points" begin
        # Compute boundary and interior points of the domain
        boundary, boundary_index = boundary_points(domain, u, N, num_boundary)
        interior = interior_points(domain, num_interior)
    end

    @timeit_debug "matrix" begin
        # Currently eigenfunctions only allow evaluation with arbs, so
        # convert λ to an arb
        λ = domain.parent(λ)

        A = Array{float(T)}(undef, num_boundary + num_interior, N)

        # TODO: Allow evaluation on ranges of indices, to be able to
        # reuse computations similar to several basis functions.
        Threads.@threads for i = 1:num_boundary
            A[i, :] = u(boundary[i], λ, 1:N, boundary = boundary_index[i])
        end
        Threads.@threads for i = 1:num_interior
            A[num_boundary+i, :] = u(interior[i], λ, 1:N)
        end
    end

    return A
end

"""
    sigma(λ, domain, u, N; num_boundary, num_interior)

Return `σ(λ)`.

Uses [`sigma_matrix`](@ref) for computing the integral which is used,
see that documentation for more details about the arguments.
"""
function sigma(
    λ::T,
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    N::Integer;
    num_boundary::Integer = 2N,
    num_interior::Integer = 2N,
) where {T}
    @timeit_debug "sigma" begin
        # Compute the matrix A
        A = sigma_matrix(
            λ,
            domain,
            u,
            N,
            num_boundary = num_boundary,
            num_interior = num_interior,
        )

        # Compute a QR factorization of A
        @timeit_debug "qr" Q = Matrix(qr(A).Q)

        # Compute the smallest singular value of the top part of Q,
        # corresponding to the boundary points
        @timeit_debug "svd" value = svdvals(Q[1:num_boundary, :])[end]
    end

    return value
end

"""
    sigma_coefficients(λ, domain, u, N; num_boundary, num_interior)

Return the coefficients corresponding to `σ(λ)`.
"""
function sigma_coefficients(
    λ::T,
    domain::AbstractDomain,
    eigenfunction::AbstractEigenfunction,
    N::Integer;
    num_boundary::Integer = 2N,
    num_interior::Integer = 2N,
) where {T}
    @timeit_debug "sigma_coefficients" begin
        # Compute the matrix A
        A = sigma_matrix(
            λ,
            domain,
            eigenfunction,
            N,
            num_boundary = num_boundary,
            num_interior = num_interior,
        )

        @timeit_debug "qr" begin
            # Compute a QR factorization of A
            q = qr(A)
            Q = Matrix(q.Q)
        end

        # Compute the right singular vector corresponding to the smallest
        # singular value

        # We create a copy of Q to be able to call svd! directly,
        # otherwise it prints a warning about alg keyword being ignored.
        @timeit_debug "svd" v = svd!(
            LinearAlgebra.copy_oftype(
                Q[1:num_boundary, :],
                LinearAlgebra.eigtype(eltype(A)),
            ),
        ).V[
            :,
            end,
        ]

        # Compute the coefficients
        coefficients = q \ (Q * v)
    end

    return coefficients
end


text:
#Wigner semicircle distribution

importall Distributions
export Semicircle

immutable Semicircle <: ContinuousUnivariateDistribution
    mean::Float64
    radius::Float64
    Semicircle(mu, r) = r > 0 ? new(float64(mu), float64(r)) :
        error("radius r must be positive") #Constructor
end

#Standardized distribution - mean 0, variance 1
Semicircle(mu) = Semicircle(mu, 2.0)
Semicircle() = Semicircle(0.0, 2.0)

# Distribution function methods
###############################

# cumulative distribution function
function cdf(X::Semicircle, x::Real)
    r, a = X.mean, X.radius
    insupport(x) ? 0.5 + (x-a)/(pi*r^2) * sqrt(r^2 - (x-a)^2) + 1/pi * asin((x-a)/r) : (x>a ? 1.0 : 0.0)
end

# probability density function
function pdf(X::Semicircle, x::Real)
    r, a = X.mean, X.radius
    insupport(x) ? 2/(pi*r^2) * sqrt(r^2 - (x-a)^2) : 0.0
end

# predicate is x in the support of the distribution?
insupport(X::Semicircle, x::Real)=abs(x)<=X.radius

#Entropy methods
################

# entropy of distribution in nats
entropy(X::Semicircle)=log(pi*X.radius) - 0.5

#Measures of central tendency methods
#####################################

# mean of distribution
mean(X::Semicircle)=X.mean

# median of distribution
median(X::Semicircle)=X.mean

# mode(s) of distribution as vector
modes(X::Semicircle)=[X.mean]

# kurtosis of the distribution
kurtosis(X::Semicircle)=2

# skewness of the distribution
skewness(X::Semicircle)=0

# standard deviation of distribution
std(X::Semicircle)=X.radius/2

# variance of distribution
var(X::Semicircle)=std(X)^2

# moment of distribution
if _HAVE_GSL
  function moment(X::Semicircle, order::Integer)
    a, r = X.mean, X.radius
    if X.mean != 0
      a^n*hypergeom([(1-n)/2, -n/2], 2, (r/a)^2)
    else
      order%2 ? (0.5*r)^(2n) * catalan(div(order,2)) : 0
    end
  end
else
  function moment(X::Semicircle, order::Integer)
    error("The moment of a semicircle distribution requires the GSL module to be loaded.")
  end
end

# cumulant of distribution
function cumulant(X::Semicircle, order::Integer)
  if X.mean != 0 error("not supported") end
  if order%2
    order==0 ? 1 : (0.5*r)^(2n) * lassalle(order/2)
  else
    0
  end
end

# free cumulant of distribution
function freecumulant(X::Semicircle, order::Integer)
  if order == 0
    return 1
  elseif order == 1
    return mean(X)
  elseif order == 2
    return var(X)
  else
    return 0
  end
end

#Generating function methods
############################

# characteristic function
function cf(X::Semicircle, t::Real)
  r = t * X.mean
  2 * besselj(1, r)/r
end

# moment generating function
function mgf(X::Semicircle, t::Real)
  r = t * X.mean
  2 * besseli(1, r)/r
end

#Sampling methods
#################

# random sampler
# Use relationship with beta distribution
function rand(X::Semicircle)
  Y = rand(Beta(1.5, 1.5))
  X.mean + 2 * X.radius * Y - X.radius
end


text:
using PAPI, Test

@testset "eventset" begin
    @testset "add/remove" begin
        evtset = EventSet()

        # add a preset
        push!(evtset, PAPI.TOT_INS)

        # add "native"
        native = event"PAPI_TOT_CYC"
        push!(evtset, native)

        @test length(evtset) == 2
        @test PAPI.TOT_INS in evtset
        @test native in evtset

        delete!(evtset, native)
        @test ! (native in evtset)

        empty!(evtset)
        @test isempty(evtset)

        append!(evtset, [PAPI.TOT_INS, native])
        @test length(evtset) == 2
        @test PAPI.TOT_INS in evtset
        @test native in evtset

        q = collect(evtset)
        @test eltype(q) == Event
        @test length(q) == 2
    end
end


text:
#
# calls.jl -
#
# Low level interface to Tcl C library.
#
# The convention is to call the functions of the Tcl C library with raw
# arguments (like pointers).  The only changes are:
#
#  - Tcl status code is returned as a `TclStatus`.
#
#  - Return values passed by reference to a C function are not needed in the
#    Julia interface and are returned as tuple of values.
#
#  - Values representing a length or an index are returned as an `Int` (not a
#    `Cint`);
#
#  - Values representing a boolean (a `Cint` in Tcl, a `Bool` in Julia) are
#    automatically converted using consistent conventions.
#
#  - Indices account for Julia convention (Julia indices start at 1, Tcl
#    indices start at 0).
#

# Fake Julia structure which reflects the layout of a Tcl_Obj and is intended
# to compute the offset of the various fields of a Tcl object.
struct __TclObj
    # Object's Reference count.  When 0 the object will be freed.
    refCount::Cint

    # Object's string representation.  This points to the first byte of the
    # object's string representation.  The array must be followed by a null
    # byte (i.e., at offset `length`) but may also contain embedded null
    # characters.  The array's storage is allocated by `ckalloc`.  NULL means
    # the string representation is invalid and must be regenerated from the
    # internal representation.  Clients should use Tcl_GetStringFromObj or
    # Tcl_GetString to get a pointer to the byte array as a readonly value.
    bytes::Ptr{Cchar}

    # The number of bytes for object's string representation, not including the
    # terminating null.
    length::Cint

    # Object's type.  Always corresponds to the type of the object's internal
    # representation.  NULL indicates the object has no specific type and no
    # internal representation.
    typePtr::Ptr{Cvoid}

    # Object's internal representation.  The value here is only valid for a
    # double precision floating-point object.  For other object types, this
    # field can be used to compute the offset of the object's internal
    # representation which is defined as an union in the C code.
    value::Cdouble
end

const __offset_of_refcount = fieldoffset(__TclObj, 1)
const __offset_of_bytes    = fieldoffset(__TclObj, 2)
const __offset_of_length   = fieldoffset(__TclObj, 3)
const __offset_of_type     = fieldoffset(__TclObj, 4)
const __offset_of_value    = fieldoffset(__TclObj, 5)

const TclObjTypePtr = fieldtype(__TclObj, :typePtr)

# Julia takes care of managing its objects so we just need to add a single
# reference for Julia for any Tcl object returned by Tcl library and make sure
# that the reference count is decremented when the Julia object is finalized.
#
# The following methods correspond to the Tcl macros which are provided to
# increment and decrement a Tcl_Obj's reference count, and to test whether an
# object is shared (i.e. has reference count > 1).
#
# The reference count of a Tcl object is an `int` which is the first member of
# the Tcl_Obj structure and we directly address it using "unsafe" operations.

@static if __offset_of_refcount != 0
    error("it is assumed that refCount comes first in Tcl_Obj structure")
end

#------------------------------------------------------------------------------
# REFERENCE COUNTING

@inline Tcl_Preserve(ptr::Ptr{T}) where {T} =
    ccall((:Tcl_Preserve, libtcl), Cvoid, (Ptr{T},), ptr)

@inline Tcl_Release(ptr::Ptr{T}) where {T} =
    ccall((:Tcl_Release, libtcl), Cvoid, (Ptr{T},), ptr)

"""
```julia
Tcl_IncrRefCount(objptr) -> objptr
```

increments the reference count of the object referenced by `objptr` and returns
this address.

"""
@inline function Tcl_IncrRefCount(objptr::TclObjPtr)
    ptr = Ptr{Cint}(objptr)
    unsafe_store!(ptr, unsafe_load(ptr) + one(Cint))
    return objptr
end

"""
```julia
Tcl_DecrRefCount(objptr) -> nothing
```

decrements the reference count of the object referenced by `objptr` and delete
it if the number of references is then smaller or equal 0.

"""
@inline function Tcl_DecrRefCount(objptr::TclObjPtr)
    ptr = Ptr{Cint}(objptr)
    newrefcount = unsafe_load(ptr) - one(Cint)
    if newrefcount ≥ 1
        unsafe_store!(ptr, newrefcount)
    else
        ccall((:TclFreeObj, libtcl), Cvoid, (TclObjPtr,), objptr)
    end
    return nothing
end

"""
```julia
Tcl_GetRefCount(objptr)
```

yields the reference count of the object referenced by `objptr`.

"""
@inline Tcl_GetRefCount(objptr::TclObjPtr) = unsafe_load(Ptr{Cint}(objptr))

"""
```julia
Tcl_IsShared(objptr)
```

yields whether the object referenced by `objptr` is shared; that is, its
reference count is greater than one.

"""
@inline Tcl_IsShared(objptr::TclObjPtr) = Tcl_GetRefCount(objptr) > 1

#------------------------------------------------------------------------------
# OBJECTS

"""
```julia
Tcl_NewBooleanObj(value) -> objptr
Tcl_NewIntObj(    value) -> objptr
Tcl_NewLongObj(   value) -> objptr
Tcl_NewWideIntObj(value) -> objptr
Tcl_NewDoubleObj( value) -> objptr
```

"""
@inline Tcl_NewBooleanObj(value::Bool) =
    ccall((:Tcl_NewBooleanObj, libtcl), TclObjPtr,
          (Cint,), (value ? one(Cint) : zero(Cint)))

for (f, Tj, Tc) in ((:Tcl_NewIntObj,     Integer, Cint),
                    (:Tcl_NewLongObj,    Integer, Clong),
                    (:Tcl_NewWideIntObj, Integer, WideInt),
                    (:Tcl_NewDoubleObj,  Real,    Cdouble))
    tup = (f, libtcl)
    @eval @inline $f(value::$Tj) = ccall($tup, TclObjPtr, ($Tc,), value)
end

@inline Tcl_NewStringObj(str::AbstractString) =
    # Use sizeof() not length() because there may be multi-byte characters
    # and use Ptr{Cchar} not Cstring because there may be embedded nulls.
    ccall((:Tcl_NewStringObj, libtcl), TclObjPtr,
          (Ptr{Cchar}, Cint), str, sizeof(str))

@inline Tcl_NewStringObj(ptr::Ptr{T}, nbytes::Integer) where {T<:Byte} =
    ccall((:Tcl_NewStringObj, libtcl), TclObjPtr, (Ptr{T}, Cint), ptr, nbytes)

@inline Tcl_NewByteArrayObj(arr::DenseArray{T}) where {T<:Byte} =
    Tcl_NewByteArrayObj(pointer(arr), sizeof(arr))

@inline Tcl_NewByteArrayObj(ptr::Ptr{T}, nbytes::Integer) where {T<:Byte} =
    ccall((:Tcl_NewByteArrayObj, libtcl), TclObjPtr,
          (Ptr{T}, Cint), ptr, nbytes)

"""
```julia
Tcl_GetBooleanFromObj(intptr, objptr) -> status::TclStatus, value::Bool
Tcl_GetIntFromObj(    intptr, objptr) -> status::TclStatus, value::Cint
Tcl_GetLongFromObj(   intptr, objptr) -> status::TclStatus, value::Clong
Tcl_GetWideIntFromObj(intptr, objptr) -> status::TclStatus, value::Tcl.WideInt
Tcl_GetDoubleFromObj( intptr, objptr) -> status::TclStatus, value::Cdouble
```

"""
@inline function Tcl_GetBooleanFromObj(intptr::TclInterpPtr, objptr::TclObjPtr)
    valref = Ref{Cint}()
    return (ccall((:Tcl_GetBooleanFromObj, libtcl), TclStatus,
                  (TclInterpPtr, TclObjPtr, Ptr{Cint}),
                  intptr, objptr, valref),
            (valref[] != zero(Cint)))
end

for (f, T) in ((:Tcl_GetIntFromObj,     Cint),
               (:Tcl_GetLongFromObj,    Clong),
               (:Tcl_GetWideIntFromObj, WideInt),
               (:Tcl_GetDoubleFromObj,  Cdouble))
    tup = (f, libtcl)
    @eval @inline function $f(intptr::TclInterpPtr, objptr::TclObjPtr)
        valref = Ref{$T}()
        return (ccall($tup, TclStatus, (TclInterpPtr, TclObjPtr, Ptr{$T}),
                      intptr, objptr, valref),
                valref[])
    end
end

"""
```julia
Tcl_GetStringFromObj(objptr) -> ptr::Ptr{Cchar}, len::Int
```

"""
@inline function Tcl_GetStringFromObj(objptr::TclObjPtr)
    lenref = Ref{Cint}()
    return (ccall((:Tcl_GetStringFromObj, libtcl), Ptr{Cchar},
                  (TclObjPtr, Ptr{Cint}), objptr, lenref),
            convert(Int, lenref[]))
end

@inline Tcl_DuplicateObj(objptr::TclObjPtr) =
    ccall((:Tcl_DuplicateObj, libtcl), TclObjPtr, (TclObjPtr,), objptr)

@inline Tcl_GetObjType(name::StringOrSymbol) =
    ccall((:Tcl_GetObjType, libtcl), TclObjTypePtr, (Cstring,), name)

#------------------------------------------------------------------------------
# INTERPRETERS AND EVALUATION OF SCRIPTS

@inline Tcl_CreateInterp() =
    ccall((:Tcl_CreateInterp, libtcl), Ptr{Cvoid}, ())

@inline Tcl_Init(intptr::TclInterpPtr) =
    ccall((:Tcl_Init, libtcl), TclStatus, (TclInterpPtr,), intptr)

@inline Tcl_InterpDeleted(intptr::TclInterpPtr) =
    (zero(Cint) != ccall((:Tcl_InterpDeleted, libtcl), Cint,
                         (TclInterpPtr,), intptr))

@inline Tcl_InterpActive(intptr::TclInterpPtr) =
    (zero(Cint) != ccall((:Tcl_InterpActive, libtcl), Cint,
                         (TclInterpPtr,), intptr))

@inline Tcl_DeleteInterp(intptr::TclInterpPtr) =
    ccall((:Tcl_DeleteInterp, libtcl), Cvoid, (TclInterpPtr,), intptr)

@inline Tcl_SetObjResult(intptr::TclInterpPtr, objptr::TclObjPtr) =
    ccall((:Tcl_SetObjResult, libtcl), Cvoid, (TclInterpPtr, TclObjPtr),
          intptr, objptr)

@inline function Tcl_SetResult(intptr::TclInterpPtr, strptr::Ptr{Cchar},
                               freeproc::Ptr{Cvoid})
    ccall((:Tcl_SetResult, libtcl), Cvoid,
          (TclInterpPtr, Ptr{Cchar}, Ptr{Cvoid}),
          intptr, strptr, freeproc)
end

@inline Tcl_GetObjResult(intptr::TclInterpPtr) =
    ccall((:Tcl_GetObjResult, libtcl), TclObjPtr, (TclInterpPtr,), intptr)

@inline function Tcl_EvalEx(intptr::TclInterpPtr, script::Ptr{Cchar},
                            nbytes::Integer, flags::Integer)
    return ccall((:Tcl_EvalEx, libtcl), TclStatus,
                 (TclInterpPtr, Ptr{Cchar}, Cint, Cint),
                 intptr, script, nbytes, flags)
end

@inline Tcl_EvalObjEx(intptr::TclInterpPtr, objptr::TclObjPtr, flags::Integer) =
    ccall((:Tcl_EvalObjEx, libtcl), TclStatus,
          (TclInterpPtr, TclObjPtr, Cint),
          intptr, objptr, flags)

@inline function Tcl_EvalObjv(intptr::TclInterpPtr, objc::Integer,
                              objv::Ptr{TclObjPtr}, flags::Integer)
    return ccall((:Tcl_EvalObjv, libtcl), TclStatus,
                 (TclInterpPtr, Cint, Ptr{TclObjPtr}, Cint),
                 intptr, objc, objv, flags)
end

"""
```julia
Tcl_DoOneEvent(flags) -> boolean
```
"""
@inline Tcl_DoOneEvent(flags::Integer) =
    (zero(Cint) != ccall((:Tcl_DoOneEvent, libtcl), Cint, (Cint,), flags))

#------------------------------------------------------------------------------
# COMMANDS

@inline function Tcl_CreateCommand(intptr::TclInterpPtr,
                                   name::StringOrSymbol,
                                   evalproc::Ptr{Cvoid},
                                   data::Ptr{Cvoid},
                                   freeproc::Ptr{Cvoid})
    return ccall((:Tcl_CreateCommand, libtcl), Ptr{Cvoid},
                 (TclInterpPtr, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                 intptr, name, evalproc, data, freeproc)
end

@inline function Tcl_CreateObjCommand(intptr::TclInterpPtr,
                                      name::StringOrSymbol,
                                      evalproc::Ptr{Cvoid},
                                      data::Ptr{Cvoid},
                                      freeproc::Ptr{Cvoid})
    return ccall((:Tcl_CreateObjCommand, libtcl), TclCommand,
                 (TclInterpPtr, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                 intptr, name, evalproc, data, freeproc)
end

@inline Tcl_DeleteCommand(intptr::TclInterpPtr, name::StringOrSymbol) =
    ccall((:Tcl_DeleteCommand, libtcl), TclStatus,
          (TclInterpPtr, Cstring), intptr, name)

@inline Tcl_DeleteCommandFromToken(intptr::TclInterpPtr, token::TclCommand) =
    ccall((:Tcl_DeleteCommandFromToken, libtcl), TclStatus,
          (TclInterpPtr, TclCommand), intptr, token)

@inline Tcl_GetCommandName(intptr::TclInterpPtr, token::TclCommand) =
    ccall((:Tcl_GetCommandName, libtcl), Cstring,
          (TclInterpPtr, TclCommand), intptr, token)

@inline function Tcl_GetCommandFullName(intptr::TclInterpPtr,
                                        token::TclCommand,
                                        objptr::TclObjPtr)
    return ccall((:Tcl_GetCommandFullName, libtcl), Cvoid,
                 (TclInterpPtr, TclCommand, TclObjPtr),
                 intptr, token, objptr)
end

@inline Tcl_GetCommandFromObj(intptr::TclInterpPtr, objptr::TclObjPtr) =
    ccall((:Tcl_GetCommandFromObj, libtcl), TclCommand,
          (TclInterpPtr, TclObjPtr), intptr, objptr)

#------------------------------------------------------------------------------
# VARIABLES

@inline function Tcl_ObjGetVar2(intptr::TclInterpPtr, name1ptr::TclObjPtr,
                                name2ptr::TclObjPtr, flags::Integer)
    return ccall((:Tcl_ObjGetVar2, libtcl), TclObjPtr,
                 (TclInterpPtr, TclObjPtr, TclObjPtr, Cint),
                 intptr, name1ptr, name2ptr, flags)
end

@inline function Tcl_ObjSetVar2(intptr::TclInterpPtr, name1ptr::TclObjPtr,
                                name2ptr::TclObjPtr, valueptr::TclObjPtr,
                                flags::Integer)
    return ccall((:Tcl_ObjSetVar2, libtcl), TclObjPtr,
                 (TclInterpPtr, TclObjPtr, TclObjPtr, TclObjPtr, Cint),
                 intptr, name1ptr, name2ptr, valueptr, flags)
end

for (Tj, Tc) in ((StringOrSymbol, Cstring),
                 (Ptr{Cchar}, Ptr{Cchar}))
    @eval begin

        @inline function Tcl_UnsetVar(intptr::TclInterpPtr, name::$Tj,
                                      flags::Integer)
            return ccall((:Tcl_UnsetVar, libtcl), TclStatus,
                         (TclInterpPtr, $Tc, Cint), intptr, name, flags)
        end

        @inline function Tcl_UnsetVar2(intptr::TclInterpPtr, name1::$Tj,
                                       name2::$Tj, flags::Integer)
            return ccall((:Tcl_UnsetVar2, libtcl), TclStatus,
                         (TclInterpPtr, $Tc, $Tc, Cint),
                         intptr, name1, name2, flags)
        end

    end
end

#------------------------------------------------------------------------------
# LISTS
#
# Note that applying a list function to any Tcl object has the side effect of
# converting the object to a list.  This may be annoying for non-temporary
# objects.

@inline Tcl_NewListObj(objc::Integer, objv::Ptr{TclObjPtr}) =
    ccall((:Tcl_NewListObj, libtcl), TclObjPtr,
          (Cint, Ptr{TclObjPtr}), objc, objv)

@inline function Tcl_SetListObj(objptr::TclObjPtr, objc::Integer,
                                objv::Ptr{TclObjPtr})
    return ccall((:Tcl_SetListObj, libtcl), Cvoid,
                 (TclObjPtr, Cint, Ptr{TclObjPtr}), objptr, objc, objv)
end

"""
```julia
Tcl_ListObjAppendList(intptr, listptr, objptr) -> status::TclStatus
```

appends each element of the list value referenced by `objptr` to to the end of
the list value referenced by `listptr`.

If `listptr` does not already point to a list value, an attempt will be made to
convert it to one.

If `objptr` is not NULL and does not already point to a list value, an attempt
will be made to convert it to one.

If an error occurs while converting a value to be a list value, an error
message is left as the result of the interpreter referenced by `intptr` unless
it is NULL.

"""
@inline function Tcl_ListObjAppendList(intptr::TclInterpPtr,
                                       listptr::TclObjPtr,
                                       objptr::TclObjPtr)
    return ccall((:Tcl_ListObjAppendList, libtcl), TclStatus,
                 (TclInterpPtr, TclObjPtr, TclObjPtr),
                 intptr, listptr, objptr)
end

"""
```julia
Tcl_ListObjAppendElement(intptr, listptr, objptr) -> status::TclStatus
```

appends the single value referenced by `objptr` to to the end of the list value
referenced by `listptr`.

The object referenced by `listptr` must not be shared (its reference count must
be ≤ 1) otherwise  Tcl will panic (and abort the program).  To avoid aborting,
an error is reported.

If `listptr` does not already point to a list value, an attempt will be made to
convert it to one.

If an error occurs while converting a value to be a list value, an error
message is left as the result of the interpreter referenced by `intptr` unless
it is NULL.

"""
@inline function Tcl_ListObjAppendElement(intptr::TclInterpPtr,
                                          listptr::TclObjPtr,
                                          objptr::TclObjPtr) :: TclStatus
    if Tcl_IsShared(listptr)
        msg = "modifying a shared Tcl list is forbidden"
        if intptr == C_NULL
            warn(msg, once=true)
        else
            Tcl_SetResult(intptr, msg)
        end
        return TCL_ERROR
    end
    return ccall((:Tcl_ListObjAppendElement, libtcl), TclStatus,
                 (TclInterpPtr, TclObjPtr, TclObjPtr),
                 intptr, listptr, objptr)
end

"""
```julia
Tcl_ListObjGetElements(intptr, listptr)
    -> status::TclStatus, objc::Int, objv::Ptr{TclObjPtr}
```

Does not touch the reference count of the list object and of its elements.

"""
@inline function Tcl_ListObjGetElements(intptr::TclInterpPtr,
                                        listptr::TclObjPtr)
    objc = Ref{Cint}()
    objv = Ref{Ptr{TclObjPtr}}()
    return (ccall((:Tcl_ListObjGetElements, libtcl), TclStatus,
                  (TclInterpPtr, TclObjPtr, Ptr{Cint}, Ptr{Ptr{TclObjPtr}}),
                  intptr, listptr, objc, objv),
            objc[], objv[])
end

"""
```julia
Tcl_ListObjLength(intptr, listptr) -> status::TclStatus, length::Int
```
"""
@inline function Tcl_ListObjLength(intptr::TclInterpPtr, listptr::TclObjPtr)
    lenref = Ref{Cint}()
    return (ccall((:Tcl_ListObjLength, libtcl), TclStatus,
                  (TclInterpPtr, TclObjPtr, Ptr{Cint}),
                  intptr, listptr, lenref),
            convert(Int, lenref[]))
end

@inline function Tcl_ListObjIndex(intptr::TclInterpPtr, listptr::TclObjPtr,
                                  index::Integer)
    objptr = Ref{TclObjPtr}()
    return (ccall((:Tcl_ListObjIndex, libtcl), TclStatus,
                  (TclInterpPtr, TclObjPtr, Cint, Ptr{TclObjPtr}),
                  intptr, listptr, index - 1, objptr),
            objptr[])
end

"""
```julia
Tcl_ListObjReplace(intptr, listptr, first, count,
                   objc, objv) -> status::TclStatus
```

`first` can be the length of the list plus one and `count` can be 0 to append
to the end of the list.

`objc = 0` and `objv = NULL` are OK to delete elements.

The object referenced by `listptr` must not be shared (its reference count must
be ≤ 1) otherwise  Tcl will panic (and abort the program).  To avoid aborting,
an error is reported.

If `listptr` does not already point to a list value, an attempt will be made to
convert it to one.

"""
@inline function Tcl_ListObjReplace(intptr::TclInterpPtr, listptr::TclObjPtr,
                                    first::Integer, count::Integer,
                                    objc::Integer, objv::Ptr{TclObjPtr})
    if Tcl_IsShared(listptr)
        msg = "modifying a shared Tcl list is forbidden"
        if intptr == C_NULL
            warn(msg, once=true)
        else
            Tcl_SetResult(intptr, msg)
        end
        return TCL_ERROR
    end
    return ccall((:Tcl_ListObjReplace, libtcl), TclStatus,
                 (TclInterpPtr, TclObjPtr, Cint, Cint, Cint, Ptr{TclObjPtr}),
                 intptr, listptr, first + 1, count, objc, objv)
end


text:
struct KeyFrame
    id::Int64
    left_pyramid::Union{Nothing, LKPyramid{Vector{Matrix{Gray{Float64}}}, LKCache}}
    right_image::Union{Nothing, Matrix{Gray{Float64}}}
end

mutable struct Mapper
    params::Params
    map_manager::MapManager
    estimator::Estimator

    current_frame::Frame
    keyframe_queue::Vector{KeyFrame}
    right_pyramid::LKPyramid
    geev_cache::GEEV4x4Cache

    exit_required::Bool
    new_kf_available::Bool

    estimator_thread
    queue_lock::ReentrantLock
end

function Mapper(params::Params, map_manager::MapManager, frame::Frame, slam_io)
    estimator = Estimator(map_manager, params, slam_io)
    estimator_thread = Threads.@spawn run!(estimator)
    @debug "[MP] Launched estimator thread."
    empty_pyr = LKPyramid(
        [Matrix{Gray{Float64}}(undef, 0, 0)],
        nothing, nothing, nothing, nothing, nothing, nothing)
    Mapper(
        params, map_manager, estimator,
        frame, KeyFrame[], empty_pyr, GEEV4x4Cache(),
        false, false, estimator_thread, ReentrantLock())
end

function run!(mapper::Mapper)
    while !mapper.exit_required
        succ, kf = get_new_kf!(mapper)
        if !succ
            sleep(1e-2)
            continue
        end

        new_keyframe = get_keyframe(mapper.map_manager, kf.id)
        new_keyframe ≡ nothing && @error "[MP] Got invalid frame $(kf.id) from Map"

        if mapper.params.stereo
            try
                t1 = time()
                if has_gradients(mapper.right_pyramid)
                    update!(mapper.right_pyramid, kf.right_image)
                else
                    mapper.right_pyramid = LKPyramid(
                        kf.right_image, mapper.params.pyramid_levels;
                        σ=mapper.params.pyramid_σ, reusable=true)
                end
                optical_flow_matching!(
                    mapper.map_manager, new_keyframe,
                    kf.left_pyramid, mapper.right_pyramid, true)
                t2 = time()
                @debug "[MP] Stereo Matching ($(t2 - t1) sec): $(new_keyframe.nb_stereo_kpts) Keypoints"
            catch e
                showerror(stdout, e)
                display(stacktrace(catch_backtrace()))
            end

            if new_keyframe.nb_stereo_kpts > 0
                lock(mapper.map_manager.map_lock)
                try
                    t1 = time()
                    triangulate_stereo!(
                        mapper.map_manager, new_keyframe,
                        mapper.params.max_reprojection_error, mapper.geev_cache)
                    t2 = time()
                    @debug "[MP] Stereo Triangulation ($(t2 - t1) sec)."
                catch e
                    showerror(stdout, e)
                    display(stacktrace(catch_backtrace()))
                finally
                    unlock(mapper.map_manager.map_lock)
                end
            end
        end

        if new_keyframe.nb_2d_kpts > 0 && new_keyframe.kfid > 0
            lock(mapper.map_manager.map_lock)
            try
                t1 = time()
                triangulate_temporal!(
                    mapper.map_manager, new_keyframe,
                    mapper.params.max_reprojection_error, mapper.geev_cache)
                t2 = time()
                @debug "[MP] Temporal Triangulation ($(t2 - t1) sec)."
            catch e
                showerror(stdout, e)
                display(stacktrace(catch_backtrace()))
            finally
                unlock(mapper.map_manager.map_lock)
            end
        end

        # Check if reset is required.
        if mapper.params.vision_initialized
            if kf.id == 1 && new_keyframe.nb_3d_kpts < 30
                @warn "[MP] Bad initialization detected. Resetting!"
                mapper.params.reset_required = true
                mapper |> reset!
                continue
            elseif kf.id < 10 && new_keyframe.nb_3d_kpts < 3
                @warn "[MP] Reset required. Nb 3D points: $(new_keyframe.nb_3d_kpts)."
                mapper.params.reset_required = true
                mapper |> reset!
                continue
            end
        end

        t1 = time()
        update_frame_covisibility!(mapper.map_manager, new_keyframe)
        t2 = time()
        @debug "[MP] Covisibility update ($(t2 - t1) sec)."

        if mapper.params.do_local_matching && kf.id > 0
            t1 = time()
            try
                match_local_map!(mapper, new_keyframe)
            catch e
                showerror(stdout, e)
                display(stacktrace(catch_backtrace()))
            end
            t2 = time()
            @debug "[MP] Local Matching ($(t2 - t1) sec)."
        end

        add_new_kf!(mapper.estimator, new_keyframe)
    end
    mapper.estimator.exit_required = true
    @debug "[MP] Exit required."
    wait(mapper.estimator_thread)
end

function triangulate_stereo!(
    map_manager::MapManager, frame::Frame, max_error, cache::GEEV4x4Cache,
)
    stereo_keypoints = get_stereo_keypoints(frame)
    if isempty(stereo_keypoints)
        @warn "[MP] No stereo keypoints to triangulate."
        return
    end

    P1 = to_4x4(frame.camera.K) * SMatrix{4, 4, Float64, 16}(I)
    P2 = to_4x4(frame.right_camera.K) * frame.right_camera.Ti0

    n_good = 0
    for kp in stereo_keypoints
        kp.is_3d && continue
        mp = get_mappoint(map_manager, kp.id)
        mp ≡ nothing &&
            (remove_mappoint_obs!(map_manager, kp.id, frame.kfid); continue)
        mp.is_3d && continue

        left_point = triangulate(
            kp.undistorted_pixel[[2, 1]], kp.right_undistorted_pixel[[2, 1]],
            P1, P2, cache)
        left_point *= 1.0 / left_point[4]
        left_point[3] < 0.1 && (remove_stereo_keypoint!(frame, kp.id); continue)

        right_point = frame.right_camera.Ti0 * left_point
        right_point[3] < 0.1 && (remove_stereo_keypoint!(frame, kp.id); continue)

        left_projection = project(frame.camera, left_point[1:3])
        lrepr = norm(kp.undistorted_pixel .- left_projection)
        lrepr > max_error && (remove_stereo_keypoint!(frame, kp.id); continue)

        right_projection = project(frame.right_camera, right_point[1:3])
        rrepr = norm(kp.right_undistorted_pixel .- right_projection)
        rrepr > max_error && (remove_stereo_keypoint!(frame, kp.id); continue)

        wpt = project_camera_to_world(frame, left_point)[1:3]
        update_mappoint!(map_manager, kp.id, wpt)
        n_good += 1
    end
end

function triangulate_temporal!(
    map_manager::MapManager, frame::Frame, max_error, cache::GEEV4x4Cache,
)
    keypoints = get_2d_keypoints(frame)
    if isempty(keypoints)
        @warn "[MP] No 2D keypoints to triangulate."
        return
    end

    K = to_4x4(frame.camera.K)
    # P1 - previous Keyframe, P2 - this `frame`.
    P1 = K * SMatrix{4, 4, Float64, 16}(I)
    P2 = K * SMatrix{4, 4, Float64, 16}(I)

    good = 0
    rel_kfid = -1
    # frame -> observer key frame.
    rel_pose::SMatrix{4, 4, Float64, 16} = SMatrix{4, 4, Float64, 16}(I)
    rel_pose_inv::SMatrix{4, 4, Float64, 16} = SMatrix{4, 4, Float64, 16}(I)

    for kp in keypoints
        @assert !kp.is_3d

        map_point = get_mappoint(map_manager, kp.id)
        map_point ≡ nothing &&
            (remove_mappoint_obs!(map_manager, kp.id, frame.kfid); continue)
        map_point.is_3d && continue

        # Get first KeyFrame id from the set of mappoint observers.
        observers = get_observers(map_point)
        length(observers) < 2 && continue
        kfid = observers[1]
        frame.kfid == kfid && continue

        observer_kf = get_keyframe(map_manager, kfid)
        if observer_kf ≡ nothing
            @error "[MP] Missing observer for triangulation."
            continue
        end

        # Compute relative motion between new KF & observer KF.
        # Don't recompute if the frame's ids don't change.
        if rel_kfid != kfid
            rel_pose = observer_kf.cw * frame.wc
            rel_pose_inv = inv(SE3, rel_pose)
            rel_kfid = kfid
            P2 = K * rel_pose_inv
        end

        observer_kp = get_keypoint(observer_kf, kp.id)
        observer_kp ≡ nothing && continue
        obup = observer_kp.undistorted_pixel
        kpup = kp.undistorted_pixel

        parallax = norm(
            obup .- project(frame.camera, rel_pose[1:3, 1:3] * kp.position))

        left_point = triangulate(obup[[2, 1]], kpup[[2, 1]], P1, P2, cache)
        left_point *= 1.0 / left_point[4]
        left_point[3] < 0.1 && parallax > 20.0 && (remove_mappoint_obs!(
            map_manager, observer_kp.id, frame.kfid); continue)

        right_point = rel_pose_inv * left_point
        right_point[3] < 0.1 && parallax > 20.0 && (remove_mappoint_obs!(
            map_manager, observer_kp.id, frame.kfid); continue)

        lrepr = norm(project(frame.camera, left_point[1:3]) .- obup)
        lrepr > max_error && parallax > 20.0 && (remove_mappoint_obs!(
            map_manager, observer_kp.id, frame.kfid); continue)

        rrepr = norm(project(frame.camera, right_point[1:3]) .- kpup)
        rrepr > max_error && parallax > 20.0 && (remove_mappoint_obs!(
            map_manager, observer_kp.id, frame.kfid); continue)

        wpt = project_camera_to_world(observer_kf, left_point)[1:3]
        update_mappoint!(map_manager, kp.id, wpt)
        good += 1
    end
end

"""
Try matching keypoints from `frame` with keypoints from frames in its
covisibility graph (aka local map).
"""
function match_local_map!(mapper::Mapper, frame::Frame)
    # Maximum number of MapPoints to track.
    max_nb_mappoints = 10 * mapper.params.max_nb_keypoints
    covisibility_map = get_covisible_map(frame)

    if length(frame.local_map_ids) < max_nb_mappoints
        # Get local map of the oldest covisible KeyFrame and add it
        # to the local map to `frame` to search for MapPoints.
        kfid = collect(keys(covisibility_map))[1]
        co_kf = get_keyframe(mapper.map_manager, kfid)
        while co_kf ≡ nothing && kfid > 0
            kfid -= 1
            co_kf = get_keyframe(mapper.map_manager, kfid)
        end

        co_kf ≢ nothing && union!(frame.local_map_ids, co_kf.local_map_ids)
        # TODO if still not enough, go for another round.
    end

    prev_new_map = do_local_map_matching(
        mapper, frame, frame.local_map_ids;
        max_projection_distance=mapper.params.max_projection_distance,
        max_descriptor_distance=mapper.params.max_descriptor_distance)

    isempty(prev_new_map) || merge_matches(mapper, prev_new_map)
end

function merge_matches(mapper::Mapper, prev_new_map::Dict{Int64, Int64})
    lock(mapper.map_manager.optimization_lock)
    lock(mapper.map_manager.map_lock)
    try
        for (prev_id, new_id) in prev_new_map
            merge_mappoints(mapper.map_manager, prev_id, new_id);
        end
    catch e
        showerror(stdout, e); println()
        display(stacktrace(catch_backtrace())); println()
    finally
        unlock(mapper.map_manager.map_lock)
        unlock(mapper.map_manager.optimization_lock)
    end
end

"""
Given a frame and its local map of Keypoints ids (triangulated),
project respective mappoints onto the frame, find surrounding keypoints (triangulated?),
match surrounding keypoints with the projection.
Best match is the new candidate for replacement.
"""
function do_local_map_matching(
    mapper::Mapper, frame::Frame, local_map::Set{Int64};
    max_projection_distance, max_descriptor_distance,
)
    prev_new_map = Dict{Int64, Int64}()
    isempty(local_map) && return prev_new_map

    # Maximum field of view.
    vfov = 0.5 * frame.camera.height / frame.camera.fy
    hfov = 0.5 * frame.camera.width / frame.camera.fx
    max_rad_fov = vfov > hfov ? atan(vfov) : atan(hfov)
    view_threshold = cos(max_rad_fov)

    # Define max distance from projection.
    frame.nb_3d_kpts < 30 && (max_projection_distance *= 2.0;)
    # matched kpid → [(local map kpid, distance)] TODO
    matches = Dict{Int64, Vector{Tuple{Int64, Float64}}}()

    # Go through all MapPoints from the local map in `frame`.
    for kpid in local_map
        is_observing_kp(frame, kpid) && continue
        mp = get_mappoint(mapper.map_manager, kpid)
        mp ≡ nothing && continue
        (!mp.is_3d || isempty(mp.descriptor)) && continue

        # Project MapPoint into KeyFrame's image plane.
        position = get_position(mp)
        camera_position = project_world_to_camera(frame, position)[1:3]
        camera_position[3] < 0.1 && continue

        view_angle = camera_position[3] / norm(camera_position)
        abs(view_angle) < view_threshold && continue

        projection = project_undistort(frame.camera, camera_position)
        in_image(frame.camera, projection) || continue

        surrounding_keypoints = get_surrounding_keypoints(frame, projection)

        # Find best match for the `mp` among `surrounding_keypoints`.
        best_id, best_distance = find_best_match(
            mapper.map_manager, frame, mp, projection, surrounding_keypoints;
            max_projection_distance, max_descriptor_distance)
        best_id == -1 && continue

        match = (kpid, best_distance)
        if best_id in keys(matches)
            push!(matches[best_id], match)
        else
            matches[best_id] = Tuple{Int64, Float64}[match]
        end
    end

    for (kpid, match) in matches
        best_distance = 1e6
        best_id = -1

        for (local_kpid, distance) in match
            if distance ≤ best_distance
                best_distance = distance
                best_id = local_kpid
            end
            best_id != -1 && (prev_new_map[kpid] = best_id;)
        end
    end
    prev_new_map
end

"""
For a given `target_mp` MapPoint, find best match among surrounding keypoints.

Given target mappoint from covisibility graph, its projection onto `frame`
and surrounding keypoints in `frame` for that projection,
find best matching keypoint (already triangulated?) in `frame`.
"""
function find_best_match(
    map_manager::MapManager, frame::Frame, target_mp::MapPoint,
    projection, surrounding_keypoints;
    max_projection_distance, max_descriptor_distance,
)
    target_mp_observers = get_observers(target_mp)
    target_mp_position = get_position(target_mp)

    # TODO parametrize descriptor size.
    min_distance = 256.0 * max_descriptor_distance
    best_distance, second_distance = min_distance, min_distance
    best_id, second_id = -1, -1

    for kp in surrounding_keypoints
        kp.id < 0 && continue
        distance = norm(projection .- kp.pixel)
        distance > max_projection_distance && continue

        mp = get_mappoint(map_manager, kp.id)
        if mp ≡ nothing
            remove_mappoint_obs!(map_manager, kp.id, frame.kfid)
            continue
        end
        isempty(mp.descriptor) && continue

        # Check that `kp` and `target_mp` are indeed candidates for matching.
        # They should have no overlap in their observers.
        mp_observers = get_observers(mp)
        isempty(intersect(target_mp_observers, mp_observers)) || continue

        avg_projection = 0.0
        n_projections = 0

        # Compute average projection distance for the `target_mp` projected
        # into each of the `mp` observers KeyFrame.
        for observer_kfid in mp_observers
            observer_kf = get_keyframe(map_manager, observer_kfid)
            observer_kf ≡ nothing && (remove_mappoint_obs!(
                map_manager, kp.id, observer_kfid); continue)

            observer_kp = get_keypoint(observer_kf, kp.id)
            observer_kp ≡ nothing && (remove_mappoint_obs!(
                map_manager, kp.id, observer_kfid); continue)

            observer_projection = project_world_to_image_distort(
                observer_kf, target_mp_position)
            avg_projection += norm(observer_kp.pixel .- observer_projection)
            n_projections += 1
        end
        avg_projection /= n_projections
        avg_projection > max_projection_distance && continue

        distance = mappoint_min_distance(target_mp, mp)
        if distance ≤ best_distance
            second_distance = best_distance
            second_id = best_id

            best_distance = distance
            best_id = kp.id
        elseif distance ≤ second_distance
            second_distance = distance
            second_id = kp.id
        end
    end

    # TODO is this necessary?
    # best_id != -1 && second_id != -1 &&
    #     0.9 * second_distance < best_distance && (best_id = -1;)

    best_id, best_distance
end

function get_new_kf!(mapper::Mapper)
    lock(mapper.queue_lock) do
        if isempty(mapper.keyframe_queue)
            mapper.new_kf_available = false
            return false, nothing
        end

        keyframe = popfirst!(mapper.keyframe_queue)
        mapper.new_kf_available = !isempty(mapper.keyframe_queue)
        true, keyframe
    end
end

function add_new_kf!(mapper::Mapper, kf::KeyFrame)
    lock(mapper.queue_lock) do
        push!(mapper.keyframe_queue, kf)
        mapper.new_kf_available = true
    end
end

function reset!(mapper::Mapper)
    lock(mapper.queue_lock) do
        mapper.new_kf_available = false
        mapper.exit_required = false
        mapper.keyframe_queue |> empty!
    end
end


text:
# Parameters:
#	N::Int
#	- the problem size
#	K::Int
#	- the problem is homogenous at N-K
# Return Values:
#	P::Array{Array{Array{Int, 1}, 1}, 1} (Partitions)
#	- P[n][p] contains the pth Partition of n that is required for the bandlimited functionality
#	- This is equivalent to the Partitions that correspond to non-zero frequencies in FFT and their immediate descendents
#	- length(P[n]) = 0 for n = 1:(N-K-1)
#	ZFI::Array{Int, 1} (Zero Frequency Information)
#	- ZFI[n] = k if, for p<=k, P[n][p] is a zero frequency partition
#	WI::Array{Int, 2} (Width Information)
#	- WI[n, w] contains the number of Partitions of n whose first element is less than or equal to w
#	- This only counts the Partitions that are in P and not all Partitions
function partitions_bl(N::Int, K::Int)
	P_K, WI_K = partitions(K)
	WI = zeros(Int, N, N)
	for n = (N - K):N
		WI[n,n] = 1
	end
	for n = (N - K):(N - 1)
		for w = (N - K - 1):(n - 1)
			WI[n,w] = WI_K[n - w, min(n - w, w)]
		end
	end
	for w = (N - K):(N - 1)
		WI[N,w] = WI_K[N - w, min(N - w, w)]
	end
	for n = (N - K):N
		for w = (N - K):n
			WI[n, w] += WI[n, w - 1]
		end
	end
	P = Array(Array{Array{Int, 1}, 1}, N)
	ZFI = Array(Int, N)
	for n = 1:(N - K - 1)
		P[n] = Array(Array{Int, 1}, 0)
		ZFI[n] = 0
	end
	for n = (N - K):(N - 1)
		Pn = Array(Array{Int, 1}, WI[n,n])
		i = 1
		for spi = 1:WI_K[n - N + K + 1, min(n - N + K + 1, N - K - 1)]
			Pn[i] = [N - K - 1; P_K[n - N + K + 1][spi]]
			i += 1
		end
		ZFI[n] = i - 1
		for w = (N - K):(n - 1)
			for spi = 1:WI_K[n - w, min(n - w, w)]
				Pn[i] = [w; P_K[n - w][spi]]
				i += 1
			end
		end
		Pn[i] = [n]
		P[n] = Pn
	end
	Pn = Array(Array{Int, 1}, WI[N, N])
	i = 1
	for w = (N - K):(N - 1)
		for spi = 1:WI_K[N - w, min(N - w, w)]
			Pn[i] = [w; P_K[N - w][spi]]
			i += 1
		end
	end
	Pn[i] = [N]
	P[N] = Pn
	ZFI[N] = 0
	return P, ZFI, WI
end


text:
in_channel = 10
out_channel = 5
num_V = 6
num_E = 7
T = Float32

adj = T[0. 1. 0. 0. 0. 0.;
       1. 0. 0. 1. 1. 1.;
       0. 0. 0. 0. 0. 1.;
       0. 1. 0. 0. 1. 0.;
       0. 1. 0. 1. 0. 1.;
       0. 1. 1. 0. 1. 0.]

struct NewGNLayer <: GraphNet
end

V = rand(T, in_channel, num_V)
E = rand(T, in_channel, 2num_E)
u = rand(T, in_channel)

@testset "gn" begin
    l = NewGNLayer()

    @testset "without aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg)

        fg = FeaturedGraph(adj, nf=V)
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (0, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    @testset "with neighbor aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg, :add)

        fg = FeaturedGraph(adj, nf=V, ef=E, gf=zeros(0))
        l = NewGNLayer()
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (in_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.update_edge(l::NewGNLayer, e, vi, vj, u) = rand(T, out_channel)
    @testset "update edge with neighbor aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg, :add)

        fg = FeaturedGraph(adj, nf=V, ef=E, gf=zeros(0))
        l = NewGNLayer()
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (out_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.update_vertex(l::NewGNLayer, ē, vi, u) = rand(T, out_channel)
    @testset "update edge/vertex with all aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg, :add, :add, :add)

        fg = FeaturedGraph(adj, nf=V, ef=E, gf=u)
        l = NewGNLayer()
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (out_channel, num_V)
        @test size(edge_feature(fg_)) == (out_channel, 2*num_E)
        @test size(global_feature(fg_)) == (in_channel,)
    end
end


text:
using Quantica: TightbindingModel, OnsiteTerm, HoppingTerm, padtotype, Selector, sublats, resolve, isinindices, isinsublats, nsites
using LinearAlgebra: norm

@testset "selector membership" begin
    @test  isinindices(1, missing)
    @test  isinindices(1, 1)
    @test !isinindices(2, 1)
    @test  isinindices(1, (1,2))
    @test !isinindices(3, (1,2))
    @test  isinindices(3, (1, 3:4))
    @test !isinindices(2, (1, 3:4))
    @test !isinindices(3, [1:2, 5:6])

    @test !isinindices(1, not(missing))
    @test !isinindices(1, not(1))
    @test  isinindices(2, not(1))
    @test !isinindices(1, not(1,2))
    @test  isinindices(3, not(1,2))
    @test !isinindices(3, not(1, 3:4))
    @test  isinindices(2, not(1, 3:4))
    @test  isinindices(3, not([1:2, 5:6]))

    @test  isinsublats(3, (not(1,2), 3))
    @test  isinsublats(3, (not(1,2), 4))
    @test !isinsublats(2, (not(1,2), 3))
    @test  isinsublats(3, (not(1,2), 3))

    @test  isinindices(1=>2, 1=>2)
    @test !isinindices(1=>2, 1=>3)
    @test  isinindices(1=>2, (1=>2, 3=>4))
    @test !isinindices(1=>2, (1=>3, 2=>1))
    @test  isinindices(1=>3, 1:2=>3:4)
    @test !isinindices(1=>2, 1:2=>3:4)
    @test  isinindices(1=>2, (1:2=>3:4, 1=>2))
    @test !isinindices(3=>4, (1:2=>3:4, 1=>2))
    @test  isinindices(1=>2, (1:2=>3:4, (1,2) .=> (2,1)))
    @test  isinindices(1=>3, (1:2=>3:4, 1:2 .=> (2,1)))
    @test !isinindices(3=>4, (1:2=>3:4, 1:2 .=> (2,1), 4=>3))

    @test !isinindices(1=>2, not(1=>2))
    @test  isinindices(1=>2, not(1=>3))
    @test !isinindices(1=>2, not(1=>2, 3=>4))
    @test  isinindices(1=>2, not(1=>3, 2=>1))
    @test !isinindices(1=>3, not(1:2=>3:4))
    @test  isinindices(1=>2, not(1:2=>3:4))
    @test !isinindices(1=>2, not(1:2=>3:4, 1=>2))
    @test  isinindices(3=>4, not(1:2=>3:4, 1=>2))
    @test !isinindices(1=>2, not(1:2=>3:4, (1,2) .=> (2,1)))
    @test !isinindices(1=>3, not(1:2=>3:4, 1:2 .=> (2,1)))
    @test  isinindices(3=>4, not(1:2=>3:4, 1:2 .=> (2,1), 4=>3))

    @test  isinindices(3=>4, not(1:2)=>3:4)
    @test !isinindices(1=>4, (not(1:2)=>3:4, 1:2 .=> (2,1), 4=>3))
    @test  isinindices(3=>4, (not(1:2)=>3:4, 1:2 .=> (2,1), 4=>3))
end

@testset "selector usage" begin
    lat = LatticePresets.honeycomb() |> unitcell(2)
    rs = (r->true, missing)
    ss = (:A, (:A, :B), missing)
    inds = (1, (1, 3), (1, 3:4), missing)
    for r in rs, s in ss, i in inds
        sel = siteselector(region = r, sublats = s, indices = i)
        rsel = resolve(sel, lat)
        @test 1 in rsel
        if s === :A || i !== missing
            @test !in(5, rsel)
        end
    end

    lat = LatticePresets.honeycomb() |> unitcell(2)
    rs =   ((r,dr)->true, missing)
    ss =   (:A=>:B, (:A=>:B, :B=>:C), (:A,:B)=>(:B,:A), ((:A,:B)=>(:B,:A), :A=>:C), ((:A,:B)=>(:B,:A), (:C,:D).=>(:C,:D)), ((:A,)=>(:A,), (:C,:D).=>(:C,:D), :C=>:C), missing)
    inds = (1=>2, (1=>2, 1=>4), 1:2=>3:4, (1:2=>3:4, 1=>2), (1:2=>3:4, (1,2) .=> (2,1)), (1:2=>3:4, 1:2 .=> (2,1), 4=>3), missing)
    for r in rs, s in ss, i in inds
        sel = hopselector(region = r, sublats = s, indices = i)
        rsel = resolve(sel, lat)
        conds = s in (:A=>:B, (:A=>:B, :B=>:C))
        condi = i in (1=>2,)
        if (conds !== missing && conds) || (condi !== missing && condi)
            @test !in(1=>4, rsel)
        else
            @test in(1=>4, rsel)
        end
    end

    h  = hamiltonian(lat, hopping(1, sublats = not(:A => :B)))
    h´ = hamiltonian(lat, hopping(1, sublats = :A => :B))
    @test !ishermitian(h) && !ishermitian(h´)
    @test bloch(h) == transpose(bloch(h´))
end

@testset "term algebra" begin
    r = SVector(0, 0)
    model = onsite(1) + hopping(2I)
    @test (t -> t(r, r)).(model.terms) == (1, 2I)
    model = onsite(1) - hopping(2I)
    @test (t -> t(r, r)).(model.terms) == (1, -2I)
    model = -onsite(@SMatrix[1 0; 1 1]) - 2hopping(2I)
    @test (t -> t(r, r)).(model.terms) == (-@SMatrix[1 0; 1 1], -4I)
end

@testset "onsite terms" begin
    rs = (r->true, missing)
    ss = (:A, (:A, :B), missing)
    inds = (1, (1, 2), (1:2, 3), missing)
    for r in rs, s in ss, i in inds
        model0 = onsite(1, region = r, sublats = s, indices = i) + hopping(1)
        model1 = onsite(1) + hopping(1)
        model2 = onsite(1, sublats = s) + hopping(1)
        model3 = onsite(1, region = r) + hopping(1)
        model4 = onsite(1, indices = i) + hopping(1)
        @test onsite(model1, region = r, sublats = s, indices = i) + hopping(model1) === model0
        @test onsite(model2, region = r, sublats = s, indices = i) + hopping(model2) === model0
        @test onsite(model3, region = r, sublats = s, indices = i) + hopping(model3) === model0
        @test onsite(model4, region = r, sublats = s, indices = i) + hopping(model3) === model0
    end
end

@testset "hopping terms" begin
    rs = ((r, dr) -> true, missing)
    ss = (:A => :B, :A => (:A,:B), (:A,:B) .=> (:A,:B), (:A,:B) => (:A,:B), missing)
    dns = ((0,1), ((0,1),(1,0)), SVector(0,1), (SVector(0,1), (0,3)), [1, 2], ([1.0,2], (0,4.0)), missing)
    inds = ((1 => 2, ), [1 => 2], (1 => 2, 3 => 4), (1:2 => 3:4, 1 => 2), (1:2 .=> 3:4, 1 => 2), missing)
    ranges = (Inf, 1)  # no missing here, because hopping range default is 1.0
    for r in rs, s in ss, dn in dns, rn in ranges
        model0 = hopping(1, region = r, sublats = s, dn = dn, range = rn, indices = inds) + onsite(1)
        model1 = hopping(1, region = r, sublats = s, dn = dn) + onsite(1)
        model2 = hopping(1, region = r, sublats = s) + onsite(1)
        model3 = hopping(1, region = r, range = rn) + onsite(1)
        model4 = hopping(1) + onsite(1)
        @test hopping(model1, region = r, sublats = s, dn = dn, range = rn, indices = inds) + onsite(model1) === model0
        @test hopping(model2, region = r, sublats = s, dn = dn, range = rn, indices = inds) + onsite(model2) === model0
        @test hopping(model3, region = r, sublats = s, dn = dn, range = rn, indices = inds) + onsite(model3) === model0
        @test hopping(model4, region = r, sublats = s, dn = dn, range = rn, indices = inds) + onsite(model4) === model0
    end
end

@testset "hopping adjoint" begin
    ts  = (1, 1im, @SMatrix[0 -im; im 0], @SMatrix[0 -im])
    ss  = (:A => :B, :A => (:A,:B), (:A,:B) .=> (:A,:B), (:A,:B) => (:B,:A), missing)
    ss´ = (:B => :A, (:A,:B) => :A, (:A,:B) .=> (:A,:B), (:B,:A) => (:A,:B), missing)
    dns = ((0,1), ((0,1),(1,0)),    SVector(0,1),  (SVector(0,1), (0,3)),    [1, 2], ([1.0,2], (0,4.0)),   missing)
    dns´= ((0,-1), ((0,-1),(-1,0)), SVector(0,-1), (SVector(0,-1), (0,-3)), -[1, 2], (-[1.0,2], (0,-4.0)), missing)
    ranges = (Inf, 1)
    for (t,t´) in zip(ts, adjoint.(ts)), (s, s´) in zip(ss, ss´), (dn, dn´) in zip(dns, dns´), rn in ranges
        hop = hopping(t, sublats = s, dn = dn, range = rn)
        hop´ = hopping(t´, sublats = s´, dn = dn´, range = rn)
        @test hop' === hop´
    end
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using NearestNeighbors
using Distances: chebyshev
using SpecialFunctions: digamma
using StatsBase: autocor
using KernelDensity

export mutinfo, mutinfo_delaycurve
#####################################################################################
#                                Mutual Information                                 #
#####################################################################################
"""
    mutinfo(k, X1, X2[, ..., Xm]) -> MI

Calculate the mutual information `MI` of the given vectors
`X1, X2, ...`, using `k` nearest-neighbors.

The method follows the second algorithm ``I^{(2)}`` outlined by Kraskov in [1].

## References
[1] : A. Kraskov *et al.*, [Phys. Rev. E **69**, pp 066138 (2004)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.69.066138)

## Performance Notes
This functin gets very slow for large `k`.

See also [`estimate_delay`](@ref) and [`mutinfo_delaycurve`](@ref).
"""
function mutinfo(k, Xm::Vararg{<:AbstractVector,M}) where M
    @assert M > 1
    @assert (size.(Xm,1) .== size(Xm[1],1)) |> prod
    k += 1
    N = size(Xm[1],1)
    invN = 1/N

    d = Dataset(Xm...)
    tree = KDTree(d.data, Chebyshev())

    n_x_m = zeros(M)

    Xm_sp = zeros(Int, N, M)
    Xm_revsp = zeros(Int, N, M)
    for m in 1:M
        Xm_sp[:,m] .= sortperm(Xm[m]; alg=QuickSort)
        Xm_revsp[:,m] .= sortperm(Xm_sp[:,m]; alg=QuickSort)
    end

    I = digamma(k) - (M-1)/k + (M-1)*digamma(N)

    nns = (x = knn(tree, d.data, k)[1]; [ind[1] for ind in x])

    I_itr = zeros(M)
    # Makes more sense computationally to loop over N rather than M
    for i in 1:N
        ϵ = abs.(d[nns[i]] - d[i])./2

        for m in 1:M # this loop takes 8% of time
            hb = lb = Xm_revsp[i,m]
            while abs(Xm[m][Xm_sp[hb,m]] - Xm[m][i]) <= ϵ[m] && hb < N
                hb += 1
            end
            while abs(Xm[m][Xm_sp[lb,m]] - Xm[m][i]) <= ϵ[m] && lb > 1
                lb -= 1
            end
            n_x_m[m] = hb - lb
        end

        I_itr .+= digamma.(n_x_m)
    end

    I_itr .*= invN

    I -= sum(I_itr)

    return max(0, I)
end

"""
    mutinfo_delaycurve(x; maxtau=100, k=1)

Return the [`mutinfo`](@ref) between `x` and itself for delays of `1:maxtau`.
"""
function mutinfo_delaycurve(X::AbstractVector; maxtau=100, k=1)
    I = zeros(maxtau)

    @views for τ in 1:maxtau
        I[τ] = mutinfo(k, X[1:end-τ],X[τ+1:end])
    end

    return I
end

### Fraser's and Swinney's algorithm ###

# Minimum number of points in the margin of a node
const minnodesize = 6

"""
    FSNode

Tree structure for 4^m partitions of a plane:
`low`, `high` are the corners of a node of the partition.
`middle` is the middle point that defines the subpartition.
`children` contains the four sub-partitions.
"""
struct FSNode{T}
    low::Tuple{T,T}
    high::Tuple{T,T}
    middle::Tuple{T,T}
    children::Array{FSNode,1}
end

# Recursive creation of partitions
FSNode(low::Int, high::Int) = FSNode((low,low), (high,high))
function FSNode(low::Tuple{T,T}, high::Tuple{T,T}) where {T}
    half = [div(high[1]-low[1], 2), div(high[2]-low[2],2)]
    middle = (low[1]+half[1], low[2]+half[2])
    node = FSNode(low, high, middle, FSNode[])
    if maximum(half) > minnodesize # or another condition
        append!(node.children,
                [FSNode((low[1],     low[2]),      (middle[1],middle[2])),
                 FSNode((middle[1]+1,low[2]),      (high[1],  middle[2])),
                 FSNode((low[1],     middle[2]+1), (middle[1],high[2])),
                 FSNode((middle[1]+1,middle[2]+1), (high[1],  high[2]))])
    end
    return node
end

"""
    mi_fs(s,τs)

Compute the mutual information between the series `s` and various
images of the same series delayed by `τ` ∈ `τs`, according to
Fraser's and Swinney's algorithm [1].

## References:
[1]: Fraser A.M. & Swinney H.L. "Independent coordinates for strange attractors
from mutual information" *Phys. Rev. A 33*(2), 1986, 1134:1140.
"""
function mi_fs(s::AbstractVector,τs)
    n = length(s)
    nτ = n-maximum(τs)       # trim the series for the maximum delay
    tree = FSNode(1, nτ)      # tree structure of the subpartitions
    perm = sortperm(s[1:nτ])  # order of the original series
    mi = zeros(length(τs))
    for (i,τ) in enumerate(τs)
        indices = sortperm(sortperm(s[τ+1:n])) # rank values of the delayed image of `s`
        mi[i] = _recursive_mi(tree, indices[perm])/nτ - log2(nτ) # eq. (19) of [1]
    end
    return mi
end

# Recursive implementation of equations (20a) and (20b) in [1]
function _recursive_mi(node::FSNode, s)
    # Get view of the `s` within the first dimension of the node
    v = @view s[node.low[1]:node.high[1]]
    if isempty(node.children) # terminal node, no possible substructure
        # Count points in the range of the second dimension
        n = _count(v, node.low[2], node.high[2])
        mi_value = (n==0) ? 0 : n*log2(n)    # eq. (20a) in [1]
    else
        # Split the view in two halves of the first dimension
        # and count point in the two subranges of the second
        n1 = node.middle[1] - node.low[1] + 1
        vhalf = @view v[1:n1]
        n11, n21 = _count(vhalf, node.low[2], node.middle[2], node.high[2])
        vhalf = @view v[n1+1:end]
        n12, n22 = _count(vhalf, node.low[2], node.middle[2], node.high[2])
        n = n11+n12+n21+n22
        if _homogeneous(n11,n12,n21,n22) # no substructure (homogeneous table)
            mi_value = n*log2(n) # eq. (20a) in [1]
        else
            mi_value = 2.0*n     # eq. (20b) in [1]
            for (i, ni) in enumerate((n11,n12,n21,n22))
                (ni > 0) && (mi_value += _recursive_mi(node.children[i], s))
            end
        end
    end
    return mi_value
end

# count values of `s` inside `low:high`
function _count(s, low, high)
    n = 0
    @inbounds for i in s
        (low <= i <= high) && (n+= 1)
    end
    return n
end
# count values of `s` inside `low:middle` and `low:high`
function _count(s, low, middle, high)
    n1 = 0
    n2 = 0
    @inbounds for i in s
        if low <= i
            if i <= middle
                n1 += 1
            elseif i <= high
                n2 += 1
            end
        end
    end
    return n1, n2
end

# homogeneity test - equation (22) in [1]
function _homogeneous(n11,n12,n21,n22)
    n = n11+n12+n21+n22
    nh = n/4
    χ2 = 4/3/n*((n11-nh)^2+(n12-nh)^2+(n21-nh)^2+(n22-nh)^2)
    return χ2 < 1.547
end

function ami_kde(x::AbstractVector{T}, τs::AbstractVector{Int}) where {T}
    n = length(x)
    np = min(2048, prevpow(2, n))
    dx = kde(x, npoints=np)
    #
    tails_width = 4.0*KernelDensity.default_bandwidth(x)
    tp = round(Int, np*tails_width/(dx.x[end]-dx.x[1]))
    #
    # dx.density .*= step(dx.x)
    px = @view dx.density[tp+1:end-tp]
    hx = -sum(px.*log2.(px))*step(dx.x)
    ami = zeros(T, length(τs))
    #
    np2 = min(256, 2^floor(Int, log2(n)))
    tp = div(tp*np2,np)
    for (i,τ) in enumerate(τs)
        hxy = zero(T)
        dxy = kde((x[1:n-τ], x[τ+1:n]), npoints=(np2,np2))
        # dxy.density .*= (step(dxy.x) * step(dxy.y))
        pxy = @view dxy.density[tp+1:end-tp, tp+1:end-tp]
        for p in pxy
            if p ≉ 0
                hxy -= p*log2(p)
            end
        end
        hxy *= (step(dxy.x)*step(dxy.y))
        mi[i] = 2hx - hxy
    end
    return ami
end


text:
# function Phi(u::Float64, v::Float64)
#     BIG = 10000000
#     if v>BIG
#         return (u,1,0)
#     else
#         sq = sqrt(u^2+v^2)
#         p = u+v-sq
#         J_u = 1 - u/sq
#         J_v = 1 - v/sq
#         return (p,J_u,J_v)
# end

# function Phi(u::Point{d},v::Point{d}) where d
#     rr = [Phi(u[i], v[i]) for i=1:d]
#     p = [r[1] for r in rr]
#     J_u = [r[2] for r in rr]
#     J_v = [r[3] for r in rr]
#     return SVector(p...), SDiagonal(J_u...), SDiagonal(J_v...)
# end

function Phi(u::Point{d},v::Point{d}) where d
    sq = sqrt.(u.^2+v.^2)
    p = u+v-sq
    J_u = 1 - u./sq
    J_v = 1 - v./sq
    return p, SDiagonal(J_u), SDiagonal(J_v)
end

function PhiPhi(f::Point{d}, x::Point{d}, a::Point{d}, b::Point{d}) where d
    y, y_f, y_x = Phi(f,x-a)
    z, z_y, z_x = Phi(-y,b-x)
    return z, (-z_y*y_f), (-z_y*y_x - z_x)
end

function PhiPhi0(f::Point{d}, x::Point{d}, a::Point{d}, b::Point{d}) where d
    y = Phi(f,x-a)[1]
    z = Phi(-y,b-x)[1]
    return z
end

function PhiPhi(f::Point{d}, D::SMatrix{d,d,Float64,q}, x::Point{d}, a::Point{d}, b::Point{d}) where d where q
    z, z_f, z_x = PhiPhi(f,x,a,b)
    z, z_f*D + z_x
end

function PhiPhi!(F::Vector{Vector{Point{d}}},
                 X::Vector{Vector{Point{d}}},
                 A::Vector{Vector{Point{d}}},
                 B::Vector{Vector{Point{d}}},
                 D::Vector{v},
                 J::Matrix{w}) where v<:AbstractVector{SMatrix{d,d,Float64,q}} where w<:AbstractVector{SMatrix{d,d,Float64,q}} where d where q

    n_m, n_M = size(J)
    N = length(F[1])
    for i=1:n_m
        for n=1:N
            f = F[i][n]
            x = X[i][n]
            a = A[i][n]
            b = B[i][n]
            z, z_f, z_x = PhiPhi(f,x,a,b)
            F[i][n] = z
            D[i][n] = z_f*D[i][n] + z_x
            for j=1:n_M
                J[i,j][n] = z_f*J[i,j][n]
            end
        end
    end
end

function PhiPhi!(F::Vector{Vector{Point{d}}},X::Vector{Vector{Point{d}}},A::Vector{Vector{Point{d}}},B::Vector{Vector{Point{d}}},D::Vector{v}) where v<:AbstractVector{SMatrix{d,d,Float64,q}} where d where q

    n_m = size(D,1)
    N = length(F[1])
    for i=1:n_m
        for n=1:N
            f = F[i][n]
            x = X[i][n]
            a = A[i][n]
            b = B[i][n]
            z, z_f, z_x = PhiPhi(f,x,a,b)
            F[i][n] = z
            D[i][n] = z_f*D[i][n] + z_x
        end
    end
end

function PhiPhi(F::Vector{Vector{Point{d}}},X::Vector{Vector{Point{d}}},A::Vector{Vector{Point{d}}},B::Vector{Vector{Point{d}}},D::Vector{Vector{SMatrix{d,d,Float64,q}}}, J::Matrix{Vector{SMatrix{d,d,Float64,q}}}) where d where q
    FF = deepcopy(F)
    DD = deepcopy(D)
    JJ = deepcopy(J)
    PhiPhi!(FF,X,A,B,DD,JJ)
    return FF,DD,JJ
end


text:
function langevin_inv(y::Real)
  return (3*y - y/5 * (6*y^2 + y^4 - 2*y^6))/(1 - y^2);
end

function free_energy_langevin(k::Real, T::Real, N::Int, γ::Real)
  linvγ = langevin_inv(γ);
  return N*k*T*( log( linvγ / (4*pi*sinh(linvγ)) ) + γ*linvγ );
end

function force_langevin(k::Real, T::Real, b::Real, γ::Real)
  return k*T/b * langevin_inv(γ);
end

function forces(As::Vector, γs::Vector, N::Int)
  n = length(As)
  grads = zeros(n);
  grads[1] = (As[2]-As[1])/(γs[2]-γs[1]);
  grads[end] = (As[n]-As[n-1])/(γs[n]-γs[n-1]);
  for i=2:(n-1)
    grads[i] = (As[i+1] - As[i-1])/(γs[i+1] - γs[i-1]);
  end
  return grads / N;
end

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
STORE_URL = nothing

"""
cloudstore_init(url)

Set cloud storage url initialize

# Example

```julia
cloudstore_init("https://firebasestorage.googleapis.com/v0/b/[PROJECT_ID].appspot.com/o")
```
"""
function cloudstore_init(url)
    global STORE_URL = url
    println("Cloud Storage set:", STORE_URL)
end

"""
cloudstore_sendfile(file2download = "")

Send a file to cloud storage to store

# Example

```julia
init("[PROJECT_SDK].json")
cloudstore_init("https://firebasestorage.googleapis.com/v0/b/[PROJECT_ID].appspot.com/o")
cloudstore_sendfile("/audio/data.mp3")
```
"""
function cloudstore_sendfile(file2upload="/test/onedance.mp3")
    # println(STORE_URL)
    string  = readchomp(`file --mime-type -b $(file2upload[2:end])`)
    file_binary = open(file2upload[2:end], "r")
    s = read(file_binary, String)
    header = Dict("Content-Type" => "$string", auth_header())
    project_id = projectid()
    spliturl= split(file2upload[2:end],"/")
    file2upload = join(spliturl, "%2F")
    final_url = "$STORE_URL$file2upload"
    body = s
    println("FINAL_URL:", final_url)
    http_response = HTTP.post(
        final_url,
        header,
        body
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)

end

# STORE_URL = "https://firebasestorage.googleapis.com/v0/b/fir-jl-457eb.appspot.com/o/"

"""
cloudstore_get(file2download = "")

Get a file cloud storage

# Example
```julia
init("[PROJECT_SDK].json")
cloudstore_init("https://firebasestorage.googleapis.com/v0/b/[PROJECT_ID].appspot.com/o")
cloudstore_get("/audio/data.mp3")
cloudstore_get("/data.mp3")
```
"""
function cloudstore_get(file2download="242160__xserra__cello-phrase.wav")
    temp = file2download
    spliturl= split(file2download[2:end],"/")
    file2download = join(spliturl, "%2F")
    final_url = "$STORE_URL$file2download?alt=media&token=878c7715-61cf-46b5-a9c0-a790300d8c43"
    println("Final URL:",final_url)
    res = HTTP.get(final_url)
    dir = "$(temp[2:end])"
    dirtocreate=join(spliturl[1:end-1],"/")
    mkpath(dirtocreate)
    # mkpath(dir)
    out =open("$dir","w")
    write(out,res.body)
    close(out)
end


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavcodec/vdpau.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_alloc_vdpaucontext,
    av_vdpau_hwaccel_get_render2,
    av_vdpau_hwaccel_set_render2,
    av_vdpau_bind_context,
    av_vdpau_alloc_context,
    av_vdpau_get_profile

function av_alloc_vdpaucontext()
    ccall((:av_alloc_vdpaucontext,libavcodec),Ptr{AVVDPAUContext},())
end

function av_vdpau_hwaccel_get_render2(arg1)
    ccall((:av_vdpau_hwaccel_get_render2,libavcodec),AVVDPAU_Render2,(Ptr{AVVDPAUContext},),arg1)
end

function av_vdpau_hwaccel_set_render2(arg1,arg2::AVVDPAU_Render2)
    ccall((:av_vdpau_hwaccel_set_render2,libavcodec),Void,(Ptr{AVVDPAUContext},AVVDPAU_Render2),arg1,arg2)
end

function av_vdpau_bind_context(avctx,device::VdpDevice,get_proc_address,flags::Integer)
    ccall((:av_vdpau_bind_context,libavcodec),Cint,(Ptr{AVCodecContext},VdpDevice,Ptr{VdpGetProcAddress},UInt32),avctx,device,get_proc_address,flags)
end

function av_vdpau_alloc_context()
    ccall((:av_vdpau_alloc_context,libavcodec),Ptr{AVVDPAUContext},())
end

function av_vdpau_get_profile(avctx,profile)
    ccall((:av_vdpau_get_profile,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{VdpDecoderProfile}),avctx,profile)
end