text:
function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::DeviceModel,
    devices::IS.FlattenIteratorWrapper{V},
) where {V <: PSY.Component}
    for ff in get_feedforwards(model)
        @debug "arguments" ff V _group = LOG_GROUP_FEEDFORWARDS_CONSTRUCTION
        add_feedforward_arguments!(container, model, devices, ff)
    end
    return
end

function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::ServiceModel,
    service::V,
) where {V <: PSY.AbstractReserve}
    for ff in get_feedforwards(model)
        @debug "arguments" ff V _group = LOG_GROUP_FEEDFORWARDS_CONSTRUCTION
        contributing_devices = get_contributing_devices(model)
        add_feedforward_arguments!(container, model, contributing_devices, ff)
    end
    return
end

function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::DeviceModel,
    devices::IS.FlattenIteratorWrapper{T},
    ff::AbstractAffectFeedforward,
) where {T <: PSY.Component}
    parameter_type = get_default_parameter_type(ff, T)
    source_key = get_optimization_container_key(ff)
    add_parameters!(container, parameter_type, source_key, model, devices)
    return
end

function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::ServiceModel{SR},
    contributing_devices::Vector{T},
    ff::AbstractAffectFeedforward,
) where {T <: PSY.Component, SR <: PSY.AbstractReserve}
    parameter_type = get_default_parameter_type(ff, SR)
    source_key = get_optimization_container_key(ff)
    add_parameters!(container, parameter_type, source_key, model, contributing_devices)
    return
end

function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::DeviceModel,
    devices::IS.FlattenIteratorWrapper{T},
    ff::SemiContinuousFeedforward,
) where {T <: PSY.Component}
    parameter_type = get_default_parameter_type(ff, T)
    source_key = get_optimization_container_key(ff)
    add_parameters!(container, parameter_type, source_key, model, devices)
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        parameter_type,
        devices,
        model,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        parameter_type,
        devices,
        model,
    )
    return
end

function add_feedforward_arguments!(
    container::OptimizationContainer,
    model::DeviceModel,
    devices::IS.FlattenIteratorWrapper{T},
    ff::EnergyTargetFeedforward,
) where {T <: PSY.Component}
    parameter_type = get_default_parameter_type(ff, T)
    source_key = get_optimization_container_key(ff)
    add_parameters!(container, parameter_type, source_key, model, devices)
    # Enabling this FF requires the addition of an extra variable
    add_variables!(container, EnergyShortageVariable, devices, get_formulation(model)())
    return
end


text:
# Autogenerated wrapper script for SFML_jll for i686-linux-musl-cxx11
export libsfml_audio, libsfml_graphics, libsfml_network, libsfml_system, libsfml_window

using Libglvnd_jll
using Ogg_jll
using FLAC_jll
using FreeType2_jll
using libvorbis_jll
using Xorg_libXrandr_jll
using Xorg_libX11_jll
using OpenAL_jll
using eudev_jll
JLLWrappers.@generate_wrapper_header("SFML")
JLLWrappers.@declare_library_product(libsfml_audio, "libsfml-audio.so.2.5")
JLLWrappers.@declare_library_product(libsfml_graphics, "libsfml-graphics.so.2.5")
JLLWrappers.@declare_library_product(libsfml_network, "libsfml-network.so.2.5")
JLLWrappers.@declare_library_product(libsfml_system, "libsfml-system.so.2.5")
JLLWrappers.@declare_library_product(libsfml_window, "libsfml-window.so.2.5")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, Ogg_jll, FLAC_jll, FreeType2_jll, libvorbis_jll, Xorg_libXrandr_jll, Xorg_libX11_jll, OpenAL_jll, eudev_jll)
    JLLWrappers.@init_library_product(
        libsfml_audio,
        "lib/libsfml-audio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_graphics,
        "lib/libsfml-graphics.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_network,
        "lib/libsfml-network.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_system,
        "lib/libsfml-system.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_window,
        "lib/libsfml-window.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
############
# Potentials
############

"""
    abstract type AbstractMaxwellPotential{T<:AbstractKernel,S}

Type that wraps around `IntegralOperator`. It is used to define
the Maxwell potentials `MaxwellEFIEPotential` and `MaxwellMFIEPotential`.
"""
abstract type AbstractMaxwellPotential{T<:AbstractKernel,S} end

# same interface as IntegralPotential
kernel(p::AbstractMaxwellPotential) = kernel(p.ip)
surface(p::AbstractMaxwellPotential) = surface(p.ip)
kernel_type(p::AbstractMaxwellPotential) = kernel_type(p.ip)
Base.getindex(p::AbstractMaxwellPotential,σ::AbstractVector) = Base.getindex(p.ip,σ)

# EFIE
struct MaxwellEFIEPotential{T,S} <: AbstractMaxwellPotential{T,S}
    ip::IntegralPotential{T,S}
end
MaxwellEFIEPotential(op::Maxwell,surf) = MaxwellEFIEPotential(SingleLayerPotential(op,surf))

# MFIE
struct MaxwellMFIEPotential{T,S} <: AbstractMaxwellPotential{T,S}
    ip::IntegralPotential{T,S}
end
MaxwellMFIEPotential(op::Maxwell,surf) = MaxwellMFIEPotential(DoubleLayerPotential(op,surf))
# overload MFIE potential evaluation
Base.getindex(p::MaxwellMFIEPotential,σ::AbstractVector) = Base.getindex(p.ip,-ncross(σ))

############
# Operators
############

function maxwell_dim(pde,X,args...;kwargs...)
    S, D = single_doublelayer_dim(pde,X,args...;kwargs...)
    N = diagonal_ncross_matrix(X) |> DiscreteOp
    T = N*S     # EFIE
    K = -N*D*N  # MFIE
    return T, K
end
maxwell_efie_dim(pde,X,args...;kwargs...) = maxwell_dim(pde,X,args...;kwargs...)[1]
maxwell_mfie_dim(pde,X,args...;kwargs...) = maxwell_dim(pde,X,args...;kwargs...)[2]

function assemble_maxwell_indirect_cfie(nmesh, α, β, K, T; exterior=true)
    σ = exterior ? 0.5 : -0.5
    N = diagonal_ncross_matrix(nmesh)  |> DiscreteOp
    J = diagonal_jacobian_matrix(nmesh) |> DiscreteOp
    dualJ = diagonal_dualjacobian_matrix(nmesh) |> DiscreteOp
    cfie = dualJ*(σ*α*I + α*K + β*T*N)*J
    return cfie
end

# TODO: delete
function assemble_dim_nystrom_matrix_Luiz(mesh, α, β, D, S)
    N, J, dualJ = diagonal_ncross_and_jacobian_matrices(mesh)
    n_qnodes = length(dofs(mesh))
    M = Matrix{SMatrix{2,2,ComplexF64,4}}(undef, n_qnodes, n_qnodes)
    M .= dualJ*(0.5*α*I + N*(-α*D + β*S)*N)*J
    return M
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
include("nnutil.jl")
include("CommonFuncs.jl")
th = 1e-2 #dm
force_scale = 5.0/th #50
tid = 200
fiber_size = 1
porder = 2
stress_scale = 1.0e3
strain_scale = 1

if length(ARGS) == 4
   global tid = parse(Int64, ARGS[1])
   global force_scale = parse(Float64, ARGS[2])
   global fiber_size = parse(Int64, ARGS[3])
   global porder = parse(Int64, ARGS[4])
end
printstyled("force_scale=$force_scale, tid=$tid, fiber_size=$fiber_size, porder=$porder\n", color=:green)

if Sys.MACHINE=="x86_64-apple-darwin18.6.0"
    matplotlib.use("macosx")
end
np = pyimport("numpy")

Δt = 2.0e-4
NT = 100
T = Δt * NT #50.0 #500.0

# Step 1: make domain and globaldata
Lx, Ly = 1.0, 0.5
# number of elements in each directions
nxc, nyc = 10,5
nx, ny =  nxc*fiber_size, nyc*fiber_size
ngp = 3
prop = Dict("name"=> "PlaneStress","rho"=> 4.5, "E"=> 1e+4, "nu"=> 0.2)
Dc = PlaneStress(prop).H

nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid, nx, ny, porder, Lx, Ly; force_scale=force_scale)
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, ngp))
    end
end

ndofs = 2
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
setGeometryPoints!(domain, npoints, node_to_point)

state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

# Step 2: assemble global data mass matrix
assembleMassMatrix!(globdat, domain)

# Step 3: compute_boundary_info to obtain time-depedent displacement and acceleration ubd, abd
ts = ExplicitSolverTime(Δt, NT)

ubd, abd = compute_boundary_info(domain, globdat, ts)

# Step 4: compute_external_force to obtain total external Fext
Fext = compute_external_force(domain, globdat, ts)

# Step 5 get initial state, velo, acceleration with SovlerInitial
d0, v0, a0 = SolverInitial(Δt, globdat, domain)

# Step 6: Solve the equation
#d, v, a =  ExplicitSolver(globdat, domain, d0, v0, a0, Δt, NT, nn_law, Fext, ubd, abd; strain_type="small")

#TODO assume no prestress/prestrain
ngpt = getNGauss(domain)
σ0, ε0 = zeros(ngpt, 3), zeros(ngpt, 3)
d, v, a, σ, ε =  ExplicitSolver(globdat, domain, d0, v0, a0, σ0, ε0, Δt, NT, nn_associated_plasticity, Fext, ubd, abd; strain_type="small")

# Step 7: Start optimization
sess = Session(); init(sess)
dobs = matread("Data/order$(porder)/data$(tid).mat")["d"]
loss = sum((d - dobs)^2*1.0e6)
@show run(sess, loss)

@info "Start"
BFGS!(sess, loss)

d_, v_, a_ = run(sess, [d,v,a])
matwrite("Data/order$(porder)/data-nn$(tid).mat", Dict("d"=>d_, "v"=>v_, "a"=>a_))


text:
module TestFind

using DocumentationOverview
using Test

function test_find()
    apis = DocumentationOverview.find(DocumentationOverview)
    allnames = sort!([api.name for api in apis])
    @test allnames == [:API, :find, :list, :list_md, :table, :table_md]
end

end  # module


text:
n = 10 # rand(1:10)
@test matrixdepot("hilb", n) == matrixdepot("hilb", n, n)
@test matrixdepot("hilb" , Float16, n) ≈ matrixdepot("hilb", Float32, n)
x = Float64[1:n;]
y = x .- 1
@test matrixdepot("hilb", n) ≈ matrixdepot("cauchy", x, y)
println("'hilb' passed test...")


text:
# Autogenerated wrapper script for SymEngine_jll for powerpc64le-linux-gnu-cxx11
export libsymengine

using GMP_jll
using MPFR_jll
using MPC_jll
JLLWrappers.@generate_wrapper_header("SymEngine")
JLLWrappers.@declare_library_product(libsymengine, "libsymengine.so.0.8")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, MPFR_jll, MPC_jll)
    JLLWrappers.@init_library_product(
        libsymengine,
        "lib/libsymengine.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ************************************************************************
# FAUST Architecture File
# Copyright (C) 2021 GRAME, Centre National de Creation Musicale
# ---------------------------------------------------------------------

# This is sample code. This file is provided as an example of minimal
# FAUST architecture file. Redistribution and use in source and binary
# forms, with or without modification, in part or in full are permitted.
# In particular you can create a derived work of this FAUST architecture
# and distribute that work under terms of your choice.

# This sample code is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# ************************************************************************

# UI architecture
abstract type UI end

# One can override the behavior by defining another set of methods
# that takes a different concrete UI type.
# Only the needed methods have to be implemented in subtypes.

# -- widget's layouts
function openTabBox!(ui_interface::UI, label::String)
end
function openHorizontalBox!(ui_interface::UI, label::String)
end
function openVerticalBox!(ui_interface::UI, label::String)
end
function closeBox!(ui_interface::UI)
end

# -- active widgets
function addButton!(ui_interface::UI, label::String, param::Symbol)
end
function addCheckButton!(ui_interface::UI, label::String, param::Symbol)
end
function addHorizontalSlider!(ui_interface::UI, label::String, param::Symbol, init::FAUSTFLOAT, min::FAUSTFLOAT, max::FAUSTFLOAT, step::FAUSTFLOAT)
end
function addVerticalSlider!(ui_interface::UI, label::String, param::Symbol, init::FAUSTFLOAT, min::FAUSTFLOAT, max::FAUSTFLOAT, step::FAUSTFLOAT)
end
function addNumEntry!(ui_interface::UI, label::String, param::Symbol, init::FAUSTFLOAT, min::FAUSTFLOAT, max::FAUSTFLOAT, step::FAUSTFLOAT)
end

# -- passive widgets
function addHorizontalBargraph!(ui_interface::UI, label::String, param::Symbol, min::FAUSTFLOAT, max::FAUSTFLOAT)
end
function addVerticalBargraph!(ui_interface::UI, label::String, param::Symbol, min::FAUSTFLOAT, max::FAUSTFLOAT)
end

# -- soundfiles
function addSoundfile!(ui_interface::UI, label::String, filename::String, soundfile::Symbol)
end

# -- metadata declarations
function declare!(ui_interface::UI, param::Symbol, key::String, val::String)
end


text:
# tests with the 5 taxon networks write to parenthetical format
# Claudia March 2015
# need to draw in dendroscope to compare, still not automatic

# types in "types.jl"
include("../src/types.jl")
include("../src/functions.jl")

# needed modules:
using Base.Collections # for updateInCycle with priority queue

tests = ["C","F","G","H","J","D","E","I"];
wrong = String[];

function whichtree(t::String)
    if(t == "tree")
        tree = "(((6:0.1,4:1.5)1:0.2,7:0.2)5:0.1,8:0.1,10:0.1);" # normal tree
    elseif(t == "C")
        tree = "((((6:0.1,4:1.5),(7:0.2)11#H1),11#H1),8:0.1,10:0.1);" # Case C: bad triangle II
    elseif(t == "F")
        tree = "(((6:0.1,(4)11#H1)1:0.2,(11#H1,7))5:0.1,8:0.1,10:0.1);" # Case F: bad diamond I
    elseif(t == "G")
        tree = "((((6:0.1,4:1.5)1:0.2,(7)11#H1)5:0.1,(11#H1,8)),10:0.1);" # Case G
    elseif(t == "H")
        tree = "((((6,4),#H1),7),(8)#H1,10);" # Case H
    elseif(t == "J")
        tree = "((((6)#H1,4),7),8,(#H1,10));" # Case J
    elseif(t == "D")
        tree = "((((6,4))#H1,(#H1,7)),8,10);" # Case D Bad triangle I
    elseif(t == "E")
        tree = "(((((8,10))#H1,7),#H1),6,4);" # Case E Bad triangle I
    elseif(t == "I")
        tree = "((((8,10))#H1,7),6,(4,#H1));" # Case I Bad diamond II
    else
        error("not a known 5 taxon network case")
    end
    return tree
end

tests = ["F","G","H","J","I"];
t="I"
for t in tests
    println("running $(t)")
    net = nothing;
    tree = whichtree(t)
    f = open("prueba_tree.txt","w")
    write(f,tree)
    close(f)
    net = readTopologyUpdate("prueba_tree.txt");
    printEdges(net)
    printNodes(net)
    written = writeTopologyLevel1(net)
    written2 = writeTopologyLevel1(net,true)
    tree
end


text:
export mouseover, mouse_selection, mouseposition, hovered_scene
export select_rectangle, select_line, select_point

"""
    mouseover(scene::SceneLike, plots::AbstractPlot...)

Returns true if the mouse currently hovers any of `plots`.
"""
function mouseover(scene::SceneLike, plots::AbstractPlot...)
    p, idx = mouse_selection(scene)
    return p in flatten_plots(plots)
end

"""
    onpick(f, scene::SceneLike, plots::AbstractPlot...)

Calls `f(idx)` whenever the mouse is over any of `plots`.
`idx` is an index, e.g. when over a scatter plot, it will be the index of the
hovered element
"""
function onpick(f, scene::SceneLike, plots::AbstractPlot...; range=1)
    fplots = flatten_plots(plots)
    args = range == 1 ? (scene,) : (scene, range)
    map_once(events(scene).mouseposition) do mp
        p, idx = mouse_selection(args...)
        (p in fplots) && f(idx)
        return Consume(false)
    end
end

"""
    mouse_selection(scene::Scene)

Returns the plot that is under the current mouse position
"""
function mouse_selection(scene::SceneLike)
    pick(scene, events(scene).mouseposition[])
end
function mouse_selection(scene::SceneLike, range)
    pick(scene, events(scene).mouseposition[], range)
end

function flatten_plots(x::Atomic, plots = AbstractPlot[])
    if isempty(x.plots)
        push!(plots, x)
    else
        flatten_plots(x.plots, plots)
    end
    plots
end

function flatten_plots(x::Combined, plots = AbstractPlot[])
    for elem in x.plots
        flatten_plots(elem, plots)
    end
    plots
end

function flatten_plots(array, plots = AbstractPlot[])
    for elem in array
        flatten_plots(elem, plots)
    end
    plots
end

"""
    mouse_in_scene(scene::Scene[, priority = 0])

Returns a new node that is true whenever the cursor is inside the given scene.

See also: [`is_mouseinside`](@ref)
"""
function mouse_in_scene(scene::SceneLike; priority = Int8(0))
    p = rootparent(scene)
    output = Node(Vec2(0.0))
    on(events(scene).mouseposition, priority = priority) do mp
        output[] = Vec(mp) .- minimum(pixelarea(scene)[])
        return Consume(false)
    end
    output
end

"""
    pick(scene, x, y)

Returns the plot under pixel position `(x, y)`.
"""
function pick(scene::SceneLike, x::Number, y::Number)
    return pick(scene, Vec{2, Float64}(x, y))
end

"""
    pick(scene::Scene, xy::VecLike)

Return the plot under pixel position xy.
"""
function pick(scene::SceneLike, xy)
    screen = getscreen(scene)
    screen === nothing && return (nothing, 0)
    pick(scene, screen, Vec{2, Float64}(xy))
end

"""
    pick(scene::Scene, xy::VecLike, range)

Return the plot closest to xy within a given range.
"""
function pick(scene::SceneLike, xy, range)
    screen = getscreen(scene)
    screen === nothing && return (nothing, 0)
    pick_closest(scene, screen, xy, range)
end

# The backend may handle this more optimally
function pick_closest(scene::SceneLike, screen, xy, range)
    w, h = widths(screen)
    ((1.0 <= xy[1] <= w) && (1.0 <= xy[2] <= h)) || return (nothing, 0)
    x0, y0 = max.(1, floor.(Int, xy .- range))
    x1, y1 = min.([w, h], floor.(Int, xy .+ range))
    dx = x1 - x0; dy = y1 - y0

    picks = pick(scene, screen, IRect2D(x0, y0, dx, dy))

    min_dist = range^2
    selected = (0, 0)
    x, y =  xy .+ 1 .- Vec2f0(x0, y0)
    for i in 1:dx, j in 1:dy
        d = (x-i)^2 + (y-j)^2
        if (d < min_dist) && (picks[i, j][1] != nothing)
            min_dist = d
            selected = (i, j)
        end
    end

    return selected == (0, 0) ? (nothing, 0) : picks[selected[1], selected[2]]
end

"""
    pick_sorted(scene::Scene, xy::VecLike, range)

Return all `(plot, index)` pairs in a `(xy .- range, xy .+ range)` region
sorted by distance to `xy`.
"""
function pick_sorted(scene::SceneLike, xy, range)
    screen = getscreen(scene)
    screen === nothing && return Tuple{AbstractPlot, Int}[]
    pick_sorted(scene, screen, xy, range)
end

function pick_sorted(scene::SceneLike, screen, xy, range)
    w, h = widths(screen)
    if !((1.0 <= xy[1] <= w) && (1.0 <= xy[2] <= h))
        return Tuple{AbstractPlot, Int}[]
    end
    x0, y0 = max.(1, floor.(Int, xy .- range))
    x1, y1 = min.([w, h], floor.(Int, xy .+ range))
    dx = x1 - x0; dy = y1 - y0

    picks = pick(scene, screen, IRect2D(x0, y0, dx, dy))

    selected = filter(x -> x[1] != nothing, unique(vec(picks)))
    distances = [range^2 for _ in selected]
    x, y =  xy .+ 1 .- Vec2f0(x0, y0)
    for i in 1:dx, j in 1:dy
        if picks[i, j][1] != nothing
            d = (x-i)^2 + (y-j)^2
            i = findfirst(isequal(picks[i, j]), selected)
            if i === nothing
                @warn "This shouldn't happen..."
            elseif distances[i] > d
                distances[i] = d
            end
        end
    end

    idxs = sortperm(distances)
    permute!(selected, idxs)
    return selected
end

"""
    pick(scene::Scene, rect::IRect2D)

Return all `(plot, index)` pairs within the given rect. The rect must be within
screen boundaries.
"""
function pick(scene::SceneLike, rect::IRect2D)
    screen = getscreen(scene)
    screen === nothing && return Tuple{AbstractPlot, Int}[]
    return pick(scene, screen, rect)
end

"""
    screen_relative(scene, pos)

Normalizes mouse position `pos` relative to the screen rectangle.
"""
function screen_relative(scene::Scene, mpos)
    return Point2f0(mpos) .- Point2f0(minimum(pixelarea(scene)[]))
end

"""
    mouseposition(scene = hovered_scene())

Return the current position of the mouse in _data coordinates_ of the
given `scene`.

By default uses the `scene` that the mouse is currently hovering over.
"""
function mouseposition(scene = hovered_scene())
    to_world(scene, mouseposition_px(scene))
end

function mouseposition_px(scene = hovered_scene())
    screen_relative(
        scene,
        events(scene).mouseposition[]
    )
end

"""
    hovered_scene()

Returns the `scene` that the mouse is currently hovering over.

Properly identifies the scene for a plot with multiple sub-plots.
"""
hovered_scene() = error("hoevered_scene is not implemented yet.")

"""
    select_rectangle(scene; kwargs...) -> rect

Interactively select a rectangle on a 2D `scene` by clicking the left mouse button,
dragging and then un-clicking. The function returns an **observable** `rect` whose
value corresponds to the selected rectangle on the scene. In addition the function
_plots_ the selected rectangle on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted rectangle disappears.

The value of the returned observable is updated **only** when the user un-clicks
(i.e. when the final value of the rectangle has been decided) and only if the
rectangle has area > 0.

The `kwargs...` are propagated into `lines!` which plots the selected rectangle.
"""
function select_rectangle(scene; blocking = false, priority = 2, strokewidth = 3.0, kwargs...)
    key = Mouse.left
    waspressed = Node(false)
    rect = Node(FRect(0, 0, 1, 1)) # plotted rectangle
    rect_ret = Node(FRect(0, 0, 1, 1)) # returned rectangle

    # Create an initially hidden rectangle
    plotted_rect = poly!(
        scene, rect, raw = true, visible = false, color = RGBAf0(0, 0, 0, 0), strokecolor = RGBAf0(0.1, 0.1, 0.8, 0.5), strokewidth = strokewidth, kwargs...,
    )

    on(events(scene).mousebutton, priority=priority) do event
        if event.button == key
            if event.action == Mouse.press && is_mouseinside(scene)
                mp = mouseposition(scene)
                waspressed[] = true
                plotted_rect[:visible] = true # start displaying
                rect[] = FRect(mp, 0.0, 0.0)
                return Consume(blocking)
            end
        end
        if !(event.button == key && event.action == Mouse.press)
            if waspressed[] # User has selected the rectangle
                waspressed[] = false
                r = absrect(rect[])
                w, h = widths(r)
                if w > 0.0 && h > 0.0 # Ensure that the rectangle has non0 size.
                    rect_ret[] = r
                end
            end
            # always hide if not the right key is pressed
            plotted_rect[:visible] = false # make the plotted rectangle invisible
            return Consume(blocking)
        end

        return Consume(false)
    end
    on(events(scene).mouseposition, priority=priority) do event
        if waspressed[]
            mp = mouseposition(scene)
            mini = minimum(rect[])
            rect[] = FRect(mini, mp - mini)
            return Consume(blocking)
        end
        return Consume(false)
    end

    return rect_ret
end

"""
    select_line(scene; kwargs...) -> line

Interactively select a line (typically an arrow) on a 2D `scene` by clicking the left mouse button,
dragging and then un-clicking. Return an **observable** whose value corresponds
to the selected line on the scene. In addition the function
_plots_ the line on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted line disappears.

The value of the returned line is updated **only** when the user un-clicks
and only if the selected line has non-zero length.

The `kwargs...` are propagated into `lines!` which plots the selected line.
"""
function select_line(scene; blocking = false, priority = 2, kwargs...)
    key = Mouse.left
    waspressed = Node(false)
    line = Node([Point2f0(0,0), Point2f0(1,1)])
    line_ret = Node([Point2f0(0,0), Point2f0(1,1)])
    # Create an initially hidden  arrow
    plotted_line = lines!(
        scene, line; visible = false, color = RGBAf0(0.1, 0.1, 0.8, 0.5),
        linewidth = 4, kwargs...,
    )

    on(events(scene).mousebutton, priority=priority) do event
        if event.button == key && is_mouseinside(scene)
            mp = mouseposition(scene)
            if event.action == Mouse.press
                waspressed[] = true
                plotted_line[:visible] = true  # start displaying
                line[][1] = mp
                line[][2] = mp
                line[] = line[]
                return Consume(blocking)
            end
        end
        if !(event.button == key && event.action == Mouse.press)
            if waspressed[] # User has selected the rectangle
                waspressed[] = false
                if line[][1] != line[][2]
                    line_ret[] = copy(line[])
                end
            end
            plotted_line[:visible] = false
            return Consume(blocking)
        end
        return Consume(false)
    end
    on(events(scene).mouseposition, priority=priority) do event
        if waspressed[]
            mp = mouseposition(scene)
            line[][2] = mp
            line[] = line[] # actually update observable
            return Consume(blocking)
        end
        return Consume(false)
    end

    return line_ret
end

"""
    select_point(scene; kwargs...) -> point

Interactively select a point on a 2D `scene` by clicking the left mouse button,
dragging and then un-clicking. Return an **observable** whose value corresponds
to the selected point on the scene. In addition the function
_plots_ the point on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted point disappears.

The value of the returned point is updated **only** when the user un-clicks.

The `kwargs...` are propagated into `scatter!` which plots the selected point.
"""
function select_point(scene; blocking = false, priority=2, kwargs...)
    key = Mouse.left
    pmarker = Circle(Point2f0(0, 0), Float32(1))
    waspressed = Node(false)
    point = Node([Point2f0(0,0)])
    point_ret = Node(Point2f0(0,0))
    # Create an initially hidden  arrow
    plotted_point = scatter!(
        scene, point; visible = false, marker = pmarker, markersize = 20px,
        color = RGBAf0(0.1, 0.1, 0.8, 0.5), kwargs...,
    )

    on(events(scene).mousebutton, priority=priority) do event
        if event.button == key && is_mouseinside(scene)
            mp = mouseposition(scene)
            if event.action == Mouse.press
                waspressed[] = true
                plotted_point[:visible] = true  # start displaying
                point[][1] = mp
                point[] = point[]
                return Consume(blocking)
            end
        end
        if !(event.button == key && event.action == Mouse.press)
            if waspressed[] # User has selected the rectangle
                waspressed[] = false
                point_ret[] = copy(point[][1])
            end
            plotted_point[:visible] = false
            return Consume(blocking)
        end
        return Consume(false)
    end
    on(events(scene).mouseposition, priority=priority) do event
        if waspressed[]
            mp = mouseposition(scene)
            point[][1] = mp
            point[] = point[] # actually update observable
            return Consume(blocking)
        end
        return Consume(false)
    end

    return point_ret
end


text:
import Pkg
Pkg.activate(joinpath(@__DIR__, "../../"))  # environment: CytofResearch

using BSON
using CSV
using CytofResearch
import DataFrames; const DF = DataFrames
using Distributions
using Flux, Flux.Tracker
import Random

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using Test

@testset "Select.jl" begin
    include("Select.jl")
end


text:
function estimate_norm(mat; tol=1e-4, itmax = 1000)
    v = rand(size(mat,2))

    v = v/norm(v)
    itermin = 3
    i = 1
    σold = 1
    σnew = 1
    while (norm(sqrt(σold)-sqrt(σnew))/norm(sqrt(σold)) > tol || i < itermin) && i < itmax
        σold = σnew
        w = mat*v
        x = adjoint(mat)*w
        σnew = norm(x)
        v = x/norm(x)
        i += 1
    end
    return sqrt(σnew)
end

function estimate_reldifference(hmat, refmat; tol=1e-4)
    #if size(hmat) != size(refmat)
    #    error("Dimensions of matrices do not match")
    #end

    v = rand(size(hmat,2))

    v = v/norm(v)
    itermin = 3
    i = 1
    σold = 1
    σnew = 1
    while norm(sqrt(σold)-sqrt(σnew))/norm(sqrt(σold)) > tol || i < itermin
        σold = σnew
        w = hmat*v - refmat*v
        x = adjoint(hmat)*w - adjoint(refmat)*w
        σnew = norm(x)
        v = x/norm(x)
        i += 1
    end

    norm_refmat = estimate_norm(refmat, tol=tol)

    return sqrt(σnew)/norm_refmat
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
# Standard deviation of a vector without Bessel's correction
function Base.std{T<:Real}(data::Vector{T}, avg::Real, n::Int)
    s = 0.0
    for i = 1:n
        @inbounds z = data[i] - avg
        s += z * z
    end
    sqrt(s / n)
end

# Per-column standard deviation of a matrix without Bessel's correction
Base.std{T<:Real}(data::Matrix{T}, avg::Vector{T}, rows::Int, cols::Int) =
    vec([std(data[:,j], avg[j], rows) for j in 1:cols])

Base.std{T<:Real}(data::Matrix{T}, avg::Matrix{T}, rows::Int, cols::Int) =
    std(data, vec(avg), rows, cols)

normalize{T<:Real}(v::Vector{T}) = vec(v) / sum(v)

normalize{T<:Real}(v::Matrix{T}) = normalize(vec(v))

function coskewness{T<:Real}(data::Matrix{T}, w::Vector{T}; standardize=false, bias=0)
    if sum(w) != 1
        w = normalize(w)
    end
    length(w) == size(data, 2) || throw(DimensionMismatch("Inconsistent array lengths."))
    first(w' * coskew(data, standardize=standardize, flatten=true, bias=bias) * kron(w, w))
end

coskewness{T<:Real}(data::Matrix{T}; standardize=false, bias=0) =
    coskewness(data, ones(size(data, 2)), standardize=standardize, bias=bias)

coskewness{T<:Real}(data::Matrix{T}, w::Matrix{T}; standardize=false, bias=0) =
    coskewness(data, vec(w), standardize=standardize, bias=bias)

function cokurtosis{T<:Real}(data::Matrix{T}, w::Vector{T}; standardize=false, bias=0)
    if sum(w) != 1
        w = normalize(w)
    end
    length(w) == size(data, 2) || throw(DimensionMismatch("Inconsistent array lengths."))
    first(w' * cokurt(data, standardize=standardize, flatten=true, bias=bias) * kron(kron(w, w), w))
end

cokurtosis{T<:Real}(data::Matrix{T}; standardize=false, bias=0) =
    cokurtosis(data, ones(size(data, 2)), standardize=standardize, bias=bias)

cokurtosis{T<:Real}(data::Matrix{T}, w::Matrix{T}; standardize=false, bias=0) =
    cokurtosis(data, vec(w), standardize=standardize, bias=bias)


text:
"""
A Trial represents a guess we've tried and the result we got back.
"""
struct Trial
    candidate::Candidate
    result::Result
end


text:
"""
    function AMat(ped)
---
"""
function AMat(ped)
end

function read_ped(file)
    @info "Reading pedigree from $file"
end

"""
    function code_ped()
---
"""
function code_ped()
end

function read_srt_ped(file)
    ped = Int32[]
    for line in eachline(file)
        (lstrip(line)[1] == '#') && continue
        append!(ped, parse.(Int32, split(line)[[1, 2]]))
    end
    reshape(ped, 2, :)'
end

function AInv(ped)
end

"""
    function A_matrix(ped)
---
Give a matrix of 2 columns:
- Sire
- Dam

With row number as ID number, this function return a **A** matrix.
`Sire`s and `Dam`s are of integers less than their offspring `ID`.
This function cannot deal very large pedigree, i.e.,
with 32G momory, this function can only calculate a pedigree with
~90k ID.
With half precision, this can be ~130k.
With double precision, this function can only handle ~65k ID.
"""
function A_matrix(ped)
    @warn "Not right, needs debug"
    n = size(ped, 1)
    A = zeros(Float32, n, n)
    for i in 1:n
        ip, im = ped[i, :]      # ID i's pa and ma
        A[i, i] = 1.            # below avoid swap
        ip >0 && im >0 && (A[i, i] += .5(A[ip, im] + A[im, ip]))
        for j in i+1:n          # upper triangular
            jp, jm = ped[j, :]
            jp > 0 && (A[i, j] += .5A[i, jp])
            jm > 0 && (A[i, j] += .5A[i, jm])
            A[j, i] = A[i, j]
        end
        i % 100 == 0 && (print("\t$i"))
    end
    A
end

"""
    function kinship(ped, i, j)
---
This function is handy if just to calculate relationship of a few (pairs of) ID.
It can also speed up by adding `Thread.@threads` before your pair loop.
"""
function kinship(ped, i, j)
    (i == 0 || j == 0) && return 0
    ipa, ima = ped[i, :]          # used both for below and the last
    i == j && (return 1 + .5kinship(ped, ipa, ima))
    if i < j
        jpa, jma = ped[j, :]
        return .5(kinship(ped, i, jpa) + kinship(ped, i, jma))
    end
    return .5(kinship(ped, j, ipa) + kinship(ped, j, ima))
end

"""
    function kinship(ped, i, j, A)
---
This function is handy if just to calculate relationship of a few (pairs of) ID.
It can also speed up by adding `Thread.@threads` before your pair loop.
"""
function kinship(ped, i, j, A)
    (i == 0 || j == 0) && return 0

    !iszero(A[i, j]) && return A[i, j] # already calculated

    ip, im = ped[i, :]          # i's pa and ma
    i == j && begin
        A[i, i] = 1 + .5kinship(ped, ip, im, A)
        return A[i, i]
    end

    i < j && begin
        jp, jm = ped[j, :]      # j's pa and ma
        A[i, j] = .5(kinship(ped, i, jp, A) + kinship(ped, i, jm, A))
        return A[i, j]
    end

    A[i, j] = .5(kinship(ped, j, ip, A) + kinship(ped, j, im, A))
    return A[i, j]
end

"""
    function sprtA(ped, idx)
---
Part `A` of pedigree `ped` using sparse arrays.
List `idx` is the relationship ID to be calculated.

## Notes:
1. `idx` is better sorted.
2. `idx`'s biggest number should not be bigger than `size(ped, 1)`.
3. `ped` is of `Array{Int32, 2}`, with `pa` and `ma` column, and row number is `id` number.
"""
function sprtA(ped, idx)
    N = size(ped, 1)
    A = spzeros(N, N)
    n = length(idx)
    for i in n:-1:1
        x = idx[i]
        for j in i:-1:1
            y = idx[j]
            print("\r",
                  lpad("Row progress: $(n-i+1)/$n", 25),
                  lpad("ID: $x", 15),
                  ", ",
                  lpad("Col progress: $(i-j+1)/$i", 25),
                  lpad("ID: $y", 15))
            A[x, y] = A[y, x] = kinship(ped, x, y, A)
        end
    end
    println()
    Matrix(A[idx, idx])
end

#=
  case 'F':
  case 'f':
    for(auto id{1}; id<=nid; ++id)
      cout << Amat(id, id, ped, mid) - 1. <<'\n';
    break;
  case 'D':
  case 'd':
  case 'T':
  case 't':
    pma[0] = -1;

    cout << nid << '\n';

    for(auto id{1}; id<=nid; ++id){
      const auto&[pa, ma] = ped[id];
      if(pma.find(pa) == pma.end()) pma[pa] = Amat(pa, pa, ped, mid) - 1.;
      if(pma.find(ma) == pma.end()) pma[ma] = Amat(ma, ma, ped, mid) - 1.;
      cout << 1./(.5 - .25*(pma[pa] + pma[ma])) << '\n';
    }

    for(auto id{1}; id<=nid; ++id){
      const auto&[pa, ma] = ped[id];
      if(pa<ma) putT(pa, ma, id);
      else      putT(ma, pa, id);
    }
    break;
=#

function assignCSC(i, j, d, k, Ti, Tj, Td)
    k += 1
    Ti[k] = i
    Tj[k] = j
    Td[k] = d
    k
end

function A_inverse(ped)
    @info "Calculating A inversed"
    N = size(ped, 1)
    Ai = spzeros(N, N)
    A = spzeros(N, N)

    prt = collect(Set(ped))
    npm = length(prt)
    ele = zeros(npm)
    @showprogress "Inbreeding coefficients of parents" for i in 1:npm
        ele[i] = kinship(ped, prt[i], prt[i], A) - 1
    end

    pm = Dict(prt .=> ele)

    @info "Computing T & D"
    D = zeros(N)
    Ti = zeros(Int, 3N)
    Tj = zeros(Int, 3N)
    Td = zeros(3N)
    k = 0
    @showprogress "Matrix D and T" for i in 1:N
        pa, ma = ped[i, :]
        D[i] = 1. / (.5 - .25(pm[pa] + pm[ma]))
        pa > 0 && (k = assignCSC(i, pa, -.5, k,  Ti, Tj, Td))
        ma > 0 && (k = assignCSC(i, ma, -.5, k,  Ti, Tj, Td))
        k = assignCSC(i, i, 1., k,  Ti, Tj, Td)
    end
    T = sparse(Ti[1:k], Tj[1:k], Td[1:k])

    @info "Constructing A^-1"
    T'Diagonal(D)T
end


text:
# JSON Microbenchmarks

using JSON
using BenchmarkTools

const suite = BenchmarkGroup()

suite["print"] = BenchmarkGroup(["serialize"])
suite["pretty-print"] = BenchmarkGroup(["serialize"])

immutable CustomListType
    x::Int
    y::Float64
    z::Union{CustomListType, Void}
end

immutable CustomTreeType
    x::String
    y::Union{CustomTreeType, Void}
    z::Union{CustomTreeType, Void}
end

list(x) = x == 0 ? nothing : CustomListType(1, 1.0, list(x - 1))
tree(x) = x == 0 ? nothing : CustomTreeType("!!!", tree(x - 1), tree(x - 1))

const micros = Dict(
    "integer" => 88,
    "float" => -88.8,
    "ascii" => "Hello World!",
    "ascii-1024" => "x" ^ 1024,
    "unicode" => "ສະບາຍດີຊາວໂລກ!",
    "unicode-1024" => "ℜ" ^ 1024,
    "bool" => true,
    "null" => nothing,
    "flat-homogenous-array-16" => collect(1:16),
    "flat-homogenous-array-1024" => collect(1:1024),
    "heterogenous-array" => [
        1, 2, 3, 7, "A", "C", "E", "N", "Q", "R", "Shuttle to Grand Central"],
    "nested-array-16^2" => [collect(1:16) for _ in 1:16],
    "nested-array-16^3" => [[collect(1:16) for _ in 1:16] for _ in 1:16],
    "small-dict" => Dict(
        :a => :b, :c => "💙💙💙💙💙💙", :e => 10, :f => Dict(:a => :b)),
    "flat-dict-128" => Dict(zip(collect(1:128), collect(1:128))),
    "date" => Date(2016, 08, 09),
    "matrix-16" => eye(16),
    "custom-list-128" => list(128),
    "custom-tree-8" => tree(8))

for (k, v) in micros
    io = IOBuffer()
    suite["print"][k] = @benchmarkable JSON.print($(IOBuffer()), $v)
    suite["pretty-print"][k] = @benchmarkable JSON.print(
        $(IOBuffer()), $v, 4)
end


text:
module NormalMeanVarianceTest

using Test
using ReactiveMP

@testset "NormalMeanVariance" begin

    @testset "Constructor" begin
        @test NormalMeanVariance <: NormalDistributionsFamily
        @test NormalMeanVariance <: UnivariateNormalDistributionsFamily

        @test NormalMeanVariance()         == NormalMeanVariance{Float64}(0.0, 1.0)
        @test NormalMeanVariance(1.0)      == NormalMeanVariance{Float64}(1.0, 1.0)
        @test NormalMeanVariance(1.0, 2.0) == NormalMeanVariance{Float64}(1.0, 2.0)
        @test NormalMeanVariance(1)        == NormalMeanVariance{Float64}(1.0, 1.0)
        @test NormalMeanVariance(1, 2)     == NormalMeanVariance{Float64}(1.0, 2.0)
        @test NormalMeanVariance(1.0, 2)   == NormalMeanVariance{Float64}(1.0, 2.0)
        @test NormalMeanVariance(1, 2.0)   == NormalMeanVariance{Float64}(1.0, 2.0)
        @test NormalMeanVariance(1f0)      == NormalMeanVariance{Float32}(1f0, 1f0)
        @test NormalMeanVariance(1f0, 2f0) == NormalMeanVariance{Float32}(1f0, 2f0)
        @test NormalMeanVariance(1f0, 2)   == NormalMeanVariance{Float32}(1f0, 2f0)
        @test NormalMeanVariance(1f0, 2.0) == NormalMeanVariance{Float64}(1.0, 2.0)

        @test eltype(NormalMeanVariance())         === Float64
        @test eltype(NormalMeanVariance(0.0))      === Float64
        @test eltype(NormalMeanVariance(0.0, 1.0)) === Float64
        @test eltype(NormalMeanVariance(0))        === Float64
        @test eltype(NormalMeanVariance(0, 1))     === Float64
        @test eltype(NormalMeanVariance(0.0, 1))   === Float64
        @test eltype(NormalMeanVariance(0, 1.0))   === Float64
        @test eltype(NormalMeanVariance(0f0))      === Float32
        @test eltype(NormalMeanVariance(0f0, 1f0)) === Float32
        @test eltype(NormalMeanVariance(0f0, 1.0)) === Float64
    end

    @testset "Stats methods" begin

        dist1 = NormalMeanVariance(0.0, 1.0)

        @test mean(dist1)         === 0.0
        @test median(dist1)       === 0.0
        @test mode(dist1)         === 0.0
        @test weightedmean(dist1) === 0.0
        @test var(dist1)          === 1.0
        @test std(dist1)          === 1.0
        @test cov(dist1)          === 1.0
        @test invcov(dist1)       === 1.0
        @test precision(dist1)    === 1.0
        @test entropy(dist1)      ≈ 1.41893853320467
        @test pdf(dist1, 1.0)     ≈ 0.24197072451914337
        @test pdf(dist1, -1.0)    ≈ 0.24197072451914337
        @test pdf(dist1, 0.0)     ≈ 0.3989422804014327
        @test logpdf(dist1, 1.0)  ≈ -1.4189385332046727
        @test logpdf(dist1, -1.0) ≈ -1.4189385332046727
        @test logpdf(dist1, 0.0)  ≈ -0.9189385332046728

        dist2 = NormalMeanVariance(1.0, 1.0)

        @test mean(dist2)         === 1.0
        @test median(dist2)       === 1.0
        @test mode(dist2)         === 1.0
        @test weightedmean(dist2) === 1.0
        @test var(dist2)          === 1.0
        @test std(dist2)          === 1.0
        @test cov(dist2)          === 1.0
        @test invcov(dist2)       === 1.0
        @test precision(dist2)    === 1.0
        @test entropy(dist2)      ≈ 1.41893853320467
        @test pdf(dist2, 1.0)     ≈ 0.3989422804014327
        @test pdf(dist2, -1.0)    ≈ 0.05399096651318806
        @test pdf(dist2, 0.0)     ≈ 0.24197072451914337
        @test logpdf(dist2, 1.0)  ≈ -0.9189385332046728
        @test logpdf(dist2, -1.0) ≈ -2.9189385332046727
        @test logpdf(dist2, 0.0)  ≈ -1.4189385332046727

        dist3 = NormalMeanVariance(1.0, 2.0)

        @test mean(dist3)         === 1.0
        @test median(dist3)       === 1.0
        @test mode(dist3)         === 1.0
        @test weightedmean(dist3) === inv(2.0)
        @test var(dist3)          === 2.0
        @test std(dist3)          === sqrt(2.0)
        @test cov(dist3)          === 2.0
        @test invcov(dist3)       === inv(2.0)
        @test precision(dist3)    === inv(2.0)
        @test entropy(dist3)      ≈ 1.7655121234846454
        @test pdf(dist3, 1.0)     ≈ 0.28209479177387814
        @test pdf(dist3, -1.0)    ≈ 0.1037768743551487
        @test pdf(dist3, 0.0)     ≈ 0.21969564473386122
        @test logpdf(dist3, 1.0)  ≈ -1.2655121234846454
        @test logpdf(dist3, -1.0) ≈ -2.2655121234846454
        @test logpdf(dist3, 0.0)  ≈ -1.5155121234846454

    end

    @testset "Base methods" begin
        @test convert(NormalMeanVariance{Float32}, NormalMeanVariance()) === NormalMeanVariance{Float32}(0f0, 1f0)
        @test convert(NormalMeanVariance{Float64}, NormalMeanVariance(0.0, 10.0)) == NormalMeanVariance{Float64}(0.0, 10.0)
        @test convert(NormalMeanVariance{Float64}, NormalMeanVariance(0.0, 0.1)) == NormalMeanVariance{Float64}(0.0, 0.1)
        @test convert(NormalMeanVariance{Float64}, 0, 1) == NormalMeanVariance{Float64}(0.0, 1.0)
        @test convert(NormalMeanVariance{Float64}, 0, 10) == NormalMeanVariance{Float64}(0.0, 10.0)
        @test convert(NormalMeanVariance{Float64}, 0, 0.1) == NormalMeanVariance{Float64}(0.0, 0.1)
        @test convert(NormalMeanVariance, 0, 1) == NormalMeanVariance{Float64}(0.0, 1.0)
        @test convert(NormalMeanVariance, 0, 10) == NormalMeanVariance{Float64}(0.0, 10.0)
        @test convert(NormalMeanVariance, 0, 0.1) == NormalMeanVariance{Float64}(0.0, 0.1)
    end

    @testset "vague" begin
        d1 = vague(NormalMeanVariance)

        @test typeof(d1) <: NormalMeanVariance
        @test mean(d1)      == 0.0
        @test var(d1)       == ReactiveMP.huge
    end

    @testset "prod" begin

        @test prod(ProdAnalytical(), NormalMeanVariance(-1, 1), NormalMeanVariance(1, 1)) ≈ NormalMeanVariance(0.0, 0.5)
        @test prod(ProdAnalytical(), NormalMeanVariance(-1, 2), NormalMeanVariance(1, 4)) ≈ NormalMeanVariance(-1/3, 4/3)
        @test prod(ProdAnalytical(), NormalMeanVariance(2, 2), NormalMeanVariance(0, 10)) ≈ NormalMeanVariance(5/3, 5/3)

    end

end

end


text:
#randomgrowth.jl
using Winston

function random_growth(M, N, q)
    G = zeros(M, N)
    T = 1000

    G[1,1] = true
    # update the possible sets
    imagesc((0,N), (M,0), G)
    for t = 1:T
        sets = next_possible_squares(G)
        ## Highlights all the possible squares
        for i = 1:length(sets)
            idx = sets[i]::(Int,Int)
            G[idx[1], idx[2]] = 0.25
        end
        imagesc((0,N), (M,0), G)
        sleep(.01)

        ## Actual growth
        for i = 1:length(sets)
            ison = 0.5 * (rand() > (1-q))
            if ison > 0
                idx = sets[i]::(Int,Int)
                G[idx[1], idx[2]] = ison
            end
        end
        imagesc((0,N), (M,0), G)
        G[G .== 0.5] = 1
        G[G .== 0.25] = 0
        sleep(.01)
    end
    G
end

function next_possible_squares(G)
    M, N = size(G)::(Int,Int)
    sets = Array((Int,Int),0)
    for ii = 1:M
        for jj = 1:N
            if G[ii, jj] == 0
                if  (ii == 1 && jj > 1 && G[ii, jj-1] == 1) ||
                    (jj == 1 && ii > 1 && G[ii-1, jj] == 1) ||
                    (ii > 1 && jj > 1 && G[ii-1, jj] == 1 && G[ii, jj-1] == 1)
                    push!(sets, (ii,jj))
                end
            end
        end
    end
    sets
end


text:
using Pkg
Pkg.activate("/home/dabajabaza/jianguoyun/Workspace/SmartParser/")
using SmartParser

function count_frequent_tokens(fns)
    stat = Dict{String,Int}()
    good_key(x) = !any(a->occursin(a,x), ["£",".", "\"", "*****","-----","%%%%%","+++++"])
    for f in fns
        if isfile(f)
            S0 = read(f,String) |> preprocess_raw_input
            MS0 = mask(S0)
            patts, code = tokenize0(MS0)
            codeinv = revert(code)
            for pt in patts
                for p in pt
                    ky = codeinv[p]
                    if p!=0 && good_key(ky)
                        increaseindex!(stat, ky)
                    end
                end
            end
        end
    end
    return stat
end

ALL_FILES = []

for RT ∈ ["/data/ReO3_phonon_configs","/data/ScF3_phonon_configs",]
    fns=vcat([[ "$RT/$fd/$ff/$fout" for ff in readdir("$RT/$fd") if isdir("$RT/$fd/$ff")
                                        for fout in readdir("$RT/$fd/$ff") if endswith(fout,".out") ]
                                            for fd ∈ readdir(RT) if isdir("$RT/$fd") ]...)
    ALL_FILES = [ALL_FILES; fns]
end

##

STATS = count_frequent_tokens(ALL_FILES)
STATS_SORTV = sort([k=>v for (k,v) ∈ STATS if v>1], by=last) ;
L = length(STATS_SORTV)

##

CODE_ALL = Dict(k=>L-i+1 for (i,(k,v)) ∈ enumerate(STATS_SORTV)) ;
CODE_ALL |> print

##

for f in ALL_FILES
    try rm("$(f).replaced.1.txt") catch ; end
    try rm("$(f).replaced.2.txt") catch ; end
    try rm("$(f).replaced.3.txt") catch ; end
    try rm("$(f).replaced.4.txt") catch ; end
end

##

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
module test_bukdu_server

using Test
using Bukdu
using HTTP
using JSON

Bukdu.start(8190)
sleep(0)
Bukdu.stop()

routes() do
    post("/") do conn::Conn
        render(JSON, conn.params.k + 1)
    end
    plug(Plug.Parsers, [:json])
end

Bukdu.start(8190, host="127.0.0.1")
r = HTTP.post("http://127.0.0.1:8190/", ["Content-Type"=>"application/json"]; body=JSON.json((k=2,)))
sleep(0)
@test HTTP.header(r, "Server") == Bukdu.server_info
@test HTTP.header(r, "Content-Type") == "application/json; charset=utf-8"
@test HTTP.header(r, "Content-Length") == "1"
@test r.body == Vector{UInt8}("3")

Plug.Loggers.config[:error_stackframes_range] = 1:2
@test_throws HTTP.ExceptionRequest.StatusError HTTP.post("http://127.0.0.1:8190/"; body=JSON.json((k=3,)))
@test_throws HTTP.ExceptionRequest.StatusError HTTP.post("http://127.0.0.1:8190/", ["Content-Type"=>"application/json"]; body=JSON.json(3))

Bukdu.stop()

Routing.empty!()

version_line = first(filter(line -> startswith(line, "version"), readlines(normpath(pathof(Bukdu), "..", "..", "Project.toml"))))
@test occursin(string(Bukdu.BUKDU_VERSION), version_line)

end # module test_bukdu_server


text:
using Pkg;Pkg.activate(".")
using Test, Plots, BenchmarkTools, Revise, VortexDistributions
gr(transpose=true)

## field with 2 points per healing length
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex

## make vortex dipole (cant use this for single vortex accuracy)
# Periodic boundary conditions
psi = Torus(copy(psi0),x,y)

xp = 1.7
yp = 2.1
vp = PointVortex(xp,yp,1)

xn = 3.24
yn = 10.21
vn = PointVortex(xn,yn,-1)

dip = ScalarVortex([vp;vn])
periodic_dipole!(psi,dip)

heatmap(x,y,angle.(psi.ψ))

vfound = findvortices(psi)
@show vfound

## best test: find a known vortex
# Periodic boundary conditions
psi = Torus(copy(psi0),x,y)

xp = 1.133
yp = 2.787
vp = PointVortex(xp,yp,1)

sp = ScalarVortex(vp)

vortex!(psi,sp)
vfound = findvortices(psi);@show vortex_array(vfound)

heatmap(x,y,angle.(psi.ψ))
## test zoom_grid

psiw,xw,yw = zoom_grid(psi.ψ,x,y,vortex_array(vp)[1:2]...)
heatmap(xw,yw,angle.(psiw))

## readme example
@btime findvortices(psi)


text:
# test clone frequency calculations
clonetype = [0,1,1]
clonefreq = [0.903,0.058,0.034]
@test isapprox(CancerSeqSim.calculateclonefreq(clonefreq, clonetype), [0.995, 0.058, 0.034])

clonefreq = [0.036,0.013,0.732,0.002]
clonetype = [0,0,1,3]
@test isapprox(CancerSeqSim.calculateclonefreq(clonefreq, clonetype), [0.77 ,0.013, 0.734, 0.002])

clonefreq = [0.053,0.02,0.002,0.924]
clonetype = [0,1,1,3]
@test isapprox(CancerSeqSim.calculateclonefreq(clonefreq, clonetype), [0.999, 0.02, 0.926, 0.924])

clonefreq = [0.036,0.003,0.91,0.008,0.004]
clonetype = [0,1,1,1,3]
@test isapprox(CancerSeqSim.calculateclonefreq(clonefreq, clonetype), [0.961, 0.003, 0.914, 0.008, 0.004])


text:
struct ReferenceTotalSplitter{F}
    f::F
end

partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::ReferenceTotalSplitter{F}, args...; kwargs...) where {Tv, Ti, F} =
    invoke(partition_stripe, Tuple{SparseMatrixCSC{Tv, Ti}, Any, DynamicTotalSplitter, Vararg}, A, K, DynamicTotalSplitter(method.f), args...; kwargs...)

struct ReferenceBottleneckSplitter{F}
    f::F
end

partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::ReferenceBottleneckSplitter{F}, args...; kwargs...) where {Tv, Ti, F} =
    invoke(partition_stripe, Tuple{SparseMatrixCSC{Tv, Ti}, Any, DynamicBottleneckSplitter, Vararg}, A, K, DynamicBottleneckSplitter(method.f), args...; kwargs...)

struct ReferenceTotalChunker{F}
    f::F
end

pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::ReferenceTotalChunker{F}, args...; kwargs...) where {Tv, Ti, F} =
    invoke(pack_stripe, Tuple{SparseMatrixCSC{Tv, Ti}, DynamicTotalChunker{<:Any}, Vararg}, A, DynamicTotalChunker(method.f), args...; kwargs...)


text:
"""
    hash(α::Cyclotomic[, h::UInt])

A basic hashing function for cyclotomic elements; Note that unlike the `Base` types hashing of `Cyclotomic`s is expensive as it necessitates reducing to minimal embeding.
This is to keep `hash`ing consistent and reliable with respect to `==`, i.e. that the equality of elements implies the equality of `hash`es.
"""
function Base.hash(α::Cyclotomic, h::UInt)
    β = reduced_embedding(α)
    return hash(coeffs(β), hash(conductor(β), hash(Cyclotomic, h)))
end

function Base.:(==)(α::Cyclotomic, β::Cyclotomic)
    coeffs(α) == coeffs(β) && return true

    if conductor(α) == conductor(β)
        normalform!(α)
        normalform!(β)
        return coeffs(α) == coeffs(β)
    else
        l = lcm(conductor(α), conductor(β))
        return embed(α, l) == embed(β, l)
    end
end

Base.:(==)(α::Cyclotomic{T}, x::R) where {T,R<:Real} = α == Cyclotomic(1, [x])
Base.:(==)(x::Real, α::Cyclotomic) = α == x

function Base.isapprox(α::Cyclotomic, a::Real; kwargs...)
    x, y = reim(Complex{Float64}(α))
    return isapprox(real(a), x; kwargs...) &&
           isapprox(one(a), 1.0 + y; kwargs...)
end

function Base.isapprox(α::Cyclotomic{T}, β::Cyclotomic{T}; kwargs...) where T<:AbstractFloat
    x, y = reim(Complex{T}(α))
    z, w = reim(Complex{T}(β))
    return isapprox(x, z; kwargs...) && isapprox(y, w; kwargs...)
end

Base.isapprox(a::Number, α::Cyclotomic; kwargs...) = isapprox(α, a; kwargs...)

Base.iszero(α::Cyclotomic) =
    all(iszero, values(α)) || (normalform!(α); all(iszero, values(α)))

Base.isreal(α::Cyclotomic) =
    α == conj(α) || conductor(reduced_embedding(α)) == 1

Base.isreal(α::Cyclotomic{T}) where T<:AbstractFloat =
    α ≈ conj(α) || conductor(reduced_embedding(α)) == 1

function Base.isone(α::Cyclotomic)
    β = reduced_embedding(α)
    conductor(β) == 1 || return false
    return isone(β[0])
end

"""
    isnormalized(α::Cyclotomic, basis)
Check if `α` is already in normal form with respect to the given basis.
"""
function isnormalized(α::Cyclotomic, basis = zumbroich_basis(conductor(α)))
    # return all(in(basis), exponents(α))
    for e in exponents(α)
        e in basis || return false
    end
    return true
end

function _all_equal(α::Cyclotomic, exps, value = α[first(exps)])
    # return all(e -> α[e] == val, exps)
    for e in exps
        α[e] == value || return false
    end
    return true
end


text:
using JetPack, Jets, Test

@testset "roughness, dot product test N=$N, dim=$dim" for (N,dim) in (
        ((256,512),1), ((16,32,64,4),1), ((16,32,64,4),2), ((16,32,64,4),3))
    A = JopRoughness(JetSpace(Float64,N),dim)

    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "roughness, linearity test N=$N, dim=$dim" for (N,dim) in (
        ((256,512),1), ((16,32,64,4),1), ((16,32,64,4),2), ((16,32,64,4),3))
    A = JopRoughness(JetSpace(Float64,N),dim)

    lhs, rhs = linearity_test(A)
    @test lhs ≈ rhs
end


text:
check_encoding(enc::TwoClassEncoding, targets) =
    all(val == enc.positives || val == enc.negatives for val in unique(targets))

ispositive(enc::TwoClassEncoding, val) = _ispositive(enc, val)
_ispositive(enc::TwoClassEncoding, val) = isone(val)

Broadcast.broadcasted(::typeof(ispositive), enc, val) =
    broadcast(_ispositive, Ref(enc), val)

isnegative(enc::TwoClassEncoding, val) = _isnegative(enc, val)
_isnegative(enc::TwoClassEncoding, val) = !_ispositive(enc, val)

Broadcast.broadcasted(::typeof(isnegative), enc, val) =
    broadcast(_isnegative, Ref(enc), val)

positive_label(enc::TwoClassEncoding) = enc.positives
negative_label(enc::TwoClassEncoding) = enc.negatives

"""
    OneZero{T<:Number} <: TwoClassEncoding{T}

Two class label encoding in which `one(T)` represents the positive class,
and `zero(T)` the negative class.
"""
struct OneZero{T<:Number} <: TwoClassEncoding{T}
    positives::T
    negatives::T

    OneZero(::Type{T} = Float64) where {T<:Number} =
        new{T}(one(T), zero(T))
end

"""
    OneMinusOne{T<:Number} <: TwoClassEncoding{T}

Two class label encoding in which `one(T)` represents the positive class,
and `-one(T)` the negative class.
"""
struct OneMinusOne{T<:Number} <: TwoClassEncoding{T}
    positives::T
    negatives::T

    OneMinusOne(::Type{T} = Float64) where {T<:Number} =
        new{T}(one(T), -one(T))
end

"""
    OneTwo{T<:Number} <: TwoClassEncoding{T}

Two class label encoding in which `one(T)` represents the positive class,
and `2*one(T)` the negative class.
"""
struct OneTwo{T<:Number} <: TwoClassEncoding{T}
    positives::T
    negatives::T

    OneTwo(::Type{T} = Float64) where {T<:Number} =
        new{T}(one(T), 2*one(T))
end

"""
    OneVsOne{T} <: TwoClassEncoding{T}

Two class label encoding ...
"""
struct OneVsOne{T} <: TwoClassEncoding{T}
    positives::T
    negatives::T

    OneVsOne(pos::P, neg::N, ::Type{T} = default_type(P, N)) where {P, N, T} =
        new{T}(T(pos), T(neg))
end

_ispositive(enc::OneVsOne, val) = val == enc.positives

"""
    OneVsRest{T} <: TwoClassEncoding{T}

Two class label encoding ...
"""
struct OneVsRest{T} <: TwoClassEncoding{T}
    positives::T
    negatives::AbstractVector{T}

    OneVsRest(pos::P, neg::AbstractVector{N}, ::Type{T} = default_type(P, N)) where {P, N, T} =
        new{T}(T(pos), T.(neg))
end

check_encoding(enc::OneVsRest, targets) =
    all(val == enc.positives || (val in enc.negatives) for val in unique(targets))

_ispositive(enc::OneVsRest, val) = val == enc.positives
negative_label(enc::OneVsRest) = enc.negatives[1]

"""
    RestVsOne{T} <: TwoClassEncoding{T}

Two class label encoding ...
"""
struct RestVsOne{T} <: TwoClassEncoding{T}
    positives::AbstractVector{T}
    negatives::T

    RestVsOne(pos::AbstractVector{P}, neg::N, ::Type{T} = default_type(P, N)) where {P, N, T} =
        new{T}(T.(pos), T(neg))
end

check_encoding(enc::RestVsOne, targets) =
    all(val in enc.positives || val == enc.negatives for val in unique(targets))

_ispositive(enc::RestVsOne, val) = val != enc.negatives
positive_label(enc::RestVsOne) = enc.positives[1]


text:
export classify_branch
export classify_solutions!

"""
$(TYPEDSIGNATURES)

Creates a solution class in `res` using the inequality `condition` (parsed into Symbolics.jl input).

The new class is labelled with `name` and stored under `res.classes[name]`.

By default, only physical (=real) solutions are classified, `false` is returned for the rest.

# Example
```julia
# solve a previously-defined problem
res = get_steady_states(problem, swept_parameters, fixed_parameters)

# classify, store in result.classes["large_amplitude"]
classify_solutions!(res, "sqrt(u1^2 + v1^2) > 1.0" , "large_amplitude")
```

"""
function classify_solutions!(res::Result, condition::String, name::String; physical=true)
    expr = Num(eval(Meta.parse(condition)))
    function cond_func(s::OrderedDict, res)
        physical && !is_physical(s, res) && return false
        s = [key => real(s[key]) for key in keys(s)] # make values real
        Bool(substitute_all(expr, s).val)
    end
    values = classify_solutions(res, cond_func)
    res.classes[name] = values
end

#   Classify solutions where for `f` which is a function accepting a solution dictionary
#   specifies all params and variables).
function classify_solutions!(res::Result, f::Function, name::String)
    values = classify_solutions(res, f)
    res.classes[name] = values
end

"Return an array of booleans classifying the solution in `res`
according to `f` (`f` takes a solution dictionary, return a boolean)"
function classify_solutions(res::Result, f::Function)
    values = similar(res.solutions, Vector{Bool})
    for (idx, soln) in enumerate(res.solutions)
        values[idx] = [f(get_single_solution(res, index=idx, branch=b), res) for b in 1:length(soln)]
    end
    values
end

"""
$(TYPEDSIGNATURES)
Returns an array of booleans classifying `branch` in the solutions in `res`
according to `class`.
"""
function classify_branch(res::Result, branch::Int64, class::String)
    branch_values = getindex.(res.classes[class], branch)
end

classify_branch(soln::Result, class::String) = [classify_branch(soln, b, class) for b in 1:length(first(soln.solutions))]

"""
$(TYPEDSIGNATURES)
Returns true if the solution `soln` of the Result `res` is physical (= real number).
`im_tol` : an absolute threshold to distinguish real/complex numbers.
"""
function is_physical(soln::StateDict, res::Result; im_tol=im_tol)
    var_values = [getindex(soln, v) for v in res.problem.variables]
    return all(abs.(imag.(var_values)).<im_tol) && any(isnan.(var_values).==false)
end

"""
$(TYPEDSIGNATURES)
Returns true if the solution `soln` of the Result `res` is stable.
Stable solutions are real and have all Jacobian eigenvalues Re[λ] < 0.
`im_tol` : an absolute threshold to distinguish real/complex numbers.
"""
function is_stable(soln::StateDict, res::Result; im_tol=im_tol)
    is_physical(soln, res ,im_tol=im_tol) || return false  # the solution is unphysical anyway
    J = real.(res.jacobian(soln))
    return all([real.(eigvals(J)) .< 0]...)
end

"""
$(TYPEDSIGNATURES)
Returns true if the solution `soln` of the problem `prob` is Hopf-unstable.
Hopf-unstable solutions are real and have exactly two Jacobian eigenvalues with +ve real parts, which
are complex conjugates of each other.
`im_tol` : an absolute threshold to distinguish real/complex numbers.
"""
function is_Hopf_unstable(soln::StateDict, res::Result; im_tol=im_tol)
    is_physical(soln, res, im_tol=im_tol) || return false  # the solution is unphysical anyway
    J = res.jacobian(soln)
    unstable = filter(x -> real(x) > 0, eigvals(J))
    (length(unstable) == 2 && abs(conj(unstable[1]) - unstable[2]) < im_tol && return true) || return false
    return all([ real.(eigvals(J)) .< 0]...)
end

"""
$(TYPEDSIGNATURES)
Create binary classification of the solutions, such that each solution points receives an identifier based
on its permutation of stable branches (allows to distinguish between different phases, which may have the same number
of stable solutions). It works by converting each bistring `[is_stable(solution_1),is_stable(solution_2),...,]` into unique labels.
"""
function classify_binaries!(res::Result)
    bin_label  = bitarr_to_int.(clean_bitstrings(res)) #mapping of binary string (stable or not) for each solution set to integer. Sensitive to ordering!
    #renormalize labels with numbers from 1 to length(unique(label))
    for (idx,el) in enumerate(unique(bin_label))
        bin_label[findall(x->x==el, bin_label)] .= idx
    end
    res.classes["binary_labels"] = bin_label
end

clean_bitstrings(res::Result) = [[el for el in bit_string[phys_string]]
for (bit_string,phys_string) in zip(res.classes["stable"],res.classes["physical"])]; #remove unphysical solutions from strings

function bitarr_to_int(arr)
    return sum(arr .* (2 .^ collect(length(arr)-1:-1:0)))
end


text:
# internal data for packed integers
type Buffer{w,T<:Unsigned}
    data::Vector{T}
    function Buffer(len::Integer, mmap::Bool=false)
        @assert w ≤ bitsof(T)
        buflen = cld(len * w, bitsof(T))
        data = mmap ? Mmap.mmap(Vector{T}, buflen) : Vector{T}(buflen)
        return new(data)
    end
end

bitsof{T}(::Type{T}) = sizeof(T) * 8

function resize!{w,T}(buffer::Buffer{w,T}, len::Integer)
    buflen = cld(len * w, bitsof(T))
    resize!(buffer.data, buflen)
    return buffer
end

@inline function mask{T}(::Type{T}, w)
    ~T(0) >> (bitsof(T) - w)
end

@inline function get_chunk_id{w}(::Buffer{w,UInt8}, i::Integer)
    j = Int(i - 1) * w
    return (j >> 3) + 1, j & 0b000111
end

@inline function get_chunk_id{w}(::Buffer{w,UInt16}, i::Integer)
    j = Int(i - 1) * w
    return (j >> 4) + 1, j & 0b001111
end

@inline function get_chunk_id{w}(::Buffer{w,UInt32}, i::Integer)
    j = Int(i - 1) * w
    return (j >> 5) + 1, j & 0b011111
end

@inline function get_chunk_id{w}(::Buffer{w,UInt64}, i::Integer)
    j = Int(i - 1) * w
    return (j >> 6) + 1, j & 0b111111
end

@inline function getindex{w,T}(buf::Buffer{w,T}, i::Integer)
    k, r = get_chunk_id(buf, i)
    W = bitsof(T)
    @inbounds begin
        a = buf.data[k] >> r
        if r + w ≤ W
            return a & mask(T, w)
        else
            b = buf.data[k+1] & mask(T, (w + r) - W)
            return a | (b << (W - r))
        end
    end
end

# these width values don't cross a boundary, therefore branching can be safely removed
for w in [1, 2, 4, 8, 16, 32, 64]
    @eval begin
        @inline function getindex{T}(buf::Buffer{$w,T}, i::Integer)
            k, r = get_chunk_id(buf, i)
            @inbounds return (buf.data[k] >> r) & mask(T, $w)
        end
    end
end

# https://graphics.stanford.edu/~seander/bithacks.html#MaskedMerge
@inline mergebits(a, b, mask) = a $ ((a $ b) & mask)

@inline function setindex!{w,T}(buf::Buffer{w,T}, x::T, i::Integer)
    k, r = get_chunk_id(buf, i)
    W = bitsof(T)
    @inbounds begin
        a = buf.data[k]
        b = x << r
        buf.data[k] = mergebits(a, b, mask(T, w) << r)
        if r + w > W
            a = buf.data[k+1]
            b = x >> (W - r)
            buf.data[k+1] = mergebits(a, b, mask(T, (w + r) - W))
        end
    end
    return
end

for w in [1, 2, 4, 8, 16, 32, 64]
    @eval begin
        @inline function setindex!{T}(buf::Buffer{$w,T}, x::T, i::Integer)
            k, r = get_chunk_id(buf, i)
            @inbounds begin
                a = buf.data[k]
                b = x << r
                buf.data[k] = mergebits(a, b, mask(T, $w) << r)
            end
            return
        end
    end
end

function fill!{w,T}(buf::Buffer{w,T}, x::T)
    x &= mask(T, w)
    W = bitsof(T)
    cycle = div(lcm(w, W), W)
    r = 0
    for j in 1:cycle
        chunk = T(0)
        if r < 0
            chunk |= (x >> (w + r)) << (W + r)
        end
        r += W
        while r > 0
            chunk >>= min(w, r)
            chunk  |= x << (W - min(w, r))
            r -= w
        end
        @inbounds for i in j:cycle:endof(buf.data)
            buf.data[i] = chunk
        end
    end
    return buf
end

for w in [1, 2, 4, 8, 16, 32, 64]
    @eval begin
        function fill!{T}(buf::Buffer{$w,T}, x::T)
            chunk = T(0)
            x &= mask(T, $w)
            for _ in 1:div(bitsof(T), $w)
                chunk = chunk << $w | x
            end
            fill!(buf.data, chunk)
            return buf
        end
    end
end

function fill0!{w,T}(buf::Buffer{w,T})
    fill!(buf.data, T(0))
    return buf
end

function fill1!{w,T}(buf::Buffer{w,T})
    fill!(buf.data, ~T(0))
    return buf
end


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ 0026826c-997f-4ed7-8c0b-c54d199010a9
using Revise

# ╔═╡ 3dd034e0-95f0-11eb-0285-7fa1f7afbcae
begin
	import PWS
	using PWS: TrajectoryCallback, Trajectory, SMCEstimate, DirectMCEstimate, marginal_configuration, MarginalEnsemble, ConditionalEnsemble, gene_expression_system, generate_configuration, mutual_information, logpdf, simulate, cooperative_chemotaxis_system, reaction_network, collect_trajectory, log_marginal
	using StaticArrays
	import Catalyst
	using DiffEqBase
	using DiffEqJump
	import ModelingToolkit
	using Plots
	using Transducers
	using Statistics
end

# ╔═╡ 11b80f90-068e-4ebd-8015-d1a92218928e
n_clusters=300

# ╔═╡ bbafdb7c-a080-4ab9-ad43-fa124be6da1c
system = cooperative_chemotaxis_system(; n_clusters, phosphorylate=0.001, dephosphorylate=0.5, dtimes=0.0:0.5:100)

# ╔═╡ 5548384b-cd4e-437b-836f-8b44715b9328
conf = generate_configuration(system)

# ╔═╡ b46de91d-350d-48f6-bdf3-d0bc9fcf4632
begin
	plot(conf.s_traj |> Map((u,t,i)::Tuple -> (u ./ 100, t, i)) |> collect_trajectory, label="Ligand concentration")
	plot!(PWS.active_receptors(conf, system) |> Map((u,t,i)::Tuple -> (u ./ n_clusters, t, i)) |> collect_trajectory, label="active fraction")
	plot!(conf.x_traj |> Map((u,t,i)::Tuple -> (u[SA[1]] ./ 10000, t, i)) |> collect_trajectory)
	plot!(fmt=:png, legend=:bottomright, dpi=300)
end

# ╔═╡ 5529e55a-9436-44b1-8268-46e054dca1d1
alg = SMCEstimate(16)

# ╔═╡ d3f75f23-4c6a-4a7d-9318-dcea0382bcf3
mi = mutual_information(system, alg, num_responses=10)

# ╔═╡ 6a41ab5b-0c2c-4acd-91d4-73f8b163221a
plot(system.dtimes, mi.MutualInformation, color=:gray, legend=false)

# ╔═╡ 6b2359a9-24b5-4a58-a10b-7f8729c3e4e9
plot(system.dtimes, mean(mi.MutualInformation), ribbon=sqrt.(var(mi.MutualInformation)./size(mi, 1)), label="SMC", ylabel="Path mutual information", xlabel="Trajectory length")

# ╔═╡ Cell order:
# ╠═0026826c-997f-4ed7-8c0b-c54d199010a9
# ╠═3dd034e0-95f0-11eb-0285-7fa1f7afbcae
# ╠═11b80f90-068e-4ebd-8015-d1a92218928e
# ╠═bbafdb7c-a080-4ab9-ad43-fa124be6da1c
# ╠═5548384b-cd4e-437b-836f-8b44715b9328
# ╠═b46de91d-350d-48f6-bdf3-d0bc9fcf4632
# ╠═d3f75f23-4c6a-4a7d-9318-dcea0382bcf3
# ╠═6a41ab5b-0c2c-4acd-91d4-73f8b163221a
# ╠═6b2359a9-24b5-4a58-a10b-7f8729c3e4e9
# ╠═5529e55a-9436-44b1-8268-46e054dca1d1


text:
function bin_dist(n, k, p, q)
    return binomial(n, k) * p^k * q^(n-k)
end

function f(p, q)
    p1 = bin_dist(6, 0, p, q)
    p2 = bin_dist(6, 1, p, q)
    p3 = bin_dist(6, 2, p, q)
    return 1 - p1 - p2 - p3
end

str = readline()
strvec = split(str, ' ')
floatvec = map(x -> parse(Float64, x), strvec)
a = floatvec[1]
b = floatvec[2]

p = a/(a+b)
q = b/(a+b)

print(round(f(p, q), 3))

max_stars_repo_path: Tutorials/10 Days of Statistics/Day 4/binomial_distribution1.jl

max_stars_repo_name: xuedong/hacker-rank

max_issues_repo_name: xuedong/hacker-rank

max_forks_repo_name: xuedong/hacker-rank


text:
using Random, Test, HDF5

import HDF5.datatype
import Base.convert
import Base.unsafe_convert

struct foo
    a::Float64
    b::String
    c::String
    d::Array{ComplexF64,2}
    e::Array{Int64,1}
end

struct foo_hdf5
    a::Float64
    b::Cstring
    c::NTuple{20, UInt8}
    d::NTuple{9, ComplexF64}
    e::HDF5.hvl_t
end

function unsafe_convert(::Type{foo_hdf5}, x::foo)
    foo_hdf5(x.a,
            Base.unsafe_convert(Cstring, x.b),
            ntuple(i -> i <= ncodeunits(x.c) ? codeunit(x.c, i) : '\0', 20),
            ntuple(i -> x.d[i], length(x.d)),
            HDF5.hvl_t(length(x.e), pointer(x.e))
            )
end

function datatype(::Type{foo_hdf5})
    dtype = HDF5.h5t_create(HDF5.H5T_COMPOUND, sizeof(foo_hdf5))
    HDF5.h5t_insert(dtype, "a", fieldoffset(foo_hdf5, 1), datatype(Float64))

    vlenstr_dtype = HDF5.h5t_copy(HDF5.H5T_C_S1)
    HDF5.h5t_set_size(vlenstr_dtype, HDF5.H5T_VARIABLE)
    HDF5.h5t_set_cset(vlenstr_dtype, HDF5.H5T_CSET_UTF8)
    HDF5.h5t_insert(dtype, "b", fieldoffset(foo_hdf5, 2), vlenstr_dtype)

    fixedstr_dtype = HDF5.h5t_copy(HDF5.H5T_C_S1)
    HDF5.h5t_set_size(fixedstr_dtype, 20 * sizeof(UInt8))
    HDF5.h5t_set_cset(fixedstr_dtype, HDF5.H5T_CSET_UTF8)
    HDF5.h5t_set_strpad(fixedstr_dtype, HDF5.H5T_STR_NULLPAD)
    HDF5.h5t_insert(dtype, "c", fieldoffset(foo_hdf5, 3), fixedstr_dtype)

    hsz = HDF5.hsize_t[3,3]
    array_dtype = HDF5.h5t_array_create(datatype(ComplexF64), 2, hsz)
    HDF5.h5t_insert(dtype, "d", fieldoffset(foo_hdf5, 4), array_dtype)

    vlen_dtype = HDF5.h5t_vlen_create(datatype(Int64))
    HDF5.h5t_insert(dtype, "e", fieldoffset(foo_hdf5, 5), vlen_dtype)

    HDF5.Datatype(dtype)
end

struct bar
    a::Vector{String}
end

struct bar_hdf5
    a::NTuple{2, NTuple{20, UInt8}}
end

function datatype(::Type{bar_hdf5})
    dtype = HDF5.h5t_create(HDF5.H5T_COMPOUND, sizeof(bar_hdf5))

    fixedstr_dtype = HDF5.h5t_copy(HDF5.H5T_C_S1)
    HDF5.h5t_set_size(fixedstr_dtype, 20 * sizeof(UInt8))
    HDF5.h5t_set_cset(fixedstr_dtype, HDF5.H5T_CSET_UTF8)

    hsz = HDF5.hsize_t[2]
    array_dtype = HDF5.h5t_array_create(fixedstr_dtype, 1, hsz)

    HDF5.h5t_insert(dtype, "a", fieldoffset(bar_hdf5, 1), array_dtype)

    HDF5.Datatype(dtype)
end

function convert(::Type{bar_hdf5}, x::bar)
    bar_hdf5(ntuple(i -> ntuple(j -> j <= ncodeunits(x.a[i]) ? codeunit(x.a[i], j) : '\0', 20), 2))
end

@testset "compound" begin
    N = 10
    v = [foo(rand(),
            randstring(rand(10:100)),
            randstring(10),
            rand(ComplexF64, 3,3),
            rand(1:10, rand(10:100))
            )
        for _ in 1:N]

    v[1] = foo(1.0,
              "uniçº∂e",
              "uniçº∂e",
              rand(ComplexF64, 3,3),
              rand(1:10, rand(10:100)))

    v_write = unsafe_convert.(foo_hdf5, v)

    w = [bar(["uniçº∂e", "hello"])]
    w_write = convert.(bar_hdf5, w)

    fn = tempname()
    h5open(fn, "w") do h5f
        foo_dtype = datatype(foo_hdf5)
        foo_space = dataspace(v_write)
        foo_dset = create_dataset(h5f, "foo", foo_dtype, foo_space)
        write_dataset(foo_dset, foo_dtype, v_write)

        bar_dtype = datatype(bar_hdf5)
        bar_space = dataspace(w_write)
        bar_dset = create_dataset(h5f, "bar", bar_dtype, bar_space)
        write_dataset(bar_dset, bar_dtype, w_write)
    end

    v_read = h5read(fn, "foo")
    for field in (:a, :b, :c, :d, :e)
        f = x -> getfield(x, field)
        @test f.(v) == f.(v_read)
    end

    w_read = h5read(fn, "bar")
    for field in (:a,)
        f = x -> getfield(x, field)
        @test f.(w) == f.(w_read)
    end

    T = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, Cstring}}
    TT = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, T}}
    TTT = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, TT}}
    TTTT = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, TTT}}

    @test HDF5.do_reclaim(TTTT) == true
    @test HDF5.do_normalize(TTTT) == true

    T = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, HDF5.FixedArray}}
    TT = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, T}}
    TTT = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, TT}}
    TTTT = NamedTuple{(:a, :b, :c, :d, :e, :f), Tuple{Int, Int, Int, Int, Int, TTT}}

    @test HDF5.do_reclaim(TTTT) == false
    @test HDF5.do_normalize(TTTT) == true
end


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ cea30482-e74a-41b0-bdc5-17c2a1c5cf91
md" ## LagrangeInterpolantGenerator.jl"

# ╔═╡ 97e64f8b-18a2-42dd-8872-2acd7983fe09
md" ##### See `https://discourse.julialang.org/t/lagrange-interpolation-function-generator` and the comments on that post. See also the BaricentricInterpolation.jl example (j-05-01b.jl)"

# ╔═╡ 19161060-226e-4469-b13c-ca41a8eace3f
function LagrangeInterpolantGenerator(xvals, yvals)
	function LagrangeInterpolant(x)
		numvalstoevaluate = length(x)
        numvalstoevaluate == 1 ? output = 0 : output = zeros(numvalstoevaluate)
        for k = 1:numvalstoevaluate
            N = length(xvals)
            LagrangePolynomials = ones(N)
            for i in 1:N
				# Surprisingly, this works even in the i=1 and i=N cases.
                for j in [1:i-1;i+1:N]
                    LagrangePolynomials[i] = LagrangePolynomials[i].*(x[k] -
						xvals[j])./(xvals[i]-xvals[j])
                end
            end
            numvalstoevaluate == 1 ? output = sum(LagrangePolynomials.*yvals) :
				output[k] = sum(LagrangePolynomials.*yvals)
        end
        return output
    end
    return LagrangeInterpolant
end

# ╔═╡ 3ce58d80-543c-42b0-9830-12ffaee6e844
begin
	x = [1.0, 3.0, 6.0, 5.0]
	y = [1.0, 5.0, 10.0, 9.0]
	xi = [2.0, 4.5]
end

# ╔═╡ e4e3e4de-2dbe-4a0d-9589-fb8456781a73
interpolantfunc = LagrangeInterpolantGenerator(x,y);

# ╔═╡ 748369b7-3df8-4bed-a125-fa812a92ed1c
interpolantfunc.(xi)

# ╔═╡ ccf3b050-7cac-41c4-baf4-c6e535590742
md" ## End of LagrageInterpolantGenerator.jl"

# ╔═╡ Cell order:
# ╟─cea30482-e74a-41b0-bdc5-17c2a1c5cf91
# ╠═97e64f8b-18a2-42dd-8872-2acd7983fe09
# ╠═19161060-226e-4469-b13c-ca41a8eace3f
# ╠═3ce58d80-543c-42b0-9830-12ffaee6e844
# ╠═e4e3e4de-2dbe-4a0d-9589-fb8456781a73
# ╠═748369b7-3df8-4bed-a125-fa812a92ed1c
# ╟─ccf3b050-7cac-41c4-baf4-c6e535590742


text:
using MLStyle
using DataStructures
using SimpleGraphs
using SimplePosets

export argtuple
argtuple(m) = arguments(m) |> astuple

astuple(x) = Expr(:tuple,x...)

function variables(expr :: Expr)
    leaf(x::Symbol) = begin
        [x]
    end
    leaf(x) = []
    branch(head, newargs) = begin
        union(newargs...)
    end
    foldall(leaf, branch)(expr)
end

variables(s::Symbol) = [s]
variables(x) = []

export arguments
arguments(m) = m.args

export stochastic
stochastic(st :: Let)        = Symbol[]
stochastic(st :: Follows)    = Symbol[st.x]
stochastic(st :: Return)     = Symbol[]
stochastic(st :: LineNumber) = Symbol[]

stochastic(m :: Model) = union(stochastic.(m.body)...)

export bound
bound(st :: Let)        = Symbol[st.x]
bound(st :: Follows)    = Symbol[]
bound(st :: Return)     = Symbol[]
bound(st :: LineNumber) = Symbol[]

bound(m :: Model) = union(bound.(m.body)...)

# """
#     parameters(m::Model)

# A _parameter_ is a stochastic node that is only assigned once
# """
export parameters
parameters(m::Model) = setdiff(
    stochastic(m),
    arguments(m) ∪ bound(m)
)

export observed
observed(m::Model) = setdiff(stochastic(m), parameters(m))

export variables
variables(m :: Model) = arguments(m) ∪ stochastic(m) ∪ bound(m)

export freeVariables
function freeVariables(m::Model)
    setdiff(arguments(m), stochastic(m))
end

export foldall
function foldall(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, args...) => branch(head, map(go, args); kwargs...)
            x                   => leaf(x; kwargs...)
        end
    end

    return go
end

export foldall1
function foldall1(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, arg1, args...) => branch(head, arg1, map(go, args); kwargs...)
            x                         => leaf(x; kwargs...)
        end
    end

    return go
end

# function condition(vs...)
#     function cond(m)
#         po = dependencies(m)

#         maybeRemove = Symbol[]
#         # Bound v ∈ vs no longer depend on other variables
#         for v ∈ (vs ∩ bound(m))
#             for x in below(po, v)
#                 delete!(po, x, v)
#                 push!(maybeRemove, x)
#             end
#         end

#         # Go upward in the poset from the stuff we're keeping
#         # Keep everything maximal from there, and below
#         keep = [vs...]
#         for m in maximals(po)
#             xs = union([m],below(po,m))
#             if !isempty(keep ∩ xs)
#                 union!(keep, xs)
#             end
#         end

#         removable = setdiff(maybeRemove, keep)

#         function proc(st::Let)
#             st.x ∈ vs && return false
#             st.x ∈ removable && return false
#             return true
#         end
#         function proc(st::Follows)
#             st.x ∈ removable && return false
#             return true
#         end
#         proc(st) = true

#         Model(m.args, filter(proc, m.body))
#     end

# end

import MacroTools: striplines, @q

# import LogDensityProblems: logdensity
# using ResumableFunctions

# export arguments, args, stochastic, observed, parameters, supports
# export paramSupport

# import MacroTools: striplines, flatten, unresolve, resyntax, @q, @capture
# using StatsFuns
using DataStructures: counter

# function arguments(model::Model)
#     model.args
# end

export digraph
function digraph(m::Model)
    g = SimpleDigraph{Symbol}()

    mvars = variables(m)
    for v in mvars
        add!(g, v)
    end

    f!(g, st::Let) =
        for v in mvars ∩ variables(st.rhs)
            add!(g, v, st.x)
        end
    f!(g, st::Follows) =
        for v in mvars ∩ variables(st.rhs)
            add!(g, v, st.x)
        end
    f!(g, st::Return)  = nothing
    f!(g, st::LineNumber) = nothing

    for st in m.body
        f!(g, st)
    end

    g
end

export poset
function poset(m::Model)
    po = SimplePoset{Symbol}()

    mvars = variables(m)
    for v in mvars
        add!(po, v)
    end

    f!(po, st::Let) =
        for v in mvars ∩ variables(st.rhs)
            add!(po, v, st.x)
        end
    f!(po, st::Follows) =
        for v in mvars ∩ variables(st.rhs)
            add!(po, v, st.x)
        end
    f!(po, st::Return)  = nothing
    f!(po, st::LineNumber) = nothing

    for st in m.body
        f!(po, st)
    end

    po
end

export dependencies
dependencies = poset

# # export paramSupport
# # function paramSupport(model)
# #     supps = Dict{Symbol, Any}()
# #     postwalk(model.body) do x
# #         if @capture(x, v_ ~ dist_(args__))
# #             if v in parameters(model)
# #                 supps[v] = support(eval(dist))
# #             end
# #         else x
# #         end
# #     end
# #     return supps
# # end

export makeLogdensity
function makeLogdensity(m :: Model)
    fpre = sourceLogdensity(m) |> eval
    f(par) = invokefrozen(fpre, Real, par)
end

export logdensity
logdensity(m::Model, par) = makeLogdensity(m)(par)

export sourceLogdensity
function sourceLogdensity(m::Model; ℓ=:ℓ, fname = gensym(:logdensity))
    proc(m, st :: Follows)    = :($ℓ += logpdf($(st.rhs), $(st.x)))
    proc(m, st :: Let)        = :($(st.x) = $(st.rhs))
    proc(m, st :: Return)     = nothing
    proc(m, st :: LineNumber) = nothing
    proc(::Nothing)        = nothing
    body = buildSource(m, proc)

    unknowns = parameters(m) ∪ arguments(m)
    unkExpr = Expr(:tuple,unknowns...)
    @gensym logdensity
    result = @q function $fname(pars)
        @unpack $(unkExpr) = pars
        $ℓ = 0.0

        $body
        return $ℓ
    end

    flatten(result)
end

allequal(xs) = all(xs[1] .== xs)

# export findsubexprs
# function findsubexprs(expr, vs)
#     intersect(getSymbols(expr), vs)
# end

export prior
function prior(m :: Model)
    po = dependencies(m)
    keep = parameters(m)
    for v in keep
        union!(keep, below(po, v))
    end
    proc(m, st::Follows) = st.x ∈ keep
    proc(m, st::Let)     = st.x ∈ keep
    proc(m, st) = true
    newbody = filter(st -> proc(m,st), m.body)
    Model([],newbody)
end

# # export priorPredictive
# # function priorPredictive(m :: Model)
# #     args = copy(m.args)
# #     body = postwalk(m.body) do x
# #         if @capture(x, v_ ~ dist_)
# #             setdiff!(args, [v])
# #             x
# #         elseif @capture(x, v_ ⩪ dist_)
# #             setdiff!(args, [v])
# #             @q ($v ~ $dist)
# #         else x
# #         end
# #     end
# #     Model(args, body)
# # end

# export likelihood
# function likelihood(m :: Model)
#     m = annotate(m)
#     args = copy(m.args)
#     body = postwalk(m.body) do x
#         if @capture(x, v_ ~ dist_)
#             union!(args, [v])
#             Nothing
#         elseif @capture(x, v_ ⩪ dist_)
#             setdiff!(args, [v])
#             @q ($v ~ $dist)
#         else x
#         end
#     end
#     Model(args, body) |> pretty
# end

# export annotate
# function annotate(m::Model)
#     newbody = postwalk(m.body) do x
#         if @capture(x, v_ ~ dist_)
#             if v ∈ observed(m)
#                 @q $v ⩪ $dist
#             else
#                 x
#             end
#         else x
#         end
#     end

#     Model(args = m.args, body=newbody, meta = m.meta)
# end

# export pretty
# pretty = stripNothing ∘ striplines ∘ flatten

# export expandSubmodels
# function expandSubmodels(m :: Model)
#     newbody = postwalk(m.body) do x
#         if @capture(x, @model expr__)
#             eval(x)
#         else x
#         end
#     end
#     Model(args=m.args, body=newbody, meta=m.meta)
# end

function unobserve(m::Model; ℓ=:ℓ)
    function proc(m, st :: Follows)
        st.x ∈ observed(m) && return nothing
        return st
    end
    proc(m, st) = st
    body = buildSource(m, proc)

    Model(freeVariables(m), body)
end

# # fold example usage:
# # ------------------
# # function leafCount(ast)
# #     leaf(x) = 1
# #     expr(head, arg1, newargs) = sum(newargs)
# #     fold(leaf, expr)(ast)
# # end

# # leaves = begin
# #     leaf(x) = [x]
# #     expr(head, arg1, newargs) = union(newargs...)
# #     fold(leaf, expr)
# # end

# # ast = :(f(x + 3y))

# # leaves(ast)

# # Example of Tamas Papp's `as` combinator:
# # julia> as((;s=as(Array, as𝕀,4), a=asℝ))(randn(5))
# # (s = [0.545324, 0.281332, 0.418541, 0.485946], a = 2.217762640580984)

function buildSource(m, proc; kwargs...)
    q = @q begin end
    for st in m.body
        ex = proc(m, st; kwargs...)
        isnothing(ex) || push!(q.args, ex)
    end
    q
end

# From https://github.com/thautwarm/MLStyle.jl/issues/66
@active LamExpr(x) begin
           @match x begin
               :($a -> begin $(bs...) end) =>
                 let exprs = filter(x -> !(x isa LineNumberNode), bs)
                   if length(exprs) == 1
                     (a, exprs[1])
                   else
                     (a, Expr(:block, bs...))
                     end
               end
                _  => nothing
           end
       end

# @inline function invokefrozen(f, rt, args...)
#     return _invokefrozen(f, rt, args)
# end

# From Chris Rackauckas: https://github.com/JuliaLang/julia/pull/32737
@inline @generated function _invokefrozen(f, ::Type{rt}, args...) where rt
    tupargs = Expr(:tuple,(a==Nothing ? Int : a for a in args)...)
    quote
        _f = $(Expr(:cfunction, Base.CFunction, :f, rt, :((Core.svec)($((a==Nothing ? Int : a for a in args)...))), :(:ccall)))
        return ccall(_f.ptr,rt,$tupargs,$((:(getindex(args,$i) === nothing ? 0 : getindex(args,$i)) for i in 1:length(args))...))
    end
end

# @cscherrer's modification of `invokelatest` does better on kwargs
export invokefrozen
@inline function invokefrozen(f, rt, args...; kwargs...)
    g(kwargs, args...) = f(args...; kwargs...)
    kwargs = (;kwargs...)
    _invokefrozen(g, rt, (;kwargs...), args...)
end

@inline function invokefrozen(f, rt, args...)
    _invokefrozen(f, rt, args...)
end

# using BenchmarkTools
# f(;kwargs...) = kwargs[:a] + kwargs[:b]

# @btime invokefrozen(f, Int; a=3,b=4)  # 3.466 ns (0 allocations: 0 bytes)
# @btime f(;a=3,b=4)                    # 1.152 ns (0 allocations: 0 bytes)

abstract type TypeLevel end
struct TLCons{Hd, Tl} <: TypeLevel end
struct TLNil <: TypeLevel end
struct TLVal{Val} <: TypeLevel end
struct TLSExp{Fn, Args} <: TypeLevel end

function typelevellist(l)
    foldr(l, init=TLNil) do each, prev
        TLCons{each, prev}
    end
end

function expr2typelevel(x)
    r = expr2typelevel
    @match x begin
        Expr(hd, tl...) =>
            let hd = r(hd),
                tl = map(r, tl) |> typelevellist,
                f = TLVal{Expr},
                args = TLCons{hd, tl}
            TLSExp{f, args}
            end
        ln :: LineNumberNode =>
            let f = TLVal{LineNumberNode},
               args = [
                    r(ln.line),
                    r(ln.file)
                ] |> typelevellist
            TLSExp{f, args}
            end
        x::QuoteNode =>
            let f = TLVal{QuoteNode},
                args = [r(x.value)] |> typelevellist

            TLSExp{f, args}
            end
        a => TLVal{a}
    end
end

function interpret(t::Type{TLNil})
    []
end

function interpret(t::Type{TLVal{Val}}) where Val
    Val
end

function interpret(t::Type{TLCons{Hd, Tl}}) where {Hd, Tl}
    tl = interpret(Tl)
    @assert tl isa Vector
    [interpret(Hd), tl...]
end

function interpret(t::Type{TLSExp{Fn, Args}}) where {Fn, Args}
    args = interpret(Args)
    @assert args isa Vector
    interpret(Fn)(args...)
end


text:
using AbInitioSoftwareBase: parentdir
using AbInitioSoftwareBase.Inputs: Setter
using Crystallography: MonkhorstPackGrid
using Dates: format, now
using QuantumESPRESSO.Commands: pw
using QuantumESPRESSO.Inputs.PWscf: KMeshCard, PWInput, VerbositySetter
using Setfield: @set!
using UnifiedPseudopotentialFormat  # To work with `download_potential`
using Unitful: ustrip, @u_str
using UnitfulAtomic

import Express.ConvergenceTestWorkflow: MakeInput, RunCmd

(::MakeInput)(template::PWInput, args...) = (customizer(args...) ∘ normalizer())(template)

struct CutoffEnergySetter <: Setter
    wfc::Number
end
function (x::CutoffEnergySetter)(template::PWInput)
    @set! template.system.ecutwfc = ustrip(u"Ry", x.wfc)
    return template
end

normalizer() = VerbositySetter("high")

struct OutdirSetter <: Setter
    timefmt::String
end
function (x::OutdirSetter)(template::PWInput)
    @set! template.control.outdir = abspath(
        joinpath(
            template.control.outdir,
            join((template.control.prefix, format(now(), x.timefmt), rand(UInt)), '_'),
        ),
    )
    mkpath(template.control.outdir)
    return template
end

struct MonkhorstPackGridSetter <: Setter
    mesh::Vector{Int}
    shift::Vector{Int}
end
function (x::MonkhorstPackGridSetter)(template::PWInput)
    @set! template.k_points = KMeshCard(MonkhorstPackGrid(x.mesh, x.shift))
    return template
end

customizer(mesh, shift, timefmt = "Y-m-d_H:M:S") =
    OutdirSetter(timefmt) ∘ MonkhorstPackGridSetter(mesh, shift)
customizer(energy::Number, timefmt::AbstractString = "Y-m-d_H:M:S") =
    OutdirSetter(timefmt) ∘ CutoffEnergySetter(energy)

(x::RunCmd)(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    pw(input, output; kwargs...)


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_gemm #
#############

function test_gemm(elty)
    A = sparse(rand(elty,m,k))
    B = sparse(rand(elty,k,n))
    C = A * B
    d_A = CudaSparseMatrixCSR(A)
    d_B = CudaSparseMatrixCSR(B)
    d_C = CUSPARSE.gemm('N','N',d_A,d_B,'O','O','O')
    r_r = to_host(d_C.rowPtr)
    r_c = to_host(d_C.colVal)
    r_v = to_host(d_C.nzVal)
    h_C = to_host(d_C)
    @test_approx_eq(C,h_C)
    @test_throws(DimensionMismatch,CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O'))
    @test_throws(DimensionMismatch,CUSPARSE.gemm('T','T',d_A,d_B,'O','O','O'))
    @test_throws(DimensionMismatch,CUSPARSE.gemm('T','N',d_A,d_B,'O','O','O'))
    @test_throws(DimensionMismatch,CUSPARSE.gemm('N','N',d_B,d_A,'O','O','O'))
    #=A = sparse(rand(elty,m,k))
    B = sparse(rand(elty,k,n))
    d_A = CudaSparseMatrixCSC(A)
    d_B = CudaSparseMatrixCSC(B)
    C = A * B
    d_C = CUSPARSE.gemm('N','N',d_A,d_B,'O','O','O')
    h_C = to_host(d_C)
    @test_approx_eq(C,h_C)
    @test_throws(DimensionMismatch,CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O'))
    @test_throws(DimensionMismatch,CUSPARSE.gemm('T','T',d_A,d_B,'O','O','O'))
    @test_throws(DimensionMismatch,CUSPARSE.gemm('T','N',d_A,d_B,'O','O','O'))
    @test_throws(DimensionMismatch,CUSPARSE.gemm('N','N',d_B,d_A,'O','O','O'))=#
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_gemm(elty)
    println("gemm took ", toq(), " for ", elty)
end


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Fisher" begin
t = HypothesisTests.FisherExactTest(1, 1, 1, 1)
@test t.ω ≈ 1.0
@test pvalue(t; tail=:left) ≈ 0.8333333333333337
@test pvalue(t; tail=:right) ≈ 0.8333333333333337
@test pvalue(t; method=:central) ≈ 1.0
@test pvalue(t; method=:minlike) ≈ 1.0
@test default_tail(t) == :both
@test_ci_approx confint(t; tail=:left) (0.0, 76.24918299781056)
@test_ci_approx confint(t; tail=:right) (0.013114894621608135, Inf)
@test_ci_approx confint(t; method=:central) (0.006400016357911029, 156.2496006379585)
#@test_approx_eq [confint(t; method=:minlike)...] [0.0131, 76.2492]
show(IOBuffer(), t)

# http://en.wikipedia.org/wiki/Fisher%27s_exact_test
t = HypothesisTests.FisherExactTest(1, 9, 11, 3)
@test t.ω ≈ 0.03720908483238119
@test pvalue(t; tail=:left) ≈ 0.0013797280926100427
@test pvalue(t; tail=:right) ≈ 0.9999663480953023
@test pvalue(t; method=:central) ≈ 0.0027594561852200853
@test pvalue(t; method=:minlike) ≈ 0.002759456185220088
@test_ci_approx confint(t; tail=:left) (0.0, 0.32600296913224913)
@test_ci_approx confint(t; tail=:right) (0.0012948958389639856, Inf)
@test_ci_approx confint(t; method=:central) (0.0006360029488751071, 0.42586647569637387)
#@test_approx_eq [confint(t; method=:minlike)...] [0.0013, 0.3567]
show(IOBuffer(), t)

# http://www.physics.csbsju.edu/stats/exact.html
t = HypothesisTests.FisherExactTest(7, 12, 0, 5)
@test t.ω ≈ Inf
@test pvalue(t; tail=:left) ≈ 1.0
@test pvalue(t; tail=:right) ≈ 0.1455862977602108
@test pvalue(t; method=:central, tail=:both) ≈ 0.2911725955204216
@test pvalue(t; method=:minlike, tail=:both) ≈ 0.2720685111989459
@test_ci_approx confint(t; tail=:left) (0.0, Inf)
@test_ci_approx confint(t; tail=:right) (0.539859556284207, Inf)
@test_ci_approx confint(t; method=:central) (0.39239937500428096, Inf)
#@test_approx_eq [confint(t; method=:minlike)...] [0.5171, Inf]
show(IOBuffer(), t)

t = HypothesisTests.FisherExactTest(12, 7, 5, 0)
@test t.ω ≈ 0.0
@test pvalue(t; tail=:left) ≈ 0.1455862977602108
@test pvalue(t; tail=:right) ≈ 1.0
@test pvalue(t; method=:central, tail=:both) ≈ 0.29117259552042146
@test pvalue(t; method=:minlike, tail=:both) ≈ 0.2720685111989459
@test_ci_approx confint(t; tail=:left) (0.0, 1.852333608546057)
@test_ci_approx confint(t; tail=:right) (0.0, Inf)
@test_ci_approx confint(t; method=:central) (0.0, 2.5484240386190433)
#@test_approx_eq [confint(t; method=:minlike)...] [0.0, 1.9338]
show(IOBuffer(), t)

t = HypothesisTests.FisherExactTest(0, 5, 7, 12)
@test t.ω ≈ 0.0
@test pvalue(t; tail=:left) ≈ 0.1455862977602108
@test pvalue(t; tail=:right) ≈ 1.0
@test pvalue(t; method=:central, tail=:both) ≈ 0.29117259552042146
@test pvalue(t; method=:minlike, tail=:both) ≈ 0.2720685111989459
@test_ci_approx confint(t; tail=:left) (0.0, 1.8523336085460567)
@test_ci_approx confint(t; tail=:right) (0.0, Inf)
@test_ci_approx confint(t; method=:central) (0.0, 2.5484240386190433)
#@test_approx_eq [confint(t; method=:minlike)...] [0.0, 1.9338]
show(IOBuffer(), t)

t = HypothesisTests.FisherExactTest(5, 0, 12, 7)
@test t.ω ≈ Inf
@test pvalue(t; tail=:left) ≈ 1.0
@test pvalue(t; tail=:right) ≈ 0.1455862977602108
@test pvalue(t; method=:central, tail=:both) ≈ 0.29117259552042146
@test pvalue(t; method=:minlike, tail=:both) ≈ 0.2720685111989459
@test_ci_approx confint(t; tail=:left) (0.0, Inf)
@test_ci_approx confint(t; tail=:right) (0.5398595562842079, Inf)
@test_ci_approx confint(t; method=:central) (0.39239937500428096, Inf)
#@test_approx_eq [confint(t; method=:minlike)...] [0.5171, Inf]
show(IOBuffer(), t)

# http://www.stata.com/support/faqs/statistics/fishers-exact-test/
t = HypothesisTests.FisherExactTest(2, 31, 136, 15532)
@test t.ω ≈ 7.3653226779913386
@test pvalue(t; tail=:left) ≈ 0.9970112864705307
@test pvalue(t; tail=:right) ≈ 0.03390271476034175
@test pvalue(t; method=:central, tail=:both) ≈ 0.06780542952068347
@test pvalue(t; method=:minlike, tail=:both) ≈ 0.03390271476034175
@test_ci_approx confint(t; tail=:left) (0.0,25.20252454804777)
@test_ci_approx confint(t; tail=:right) (1.2436262312601785, Inf)
@test_ci_approx confint(t; method=:central) (0.8458141614657836, 29.44434308524672)
#@test_approx_eq [confint(t; method=:minlike)...] [1.2436, 28.3557]
show(IOBuffer(), t)

# http://www.utstat.toronto.edu/~brunner/oldclass/312f12/lectures/312f12FisherWithR.pdf
# http://vassarstats.net/odds2x2.html
t = HypothesisTests.FisherExactTest(4, 1, 20, 1)
@test t.ω ≈ 0.21821789023599236
@test pvalue(t; tail=:left) ≈ 0.353846153846154
@test pvalue(t; tail=:right) ≈ 0.9692307692307692
@test pvalue(t; method=:central) ≈ 0.7076923076923074
@test pvalue(t; method=:minlike) ≈ 0.353846153846154
@test_ci_approx confint(t; tail=:left) (0.0, 9.594302003876502)
@test_ci_approx confint(t; tail=:right) (0.004963263361921223, Inf)
@test_ci_approx confint(t; method=:central) (0.002430190787475382, 19.59477744071154)
#@test_approx_eq [confint(t; method=:minlike)...] [0.005, 9.5943]
show(IOBuffer(), t)

t = HypothesisTests.FisherExactTest(1, 1, 1, 1)
@test HypothesisTests.pvalue(t, tail=:both) <= 1
show(IOBuffer(), t)
end


text:
using StipplePlotlyExport
using Test

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


text:
using LearnBenchmark
using Test

@testset "LearnBenchmark.jl" begin

    @testset "countpoints" begin
        @testset "1D" begin
            dataset = zeros(1,10)
            dataset[1,:] .= 1:10

            @test LearnBenchmark.countpoints([1],1.1,dataset) == 2
            @test LearnBenchmark.countpoints([2],1.1,dataset) == 3
            @test LearnBenchmark.countpoints([3],4.1,dataset) == 5
        end

        @testset "2D" begin
            dataset = zeros(2,100)
            dataset[1,:] .= repeat(1:10,10)
            dataset[2,:] .= repeat(1:10,inner = 10)

            @test LearnBenchmark.countpoints([2;2],1.1,dataset) == 5

        end
    end

    @testset "density_ratio_estimator!" begin
        @testset "1D" begin
            data1 = zeros(1,10)
            data1[1:10] .= 1:10
            data2 = zeros(1,10)
            data2[1:10] = 1:10
            dataset = [data1, data2]

            ratios = zeros(2)
            center = [3]
            LearnBenchmark.density_ratio_estimator!(ratios,center, 1.5, dataset)
            @test ratios[1] ≈ 1

            data2[1:10] = 1:.5:5.5
            dataset = [data1, data2]
            ratios = zeros(2)
            LearnBenchmark.density_ratio_estimator!(ratios,center, 1.6^2, dataset)
            @test ratios[1] ≈ 3/7
        end
        @testset "2D" begin
            data1 = zeros(2,100)
            data2 = zeros(2,100)
            data1[1,:] = repeat(1:10,10)
            data1[2,:] = repeat(1:10,inner=10)
            data2[1,:] = repeat(.5:.5:5,10)
            data2[2,:] = repeat(.5:.5:5,inner=10)
            dataset = [data1, data2]

            ratios = zeros(2)
            center = [2,2]
            ϵ = 1.1
            LearnBenchmark.density_ratio_estimator!(ratios,center,ϵ^2,dataset)
            @test ratios[1] ≈ 5/13

        end
    end

    @testset "Chebyshev Polynomials" begin
        @testset "chebyshev_poly" begin
            # First polynomial
            @test LearnBenchmark.chebyshev_poly(1, 0) == 1

            # Second
            @test LearnBenchmark.chebyshev_poly(.5, 1) == .5
            @test LearnBenchmark.chebyshev_poly(1, 1) == 1

            # Third f(x) = 2x^2 - 1
            @test LearnBenchmark.chebyshev_poly(.5, 2) == -.5
            @test LearnBenchmark.chebyshev_poly(1, 2) == 1
        end

        @testset "chebyshev_roots" begin
            @test all(LearnBenchmark.chebyshev_roots(1,2) .≈ [1])
            @test all(LearnBenchmark.chebyshev_roots(2,2) .≈ [1+sqrt(2)/2,1-sqrt(2)/2])
        end

        @testset "chebyshev_weights" begin
            @testset "Constraints" begin
                @test sum(LearnBenchmark.chebyshev_weights(5,2)) ≈ 1

                d=4
                L=6
                weight_test = LearnBenchmark.chebyshev_weights(L,d)
                root_test = LearnBenchmark.chebyshev_roots(L)
                for i = 1:d
                    @test sum(weight_test .* root_test.^i) ≈ 0 atol=1e-14
                end
            end
        end
    end

    @testset "tₖ" begin
        @testset "Even Prior 2 variables" begin
            prior = [0.5, 0.5]
            data =  [0.5, 10]
            @test LearnBenchmark.tₖ(prior,data) == (1 - data[1])/2
        end
    end

    @testset "hashmap" begin
        @testset "hashblock" begin
            x = 1
            b = 0
            ϵ = 1
            @test LearnBenchmark.hashblock(x,b,ϵ) == x

            x = 1
            b = 0
            ϵ = 0.1
            @test LearnBenchmark.hashblock(x,b,ϵ) == 10

            x = 1
            b = 1
            ϵ = 0.099
            @test LearnBenchmark.hashblock(x,b,ϵ) == 20
        end
    end

end