text:
"""
Collection of anisotropy models.

* [`IsotropicModel`](@ref): isotropic model
* [`ConstantBetaModel`](@ref): constant beta model
* [`RSBetaModel`](@ref): flexible beta model from Read & Steger
* [`beta`](@ref): orbital anisotropy profile
"""
module Anisotropy

import Slomo.Utils: B_inc, gamma
import Slomo.Models: Model
export IsotropicModel, ConstantBetaModel, RSBetaModel, beta

abstract type AnisotropyModel <: Model end

"""
    beta(model::AnisotropyModel, r)

Orbital anisotropy as a function of radius.  The orbital anisotropy is defined as

```math
\\beta = 1 - \\sigma_\\theta^2 / \\sigma_r^2
```
"""
function beta(model::AnisotropyModel, r)
    throw(NotImplemented("no defined beta profile"))
end

"""
    g_jean(model::AnisotropyModel, r)

Integrating factor when solving for \$\\sigma_r^2\$
"""
function g_jeans(model::AnisotropyModel, r)
    exp.(2.0 * integrate(x -> beta(model, x) / x, 1.0, r))
end

"""
    K_jeans(model::AnisotropyModel, r, R)

Jeans projection kernel for an anisotropy model.
"""
function K_jeans(model::AnisotropyModel, r, R)
    throw(NotImplemented("no defined Jeans Kernel"))
end

"""
    IsotropicModel()

Parameter-less model representing an isotropic system (i.e. beta = 0).
"""
struct IsotropicModel <: AnisotropyModel end

beta(model::IsotropicModel, r) = zeros(size(r))

K_jeans(model::IsotropicModel, r, R) = begin
    @. sqrt(1.0 - 1.0 / (r / R) ^ 2)
end

g_jeans(model::IsotropicModel, r) = ones(size(r))

"""
    ConstantBetaModel(beta)

Model that assumes beta is constant with galactocentric radius.
"""
struct ConstantBetaModel <: AnisotropyModel
    beta::Float64
end
ConstantBetaModel() = ConstantBetaModel(0.0)

beta(model::ConstantBetaModel, r) = model.beta .* ones(size(r))

g_jeans(model::ConstantBetaModel, r) = r .^ (2 * model.beta)

K_jeans(model::ConstantBetaModel, r, R) = begin
    Γ = gamma
    β = model.beta
    # nudge half-integral beta
    if 2.0 * β - floor(2.0 * β) == 0.0
        β = β + rand([1, -1]) * eps()^0.5
    end
    u = r ./ R
    um2 = u .^ -2.0
    term1 = (1.5 - β) * √(π) * Γ(β - 0.5) / Γ(β)
    term2 = β * B_inc(β + 0.5, 0.5, um2)
    term3 = -B_inc(β - 0.5, 0.5, um2)
    return 0.5 * u .^ (2.0 * β - 1.0) .* (term1 .+ term2 .+ term3)
end

"""
    RSBetaModel(beta0, betaInf, rbeta, nbeta)

Flexible velocity anisotropy model from Read & Steger 2017 (eq. 9)

```math
\\beta(r) = \\beta_0 + (\\beta_\\infty - \\beta_0) / (1 + (r_\\beta / r))^{n_\\beta}
```

* `beta0`: inner asymptotic anisotropy
* `betaInf`: outer asymptotic anisotropy
* `rbeta`: transition radius
* `nbeta`: transition sharpness, higher is a faster transition

`beta0, betaInf, nbeta = 0, 1, 2` corresponds to the Osipkov-Merritt profile

`beta0, betaInf, nbeta = 0, 0.5, 1` corresponds to the Mamon-Lokas profile
"""
struct RSBetaModel <: AnisotropyModel
    beta0::Float64
    betaInf::Float64
    rbeta::Float64
    nbeta::Float64
end

RSBetaModel() = RSBetaModel(0.0, 0.5, 10.0, 1.0)

beta(model::RSBetaModel, r) = begin
    @. model.beta0 + (model.betaInf - model.beta0) / (1.0 + (model.rbeta / r) ^ model.nbeta)
end

g_jeans(model::RSBetaModel, r) = begin
    n = model.nbeta
    r_β = model.rbeta
    β_inf = model.betaInf
    Δβ = model.betaInf - model.beta0
    return @. r ^ (2.0 * β_inf) * ((r_β / r)^n + 1.0) ^ (2.0 * Δβ / n)
end

end


text:
include("libswscale_h.jl")

include("swscale.jl")
#include("version.jl")

ext: jl


text:
abstract AbstractRunner

### Factory function
runner(p::AbstractPolicy,args...;keyargs...) = _runner(p,args...;keyargs...)


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    @test f(0,0) == 0
    @test f(1,0) == 2
    @test f(0,1) == 3
    @test f(1,1) == 5
end


text:
@everywhere using DifferentialEquations

prob = ODEProblem((u,p,t)->1.01u,0.5,(0.0,1.0))

@everywhere function prob_func(prob,i,repeat)
  ODEProblem(prob.f,rand()*prob.u0,prob.tspan)
end

monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
sim = solve(monte_prob,Tsit5(),trajectories=100)

@info "Simulation finished!"


text:
#!/usr/bin/env julia

for c in ARGS[1] print(c == '\n' ? c : 'z' - (c - 'a')) end


text:
# # [Modeling with JuMP](@id modeling_with_jump)
#
# This page shows a minimal example of `PowerSystems.jl` used to develop and Economic Dispatch
# model. The code  shows the stages to develop modeling code
#
# 1. Make the data set from power flow and time series data,
# 2. Serialize the data,
# 3. Pass the data and algorithm to the model.
#
# One of the main uses of ``PowerSystems.jl` is not having re-run the data generation for every
# model execution. The model code shows an example of populating the constraints and cost
# functions using accessor functions inside the model function. The example concludes by
# reading the data created earlier and passing the algorithm with the data.

using PowerSystems
const PSY = PowerSystems
using JuMP
using Ipopt

DATA_DIR = "../../../data" #hide
system_data = System(joinpath(DATA_DIR, "matpower/case5_re.m"))
add_time_series!(system_data, joinpath(DATA_DIR,"forecasts/5bus_ts/timeseries_pointers_da.json"))
to_json(system_data, "system_data.json")

function ed_model(system::System, optimizer)
    ed_m = Model(optimizer)
    time_periods = 24
    thermal_gens_names = get_name.(get_components(ThermalStandard, system))
    @variable(ed_m, pg[g in thermal_gens_names, t in time_periods] >= 0)

    for g in get_components(ThermalStandard, system), t in time_periods
        name = get_name(g)
        @constraint(ed_m, pg[name, t] >= get_active_power_limits(g).min)
        @constraint(ed_m, pg[name, t] <= get_active_power_limits(g).max)
    end

    net_load = zeros(time_periods)
    for g in get_components(RenewableGen, system)
        net_load -= get_time_series_values(SingleTimeSeries, g, "max_active_power")
    end

    for g in get_components(StaticLoad, system)
        net_load += get_time_series_values(SingleTimeSeries, g, "max_active_power")
    end

    for t in time_periods
        @constraint(ed_m, sum(pg[g, t] for g in thermal_gens_names) == net_load[t])
    end

    @objective(ed_m, Min, sum(
            pg[get_name(g), t]^2*get_cost(get_variable(get_operation_cost(g)))[1] +
            pg[get_name(g), t]*get_cost(get_variable(get_operation_cost(g)))[2]
            for g in get_components(ThermalGen, system), t in time_periods
                )
            )

    return optimize!(ed_m)
end

system_data = System("system_data.json")
results = ed_model(system_data, Ipopt.Optimizer)


text:
import MultifidelityTrajectories
MT = MultifidelityTrajectories
# const MT = MultifidelityTrajectories
import PyPlot
PP = PyPlot

##### TEST #####
nwings = 2
gamma = 0.0 * pi/180

# include("test_aerodynamics.jl")
# include("test_dynamics.jl")
# include("test_integration.jl")
include("test_trim.jl")

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
module GGPlot

export
    ggplot,
    aes,
    geom_point,
    geom_line

using DataFrames
using GR

abstract AbstractLayer
abstract AbstractGeom <: AbstractLayer

include("aes.jl")
include("plot.jl")
include("geom_point.jl")
include("geom_line.jl")
include("gr_backend.jl")

end # module


text:
abstract type InitialMassDistribution{N} end

struct MassDistributionData{N}
    "Mass of the fluid particles."
    m::Float
    "Positions of the fluid particles."
    positions::Vector{SVectorF{N}}
    "Fluid particle mass densities."
    mass_densities::Vector{Float}

    function MassDistributionData(
        m::Number,
        positions::AbstractVector{SVectorF{N}},
        mass_densities::AbstractVector{Float},
    ) where {N}
        new{N}(m, positions, mass_densities)
    end
end

"""
    (::InitialMassDistribution{N})() -> MassDistributionData{N}

Compute the mass distribution data corresponding to this initial mass
distribition.
"""
function (::InitialMassDistribution) end

include("mass/empty.jl")
include("mass/uniform.jl")
include("mass/piecewise_uniform.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
function plane_strain_transformation_stress(lambda, mu, theta)
    pt = (lambda + 2mu / 3) * theta
    transfstress = [pt, pt, 0.0]
    return transfstress
end

function bulk_transformation_rhs(
    transfstress,
    basis,
    quad,
    dim,
    jac,
    detjac,
    vectosymmconverter,
)

    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)

    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jac)
        NK = sum([
            make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim
        ])
        rhs .+= NK' * transfstress * detjac * w
    end
    return rhs
end

function coherent_interface_transformation_linear_form(
    transfstress,
    basis,
    quad,
    normals,
    dim,
    scalearea,
    vectosymmconverter,
)

    numqp = length(quad)
    @assert size(normals) == (dim, numqp)
    @assert length(scalearea) == numqp

    nf = number_of_basis_functions(basis)
    totalndofs = dim * nf
    rhs = zeros(totalndofs)

    for qpidx = 1:numqp
        p, w = quad[qpidx]
        vals = basis(p)
        NI = interpolation_matrix(vals, dim)

        normal = normals[:, qpidx]
        NK = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])

        rhs .+= NI' * NK * transfstress * scalearea[qpidx] * w
    end
    return rhs
end

function component_interface_transformation_linear_form(
    transfstress,
    basis,
    quad,
    components,
    normals,
    dim,
    scalearea,
    vectosymmconverter,
)

    numqp = length(quad)
    @assert size(components) == (dim, numqp)
    @assert size(normals) == (dim, numqp)
    @assert length(scalearea) == numqp
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)

    for qpidx = 1:numqp
        p, w = quad[qpidx]

        vals = basis(p)
        normal = normals[:, qpidx]
        component = components[:, qpidx]

        projector = component * component'
        NK = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])
        NI = interpolation_matrix(vals, dim)

        rhs .+= NI' * projector * NK * transfstress * scalearea[qpidx] * w
    end
    return rhs
end


text:
"""
    init(ctx::Context)
Initialize GLFW backend. Do initialize a IMGUI context before calling this function.
"""
function init(ctx::Context)
    ctx.Time = 0.0

    # setup back-end capabilities flags
    io::Ptr{ImGuiIO} = igGetIO()
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_HasMouseCursors
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_HasSetMousePos
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_PlatformHasViewports
    # TODO: pending glfw 3.4
    # io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_HasMouseHoveredViewport
    io.BackendPlatformName = pointer(GLFW_BACKEND_PLATFORM_NAME)

    # store the contextual object reference to IMGUI
    io.UserData = pointer_from_objref(ctx)

    # keyboard mapping
    c_set!(io.KeyMap, ImGuiKey_Tab, GLFW_KEY_TAB)
    c_set!(io.KeyMap, ImGuiKey_LeftArrow, GLFW_KEY_LEFT)
    c_set!(io.KeyMap, ImGuiKey_RightArrow, GLFW_KEY_RIGHT)
    c_set!(io.KeyMap, ImGuiKey_UpArrow, GLFW_KEY_UP)
    c_set!(io.KeyMap, ImGuiKey_DownArrow, GLFW_KEY_DOWN)
    c_set!(io.KeyMap, ImGuiKey_PageUp, GLFW_KEY_PAGE_UP)
    c_set!(io.KeyMap, ImGuiKey_PageDown, GLFW_KEY_PAGE_DOWN)
    c_set!(io.KeyMap, ImGuiKey_Home, GLFW_KEY_HOME)
    c_set!(io.KeyMap, ImGuiKey_End, GLFW_KEY_END)
    c_set!(io.KeyMap, ImGuiKey_Insert, GLFW_KEY_INSERT)
    c_set!(io.KeyMap, ImGuiKey_Delete, GLFW_KEY_DELETE)
    c_set!(io.KeyMap, ImGuiKey_Backspace, GLFW_KEY_BACKSPACE)
    c_set!(io.KeyMap, ImGuiKey_Space, GLFW_KEY_SPACE)
    c_set!(io.KeyMap, ImGuiKey_Enter, GLFW_KEY_ENTER)
    c_set!(io.KeyMap, ImGuiKey_Escape, GLFW_KEY_ESCAPE)
    c_set!(io.KeyMap, ImGuiKey_KeyPadEnter, GLFW_KEY_KP_ENTER)
    c_set!(io.KeyMap, ImGuiKey_A, GLFW_KEY_A)
    c_set!(io.KeyMap, ImGuiKey_C, GLFW_KEY_C)
    c_set!(io.KeyMap, ImGuiKey_V, GLFW_KEY_V)
    c_set!(io.KeyMap, ImGuiKey_X, GLFW_KEY_X)
    c_set!(io.KeyMap, ImGuiKey_Y, GLFW_KEY_Y)
    c_set!(io.KeyMap, ImGuiKey_Z, GLFW_KEY_Z)

    # set clipboard
    io.GetClipboardTextFn = GLFW_GET_CLIPBOARD_TEXT_FUNCPTR[]
    io.SetClipboardTextFn = GLFW_SET_CLIPBOARD_TEXT_FUNCPTR[]
    io.ClipboardUserData = Ptr{Cvoid}(ctx.Window)

    # create mouse cursors
    ctx.MouseCursors[ImGuiMouseCursor_Arrow+1] = glfwCreateStandardCursor(GLFW_ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_TextInput+1] = glfwCreateStandardCursor(GLFW_IBEAM_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeNS+1] = glfwCreateStandardCursor(GLFW_VRESIZE_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeEW+1] = glfwCreateStandardCursor(GLFW_HRESIZE_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_Hand+1] = glfwCreateStandardCursor(GLFW_HAND_CURSOR)

    # prepare for GLFW 3.4+
    # ctx.MouseCursors[ImGuiMouseCursor_ResizeAll+1] = glfwCreateStandardCursor(GLFW_RESIZE_ALL_CURSOR)
    # ctx.MouseCursors[ImGuiMouseCursor_ResizeNESW+1] = glfwCreateStandardCursor(GLFW_RESIZE_NESW_CURSOR)
    # ctx.MouseCursors[ImGuiMouseCursor_ResizeNWSE+1] = glfwCreateStandardCursor(GLFW_RESIZE_NWSE_CURSOR)
    # ctx.MouseCursors[ImGuiMouseCursor_NotAllowed+1] = glfwCreateStandardCursor(GLFW_NOT_ALLOWED_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeAll+1] = glfwCreateStandardCursor(GLFW_ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeNESW+1] = glfwCreateStandardCursor(GLFW_ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeNWSE+1] = glfwCreateStandardCursor(GLFW_ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_NotAllowed+1] = glfwCreateStandardCursor(GLFW_ARROW_CURSOR)

    # chain GLFW callbacks
    ctx.PrevUserCallbackMousebutton = C_NULL
    ctx.PrevUserCallbackScroll = C_NULL
    ctx.PrevUserCallbackKey = C_NULL
    ctx.PrevUserCallbackChar = C_NULL
    ctx.PrevUserCallbackMonitor = C_NULL
    if ctx.InstalledCallbacks
        ctx.PrevUserCallbackMousebutton = glfwSetMouseButtonCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_MouseButtonCallback, Cvoid, (Ptr{GLFWwindow}, Cint, Cint, Cint)))
        ctx.PrevUserCallbackScroll = glfwSetScrollCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_ScrollCallback, Cvoid, (Ptr{GLFWwindow}, Cdouble, Cdouble)))
        ctx.PrevUserCallbackKey = glfwSetKeyCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_KeyCallback, Cvoid, (Ptr{GLFWwindow}, Cint, Cint, Cint, Cint)))
        ctx.PrevUserCallbackChar = glfwSetCharCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_CharCallback, Cvoid, (Ptr{GLFWwindow}, Cuint)))
        ctx.PrevUserCallbackMonitor = glfwSetMonitorCallback(@cfunction(ImGui_ImplGlfw_MonitorCallback, Cvoid, (Ptr{GLFWmonitor}, Cint)))
    end

    # update monitors the first time
    ImGui_ImplGlfw_UpdateMonitors(ctx)
    glfwSetMonitorCallback(@cfunction(ImGui_ImplGlfw_MonitorCallback, Cvoid, (Ptr{GLFWmonitor}, Cint)))

    # the mouse update function expect PlatformHandle to be filled for the main viewport
    main_viewport::Ptr{ImGuiViewport} = igGetMainViewport()
    main_viewport.PlatformHandle = Ptr{Cvoid}(ctx.Window)
    if Sys.iswindows()
        main_viewport.PlatformHandleRaw = ccall((:glfwGetWin32Window, libglfw), Ptr{Cvoid}, (Ptr{Cvoid},), ctx.Window)
    end

    if unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
        ImGui_ImplGlfw_InitPlatformInterface(ctx)
    end

    return true
end

"""
    shutdown((ctx::Context)
Clean up resources.
"""
function shutdown(ctx::Context)
    ImGui_ImplGlfw_ShutdownPlatformInterface(ctx)

    if ctx.InstalledCallbacks
        glfwSetMouseButtonCallback(ctx.Window, ctx.PrevUserCallbackMousebutton)
        glfwSetScrollCallback(ctx.Window, ctx.PrevUserCallbackScroll)
        glfwSetKeyCallback(ctx.Window, ctx.PrevUserCallbackKey)
        glfwSetCharCallback(ctx.Window, ctx.PrevUserCallbackChar)
        ctx.InstalledCallbacks = false
    end

    for i = 1:length(ctx.MouseCursors)
        glfwDestroyCursor(ctx.MouseCursors[i])
        ctx.MouseCursors[i] = C_NULL
    end

    ctx.ClientApi = GlfwClientApi_Unknown

    return true
end

function ImGui_ImplGlfw_UpdateMousePosAndButtons(ctx::Context)
    # update buttons
    io::Ptr{ImGuiIO} = igGetIO()
    for n = 0:length(ctx.MouseJustPressed)-1
        # if a mouse press event came, always pass it as "mouse held this frame",
        # so we don't miss click-release events that are shorter than 1 frame.
        is_down = ctx.MouseJustPressed[n+1] || glfwGetMouseButton(ctx.Window, n) == GLFW_PRESS
        c_set!(io.MouseDown, n, is_down)
        ctx.MouseJustPressed[n+1] = false
    end

    # update mouse position
    mouse_pos_backup = unsafe_load(io.MousePos)
    FLT_MAX = igGET_FLT_MAX()
    io.MousePos = ImVec2(-FLT_MAX, -FLT_MAX)
    io.MouseHoveredViewport = 0
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    vp = unsafe_load(platform_io.Viewports)
    viewport_ptrs = unsafe_wrap(Vector{Ptr{ImGuiViewport}}, vp.Data, vp.Size)
    for viewport in viewport_ptrs
        window = unsafe_load(viewport.PlatformHandle)
        @assert window != C_NULL

        if glfwGetWindowAttrib(window, GLFW_FOCUSED) != 0
            if unsafe_load(io.WantSetMousePos)
                x = mouse_pos_backup.x - unsafe_load(viewport.Pos.x)
                y = mouse_pos_backup.y - unsafe_load(viewport.Pos.y)
                glfwSetCursorPos(window, Cdouble(x), Cdouble(y))
            else
                mouse_x, mouse_y = Cdouble(0), Cdouble(0)
                @c glfwGetCursorPos(window, &mouse_x, &mouse_y)
                if unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
                    # Multi-viewport mode: mouse position in OS absolute coordinates (io.MousePos is (0,0) when the mouse is on the upper-left of the primary monitor)
                    window_x, window_y = Cint(0), Cint(0)
                    @c glfwGetWindowPos(window, &window_x, &window_y)
                    io.MousePos = ImVec2(Cfloat(mouse_x + window_x), Cfloat(mouse_y + window_y))
                else
                    # Single viewport mode: mouse position in client window coordinates (io.MousePos is (0,0) when the mouse is on the upper-left corner of the app window)
                    io.MousePos = ImVec2(Cfloat(mouse_x), Cfloat(mouse_y))
                end
            end

        end
        for n = 0:length(ctx.MouseJustPressed)-1
            c_set!(io.MouseDown, n, c_get(io.MouseDown, n) || glfwGetMouseButton(window, n) == GLFW_PRESS)
        end
    end
    # TODO: pending glfw 3.4 GLFW_HAS_MOUSE_PASSTHROUGH
    # TODO: _WIN32
end

function ImGui_ImplGlfw_UpdateMouseCursor(ctx::Context)
    io::Ptr{ImGuiIO} = igGetIO()
    if (unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_NoMouseCursorChange == ImGuiConfigFlags_NoMouseCursorChange) ||
        glfwGetInputMode(ctx.Window, GLFW_CURSOR) == GLFW_CURSOR_DISABLED
        return nothing
    end

    imgui_cursor = igGetMouseCursor()
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    vp = unsafe_load(platform_io.Viewports)
    viewport_ptrs = unsafe_wrap(Vector{Ptr{ImGuiViewport}}, vp.Data, vp.Size)
    for viewport in viewport_ptrs
        window = unsafe_load(viewport.PlatformHandle)
        @assert window != C_NULL
        if imgui_cursor == ImGuiMouseCursor_None || unsafe_load(io.MouseDrawCursor)
            # hide OS mouse cursor if imgui is drawing it or if it wants no cursor
            glfwSetInputMode(window, GLFW_CURSOR, GLFW_CURSOR_HIDDEN)
        else
            # show OS mouse cursor
            cursor = ctx.MouseCursors[imgui_cursor+1]
            glfwSetCursor(window, cursor != C_NULL ? ctx.MouseCursors[imgui_cursor+1] : ctx.MouseCursors[ImGuiMouseCursor_Arrow+1])
            glfwSetInputMode(window, GLFW_CURSOR, GLFW_CURSOR_NORMAL)
        end
    end

    return nothing
end

function ImGui_ImplGlfw_UpdateMonitors(ctx::Context)
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    monitors_count = Cint(0)
	ptr = @c glfwGetMonitors(&monitors_count)
	glfw_monitors = unsafe_wrap(Array, ptr, monitors_count)
    monitors_ptr::Ptr{ImGuiPlatformMonitor} = Libc.malloc(monitors_count * sizeof(ImGuiPlatformMonitor))
    for i = 1:monitors_count
        glfw_monitor = glfw_monitors[i]
        mptr::Ptr{ImGuiPlatformMonitor} = monitors_ptr + (i-1) * sizeof(ImGuiPlatformMonitor)

        x, y = Cint(0), Cint(0)
        @c glfwGetMonitorPos(glfw_monitor, &x, &y)
        vid_mode = unsafe_load(glfwGetVideoMode(glfw_monitor))
        mptr.MainPos = ImVec2(x, y)
        mptr.MainSize = ImVec2(vid_mode.width, vid_mode.height)
        mptr.WorkPos = ImVec2(x, y)
        mptr.WorkSize = ImVec2(vid_mode.width, vid_mode.height)

        w, h = Cint(0), Cint(0)
        @c glfwGetMonitorWorkarea(glfw_monitors[i], &x, &y, &w, &h)
        if w > 0 && h > 0
            mptr.WorkPos = ImVec2(Cfloat(x), Cfloat(y))
            mptr.WorkSize = ImVec2(Cfloat(w), Cfloat(h))
        end

        x_scale, y_scale = Cfloat(0), Cfloat(0)
        @c glfwGetMonitorContentScale(glfw_monitor, &x_scale, &y_scale)
        mptr.DpiScale = x_scale
    end

    platform_io.Monitors = ImVector_ImGuiPlatformMonitor(monitors_count, monitors_count, monitors_ptr)
    ctx.WantUpdateMonitors = false

    return nothing
end

function new_frame(ctx::Context)
    io::Ptr{ImGuiIO} = igGetIO()
    @assert ImFontAtlas_IsBuilt(unsafe_load(io.Fonts)) "Font atlas not built! It is generally built by the renderer back-end. Missing call to renderer _NewFrame() function? e.g. ImGui_ImplOpenGL3_NewFrame()."

    # setup display size (every frame to accommodate for window resizing)
    w, h = Cint(0), Cint(0)
    display_w, display_h = Cint(0), Cint(0)
    @c glfwGetWindowSize(ctx.Window, &w, &h)
    @c glfwGetFramebufferSize(ctx.Window, &display_w, &display_h)
    io.DisplaySize = ImVec2(Cfloat(w), Cfloat(h))
    if w > 0 && h > 0
        w_scale = Cfloat(display_w / w)
        h_scale = Cfloat(display_h / h)
        io.DisplayFramebufferScale = ImVec2(w_scale, h_scale)
    end

    ctx.WantUpdateMonitors && ImGui_ImplGlfw_UpdateMonitors(ctx)

    # setup time step
    current_time = glfwGetTime()
    io.DeltaTime = ctx.Time > 0.0 ? Cfloat(current_time - ctx.Time) : Cfloat(1.0/60.0)
    ctx.Time = current_time

    ImGui_ImplGlfw_UpdateMousePosAndButtons(ctx)
    ImGui_ImplGlfw_UpdateMouseCursor(ctx)

    # TODO: Gamepad navigation mapping
    # ImGui_ImplGlfw_UpdateGamepads()
end


text:
@inline function check_positive_maturity(t::AbstractFloat)
    @assert t >= 0 "The effective rate factor for negative maturity is undefined."
end

"""
    erf(compounding, rate, yearfraction)

The *Effective Rate Factor* is the result of an investment
of 1 unit of a currency for a fixed interest rate `rate`
and period `yearfraction`.

Given an Effective Rate Factor `erf`, the associated discount
factor can be calculated as `1.0/erf`.
"""
function erf
end

function erf(::Continuous, r::AbstractFloat, t::AbstractFloat)
    check_positive_maturity(t)
    return t == 0.0 ? 1.0 : exp(r*t)
end

function erf(::Simple, r::AbstractFloat, t::AbstractFloat)
    check_positive_maturity(t)
    return t == 0.0 ? 1.0 : 1.0 + r*t
end

function erf(::Exponential, r::AbstractFloat, t::AbstractFloat)
    check_positive_maturity(t)
    return t == 0.0 ? 1.0 : (1.0 + r)^t
end

"""
    discountfactor(compounding, rate, yearfraction)

The *Discount Factor* is the amount of currency one must
hold now to receive 1 unit of currency after a period
given by `yearfraction`, when investing using
and interest rate given by `rate`.
"""
function discountfactor
end

discountfactor(c::Compounding, r::AbstractFloat, t::AbstractFloat) = 1.0 / erf(c, r, t)

"""
    erf_to_rate(compounding, erf, yearfraction)

Converts an Effective Rate Factor `erf` to an Interest Rate.
"""
function erf_to_rate
end

function erf_to_rate(::Continuous, factor::AbstractFloat, t::AbstractFloat)
    check_positive_maturity(t)
    return log(factor) / t
end

function erf_to_rate(::Simple, factor::AbstractFloat, t::AbstractFloat)
    check_positive_maturity(t)
    return (factor-1.0)/t
end

function erf_to_rate(::Exponential, factor::AbstractFloat, t::AbstractFloat)
    check_positive_maturity(t)
    return factor^(1.0/t) - 1.0
end

"""
    discountfactor_to_rate(compounding, discountfactor, yearfraction)

Converts a Discount Factor `discountfactor` to an Interest Rate.
"""
function discountfactor_to_rate(c::Compounding, discountfactor::AbstractFloat, t::AbstractFloat)
    return erf_to_rate(c, 1.0 / discountfactor, t)
end

# allows `erf(Exponential, 0.1, 1.0)` syntax.
for fun in (:erf, :discountfactor, :erf_to_rate, :discountfactor_to_rate)
    @eval begin
        function ($fun)(::Type{T}, a::AbstractFloat, b::AbstractFloat) where {T<:Compounding}
            return ($fun)(T(), a, b)
        end
    end
end


text:
using Distributions
using RmathDist

distlist = [
            Weibull(1,1),
            Weibull(1,1e10),
            Weibull(1,1e-10),
            Weibull(1e3,1),
            Weibull(1e3,1e10),
            Weibull(1e3,1e-10),
            Weibull(1e-2,1),
            Weibull(1e-2,1e10),
            # Weibull(1e-2,1e-10) rmath not so good here
            ]

for d in distlist
    for x in [logspace(-300,300,1_000_000)]

        if d.shape != 1e3 || (0.5 < x/d.scale < 1.8) # rmath overflow/NaN
            test_rmath(pdf,d,x,1e4)
            if x < 1e290 # dist.jl overflow
                test_rmath(logpdf,d,x,1e5)
            end
        end

        test_rmath(cdf,d,x)
        test_rmath(ccdf,d,x)
        test_rmath(logcdf,d,x)
        test_rmath(logccdf,d,x)
    end

    for p in logspace(-300,0,1_000_000)
        test_rmath(quantile,d,p,1e4)
        test_rmath(cquantile,d,p,1e4)
    end

    for lp in -logspace(-300,300,1_000_000)
        test_rmath(invlogcdf,d,lp,1e4)
        test_rmath(invlogccdf,d,lp,1e4)
    end
end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  x1'(t) = -(k3 + k7) * x1(t) + k4 * x2(t),
  x2'(t) = k3 * x1(t) - (k4 + a(t) * k5 + b(t) * d(t) * k5) * x2(t) + k6 * x3(t) + k6 * x4(t) + k5 * x2(t) * x3(t) + k5 * x2(t) * x4(t),
  x3'(t) = a(t) * k5 * x2(t) - k6 * x3(t) - k5 * x2(t) * x3(t),
  x4'(t) = b(t) * d(t) * k5 * x2(t) - k6 * x4(t) - k5 * x2(t) * x4(t),
  x5'(t) = k7 * x1(t),
  a'(t) = 0,
  b'(t) = 0,
  d'(t) = 0,
  y1(t) = x5(t),
  y2(t) = a(t),
  y3(t) = b(t),
  y4(t) = d(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
module NewOptimizer

include("ir.jl")
include("domtree.jl")
include("passes.jl")
include("slot2ssa.jl")
include("new_ir.jl")
include("legacy.jl")
include("queries.jl")
include("verify.jl")

end # module


text:
using Gimbal
using Test
using BenchmarkTools

include("examples.jl")

@testset "pid-example" begin
    sol = solve(ex_pid())#, ImplicitEuler(),adaptive=false, dt=4e-4)
    #plot(sol, vars = [p.x])#, pid1.e, pid1.u])
    steady_state = last(sol.u)[3]

    @test isapprox(steady_state, 1.0 ,atol=1e-4)
    time = @belapsed solve(ex_pid())
    @test time < 10 #why?
end

@testset "lead lag-example" begin
    sol = solve(ex_lead_lag())
    steady_state = last(sol.u)[1]
    @test isapprox(steady_state, 1.0 ,rtol=0.02) #%2 steady state error
end


text:
using Packing, GeometryTypes
using Test

# write your own tests here
include("test_rectangle.jl")


text:
# Compare with the sample dataset from E2LSH

using LSH
using Distances
using Base.Test

points = readdlm("mnist1k.dts")
querypoints = readdlm("mnist1k.q")

# d = 784
# w = 4
# Parameters computed by E2LSH was
# k = 14
# L = 153

T = LSHTable(0.6,LSH.createUHashesAM04(784,4.0,14,153,0.6,18),
    Vector{Float64}[ vec(points[i,:]) for i=1:size(points,1) ])

qs = Vector{Float64}[ vec(querypoints[i,:]) for i=1:size(querypoints,1) ]
for (i,q) in enumerate(qs)
    println("=== Querying point $i ===")
    @time found_points = T[q]
    println("Found $(length(found_points)) points. Distances are:")
    for p in found_points
        println("  $(euclidean(p,q))")
    end
end


text:
### wiener-cameraman.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: deconvolution, signal processing
#
# This file is a part of Deconvolution.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file produces an example image used in Deconvolution.jl documentation.
#
### Code:

using Images, TestImages, Deconvolution, Random, FFTW

Random.seed!(42)

# Open the test image
img = channelview(testimage("cameraman"))
# Create the blurring kernel in frequency domain
x = hcat(ntuple(x -> collect((1:512) .- 257), 512)...)
k = 0.001
blurring_ft = @. exp(-k*(x ^ 2 + x ^ 2)^(5//6))
# Create additive noise
noise = rand(Float64, size(img))
# Fourier transform of the blurred image, with additive noise
blurred_img_ft = fftshift(blurring_ft) .* fft(img) .+ fft(noise)
# Get the blurred image from its Fourier transform
blurred_img = real(ifft(blurred_img_ft))
# Get the blurring kernel in the space domain
blurring = ifft(fftshift(blurring_ft))
# Polish the image with Deconvolution deconvolution
polished = wiener(blurred_img, img, noise, blurring)

# Wiener deconvolution works also when you don't have the real image and noise,
# that is the most common and useful case.  This happens because the Wiener
# filter only cares about the power spectrum of the signal and the noise, so you
# don't need to have the exact signal and noise but something with a similar
# power spectrum.
img2 = channelview(testimage("livingroom")) # Load another image
noise2 = rand(Float64, size(img)) # Create another additive noise
# Polish the image with Deconvolution deconvolution
polished2 = wiener(blurred_img, img2, noise2, blurring)
# Wiener also works using a real number instead of a noise array
polished3 = 0.5 .* wiener(blurred_img, img2, 1000, blurring)
polished4 = 0.5 .* wiener(blurred_img, img2, 10000, blurring)

save("original.jpg", Images.clamp01nan.(img))
save("blurred.jpg", Images.clamp01nan.(blurred_img))
save("polished.jpg", Images.clamp01nan.(polished))
save("polished2.jpg", Images.clamp01nan.(polished2))
save("polished3.jpg", Images.clamp01nan.(polished3))
save("polished4.jpg", Images.clamp01nan.(polished4))


text:
"""
    clebsch_gordan_condon_shortley(j₁, m₁, j₂, m₂, j₃, m₃=m₁+m₂)

Calculate the vector coupling coefficient `⟨j₁ m₁, j₂ m₂|j₃ m₃⟩`
according to the Condon–Shortley phase convention and using the
definition of Eq. (8.1.12) in Varshalovich (1988).
"""
clebsch_gordan_condon_shortley(j₁, m₁, j₂, m₂, j₃, m₃=m₁+m₂) =
    powneg1(Int(j₁ - j₂ + m₃))*∏(j₃)*wigner3j(j₁, j₂, j₃,
                                              m₁, m₂, -m₃)

export clebsch_gordan_condon_shortley


text:
using MathOptInterface

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIT = MOI.Test
const MOIB = MOI.Bridges

const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

MOIU.@model(LPIPModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
            (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives),
            (),
            (),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

const optimizer = MOIU.CachingOptimizer(LPIPModelData{Float64}(), LPIP.Optimizer())
const config = MOIT.TestConfig(atol=1e-5, rtol=1e-5, infeas_certificates = false)

@testset "MOI tests" begin
    @testset "SolverName" begin
        @test MOI.get(optimizer, MOI.SolverName()) == "LPIP"
    end
    @testset "Unit" begin
        MOIT.unittest(MOIB.full_bridge_optimizer(optimizer, Float64), config,[
            # Quadratic functions are not supported
            "solve_qcp_edge_cases", "solve_qp_edge_cases",
            # Integer and ZeroOne sets are not supported
            "solve_integer_edge_cases", "solve_objbound_edge_cases"
            ]
        )
    end
    @testset "MOI Continuous Linear" begin
        MOIT.contlineartest(MOIB.full_bridge_optimizer(optimizer, Float64), config, [
            # infeasible/unbounded
            # "linear8a", "linear8b", "linear8c", "linear12",
            # linear10 is poorly conditioned
            "linear10",
            # linear9 is requires precision
            "linear9",
            # primalstart not accepted
            "partial_start"
            ]
        )
        # MOIT.linear9test(MOIB.SplitInterval{Float64}(optimizer_lin_hd), config)
    end
end


text:
using Fermi.DIIS

using TensorOperations
using LinearAlgebra
using Formatting
import Base: show

export RHF

"""
    Fermi.HartreeFock.RHFAlgorithm

Abstract type for RHF implementations.
"""
abstract type RHFAlgorithm end

"""
    Fermi.HartreeFock.RHF

Wave function object for Restricted Hartree-Fock methods

# High Level Interface
Run a RHF computation and return the RHF object:
```
julia> @energy rhf
```
Equivalent to
```
julia> Fermi.HartreeFock.RHF()
```
Computes RHF using information from `Fermi.Options.Current`

# Fields

| Name   |   Description     |
|--------|---------------------|
| `molecule` |   Molecule object |
| `energy`   |   RHF Energy      |
| `ndocc`    | Number of doubly occupied spatial orbitals |
| `nvir`  | Number of virtual spatial orbitals |
| `orbitals` |    RHF Orbitals object      |
| `e_conv`   | ΔE from the last iteration  |
| `d_conv`   |  Orbitals RMS change from the last iteration|

# Relevant options

These options can be set with `@set <option> <value>`

| Option         | What it does                      | Type      | choices [default]     |
|----------------|-----------------------------------|-----------|-----------------------|
| `rhf_alg`      | Picks RHF algorithm               | `Int`     | [1]                   |
| `scf_max_rms`  | RMS density convergence criterion | `Float64` | [10^-9]               |
| `scf_max_iter` | Max number of iterations          | `Int`     | [50]                  |
| `scf_e_conv`   | Energy convergence criterion      | `Float64` | [10^-10]              |
| `basis`        | What basis set to use             | `String`  | ["sto-3g"]            |
| `df`           | Whether to use density fitting    | `Bool`    | `true` [`false`]      |
| `jkfit`        | What aux. basis set to use for JK | `String`  | ["auto"]              |
| `diis`         | Whether to use DIIS               | `Bool`    | [`true`] `false`      |
| `oda`          | Whether to use ODA                | `Bool`    | [`true`] `false`      |
| `oda_cutoff`   | When to turn ODA off (RMS)        | `Float64` | [1E-1]                |
| `oda_shutoff`  | When to turn ODA off (iter)       | `Int`     | [20]                  |
| `scf_guess`    | Which guess density to use        | `String`  | "core" ["gwh"]        |

# Struct tree

**RHF** <: AbstractHFWavefunction <: AbstractWavefunction
"""
struct RHF <: AbstractHFWavefunction
    molecule::Molecule
    energy::Float64
    ndocc::Int
    nvir::Int
    orbitals::RHFOrbitals
    e_conv::Float64
    d_conv::Float64
end

function RHF(x...)
    if !any(i-> i isa RHFAlgorithm, x)
        RHF(x..., get_rhf_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RHF method: $args"))
    end
end

"""
    Fermi.HartreeFock.get_rhf_alg()

Returns a singleton type corresponding to a RHF implementation.
"""
function get_rhf_alg(N::Int = Options.get("rhf_alg"))
    try
        return get_rhf_alg(Val(N))
    catch MethodError
        throw(FermiException("implementation number $N not available for RHF."))
    end
end

# Actual RHF routine is in here
# For each implementation a singleton type must be create
struct RHFa <: RHFAlgorithm end
include("RHFa.jl")
include("RHFHelper.jl")
# And a number is assigned to the implementation
get_rhf_alg(x::Val{1}) = RHFa()

struct Direct <: RHFAlgorithm end
include("Direct.jl")
include("DirectHelper.jl")
get_rhf_alg(x::Val{2}) = Direct()

# Gradient methods
include("Gradients/RHFgrad.jl")

### MISCELLANEOUS
# Pretty printing
function string_repr(X::RHF)
    out = ""
    out = out*" ⇒ Fermi Restricted Hartree--Fock Wave function\n"
    out = out*" ⋅ Basis:                  $(X.orbitals.basis)\n"
    out = out*" ⋅ Energy:                 $(X.energy)\n"
    out = out*" ⋅ Occ. Spatial Orbitals:  $(X.ndocc)\n"
    out = out*" ⋅ Vir. Spatial Orbitals:  $(X.nvir)\n"
    out = out*"Convergence: "
    out = out*"ΔE => $(format("{:1.2e}",abs(X.e_conv)))"
    out = out*" Dᵣₘₛ => $(format("{:1.2e}",abs(X.d_conv)))"
    return out
end

function show(io::IO, ::MIME"text/plain", X::RHF)
    print(io, string_repr(X))
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using Scratch, NPZ

download_cache = ""

function __init__()
    global download_cache = @get_scratch!("downloaded_files")
end

function download_data(url)
    fname = joinpath(download_cache, basename(url))
    if !isfile(fname)
        download(url, fname)
    end
    return fname
end

function get_data()
    fname = download_data("https://people.eecs.berkeley.edu/~bmild/nerf/tiny_nerf_data.npz")
    npzread(fname)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Uniform random generation

# This file contains the creation of Sampler objects and the associated generation of
# random values from them. More specifically, given the specification S of a set
# of values to pick from (e.g. 1:10, or "a string"), we define
#
# 1) Sampler(rng, S, ::Repetition) -> sampler
# 2) rand(rng, sampler) -> random value
#
# Note that the 1) is automated when the sampler is not intended to carry information,
# i.e. the default fall-backs SamplerType and SamplerTrivial are used.

## from types: rand(::Type, [dims...])

### random floats

Sampler(rng::AbstractRNG, ::Type{T}, n::Repetition) where {T<:AbstractFloat} =
    Sampler(rng, CloseOpen(T), n)

# generic random generation function which can be used by RNG implementors
# it is not defined as a fallback rand method as this could create ambiguities

rand_generic(r::AbstractRNG, ::CloseOpen{Float16}) =
    Float16(reinterpret(Float32,
                        (rand_ui10_raw(r) % UInt32 << 13) & 0x007fe000 | 0x3f800000) - 1)

rand_generic(r::AbstractRNG, ::CloseOpen{Float32}) =
    reinterpret(Float32, rand_ui23_raw(r) % UInt32 & 0x007fffff | 0x3f800000) - 1

rand_generic(r::AbstractRNG, ::Close1Open2_64) =
    reinterpret(Float64, 0x3ff0000000000000 | rand(r, UInt64) & 0x000fffffffffffff)

rand_generic(r::AbstractRNG, ::CloseOpen_64) = rand(r, Close1Open2()) - 1.0

#### BigFloat

const bits_in_Limb = sizeof(Limb) << 3
const Limb_high_bit = one(Limb) << (bits_in_Limb-1)

struct SamplerBigFloat{I<:FloatInterval{BigFloat}} <: Sampler
    prec::Int
    nlimbs::Int
    limbs::Vector{Limb}
    shift::UInt

    function SamplerBigFloat{I}(prec::Int) where I<:FloatInterval{BigFloat}
        nlimbs = (prec-1) ÷ bits_in_Limb + 1
        limbs = Vector{Limb}(uninitialized, nlimbs)
        shift = nlimbs * bits_in_Limb - prec
        new(prec, nlimbs, limbs, shift)
    end
end

Sampler(::AbstractRNG, I::FloatInterval{BigFloat}, ::Repetition) =
    SamplerBigFloat{typeof(I)}(precision(BigFloat))

function _rand(rng::AbstractRNG, sp::SamplerBigFloat)
    z = BigFloat()
    limbs = sp.limbs
    rand!(rng, limbs)
    @inbounds begin
        limbs[1] <<= sp.shift
        randbool = iszero(limbs[end] & Limb_high_bit)
        limbs[end] |= Limb_high_bit
    end
    z.sign = 1
    Base.@gc_preserve limbs unsafe_copy!(z.d, pointer(limbs), sp.nlimbs)
    (z, randbool)
end

function _rand(rng::AbstractRNG, sp::SamplerBigFloat, ::Close1Open2{BigFloat})
    z = _rand(rng, sp)[1]
    z.exp = 1
    z
end

function _rand(rng::AbstractRNG, sp::SamplerBigFloat, ::CloseOpen{BigFloat})
    z, randbool = _rand(rng, sp)
    z.exp = 0
    randbool &&
        ccall((:mpfr_sub_d, :libmpfr), Int32,
              (Ref{BigFloat}, Ref{BigFloat}, Cdouble, Int32),
              z, z, 0.5, Base.MPFR.ROUNDING_MODE[])
    z
end

# alternative, with 1 bit less of precision
# TODO: make an API for requesting full or not-full precision
function _rand(rng::AbstractRNG, sp::SamplerBigFloat, ::CloseOpen{BigFloat}, ::Void)
    z = _rand(rng, sp, Close1Open2(BigFloat))
    ccall((:mpfr_sub_ui, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Culong, Int32),
          z, z, 1, Base.MPFR.ROUNDING_MODE[])
    z
end

rand(rng::AbstractRNG, sp::SamplerBigFloat{T}) where {T<:FloatInterval{BigFloat}} =
    _rand(rng, sp, T())

### random integers

rand_ui10_raw(r::AbstractRNG) = rand(r, UInt16)
rand_ui23_raw(r::AbstractRNG) = rand(r, UInt32)

rand_ui52_raw(r::AbstractRNG) = reinterpret(UInt64, rand(r, Close1Open2()))
rand_ui52(r::AbstractRNG) = rand_ui52_raw(r) & 0x000fffffffffffff

### random complex numbers

rand(r::AbstractRNG, ::SamplerType{Complex{T}}) where {T<:Real} =
    complex(rand(r, T), rand(r, T))

### random characters

# returns a random valid Unicode scalar value (i.e. 0 - 0xd7ff, 0xe000 - # 0x10ffff)
function rand(r::AbstractRNG, ::SamplerType{Char})
    c = rand(r, 0x00000000:0x0010f7ff)
    (c < 0xd800) ? Char(c) : Char(c+0x800)
end

## Generate random integer within a range

### BitInteger

# remainder function according to Knuth, where rem_knuth(a, 0) = a
rem_knuth(a::UInt, b::UInt) = a % (b + (b == 0)) + a * (b == 0)
rem_knuth(a::T, b::T) where {T<:Unsigned} = b != 0 ? a % b : a

# maximum multiple of k <= 2^bits(T) decremented by one,
# that is 0xFFFF...FFFF if k = typemax(T) - typemin(T) with intentional underflow
# see http://stackoverflow.com/questions/29182036/integer-arithmetic-add-1-to-uint-max-and-divide-by-n-without-overflow
maxmultiple(k::T) where {T<:Unsigned} =
    (div(typemax(T) - k + oneunit(k), k + (k == 0))*k + k - oneunit(k))::T

# maximum multiple of k within 1:2^32 or 1:2^64 decremented by one, depending on size
maxmultiplemix(k::UInt64) = k >> 32 != 0 ?
    maxmultiple(k) :
    (div(0x0000000100000000, k + (k == 0))*k - oneunit(k))::UInt64

struct SamplerRangeInt{T<:Integer,U<:Unsigned} <: Sampler
    a::T   # first element of the range
    k::U   # range length or zero for full range
    u::U   # rejection threshold
end

# generators with 32, 128 bits entropy
SamplerRangeInt(a::T, k::U) where {T,U<:Union{UInt32,UInt128}} =
    SamplerRangeInt{T,U}(a, k, maxmultiple(k))

# mixed 32/64 bits entropy generator
SamplerRangeInt(a::T, k::UInt64) where {T} =
    SamplerRangeInt{T,UInt64}(a, k, maxmultiplemix(k))

function Sampler(::AbstractRNG, r::AbstractUnitRange{T}, ::Repetition) where T<:Unsigned
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    SamplerRangeInt(first(r), last(r) - first(r) + oneunit(T))
end

for (T, U) in [(UInt8, UInt32), (UInt16, UInt32),
               (Int8, UInt32), (Int16, UInt32), (Int32, UInt32),
               (Int64, UInt64), (Int128, UInt128), (Bool, UInt32)]

    @eval Sampler(::AbstractRNG, r::AbstractUnitRange{$T}, ::Repetition) = begin
        isempty(r) && throw(ArgumentError("range must be non-empty"))
        # overflow ok:
        SamplerRangeInt(first(r), convert($U, unsigned(last(r) - first(r)) + one($U)))
    end
end

# this function uses 32 bit entropy for small ranges of length <= typemax(UInt32) + 1
# SamplerRangeInt is responsible for providing the right value of k
function rand(rng::AbstractRNG, sp::SamplerRangeInt{T,UInt64}) where T<:Union{UInt64,Int64}
    local x::UInt64
    if (sp.k - 1) >> 32 == 0
        x = rand(rng, UInt32)
        while x > sp.u
            x = rand(rng, UInt32)
        end
    else
        x = rand(rng, UInt64)
        while x > sp.u
            x = rand(rng, UInt64)
        end
    end
    return reinterpret(T, reinterpret(UInt64, sp.a) + rem_knuth(x, sp.k))
end

function rand(rng::AbstractRNG, sp::SamplerRangeInt{T,U}) where {T<:Integer,U<:Unsigned}
    x = rand(rng, U)
    while x > sp.u
        x = rand(rng, U)
    end
    (unsigned(sp.a) + rem_knuth(x, sp.k)) % T
end

### BigInt

struct SamplerBigInt <: Sampler
    a::BigInt         # first
    m::BigInt         # range length - 1
    nlimbs::Int       # number of limbs in generated BigInt's (z ∈ [0, m])
    nlimbsmax::Int    # max number of limbs for z+a
    mask::Limb        # applied to the highest limb
end

function Sampler(::AbstractRNG, r::AbstractUnitRange{BigInt}, ::Repetition)
    m = last(r) - first(r)
    m < 0 && throw(ArgumentError("range must be non-empty"))
    nd = ndigits(m, 2)
    nlimbs, highbits = divrem(nd, 8*sizeof(Limb))
    highbits > 0 && (nlimbs += 1)
    mask = highbits == 0 ? ~zero(Limb) : one(Limb)<<highbits - one(Limb)
    nlimbsmax = max(nlimbs, abs(last(r).size), abs(first(r).size))
    return SamplerBigInt(first(r), m, nlimbs, nlimbsmax, mask)
end

function rand(rng::AbstractRNG, sp::SamplerBigInt)
    x = MPZ.realloc2(sp.nlimbsmax*8*sizeof(Limb))
    limbs = unsafe_wrap(Array, x.d, sp.nlimbs)
    while true
        rand!(rng, limbs)
        @inbounds limbs[end] &= sp.mask
        MPZ.mpn_cmp(x, sp.m, sp.nlimbs) <= 0 && break
    end
    # adjust x.size (normally done by mpz_limbs_finish, in GMP version >= 6)
    x.size = sp.nlimbs
    while x.size > 0
        @inbounds limbs[x.size] != 0 && break
        x.size -= 1
    end
    MPZ.add!(x, sp.a)
end

## random values from AbstractArray

Sampler(rng::AbstractRNG, r::AbstractArray, n::Repetition) =
    SamplerSimple(r, Sampler(rng, linearindices(r), n))

rand(rng::AbstractRNG, sp::SamplerSimple{<:AbstractArray,<:Sampler}) =
    @inbounds return sp[][rand(rng, sp.data)]

## random values from Dict

function Sampler(rng::AbstractRNG, t::Dict, ::Repetition)
    isempty(t) && throw(ArgumentError("collection must be non-empty"))
    # we use Val(Inf) below as rand is called repeatedly internally
    # even for generating only one random value from t
    SamplerSimple(t, Sampler(rng, linearindices(t.slots), Val(Inf)))
end

function rand(rng::AbstractRNG, sp::SamplerSimple{<:Dict,<:Sampler})
    while true
        i = rand(rng, sp.data)
        Base.isslotfilled(sp[], i) && @inbounds return (sp[].keys[i] => sp[].vals[i])
    end
end

## random values from Set

Sampler(rng::AbstractRNG, t::Set, n::Repetition) = SamplerTag{Set}(Sampler(rng, t.dict, n))

rand(rng::AbstractRNG, sp::SamplerTag{Set,<:Sampler}) = rand(rng, sp.data).first

## random values from BitSet

function Sampler(rng::AbstractRNG, t::BitSet, n::Repetition)
    isempty(t) && throw(ArgumentError("collection must be non-empty"))
    SamplerSimple(t, Sampler(rng, linearindices(t.bits), Val(Inf)))
end

function rand(rng::AbstractRNG, sp::SamplerSimple{BitSet,<:Sampler})
    while true
        n = rand(rng, sp.data)
        @inbounds b = sp[].bits[n]
        b && return n
    end
end

## random values from Associative/AbstractSet

# we defer to _Sampler to avoid ambiguities with a call like Sampler(rng, Set(1), Val(1))
Sampler(rng::AbstractRNG, t::Union{Associative,AbstractSet}, n::Repetition) =
    _Sampler(rng, t, n)

# avoid linear complexity for repeated calls
_Sampler(rng::AbstractRNG, t::Union{Associative,AbstractSet}, n::Val{Inf}) =
    Sampler(rng, collect(t), n)

# when generating only one element, avoid the call to collect
_Sampler(::AbstractRNG, t::Union{Associative,AbstractSet}, ::Val{1}) =
    SamplerTrivial(t)

function nth(iter, n::Integer)::eltype(iter)
    for (i, x) in enumerate(iter)
        i == n && return x
    end
end

rand(rng::AbstractRNG, sp::SamplerTrivial{<:Union{Associative,AbstractSet}}) =
    nth(sp[], rand(rng, 1:length(sp[])))

## random characters from a string

# we use collect(str), which is most of the time more efficient than specialized methods
# (except maybe for very small arrays)
Sampler(rng::AbstractRNG, str::AbstractString, n::Val{Inf}) = Sampler(rng, collect(str), n)

# when generating only one char from a string, the specialized method below
# is usually more efficient
Sampler(rng::AbstractRNG, str::AbstractString, ::Val{1}) =
    SamplerSimple(str, Sampler(rng, 1:_endof(str), Val(Inf)))

isvalid_unsafe(s::String, i) = !Base.is_valid_continuation(Base.@gc_preserve s unsafe_load(pointer(s), i))
isvalid_unsafe(s::AbstractString, i) = isvalid(s, i)
_endof(s::String) = sizeof(s)
_endof(s::AbstractString) = endof(s)

function rand(rng::AbstractRNG, sp::SamplerSimple{<:AbstractString,<:Sampler})::Char
    str = sp[]
    while true
        pos = rand(rng, sp.data)
        isvalid_unsafe(str, pos) && return str[pos]
    end
end

max_stars_repo_name: JuliaLangEs/julieta

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04dy(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    LDA = NMAX
    LDQG = NMAX
    LDWORK = NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    QG = Array{Float64,2}(undef, LDQG,NMAX+1)
    D = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    JOBSCL =  vs[2][1]
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*(_jsz+1)÷2
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    _i0 = 0
    for j in 1:_jsz
       QG[j,j+1:_jsz+1] .= parsex.(Float64, vs[_i0+1:_i0+_jsz-j+1])
       _i0 += _jsz-j+1
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*(_jsz+1)÷2
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    _i0 = 0
    for j in 1:_jsz
       QG[j:_jsz,j] .= parsex.(Float64, vs[_i0+1:_i0+_jsz-j+1])
       _i0 += _jsz-j+1
    end
    close(f)

    # interp call 1

    INFO = SLICOT.mb04dy!(JOBSCL, N, A, QG, D)
    @test INFO == 0
    INFO != 0 && return

# interp output 1
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "QG:")
    _nc = N+1
    _nr = N
    show(io, "text/plain", QG[1:_nr,1:_nc])
    println(io)

    if ( LSAME( JOBSCL, 'S' ) )
        # interp output 3
        println(io, "D:")
        _nr = N
        show(io, "text/plain", D[1:_nr])
        println(io)

    elseif ( LSAME( JOBSCL, '1' ) || LSAME( JOBSCL, 'O' ) )
        println(io, "tau = ",D[1])
    end # if

end # run_mb04dy()


text:
using PyCall
PyCall.pyimport("sys")[:executable]

has_numpy = try
    PyCall.pyimport("numpy")
    true
catch
    false
end
if has_numpy
    # Invoke numpy support (which mutates various states in PyCall) to
    # test that it does not introduce any run-time bugs:
    PyCall.npyinitialize()
end


text:
nq = 3 + 4 + 3 # position + quaternion + leg joints

function quaternion_rotation_matrix(q)
	r, i, j, k  = q

	r11 = 1.0 - 2.0 * (j^2.0 + k^2.0)
	r12 = 2.0 * (i * j - k * r)
	r13 = 2.0 * (i * k + j * r)

	r21 = 2.0 * (i * j + k * r)
	r22 = 1.0 - 2.0 * (i^2.0 + k^2.0)
	r23 = 2.0 * (j * k - i * r)

	r31 = 2.0 * (i * k - j * r)
	r32 = 2.0 * (j * k + i * r)
	r33 = 1.0 - 2.0 * (i^2.0 + j^2.0)

	[r11 r12 r13;
     r21 r22 r23;
     r31 r32 r33]
end

function rotation_matrix_x(a)
    [1.0 0.0    0.0
     0.0 cos(a) -sin(a)
     0.0 sin(a)  cos(a)]
end

leg1_offset = [1.0; 1.0; 0.0]
l_shoulder = 1.0
l_thigh = 1.0
l_calf = 1.0
function kinematics_thigh_calf(a1, a2)
    [
     l_thigh * sin(a1) + l_calf * sin(a2);
     0.0;
    -l_thigh * cos(a1) - l_calf * cos(a2);
    ]
end

function foot_kinematics(q)
    pos = q[1:3]
    quat = q[4:7]
    joint = q[8:10]

    k_thigh_calf = kinematics_thigh_calf(q[9], q[10])
    k_shoulder_thigh_calf = rotation_matrix_x(q[8]) * (k_thigh_calf + [0.0; l_shoulder; 0.0])
    k_offset_shoulder_thigh_calf = k_shoulder_thigh_calf + leg1_offset
    k_body_offset_shoulder_thigh_calf = quaternion_rotation_matrix(quat) * k_offset_shoulder_thigh_calf
    k = k_body_offset_shoulder_thigh_calf + pos

    return k
end

q0 = [0.0; 0.0; 0.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
foot_kinematics(q0)

@variables q[1:nq] q̇[1:nq]

k = foot_kinematics(q)
j = Symbolics.jacobian(k, q)

k_func = eval(Symbolics.build_function(k, q)[1])
j_func = eval(Symbolics.build_function(j, q)[1])

function lagrangian(q, q̇)
    L = 0.0
    mass = 1.0
    inertia = 1.0
    gravity = 9.81
    k_foot = k_func(q)
    v_foot = j_func(q) * q̇
    L += mass * transpose(v_foot) * v_foot
    L -= mass * gravity * k_foot[3]

    return L
end

L = lagrangian(q, q̇)

ddL = Symbolics.hessian(L, [q; q̇])#, simplify=true)
dLq = Symbolics.gradient(L, q)#, simplify=true)
ddLq̇q = ddL[nq .+ (1:nq), 1:nq]

M = ddL[nq .+ (1:nq), nq .+ (1:nq)]
C = ddLq̇q * q̇ - dLq

mass_matrix_sym = eval(Symbolics.build_function(M, q)[1])
dynamics_bias_sym = eval(Symbolics.build_function(C, q, q̇)[1])

function lagrangian_derivatives(q, v)
	D1L = -1.0 * dynamics_bias_sym(q, v)
    D2L = mass_matrix_sym(q) * v
	return D1L, D2L
end

function dynamics(h, q0, q1, u1, w1, λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(qm2, vm2)

	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2 + λ1)
end

@variables h[1:1] q0[1:nq] q1[1:nq] u1[1:0] w1[1:0] λ1[1:nq] q2[1:nq]
d = dynamics(h, q0, q1, u1, w1, λ1, q2)


text:
# ---
# title: 1385. Find the Distance Value Between Two Arrays
# id: problem1385
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/find-the-distance-value-between-two-arrays/description/>
# hidden: true
# ---
#
# Given two integer arrays `arr1` and `arr2`, and the integer `d`, _return the
# distance value between the two  arrays_.
#
# The distance value is defined as the number of elements `arr1[i]` such that
# there is not any element `arr2[j]` where `|arr1[i]-arr2[j]| <= d`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr1 = [4,5,8], arr2 = [10,9,1,8], d = 2
#     Output: 2
#     Explanation:
#     For arr1[0]=4 we have:
#     |4-10|=6 > d=2
#     |4-9|=5 > d=2
#     |4-1|=3 > d=2
#     |4-8|=4 > d=2
#     For arr1[1]=5 we have:
#     |5-10|=5 > d=2
#     |5-9|=4 > d=2
#     |5-1|=4 > d=2
#     |5-8|=3 > d=2
#     For arr1[2]=8 we have:
#     **|8-10|=2 <= d=2**
#     **|8-9|=1 <= d=2**
#     |8-1|=7 > d=2
#     **|8-8|=0 <= d=2**
#
#
# **Example 2:**
#
#
#
#     Input: arr1 = [1,4,2,3], arr2 = [-4,-3,6,10,20,30], d = 3
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: arr1 = [2,1,100,3], arr2 = [-5,-2,10,-3,7], d = 6
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr1.length, arr2.length <= 500`
#   * `-10^3 <= arr1[i], arr2[j] <= 10^3`
#   * `0 <= d <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ---
# title: 709. To Lower Case
# id: problem709
# author: Tian Jun
# date: 2020-11-14
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/to-lower-case/description/>
# hidden: true
# ---
#
# Implement function ToLowerCase() that has a string parameter str, and returns
# the same string in lowercase.
#
#
#
# **Example 1:**
#
#
#
#     Input: "Hello"
#     Output: "hello"
#
#
# **Example 2:**
#
#
#
#     Input: "here"
#     Output: "here"
#
#
# **Example 3:**
#
#
#
#     Input: "LOVELY"
#     Output: "lovely"
#
#
#
## @lc code=start
using LeetCode

function to_lower_case(s::String)::String
    return lowercase(s)
end
## @lc code=end


text:
julia> """
       # Y combinator

       * `  `
       """
       Y = f -> (x -> x(x))(y ->
           f((t...) -> y(y)(t...)))

Usage:

julia> "# Factorial"
       fac = f -> (n -> n < 2 ? 1 : n * f(n - 1))

julia> "# Fibonacci"
       fib = f -> (n -> n == 0 ? 0 : (n == 1 ? 1 : f(n - 1) + f(n - 2)))
 countList = f -> (n -> n == )
julia> [Y(fac)(i) for i = 1:10]
10-element Array{Any,1}:
       1
       2
       6
      24
     120
     720
    5040
   40320
  362880
 3628800

julia> [Y(fib)(i) for i = 1:10]
10-element Array{Any,1}:
  1
  1
  2
  3
  5
  8
 13
 21
 34
 55

 YFact = Y(fac)
 [YFact(x) for x in 1:10]

max_stars_repo_name: HerbM/WindowsPowerShell

max_issues_repo_name: HerbM/WindowsPowerShell

max_forks_repo_name: HerbM/WindowsPowerShell


text:
using QuasiGeostrophy, LinearAlgebra, Test, FFTW, BenchmarkTools, Plots

a = 0; b = 2π
Ω = S¹(a,b) # Domain
Nx = 2^7;
fourier_grid = FourierGrid(Nx, Ω) # Grid
x = fourier_grid.grid[1]
kx = fourier_grid.wavenumbers[1]
f = @. exp(-2 * (b-a) / 3 * (x - (b-a)/2)^2)*(1+0im)
f̂ = copy(f)
FFTW.set_num_threads(Threads.nthreads())
P = plan_fft(f) # typeof(iP) <: AbstractFFTs.ScaledPlan
iP = plan_ifft(f)
mul!(f̂, P, f)
∂x = FourierOperator(im .* kx)
ν = 0.01

function rhs(f̂, ∂x, P, iP, ν)
    ff = filter_convolve(f̂, f̂, P, iP)
    return -∂x(ff) + ν * ∂x(∂x(f̂))
end
function closure_rhs(∂x, P, iP, ν)
    function closed_rhs(f̂)
        ff = filter_convolve(f̂, f̂, P, iP)
        #ff = convolve(f̂, f̂, P, iP)
        return -∂x(ff) + ν * ∂x(∂x(f̂))
    end
end
old_rhs(f̂) = rhs(f̂, ∂x, P, iP, ν)
new_rhs = closure_rhs(∂x, P, iP, ν)
plot(x, real.(f))
plot(x, real.(iP * ∂x(f̂)))
plot(x, real.(iP * rhs(f̂, ∂x, P, iP, ν)))

function step!(f̂, rhs, Δt)
    rhsf = rhs(f̂)
    @. f̂ += Δt * rhsf
end
function imex_step!(f̂, rhs, Δt)
    ff = convolve(f̂, f̂, P, iP)
    rhsf = -∂x(ff)
    @. f̂ += Δt * rhsf
    lhsf = 1 .- Δt .* ν .* (∂x^2).op
    @. f̂ = f̂ / lhsf
end
##
f = @. exp(-2 * (b-a) / 3 * (x - (b-a)/2)^2)*(1+0im)
f̂ = copy(f)
mul!(f̂, P, f)
for i in 1:2π*200
    local Δt = 0.001
    # step!(f̂, new_rhs, Δt)
    imex_step!(f̂, new_rhs, Δt)
    if (i%100) == 0
        mul!(f, iP, f̂)
        p1 = plot(x, real.(f), ylims = (0,1), label = "real space", xlabel = "position", ylabel  = "amplitude")
        p2 = scatter(log.(abs.(f̂))[1:div(length(f̂),2)+1], label = "spectrum", ylims = (-7, 3), xlabel = "wavenumber", ylabel = "spectral amplitude")
        display(plot(p1,p2))
        sleep(0.1)
    end
end
##
gr(size = (250,400))
p1 = plot(x, real.(f), ylims = (0,1), label = "real space", xlabel = "position", ylabel  = "amplitude")
p2 = scatter(log.(abs.(f̂))[1:div(length(f̂),2)+1], label = "spectrum", ylims = (-7, 3), xlabel = "wavenumber", ylabel = "spectral amplitude")

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using NLSolvers
function f∇f!(∇f, x)
   if !(∇f==nothing)
       if ( x[1]^2 + x[2]^2 == 0 )
           dtdx1 = 0;
           dtdx2 = 0;
       else
           dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
           dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
       end
       ∇f[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
           200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
       ∇f[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
           200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
       ∇f[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
   end

   fx = f(x)
   return ∇f==nothing ? fx : (fx, ∇f)
end

f(x) = 100.0 * ((x[3] - 10.0 * theta(x))^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
function theta(x)
  if x[1] > 0
      return atan(x[2] / x[1]) / (2.0 * pi)
  else
      return (pi + atan(x[2] / x[1])) / (2.0 * pi)
  end
end

obj! = OnceDiffed(f∇f!)

minimize!(obj!, -rand(3)*9 .- 3, NLSolvers.NelderMead())
V = [[1.0,1.0,1.0], [0.0,1.0,1.0],[0.40,0.0,0.0],[-1.0,2.0,.03]]
F = f∇f!.(nothing, V)
splx = NLSolvers.ValuedSimplex(V, F)
minimize!(obj!, splx, NLSolvers.NelderMead())
@profiler minimize!(obj!, splx, NLSolvers.NelderMead(); itermax=1500)

function powell(∇f, x)
    fx = (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
        (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4

    if !isa(∇f, Nothing)
        ∇f[1] = 2.0 * (x[1] + 10.0 * x[2]) + 40.0 * (x[1] - x[4])^3
        ∇f[2] = 20.0 * (x[1] + 10.0 * x[2]) + 4.0 * (x[2] - 2.0 * x[3])^3
        ∇f[3] = 10.0 * (x[3] - x[4]) - 8.0 * (x[2] - 2.0 * x[3])^3
        ∇f[4] = -10.0 * (x[3] - x[4]) - 40.0 * (x[1] - x[4])^3
        return ∇f, fx
    end
    return fx
end

function powell_optim(F, ∇f, x)
    if !isa(F, Nothing)
        fx = (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
            (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4
    end
    if !isa(∇f, Nothing)
        ∇f[1] = 2.0 * (x[1] + 10.0 * x[2]) + 40.0 * (x[1] - x[4])^3
        ∇f[2] = 20.0 * (x[1] + 10.0 * x[2]) + 4.0 * (x[2] - 2.0 * x[3])^3
        ∇f[3] = 10.0 * (x[3] - x[4]) - 8.0 * (x[2] - 2.0 * x[3])^3
        ∇f[4] = -10.0 * (x[3] - x[4]) - 40.0 * (x[1] - x[4])^3
        if !isa(F, Nothing)
            return fx
        else
            return Nothing
        end
    end
    return fx
end
obj_powell = OnceDiffed(powell)
# Define vertices in V
x0 = [1.0,1.0,1.0,1.0]
V = [copy(x0)]
for i = 1:4
    push!(V, x0+39*rand(4))
end
V
F = obj_powell.(V)
splx = NLSolvers.ValuedSimplex(V, F)
@time minimize!(obj_powell, splx, NLSolvers.NelderMead(); itermax=3000)
@time minimize!(obj_powell, copy(x0), NLSolvers.NelderMead(); itermax=3000)
@allocated minimize!(obj_powell, splx, NLSolvers.NelderMead(); itermax=3000)
@profiler minimize!(obj_powell, copy(x0), NLSolvers.NelderMead(); itermax=3000)

using Optim
res = optimize(Optim.only_fg!(powell_optim), copy(x0), NelderMead())
res = optimize(Optim.only_fg!(powell_optim), copy(x0), Optim.GradientDescent())
@time optimize(Optim.only_fg!(powell_optim), x0, NelderMead(), Optim.Options(g_tol=0, iterations =3000))

function extros!(storage, x::AbstractArray)
   n = length(x)
   jodd = 1:2:n-1
   jeven = 2:2:n
   xt = similar(x)
   @. xt[jodd] = 10.0 * (x[jeven] - x[jodd]^2)
   @. xt[jeven] = 1.0 - x[jodd]

   if !isa(storage, Nothing)
       @. storage[jodd] = -20.0 * x[jodd] * xt[jodd] - xt[jeven]
       @. storage[jeven] = 10.0 * xt[jodd]
       return 0.5*sum(abs2, xt), storage
   end
   return 0.5*sum(abs2, xt)
end
x0=rand(300)
V = [copy(x0)]
for i = 1:300
    push!(V, x0+39*rand(length(x0)))
end
V
extros_obj! = OnceDiffed(extros!)
F = extros_obj!.(V)
splx = NLSolvers.ValuedSimplex(V, F)
@time minimize!(extros_obj!, splx, NLSolvers.NelderMead(); itermax=3000)


text:
__precompile__()

module MelGeneralizedCepstrums

export
    # Spectral parameter types
    SpectralParam,
    MelGeneralizedCepstrum,
    GeneralizedCepstrum,
    MelCepstrum,
    LinearCepstrum,
    AllPoleCepstrum,
    MelAllPoleCepstrum,
    LinearPredictionCoefVariants,
    LinearPredictionCoef,
    LineSpectralPair,
    PartialAutoCorrelation,

    # State, which keeps actual spectral parameter values
    SpectralParamState, # subtype of AbstractArray

    # typealiases for SpectralParamState
    SpectralParamStateVector,
    SpectralParamStateMatrix,

    # Generic interfaces
    estimate,        # estimate spectral parameter values
    params,          # returns parameters of a spectral parameter
    param_order,     # returns order of a spectral parameter

    # Mel-generalized cepstrum properties
    allpass_alpha,   # all-pass constant (α)
    glog_gamma,      # parameter of generalized log function (γ)

    # Spectral state properties
    paramdef,        # returns parameter definition from state
    rawdata,         # returns rawdata in state
    has_loggain,
    gain_normalized,
    ready_to_filt,

    loggain!,
    unloggain!,

    # Spectral parameter estimation
    gcep,            # generalized cepstrum analysis
    mcep,            # mel-cepstrum analysis
    mgcep,           # mel-generalized cepstrum analysis
    lpc,             # linear prediction (pipe to mgcep for now)
    periodogram2mcep,# periodogram to mel-cepstrum

    mcepalpha,       # determine appropriate all pass constant in mel-cesptum

    # LPC, LSP and PARCOR conversions
    lpc2c,           # LPC to cepstrum
    lpc2lsp,         # LPC -> LSP
    lpc2par,         # LPC -> PARCOR
    par2lpc,         # PARCOR -> LPC

    # Mel-generalized cepstrum conversions
    mgc2b,           # mel-generalized cepstrum -> MGLSADF coef.
    mgc2b!,          #
    mc2b,            # mel-cepstrum -> MLSADF coef.
    mc2b!,           #
    b2mc,            # MLSADF coef. -> mel-cepstrum
    b2mc!,           #
    c2ir,            # cepstrum -> impulse response
    gc2gc,           # generalized cepstrum -> generalized cepstrum
    gnorm,           # gain normalization
    gnorm!,          #
    ignorm,          # inverse gain normalization
    ignorm!,         #
    freqt,           # frequency transform
    mgc2mgc,         # mel-generalized cepstrum -> mel-generalized cepstrum
    mgc2sp           # mel-generalized cepstrum -> spectrum envelope

for fname in [
              "common",
              "gcep",
              "mcep",
              "mgcep",
              "lpc",
              "lpc2c",
              "lpc2lsp",
              "lpc2par",
              "par2lpc",
              "mcepalpha",
              "mgc2b",
              "b2c",
              "mc2e",
              "frqtr",
              "mc2b",
              "b2mc",
              "c2ir",
              "gc2gc",
              "gnorm",
              "ignorm",
              "freqt",
              "mgc2mgc",
              "mgc2sp",
              "extend"
    ]
    include(string(fname, ".jl"))
end

end # module


text:
module FusedLassoMod
using StatsBase
import Base: +, -, *
export FusedLasso

# Implements the algorithm described in Johnson, N. A. (2013). A
# Dynamic Programming Algorithm for the Fused Lasso and
# L0-Segmentation. Journal of Computational and Graphical Statistics,
# 22(2), 246–260. doi:10.1080/10618600.2012.681238

struct NormalCoefs{T}
    lin::T
    quad::T

    NormalCoefs{T}(lin::Real) where {T} = new(lin, 0)
    NormalCoefs{T}(lin::Real, quad::Real) where {T} = new(lin, quad)
end
+(a::NormalCoefs{T}, b::NormalCoefs{T}) where {T} = NormalCoefs{T}(a.lin+b.lin, a.quad+b.quad)
-(a::NormalCoefs{T}, b::NormalCoefs{T}) where {T} = NormalCoefs{T}(a.lin-b.lin, a.quad-b.quad)
+(a::NormalCoefs{T}, b::Real) where {T} = NormalCoefs{T}(a.lin+b, a.quad)
-(a::NormalCoefs{T}, b::Real) where {T} = NormalCoefs{T}(a.lin-b, a.quad)
*(a::Real, b::NormalCoefs{T}) where {T} = NormalCoefs{T}(a*b.lin, a*b.quad)

# Implements Algorithm 2 lines 8 and 19
solveforbtilde(a::NormalCoefs{T}, lhs::Real) where {T} = (lhs - a.lin)/(2 * a.quad)

# These are marginally faster than computing btilde explicitly because
# they avoid division
btilde_lt(a::NormalCoefs{T}, lhs::Real, x::Real) where {T} = lhs - a.lin > 2 * a.quad * x
btilde_gt(a::NormalCoefs{T}, lhs::Real, x::Real) where {T} = lhs - a.lin < 2 * a.quad * x

struct Knot{T,S}
    pos::T
    coefs::S
    sign::Int8
end

struct FusedLasso{T,S} <: RegressionModel
    β::Vector{T}              # Coefficients
    knots::Vector{Knot{T,S}}  # Active knots
    bp::Matrix{T}             # Backpointers
end

function StatsBase.fit(::Type{FusedLasso}, y::AbstractVector{T}, λ::Real; dofit::Bool=true) where T
    S = NormalCoefs{T}
    flsa = FusedLasso{T,S}(Array{T}(undef, length(y)), Array{Knot{T,S}}(undef, 2), Array{T}(undef, 2, length(y)-1))
    dofit && fit!(flsa, y, λ)
    flsa
end

function StatsBase.fit!(flsa::FusedLasso{T,S}, y::AbstractVector{T}, λ::Real) where {T,S}
    β = flsa.β
    knots = flsa.knots
    bp = flsa.bp

    length(y) == length(β) || throw(ArgumentError("input size $(length(y)) does not match model size $(length(β))"))

    resize!(knots, 2)
    knots[1] = Knot{T,S}(-Inf, S(0), 1)
    knots[2] = Knot{T,S}(Inf, S(0), -1)

    # Algorithm 1 lines 2-5
    @inbounds for k = 1:length(y)-1
        t1 = 0
        t2 = 0
        aminus = NormalCoefs{T}(y[k], -0.5)                # Algorithm 2 line 4
        for outer t1 = 1:length(knots)-1                         # Algorithm 2 line 5
            knot = knots[t1]
            aminus += knot.sign*knot.coefs                 # Algorithm 2 line 6
            btilde_lt(aminus, λ, knots[t1+1].pos) && break # Algorithm 2 line 7-8
        end
        bminus = solveforbtilde(aminus, λ)

        aplus = NormalCoefs{T}(y[k], -0.5)                 # Algorithm 2 line 15
        t2 = length(knots)
        while t2 >= 2                                      # Algorithm 2 line 16
            knot = knots[t2]
            aplus -= knot.sign*knot.coefs                  # Algorithm 2 line 17
            btilde_gt(aplus, -λ, knots[t2-1].pos) && break # Algorithm 2 line 18-19
            t2 -= 1
        end
        bplus = solveforbtilde(aplus, -λ)

        # Resize knots so that we have only knots[t1+1:t2-1] and 2
        # elements at either end. It would be better to use a different
        # data structure here.
        estlen = t2 - t1 + 3
        if estlen == 4
            resize!(knots, 4)
        else
            if t2 == length(knots)
                resize!(knots, t2+1)
            else
                deleteat!(knots, t2+2:length(knots))
            end
            if t1 == 1
                pushfirst!(knots, Knot{T,S}(-Inf, S(0), 1))
            else
                deleteat!(knots, 1:t1-2)
            end
        end
        knots[1] = Knot{T,S}(-Inf, S(λ), 1)                # Algorithm 2 line 28
        knots[2] = Knot{T,S}(bminus, aminus-λ, 1)          # Algorithm 2 line 29
        knots[end-1] = Knot{T,S}(bplus, aplus+λ, -1)       # Algorithm 2 line 20
        knots[end] = Knot{T,S}(Inf, S(-λ), -1)             # Algorithm 2 line 31
        bp[1, k] = bminus
        bp[2, k] = bplus
    end

    # Algorithm 1 line 6
    aminus = NormalCoefs{T}(y[end], -0.5)
    for t1 = 1:length(knots)
        knot = knots[t1]
        aminus += knot.sign*knot.coefs
        btilde_lt(aminus, 0, knots[t1+1].pos) && break
    end
    β[end] = solveforbtilde(aminus, 0)

    # Backtrace
    for k = length(y)-1:-1:1                        # Algorithm 1 line 6
        β[k] = min(bp[2, k], max(β[k+1], bp[1, k])) # Algorithm 1 line 7
    end
    flsa
end

StatsBase.coef(flsa::FusedLasso) = flsa.β

end


text:
"""
Tools that help using the LogNormal distribution.

- Fitting to various aggregate statistics
- Sum of correlated lognormal random variables
"""
module LogNormals

export AbstractMoments, Moments, n_moments, moments,
    QuantilePoint,
    fit_mean_quantile, fit_mode_quantile, fit_median_quantile,
    @qp, @qp_ll, @qp_l, @qp_m, @qp_u, @qp_uu,
    @qs_cf90, @qs_cf95,
    Σstar, σstar,
    autocor_effective, effective_n_cor, count_for_lag, count_forlags, sem_cor, var_cor

import Random: rand, rand!
import Base: sum, size
import Statistics: mean
import StatsBase: coef, autocor

using StatsBase, Distributions, Missings, MissingStrategies
using BandedMatrices, MappedArrays, RecursiveArrayTools, FillArrays, OffsetArrays
using StaticArrays, LinearAlgebra
using Random
using SimpleTraits

# general fitting to statistics
include("fitstats.jl")

# lognormal fitting
include("lognormal.jl")

# normal fitting
include("normal.jl")

# logitnormal fitting
include("logitnormal.jl")

# standard error of the mean of correlated normal variables
include("semcor.jl")

end # module


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using RigidBodySim, Documenter, RigidBodyDynamics
import DiffEqBase, MeshCatMechanisms

makedocs(
    doctest = true,
    modules = [RigidBodySim],
    checkdocs = :none, # because of the reexported functionality
    linkcheck = true,
    format = :html,
    root = @__DIR__,
    sitename ="RigidBodySim.jl",
    authors = "Twan Koolen and contributors.",
    pages = [
        "Home" => "index.md",
        "Quick start guide" => "quickstart.md",
        "Details" => "details.md"
    ],
    html_prettyurls = parse(Bool, get(ENV, "CI", "false")),
)

deploydocs(
    repo = "github.com/JuliaRobotics/RigidBodySim.jl.git"
)


text:
module GeodesyGPX

import Geodesy
import GPX

"""
    LatLon(point::GPXPoint)

Latitude and longitude co-ordinates, constructed from a GPX point. Elevation is simply ignored.
"""
Geodesy.LatLon(point::GPX.GPXPoint) = Geodesy.LatLon(point.lat, point.lon)

"""
    LLA(point::GPXPoint)

Latitude, longitude, and altitude co-ordinates, constructed from a GPX point.
"""
Geodesy.LLA(point::GPX.GPXPoint) = Geodesy.LLA(point.lat, point.lon, point.ele)

"""
	ECEF(point::GPXPoint, [datum = wgs84])

Earth-Centered-Earth-Fixed (ECEF) coordinates, constructed from a GPX point.
"""
Geodesy.ECEF(point::GPX.GPXPoint, datum=wgs84) = Geodesy.ECEF(Geodesy.LLA(point), datum)

"""
    euclidean_distance(a::GPXPoint, b::GPXPoint, [datum = wgs84])

The straight-line distance between GPX points `a` and `b`.
"""
Geodesy.euclidean_distance(a::GPX.GPXPoint, b::GPX.GPXPoint, datum=wgs84) = Geodesy.euclidean_distance(Geodesy.ECEF.((a, b))..., datum)
end


text:
#///////////////////////////////////////
#// File Name: cost_gpu.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: GPU (CUDA) Cost model for Risk Sensitive Stochastic SAC
#///////////////////////////////////////

using CUDA
using LinearAlgebra

# Instantaneous Cost
# # Position
function kernel_instant_position_cost!(out::AbstractArray{Float32, 2},
                                       ex_array::AbstractArray{Float32, 3},
                                       target_pos_array::AbstractArray{Float32, 2},
                                       Cep::AbstractArray{Float32, 2})
    # out : (num_controls, T)
    # ex_array : (num_controls, T, 4)
    # target_pos_array : (T, 2)
    # Cep : (2, 2)

    ii = (blockIdx().x - 1)*blockDim().x + threadIdx().x; # dimension for nominal control candidates
    jj = (blockIdx().y - 1)*blockDim().y + threadIdx().y; # dimension for timesteps

    if (ii <= size(out, 1)) && (jj <= size(out, 2))
        position_error_x = ex_array[ii, jj, 1] - target_pos_array[jj, 1];
        position_error_y = ex_array[ii, jj, 2] - target_pos_array[jj, 2];
        out[ii, jj] = 0.5f0*(position_error_x*Cep[1, 1]*position_error_x +
                             position_error_x*Cep[1, 2]*position_error_y +
                             position_error_y*Cep[2, 1]*position_error_x +
                             position_error_y*Cep[2, 2]*position_error_y);
    end
    return nothing
end

function instant_position_cost(ex_array::AbstractArray{Float32, 3},
                               target_pos_array::AbstractArray{Float32, 2},
                               param::CostParameter;
                               threads::NTuple{2, Int}=(8, 32))

    # out : (num_controls, total_timesteps - 1)
    # ex_array : (num_controls, total_timesteps, 4)
    # target_pos_array : (total_timesteps, 2)

    Cep = cu(param.Cep);
    out = CuArray{Float32, 2}(undef, size(ex_array, 1), size(ex_array, 2) - 1)
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    blocks = (numblocks_x, numblocks_y)
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_instant_position_cost!(out, ex_array[:, 1:end-1, :], target_pos_array[1:end-1, :], Cep)
    end
    return out
end

# # Control
function kernel_instant_control_cost!(out::AbstractArray{Float32, 2},
                                      u_array::AbstractArray{Float32, 3},
                                      Cu::AbstractArray{Float32, 2})
    # out : (num_controls, T)
    # u_array : (num_controls, T, 2)
    # Cu : (2, 2)

    ii = (blockIdx().x - 1)*blockDim().x + threadIdx().x; # dimension for nominal control candidates
    jj = (blockIdx().y - 1)*blockDim().y + threadIdx().y; # dimension for timesteps

    if (ii <= size(out, 1)) && (jj <= size(out, 2))
        out[ii, jj] = 0.5f0*(u_array[ii, jj, 1]*Cu[1, 1]*u_array[ii, jj, 1] +
                             u_array[ii, jj, 1]*Cu[1, 2]*u_array[ii, jj, 2] +
                             u_array[ii, jj, 2]*Cu[2, 1]*u_array[ii, jj, 1] +
                             u_array[ii, jj, 2]*Cu[2, 2]*u_array[ii, jj, 2])
    end
    return nothing
end

function instant_control_cost(u_array::AbstractArray{Float32, 3},
                              param::CostParameter;
                              threads::NTuple{2, Int}=(8, 32))

    # u_array : (num_controls, total_timesteps - 1)
    # out : (num_controls, total_timesteps - 1)

    Cu = cu(param.Cu);
    out = CuArray{Float32, 2}(undef, size(u_array, 1), size(u_array, 2))
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    blocks = (numblocks_x, numblocks_y)
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_instant_control_cost!(out, u_array, Cu)
    end
    return out
end

# # Collision
function kernel_instant_collision_cost!(out::AbstractArray{Float32, 3},
                                        ex_array::AbstractArray{Float32, 3},
                                        ap_array::AbstractArray{Float32, 4},
                                        time_idx_ap_array::AbstractArray{Int32, 1},
                                        control_idx_ex_array::AbstractArray{Int32, 1},
                                        α_col::Float32,
                                        λ_col::Float32)

    # out : (num_samples*num_controls, T, num_ado_agents)
    # ex_array : (num_controls, T, 4)
    # ap_array : (num_samples*num_controls, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (T)
    # control_idx_ex_array : (num_samples*num_controls)

    ii = (blockIdx().x - 1)*blockDim().x + threadIdx().x; # dimension for samples (and nominal control candidates)
    jj = (blockIdx().y - 1)*blockDim().y + threadIdx().y; # dimension for timesteps
    kk = (blockIdx().z - 1)*blockDim().z + threadIdx().z; # dimension for ado agents

    if (ii <= size(out, 1)) && (jj <= size(out, 2)) && (kk <= size(out, 3))
        ll = control_idx_ex_array[ii];
        mm = time_idx_ap_array[jj];

        squared_dist = (ex_array[ll, jj, 1] - ap_array[ii, mm, kk, 1])^2 +
                       (ex_array[ll, jj, 2] - ap_array[ii, mm, kk, 2])^2;
        out[ii, jj, kk] = α_col*CUDA.exp(-1/(2*λ_col)*squared_dist);
    end
    return nothing
end

function instant_collision_cost(ex_array::AbstractArray{Float32, 3},
                                ap_array::AbstractArray{Float32, 4},
                                time_idx_ap_array::AbstractArray{Int32, 1},
                                control_idx_ex_array::AbstractArray{Int32, 1},
                                param::CostParameter;
                                threads::NTuple{3, Int}=(32, 4, 2))

    # out : (num_samples*num_controls, total_timesteps - 1, num_ado_agents)
    # ex_array : (num_controls, total_timesteps, 4)
    # ap_array : (num_samples*num_controls, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (total_timesteps)
    # control_idx_ex_array : (num_samples*num_controls)

    α_col = Float32(param.α_col);
    λ_col = Float32(param.λ_col);
    out = CuArray{Float32, 3}(undef, size(ap_array, 1), size(ex_array, 2) - 1, size(ap_array, 3));
    if size(out, 3) == 0 # If no ado agents exist, return all zeros.
        return cu(zeros(size(out)));
    end
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    numblocks_z = ceil(Int, size(out, 3)/threads[3]);
    blocks = (numblocks_x, numblocks_y, numblocks_z);
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_instant_collision_cost!(out, ex_array[:, 1:end-1, :],
                                                                           ap_array,
                                                                           time_idx_ap_array[1:end-1],
                                                                           control_idx_ex_array,
                                                                           α_col, λ_col);
    end
    return out
end

# Terminal Cost
# # Position
function kernel_terminal_position_cost!(out::AbstractArray{Float32, 2},
                                        ex_array::AbstractArray{Float32, 3},
                                        target_pos_array::AbstractArray{Float32, 2},
                                        Cep::AbstractArray{Float32, 2})

    # out : (num_controls, 1)
    # ex_array : (num_controls, 1, 4)
    # target_pos_array : (1, 2)
    # Cep : (2, 2)

    kernel_instant_position_cost!(out, ex_array, target_pos_array, Cep)
    return nothing
end

function terminal_position_cost(ex_array::AbstractArray{Float32, 3},
                                target_pos_array::AbstractArray{Float32, 2},
                                param::CostParameter;
                                threads::NTuple{2, Int}=(256, 1))

    # out : (num_controls, 1)
    # ex_array : (num_controls, total_timesteps, 4)
    # target_pos_array : (total_timesteps, 2)

    Cep = cu(param.β_pos*param.Cep);
    out = CuArray{Float32, 2}(undef, size(ex_array, 1), 1)
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    blocks = (numblocks_x, numblocks_y)
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_terminal_position_cost!(out, ex_array[:, end:end, :], target_pos_array[end:end, :], Cep)
    end
    return out
end

# # Collision
function kernel_terminal_collision_cost!(out::AbstractArray{Float32, 3},
                                         ex_array::AbstractArray{Float32, 3},
                                         ap_array::AbstractArray{Float32, 4},
                                         time_idx_ap_array::AbstractArray{Int32, 1},
                                         control_idx_ex_array::AbstractArray{Int32, 1},
                                         α_col::Float32,
                                         λ_col::Float32)

    # out : (num_samples*num_controls, 1, num_ado_agents)
    # ex_array : (num_controls, 1, 4)
    # ap_array : (num_samples*num_controls, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (1)
    # control_idx_ex_array : (num_samples*num_controls)

    kernel_instant_collision_cost!(out, ex_array, ap_array,
                                   time_idx_ap_array, control_idx_ex_array,
                                   α_col, λ_col)
    return nothing
end

function terminal_collision_cost(ex_array::AbstractArray{Float32, 3},
                                 ap_array::AbstractArray{Float32, 4},
                                 time_idx_ap_array::AbstractArray{Int32, 1},
                                 control_idx_ex_array::AbstractArray{Int32, 1},
                                 param::CostParameter;
                                 threads::NTuple{3, Int}=(64, 1, 4))

    # out : (num_samples*num_controls, 1, num_ado_agents)
    # ex_array : (num_controls, 1, 4)
    # ap_array : (num_samples*num_controls, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (1)
    # control_idx_ex_array : (num_samples*num_controls)

    α_col = Float32(param.β_col*param.α_col);
    λ_col = Float32(param.λ_col);
    out = CuArray{Float32, 3}(undef, size(ap_array, 1), 1, size(ap_array, 3));
    if size(out, 3) == 0 # If no ado agents exist, return all zeros.
        return cu(zeros(size(out)));
    end
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    numblocks_z = ceil(Int, size(out, 3)/threads[3]);
    blocks = (numblocks_x, numblocks_y, numblocks_z);
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_terminal_collision_cost!(out, ex_array[:, end:end, :],
                                                                            ap_array,
                                                                            time_idx_ap_array[end:end],
                                                                            control_idx_ex_array,
                                                                            α_col, λ_col)
    end
    return out
end

# Instantaneous Cost Gradient (with respect to ego robot state)
# # Position
function kernel_instant_position_cost_gradient!(out::AbstractArray{Float32, 2},
                                                ex_array::AbstractArray{Float32, 2},
                                                target_pos_array::AbstractArray{Float32, 2},
                                                Cep::AbstractArray{Float32, 2})
    # out : (T, 4)
    # ex_array : (T, 4)
    # target_pos_array : (T, 2)
    # Cep : (2, 2)

    ii = (blockIdx().x - 1)*blockDim().x + threadIdx().x; # dimension for timesteps

    if (ii <= size(out, 1))
        position_error_x = ex_array[ii, 1] - target_pos_array[ii, 1];
        position_error_y = ex_array[ii, 2] - target_pos_array[ii, 2];
        out[ii, 1] = Cep[1, 1]*position_error_x + Cep[1, 2]*position_error_y;
        out[ii, 2] = Cep[2, 1]*position_error_x + Cep[2, 2]*position_error_y;
        out[ii, 3] = 0.0f0;
        out[ii, 4] = 0.0f0;
    end
    return nothing
end

function instant_position_cost_gradient(ex_array::AbstractArray{Float32, 2},
                                        target_pos_array::AbstractArray{Float32, 2},
                                        param::CostParameter;
                                        threads::Int=64)

    # out : (total_timesteps - 1, 4)
    # ex_array : (total_timesteps, 4)
    # target_pos_array : (total_timesteps, 2)

    Cep = cu(param.Cep);
    out = CuArray{Float32, 2}(undef, size(ex_array, 1) - 1, 4)
    threads = threads;
    blocks = ceil(Int, size(out, 1)/threads);
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_instant_position_cost_gradient!(out, ex_array[1:end-1, :], target_pos_array[1:end-1, :], Cep)
    end
    return out
end

## Collision
function kernel_instant_collision_cost_gradient!(out::AbstractArray{Float32, 4},
                                                 ex_array::AbstractArray{Float32, 2},
                                                 ap_array::AbstractArray{Float32, 4},
                                                 time_idx_ap_array::AbstractArray{Int32, 1},
                                                 α_col::Float32,
                                                 λ_col::Float32)

    # out : (num_samples, T, num_ado_agents, 4)
    # ex_array : (T, 4)
    # ap_array : (num_samples, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (T)

    ii = (blockIdx().x - 1)*blockDim().x + threadIdx().x; # dimension for samples
    jj = (blockIdx().y - 1)*blockDim().y + threadIdx().y; # dimension for timesteps
    kk = (blockIdx().z - 1)*blockDim().z + threadIdx().z; # dimension for ado agents

    if (ii <= size(out, 1)) && (jj <= size(out, 2)) && (kk <= size(out, 3))
        mm = time_idx_ap_array[jj];

        squared_dist = (ex_array[jj, 1] - ap_array[ii, mm, kk, 1])^2 +
                       (ex_array[jj, 2] - ap_array[ii, mm, kk, 2])^2;
        collision_cost = α_col*CUDA.exp(-1/(2*λ_col)*squared_dist);
        out[ii, jj, kk, 1] =
                    collision_cost/λ_col*(ap_array[ii, mm, kk, 1] - ex_array[jj, 1]);
        out[ii, jj, kk, 2] =
                    collision_cost/λ_col*(ap_array[ii, mm, kk, 2] - ex_array[jj, 2]);
        out[ii, jj, kk, 3] = 0.0f0;
        out[ii, jj, kk, 4] = 0.0f0;
    end
    return nothing
end

function instant_collision_cost_gradient(ex_array::AbstractArray{Float32, 2},
                                         ap_array::AbstractArray{Float32, 4},
                                         time_idx_ap_array::AbstractArray{Int32, 1},
                                         param::CostParameter;
                                         threads::NTuple{3, Int}=(32, 4, 2))

    # out : (num_samples, total_timesteps - 1, num_ado_agents, 4)
    # ex_array : (total_timesteps, 4)
    # ap_array : (num_samples, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (total_timesteps)

    α_col = Float32(param.α_col);
    λ_col = Float32(param.λ_col);
    out = CuArray{Float32, 4}(undef, size(ap_array, 1), size(ex_array, 1) - 1, size(ap_array, 3), 4);
    if size(out, 3) == 0 # If no ado agents exist, return all zeros.
        return cu(zeros(size(out)));
    end
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    numblocks_z = ceil(Int, size(out, 3)/threads[3]);
    blocks = (numblocks_x, numblocks_y, numblocks_z);
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_instant_collision_cost_gradient!(out, ex_array[1:end-1, :],
                                                                                ap_array,
                                                                                time_idx_ap_array[1:end-1],
                                                                                α_col, λ_col);
    end
    return out
end

# Terminal Cost Gradient (with respect to ego robot state)
# # Position
function kernel_terminal_position_cost_gradient!(out::AbstractArray{Float32, 2},
                                                 ex_array::AbstractArray{Float32, 2},
                                                 target_pos_array::AbstractArray{Float32, 2},
                                                 Cep::AbstractArray{Float32, 2})
    # out : (1, 4)
    # ex_array : (1, 4)
    # target_pos_array : (1, 2)
    # Cep : (2, 2)

    kernel_instant_position_cost_gradient!(out, ex_array, target_pos_array, Cep)
    return nothing
end

function terminal_position_cost_gradient(ex_array::AbstractArray{Float32, 2},
                                         target_pos_array::AbstractArray{Float32, 2},
                                         param::CostParameter;
                                         threads::Int=1)

    # out : (1, 4)
    # ex_array : (1, 4)
    # target_pos_array : (1, 2)

    Cep = cu(param.β_pos*param.Cep);
    out = CuArray{Float32, 2}(undef, 1, 4)
    threads = threads;
    blocks = ceil(Int, size(out, 1)/threads);
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_terminal_position_cost_gradient!(out, ex_array[end:end, :], target_pos_array[end:end, :], Cep)
    end
    return out
end

# # Collision
function kernel_terminal_collision_cost_gradient!(out::AbstractArray{Float32, 4},
                                                  ex_array::AbstractArray{Float32, 2},
                                                  ap_array::AbstractArray{Float32, 4},
                                                  time_idx_ap_array::AbstractArray{Int32, 1},
                                                  α_col::Float32,
                                                  λ_col::Float32)

    # out : (num_samples, 1, num_ado_agents, 4)
    # ex_array : (1, 4)
    # ap_array : (num_samples, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (1)

    kernel_instant_collision_cost_gradient!(out, ex_array, ap_array, time_idx_ap_array, α_col, λ_col)
    return nothing
end

function terminal_collision_cost_gradient(ex_array::AbstractArray{Float32, 2},
                                          ap_array::AbstractArray{Float32, 4},
                                          time_idx_ap_array::AbstractArray{Int32, 1},
                                          param::CostParameter;
                                          threads::NTuple{3, Int}=(64, 1, 4))

    # out : (num_samples, 1, num_ado_agents, 4)
    # ex_array : (1, 4)
    # ap_array : (num_samples, prediction_steps + 1, num_ado_agents, 2)
    # time_idx_ap_array : (1)

    α_col = Float32(param.β_col*param.α_col);
    λ_col = Float32(param.λ_col);
    out = CuArray{Float32, 4}(undef, size(ap_array, 1), 1, size(ap_array, 3), 4);
    if size(out, 3) == 0 # If no ado agents exist, return all zeros.
        return cu(zeros(size(out)));
    end
    threads = threads;
    numblocks_x = ceil(Int, size(out, 1)/threads[1]);
    numblocks_y = ceil(Int, size(out, 2)/threads[2]);
    numblocks_z = ceil(Int, size(out, 3)/threads[3]);
    blocks = (numblocks_x, numblocks_y, numblocks_z);
    CUDA.@sync begin
        @cuda threads=threads blocks=blocks kernel_terminal_collision_cost_gradient!(out, ex_array[end:end, :],
                                                                                ap_array,
                                                                                time_idx_ap_array[end:end],
                                                                                α_col, λ_col);
    end
    return out
end


text:
using Demo
using Test

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


text:
# Recursively create benchmark groups
function mkgrp!(base, path)
    path = map(string, path)
    curr = base
    for k in path
        if !(k in keys(curr))
            curr[k] = BenchmarkGroup()
        end
        curr = curr[k]
    end
    curr
end

function mkbench!(f, base, path)
    grp = mkgrp!(base, path[1:end-1])
    grp[path[end]] = f()
end


text:
## Conflict-Based Allocation
const TaskUtil = NamedTuple{(:task, :util), Tuple{String,Float64}}

@with_kw mutable struct SCoBAHighLevelNode
    task_allocation::Dict{String,TaskUtil}      = Dict{String,TaskUtil}()    # Agent name |--> task name
    considered_tasks::Dict{String,Set{String}}  = Dict{String,Set{String}}()  # Agent name |--> tasks to consider
    util::Float64                               = 0.0
    id::Int64                                   = 0
end
Base.isless(hln1::SCoBAHighLevelNode, hln2::SCoBAHighLevelNode) = hln1.util < hln2.util

@with_kw mutable struct SCoBAAlgorithm
    allocation::GenericAllocation
    heap::MutableBinaryMaxHeap{SCoBAHighLevelNode}    = MutableBinaryMaxHeap{SCoBAHighLevelNode}()
    num_total_conflicts::Int64                        = 0
end

# Inspired by CBS but with allocation instead of path search
function coordinate_allocation!(solver::SCoBAAlgorithm, initial_task_allocation::Dict{String,TaskUtil},
                                initial_considered_tasks::Dict{String,Set{String}}, init_util::Float64,
                                success_prob_fn::Function, util_val_fn::Function,
                                downtime::Float64)

    # Create start node with initial info
    start = SCoBAHighLevelNode(task_allocation=initial_task_allocation,
                                considered_tasks=initial_considered_tasks,
                                util=init_util)

    push!(solver.heap, start)
    id = 1

    # Iterative case: Pop top of heap
    # Check if any conflict in allocation
    # If no, then return solution. Otherwise, generate child nodes that re-assign tasks accordingly
    while ~(isempty(solver.heap))

        P = pop!(solver.heap)

        # Check for conflicting allocations
        task_to_agent = Dict{String,Vector{String}}()
        for (agent, task_util) in P.task_allocation

            task = task_util.task

            if ~(haskey(task_to_agent, task))
                task_to_agent[task] = [agent]
            else
                push!(task_to_agent[task], agent)
            end
        end # for (a,t) in task_alloc

        # Now look for tasks assigned to multiple agents
        num_conflicts = 0

        tasks_to_ignore = Set{String}(collect(keys(task_to_agent)))

        for (task, alloc_agents) in task_to_agent

            # If alloc_agents has more than one, then we need to generate constraints etc.
            if length(alloc_agents) > 1

                # For each of the assigned agents, generate a new high-level node
                # where the over-assigned task is removed from each of the other assigned agents
                for agt in alloc_agents

                    new_node = deepcopy(P)
                    new_node.id = id

                    # Create a set of all agents other than agt
                    other_agents_set = Set{String}(alloc_agents)
                    setdiff!(other_agents_set, [agt])

                    # For each of other agents, remove task from considered_tasks and recompute low level tree
                    for other_agt in other_agents_set

                        # We KNOW that other_agt will not be assigned to task again
                        # Reduce util and remove task allocation
                        new_node.util -= new_node.task_allocation[other_agt].util
                        delete!(new_node.task_allocation, other_agt)

                        # Generate new tree for other agent
                        setdiff!(new_node.considered_tasks[other_agt], tasks_to_ignore)

                        generate_search_tree!(solver.allocation, other_agt,
                                              new_node.considered_tasks[other_agt],
                                              success_prob_fn, util_val_fn, downtime)

                        tree = solver.allocation.agent_prop_set[other_agt].tree
                        dec_idx = get_next_attempt_idx(tree)

                        # If agent is assigned to a new task, update hlnode
                        if dec_idx != -1
                            dec_node = tree.nodes[dec_idx]

                            # Update util and allocation of new node
                            new_node.util += dec_node.util
                            new_node.task_allocation[other_agt] = (task=dec_node.task_name, util=dec_node.util)

                        end # if dec_idx
                    end # for other_agt

                    push!(solver.heap, new_node)
                    id += 1

                end # for agt in alloc_agents

                num_conflicts += 1
                solver.num_total_conflicts += 1

                if solver.num_total_conflicts > solver.allocation.conflict_threshold
                    @warn "High level conflict thresh exceeded!"
                    return P.task_allocation
                end

            end # if length > 1
        end # for (t,a) in task_to_agent

        if num_conflicts == 0 # This should be guaranteed to happen?
            @debug "No more conflicts!"
            return P.task_allocation
        end
    end # while ~isempty(heap)

    # Return empty allocation?
    @info "No coordinated allocation!"
    return Dict{String,TaskUtil}()

end # function coordinate_allocation