text:
"""
contour(x, y, z)

Creates a contour plot of the plane spanning x::Vector, y::Vector, z::Matrix

## Attributes
$(ATTRIBUTES)
"""
@recipe(Contour) do scene
    default = default_theme(scene)
    pop!(default, :color)
    Attributes(;
        default...,
        color = nothing,
        colormap = theme(scene, :colormap),
        colorrange = AbstractPlotting.automatic,
        levels = 5,
        linewidth = 1.0,
        alpha = 1.0,
        fillrange = false,
    )
end

"""
contour3d(x, y, z)

Creates a 3D contour plot of the plane spanning x::Vector, y::Vector, z::Matrix,
with z-elevation for each level.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Contour3d) do scene
    default_theme(scene, Contour)
end

function contourlines(::Type{<: Contour}, contours, cols)
    result = Point2f0[]
    colors = RGBA{Float32}[]
    for (color, c) in zip(cols, Contours.levels(contours))
        for elem in Contours.lines(c)
            append!(result, elem.vertices)
            push!(result, Point2f0(NaN32))
            append!(colors, fill(color, length(elem.vertices) + 1))
        end
    end
    result, colors
end

function contourlines(::Type{<: Contour3d}, contours, cols)
    result = Point3f0[]
    colors = RGBA{Float32}[]
    for (color, c) in zip(cols, Contours.levels(contours))
        for elem in Contours.lines(c)
            for p in elem.vertices
                push!(result, Point3f0(p[1], p[2], c.level))
            end
            push!(result, Point3f0(NaN32))
            append!(colors, fill(color, length(elem.vertices) + 1))
        end
    end
    result, colors
end

to_levels(x::AbstractVector{<: Number}, cnorm) = x

function to_levels(n::Integer, cnorm)
    zmin, zmax = cnorm
    dz = (zmax - zmin) / (n + 1)
    range(zmin + dz; step = dz, length = n)
end

conversion_trait(::Type{<: Contour3d}) = SurfaceLike()
conversion_trait(::Type{<: Contour}) = SurfaceLike()
conversion_trait(::Type{<: Contour{<: Tuple{X, Y, Z, Vol}}}) where {X, Y, Z, Vol} = VolumeLike()
conversion_trait(::Type{<: Contour{<: Tuple{<: AbstractArray{T, 3}}}}) where T = VolumeLike()

function plot!(plot::Contour{<: Tuple{X, Y, Z, Vol}}) where {X, Y, Z, Vol}
    x, y, z, volume = plot[1:4]
    @extract plot (colormap, levels, linewidth, alpha)
    valuerange = lift(nan_extrema, volume)
    cliprange = replace_automatic!(plot, :colorrange) do
        valuerange
    end
    cmap = lift(colormap, levels, linewidth, alpha, cliprange, valuerange) do _cmap, l, lw, alpha, cliprange, vrange
        levels = to_levels(l, vrange)
        nlevels = length(levels)
        N = nlevels * 50
        iso_eps = nlevels * ((vrange[2] - vrange[1]) / N) # TODO calculate this
        cmap = to_colormap(_cmap)
        v_interval = cliprange[1] .. cliprange[2]
        # resample colormap and make the empty area between iso surfaces transparent
        map(1:N) do i
            i01 = (i-1) / (N - 1)
            c = AbstractPlotting.interpolated_getindex(cmap, i01)
            isoval = vrange[1] + (i01 * (vrange[2] - vrange[1]))
            line = reduce(levels, init = false) do v0, level
                (isoval in v_interval) || return false
                v0 || (abs(level - isoval) <= iso_eps)
            end
            RGBAf0(Colors.color(c), line ? alpha : 0.0)
        end
    end
    volume!(
        plot, x, y, z, volume, colormap = cmap, colorrange = cliprange, algorithm = 7,
        transparency = plot.transparency, overdraw = plot.overdraw,
        ambient = plot.ambient, diffuse = plot.diffuse, lightposition = plot.lightposition,
        shininess = plot.shininess, specular = plot.specular
    )
end

function color_per_level(color, colormap, colorrange, alpha, levels)
    color_per_level(to_color(color), colormap, colorrange, alpha, levels)
end

function color_per_level(color::Colorant, colormap, colorrange, alpha, levels)
    fill(color, length(levels))
end

function color_per_level(colors::AbstractVector, colormap, colorrange, alpha, levels)
    color_per_level(to_colormap(colors), colormap, colorrange, alpha, levels)
end

function color_per_level(colors::AbstractVector{<: Colorant}, colormap, colorrange, alpha, levels)
    if length(levels) == length(colors)
        return colors
    else
        # TODO resample?!
        error("For a contour plot, `color` with an array of colors needs to
        have the same length as `levels`.
        Found $(length(colors)) colors, but $(length(levels)) levels")
    end
end

function color_per_level(::Nothing, colormap, colorrange, a, levels)
    cmap = to_colormap(colormap)
    map(levels) do level
        c = interpolated_getindex(cmap, level, colorrange)
        RGBAf0(color(c), alpha(c) * a)
    end
end

function plot!(plot::T) where T <: Union{Contour, Contour3d}
    x, y, z = plot[1:3]
    if to_value(plot[:fillrange])
        plot[:interpolate] = true
        # TODO normalize linewidth for heatmap
        plot[:linewidth] = map(x-> x ./ 10f0, plot[:linewidth])
        heatmap!(plot, Attributes(plot), x, y, z)
    else
        zrange = lift(nan_extrema, z)
        levels = lift(plot[:levels], zrange) do levels, zrange
            if levels isa AbstractVector{<: Number}
                return levels
            elseif levels isa Integer
                to_levels(levels, zrange)
            else
                error("Level needs to be Vector of iso values, or a single integer to for a number of automatic levels")
            end
        end
        replace_automatic!(plot, :colorrange) do
            lift(nan_extrema, levels)
        end
        args = @extract plot (color, colormap, colorrange, alpha)
        level_colors = lift(color_per_level, args..., levels)
        result = lift(x, y, z, levels, level_colors) do x, y, z, levels, level_colors
            t = eltype(z)
            # Compute contours
            xv, yv = to_vector(x, size(z,1), t), to_vector(y, size(z,2), t)
            contours = Contours.contours(xv, yv, z,  convert(Vector{eltype(z)}, levels))
            contourlines(T, contours, level_colors)
        end
        lines!(
            plot, lift(first, result);
            color = lift(last, result), linewidth = plot[:linewidth]
        )
    end
    plot
end

function data_limits(x::Contour{<: Tuple{X, Y, Z}}) where {X, Y, Z}
    return xyz_boundingbox(transform_func(x), to_value.((x[1], x[2]))...)
end


text:
module NCU
using QuantExQASM.GateOps
using QuantExQASM.Circuit

ops_mat_cache = Dict{GateOps.GateLabel, Matrix{<:Number}}()

# =========================================================================== #
#
# =========================================================================== #

"""
    init_intermed_gates(circ::Circuit.Circ, num_ctrl::Union{Nothing, Int})

Generates the required intermediate gates for implementing the NCU using the given Circuit gate-set. Calling this routine is essential before running the NCU.apply_ncu! function.
"""
function init_intermed_gates(circ::Circuit.Circ, num_ctrl::Union{Nothing, Int})
    for k in circ.gate_set
        gen_intermed_gates(num_ctrl == nothing ? 8 : num_ctrl, k)
    end
end

"""
    register_gate(circ::Circuit.Circ, U::GateOps.GateLabel, gate::Matrix{<:Number})

Adds the user-defined gate to the cache with the pairing U=>gate
"""
function register_gate(circ::Circuit.Circ, U::GateOps.GateLabel, gate::Matrix{<:Number})
    # Some characters are reserved for internal gates
    @assert !(String(U.label) in ["x","y","z","h","s","t","a","c"])
    Circuit.gate_cache[U] = gate
end

"""
    gen_intermed_gates(ctrl_depth::Int, U::GateOps.GateLabel)

Generates all intermediate n-th root gates upto the given control depth.
"""
function gen_intermed_gates(ctrl_depth::Int, U::GateOps.GateLabel)
    su,asu = get_intermed_gate(U)
    for i in 2:ctrl_depth-1
        su,asu = get_intermed_gate(su)
    end
end

"""
    get_intermed_gate(U::GateOps.GateLabel)

Generate sqrt U and adjoint sqrt U, given a previously cached gate U.
"""
function get_intermed_gate(U::GateOps.GateLabel)
    su = GateOps.GateLabel( Symbol("s" * String(U.label)) )
    asu = GateOps.GateLabel( Symbol("as" * String(U.label)) )

    if haskey(Circuit.gate_cache, su)
        SU = Circuit.gate_cache[su]
        ASU = Circuit.gate_cache[asu]
    else
        SU = sqrt(Circuit.gate_cache[U])
        ASU = collect(adjoint(sqrt(Circuit.gate_cache[U])))

        # Cache the gates
        Circuit.gate_cache[su] = SU
        Circuit.gate_cache[asu] = ASU

    end

    return su,asu
end

# =========================================================================== #
#
# =========================================================================== #

"""
    apply_cx!(c::Circuit.Circ, ctrl, tgt, reg)

Module-specific CX gate. Defaults to using the implementation from GateOps.
Override for custom functionality.
"""
function apply_cx!(c::Circuit.Circ, ctrl, tgt, reg)
    Circuit.add_gatecall!(c, GateOps.c_pauli_x(ctrl, tgt, reg))
end

"""
    apply_cu!(c::Circuit.Circ, ctrl, tgt, reg, gl::GateOps.GateLabel)

Module-specific CU gate. Defaults to using the implementation from GateOps.
Override for custom functionality.
"""
function apply_cu!(c::Circuit.Circ, ctrl, tgt, reg, gl::GateOps.GateLabel)
    if String(gl.label)[end] == 'x'
        glz_s = String(gl.label)[1:end-1] * 'z'
        gl_z = GateOps.GateLabel(Symbol(glz_s))

        push!(c.gate_set, gl_z)

        Circuit.add_gatecall!(c, GateOps.hadamard(tgt, reg))
        Circuit.add_gatecall!(c, GateOps.c_u(gl_z, ctrl, tgt, GateOps.u( GateOps.GateLabel(Symbol(split(glz_s, "c_")[end])), tgt, reg), reg))
        Circuit.add_gatecall!(c, GateOps.hadamard(tgt, reg))

    elseif String(gl.label)[end] == 'z'
        Circuit.add_gatecall!(c, GateOps.c_u(gl, ctrl, tgt, GateOps.pauli_z(tgt, reg), reg))
    else
        error("Currently only PauliX and PauliZ decomposed gates supported")
    end
end

# =========================================================================== #
#
# =========================================================================== #

"""
    apply_ncu!(circuit::Circuit.Circ, q_ctrl::Vector, q_aux::Vector, q_tgt, U::GateOps.GateLabel)

Apply an n-qubit controlled gate operation on the given target qubit.
Ensure the gate corresponding with symbol U is registered with g_map before use.
Appends the GateCall operations to circuit

# Arguments
- `circuit::Circuit.Circ`
- `ctrls::Vector`:
- `aux::Vector`:
- `tgt::Int`:
- `U::GateOps.GateLabel`:
"""
function apply_ncu!(circuit::Circuit.Circ, q_ctrl::Vector, q_aux::Vector, q_tgt, U::GateOps.GateLabel)

    # Check global circuit cache for gate
    if ~haskey(Circuit.gate_cache, U )
        error("Gate $(U) does not exist")
    end

    su,asu = get_intermed_gate(U)
    push!(circuit.gate_set, su)
    push!(circuit.gate_set, asu)

    if length(q_ctrl) == 2

        apply_cu!(circuit, q_ctrl[2], q_tgt, nothing, su)
        apply_cx!(circuit, q_ctrl[1], q_ctrl[2], nothing)
        apply_cu!(circuit, q_ctrl[2], q_tgt,  nothing, asu)
        apply_cx!(circuit, q_ctrl[1], q_ctrl[2], nothing)
        apply_cu!(circuit, q_ctrl[1], q_tgt,  nothing, su)

    elseif length(q_ctrl) == 3

        ssu,assu = get_intermed_gate(su)
        push!(circuit.gate_set, ssu)
        push!(circuit.gate_set, assu)

        apply_cu!(circuit, q_ctrl[1], q_tgt, nothing, ssu)
        apply_cx!(circuit, q_ctrl[1], q_ctrl[2], nothing)
        apply_cu!(circuit, q_ctrl[2], q_tgt, nothing, assu)
        apply_cx!(circuit, q_ctrl[1], q_ctrl[2], nothing)
        apply_cu!(circuit, q_ctrl[2], q_tgt, nothing, ssu)
        apply_cx!(circuit, q_ctrl[2], q_ctrl[3], nothing)
        apply_cu!(circuit, q_ctrl[3], q_tgt, nothing, assu)
        apply_cx!(circuit ,q_ctrl[1], q_ctrl[3], nothing)
        apply_cu!(circuit, q_ctrl[3], q_tgt, nothing, ssu)
        apply_cx!(circuit, q_ctrl[2], q_ctrl[3], nothing)
        apply_cu!(circuit, q_ctrl[3], q_tgt, nothing, assu)
        apply_cx!(circuit, q_ctrl[1], q_ctrl[3], nothing)
        apply_cu!(circuit, q_ctrl[3], q_tgt, nothing, ssu)

    elseif (length(q_ctrl)>=5) && (length(q_aux)>=length(q_ctrl)-2) && (U == GateOps.GateLabel(:x))

        apply_ncu!(circuit, [q_ctrl[end], q_aux[ 1 + length(q_ctrl)-3 ]], Int[], q_tgt, U)
        for i in reverse(2:length(q_ctrl)-2)
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], U)
        end

        apply_ncu!(circuit, [q_ctrl[1], q_ctrl[2]], Int[], q_aux[1], U)
        for i in 2:length(q_ctrl)-2
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], U)
        end

        apply_ncu!(circuit, [q_ctrl[end], q_aux[1 + length(q_ctrl) - 3]], Int[], q_tgt, U)
        for i in reverse(2:length(q_ctrl)-2)
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], U)
        end

        apply_ncu!(circuit, [q_ctrl[1], q_ctrl[2]], Int[], q_aux[1], U)
        for i in 2:length(q_ctrl)-2
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], U)
        end

    elseif (length(q_ctrl)>=5) && (length(q_aux)>=length(q_ctrl)-2) && (U == GateOps.GateLabel(:z))

        gl = GateOps.GateLabel(:x)
        Circuit.add_gatecall!(circuit, GateOps.hadamard(q_tgt, nothing))

        apply_ncu!(circuit, [q_ctrl[end], q_aux[ 1 + length(q_ctrl)-3 ]], Int[], q_tgt, gl)
        for i in reverse(2:length(q_ctrl)-2)
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], gl)
        end

        apply_ncu!(circuit, [q_ctrl[1], q_ctrl[2]], Int[], q_aux[1], gl)
        for i in 2:length(q_ctrl)-2
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], gl)
        end

        apply_ncu!(circuit, [q_ctrl[end], q_aux[1 + length(q_ctrl) - 3]], Int[], q_tgt, gl)
        for i in reverse(2:length(q_ctrl)-2)
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], gl)
        end

        apply_ncu!(circuit, [q_ctrl[1], q_ctrl[2]], Int[], q_aux[1], gl)
        for i in 2:length(q_ctrl)-2
            apply_ncu!(circuit, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], gl)
        end

        Circuit.add_gatecall!(circuit, GateOps.hadamard(q_tgt, nothing))

    else

        apply_cu!(circuit, q_ctrl[end], q_tgt, nothing, su)
        apply_ncu!(circuit, q_ctrl[1:end-1], q_aux, q_ctrl[end], GateOps.GateLabel(:x))
        apply_cu!(circuit, q_ctrl[end], q_tgt, nothing, asu)
        apply_ncu!(circuit, q_ctrl[1:end-1], q_aux, q_ctrl[end], GateOps.GateLabel(:x))
        apply_ncu!(circuit, q_ctrl[1:end-1], q_aux, q_tgt, su)

    end

    return circuit

end

end

ext: jl


text:
using Flipping
using LsqFit
using StatsPlots

##colors
c9090 = RGBA(108/255,218/255,165/255)#RGBA(181/255,94/255,50/255)
c9030 = RGBA(23/255,116/255,206/255)#RGBA(181/255,73/255,144/255)
c3030 = RGBA(181/255,94/255,51/255)#RGBA(111/255,108/255,166/255)
gr(grid=false,background_color = RGBA(1,1,1,1))
##
#probability that the state has change after n_fail consecutive failures
function rec_lr(n_fail, Prew, Psw; E1=0)
    #= estimate P(depleted) / P(rich); as a ratio not conditioned probability
    n_fail = number of pokes
    Prew = probability of reward
    Psw = probability of state transition
    E1 = first evidence need to be given been calculated recursively
    after a reward the evidence of P(depleted) / P(rich) = 0
    =#
    Es = zeros(n_fail)
    Es[1] = E1
    for i = 2:n_fail
        Es[i] = (Es[i-1]+Psw)/((1-Prew)*(1-Psw))
    end
    return Es
end
function rec_lr_with_uncertainty(n_fail, Prew, Psw, α; E1 = 0)
    #n_fail = numero di poke
    #Prew = probabilitá reward
    #Psw = probabilitá di transizione
    #E1= la prima evidenza deve essere assegnata perché la probabilitaá è calcolata in maniera ricorsiva
    #α = uncertainty
    Es = zeros(n_fail)
    Es[1] = E1
    for i = 2:n_fail
        Es[i] = (Es[i-1]+Psw)/((1-Prew)*(1-Psw))*α + (1-α)
    end
    return Es
end
#Probability of being in the correct side after n_fail consecutive failures, requires to compute rec_lr
function Pwrong(evidence_accumulation)
    # x = pdpl / prich; x è la stima ottenuta dalla funzione rec_lr
    # x*(1-pdpl)  = pdpl
    # x - x*pdpl = pdpl
    # x = pdpl + x*pdpl
    # x = pdpl(1+x)
    # x/(1+x) = pdpl
    evidence_accumulation ./ (1 .+ evidence_accumulation)
end

function Pwrong_with_uncertainty(n_fail, Prew, Psw, α; E1 = 0)
    evidence_accumulation = rec_lr(n_fail, Prew, Psw, E1, α)
    evidence_accumulation ./ (1 .+ evidence_accumulation)
end

function Pwrong(n_fail, Prew, Psw; E1 = 0)
    evidence_accumulation = rec_lr(n_fail, Prew, Psw; E1 = E1)
    evidence_accumulation ./ (1 .+ evidence_accumulation)
end

function Pcorrect(n_fail, Prew, Psw; E1 = 0)
    1 .- Pwrong(n_fail, Prew, Psw; E1 = E1)
end
@. expon(x,p) = p[1]*exp(x*p[2]) ## model for an exponent for fitting
inverse_exp(y,p) = (log(y/p[1]))*1/p[2] ## model for an inverse exponent

function fit_protocol(prot)
    xdata = 0:size(prot,1)-1
    ydata = prot #evidence accumulation of n_fail consecutive omissions ::array
    p0 = [0.5, 0.5] # initialization of parameters
    fit = curve_fit(expon, xdata, ydata, p0)
end

function fit_protocol(n_fail, Prew, Psw; E1 = 0, step = 0.1)
    ydata = Pcorrect(n_fail, Prew, Psw; E1 = E1) #evidence accumulation of n_fail consecutive omissions ::array
    xdata = 0:size(ydata,1)-1
    p0 = [0.5, 0.5] # initialization of parameters
    fit = curve_fit(expon, xdata, ydata, p0)
    x = 0.0:step:n_fail
    (collect(x),expon(x,fit.param))
end
##
c = fit_protocol(4,0.8,0.4)
c
##
plot(fit_protocol(4,0.8,0.4))
##
colorscheme = Dict("80/40" =>RGB(0/255, 158/255, 115/255),
    "60/30" => RGB(86/255, 180/255, 233/255),
    "40/20" => RGB(230/255, 159/255, 0/255))
##
gr(grid=false,background_color = RGBA(1,1,1,1), linewidth = 3)
plt = plot(0:1:9,Pcorrect(10, 0.4, 0.2),#plot(fit_protocol(5,0.8,0.4; step = 1),
    xlabel = "Consecutive failures",
    ylabel =  "Probability of current side high",
    label = "40/20",
    color = colorscheme["40/20"])

plot!(plt, 0:1:6,Pcorrect(7, 0.6, 0.3),
    label= "60/30",
    color = colorscheme["60/30"])

plot!(plt, 0:1:4,Pcorrect(5, 0.8, 0.4),
    label = "80/40",
    color = colorscheme["80/40"])

Plots.abline!(plt,0,0.15,
    label = "example threshold 15%",
    color = :black)

Plots.abline!(plt,0,0.05,label = "example threshold 5%",
    linestyle = :dash, color = :black)
xticks!(plt,0:1:9)

##
dir = "/Volumes/GoogleDrive/My Drive/Flipping/shared docs pietro dario/lab meeting 3-02-2020/inference qualitative observations"
f = joinpath(dir,"updated_inference.pdf")
savefig(plt,f)
##
plt = plot(fit_protocol(10,0.3,0.15),
    xticks = 0:10,
    title  = "Evidence accumulation",
    xlabel = "Evidence",
    ylabel =  "Latent state accumulation P(correct)",
    color = :black,
    legend = false,
    linewidth = 3,
    tickfont = font(:Bookman,12),
    background_color = RGBA(1,1,1,0)
    )
##

##
Plots.abline!(plt,0,0.15,linewidth = 3,
    color = :black,
    annotations=(5, 0.18, text("Latent threshold", :left)))
savefig(plt,"/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/example2.png")
##
Plots.abline!(plt,0,0.05,linewidth = 3,
    color = :blue)
savefig(plt,"/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/example3.png")
##
plot!(fit_protocol(10,0.2,0.1),
    linewidth = 3,
    color = :blue)
savefig(plt,"/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/example4.png")


text:
using Oceananigans.BoundaryConditions: PBC, ZFBC, NFBC, BoundaryFunction, ParameterizedDiscreteBoundaryFunction

using Oceananigans.Fields: Face, Cell

function instantiate_boundary_function(B, X1, X2, func)
    boundary_function = BoundaryFunction{B, X1, X2}(func)
    return true
end

bc_location(bc::BoundaryFunction{B, X, Y}) where {B, X, Y} = (B, X, Y)

simple_bc(ξ, η, t) = exp(ξ) * cos(η) * sin(t)
simple_parameterized_bc(ξ, η, t, p) = p.a * exp(ξ) * cos(η) * sin(t)

complicated_bc(i, j, grid, clock, state) = rand()
complicated_parameterized_bc(i, j, grid, clock, state, p) = p.a * rand()

@testset "Boundary conditions" begin
    @info "Testing boundary conditions..."

    @testset "Boundary functions" begin
        @info "  Testing boundary functions..."

        for B in (:x, :y, :z)
            for X1 in (:Face, :Cell)
                @test instantiate_boundary_function(B, X1, Cell, simple_bc)
            end
        end

        bc = BoundaryCondition(Value, simple_bc)
        @test typeof(bc.condition) <: BoundaryFunction
        @test bc.condition.func === simple_bc

        bc = BoundaryCondition(Value, complicated_bc; discrete_form=true)
        @test bc.condition === complicated_bc

        bc = BoundaryCondition(Value, simple_parameterized_bc; parameters=(a=π,))
        @test bc.condition.parameters.a == π

        bc = BoundaryCondition(Value, complicated_parameterized_bc; parameters=(a=π,), discrete_form=true)
        @test typeof(bc.condition) <: ParameterizedDiscreteBoundaryFunction
        @test bc.condition.func === complicated_parameterized_bc
    end

    @testset "Field boundary conditions" begin
        @info "  Testing field boundary functions..."

        # Triply periodic
        ppp_topology = (Periodic, Periodic, Periodic)
        ppp_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=ppp_topology)

        u_bcs = UVelocityBoundaryConditions(ppp_grid)
        v_bcs = VVelocityBoundaryConditions(ppp_grid)
        w_bcs = WVelocityBoundaryConditions(ppp_grid)
        T_bcs = TracerBoundaryConditions(ppp_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa PBC
        @test u_bcs.y.right isa PBC
        @test u_bcs.z.left  isa PBC
        @test u_bcs.z.right isa PBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa PBC
        @test v_bcs.y.right isa PBC
        @test v_bcs.z.left  isa PBC
        @test v_bcs.z.right isa PBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa PBC
        @test w_bcs.y.right isa PBC
        @test w_bcs.z.left  isa PBC
        @test w_bcs.z.right isa PBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa PBC
        @test T_bcs.y.right isa PBC
        @test T_bcs.z.left  isa PBC
        @test T_bcs.z.right isa PBC

        # Doubly periodic. Engineers call this a "Channel geometry".
        ppb_topology = (Periodic, Periodic, Bounded)
        ppb_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=ppb_topology)

        u_bcs = UVelocityBoundaryConditions(ppb_grid)
        v_bcs = VVelocityBoundaryConditions(ppb_grid)
        w_bcs = WVelocityBoundaryConditions(ppb_grid)
        T_bcs = TracerBoundaryConditions(ppb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa PBC
        @test u_bcs.y.right isa PBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa PBC
        @test v_bcs.y.right isa PBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa PBC
        @test w_bcs.y.right isa PBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa PBC
        @test T_bcs.y.right isa PBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        # Singly periodic. Oceanographers call this a "Channel", engineers call it a "Pipe"
        pbb_topology = (Periodic, Bounded, Bounded)
        pbb_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=pbb_topology)

        u_bcs = UVelocityBoundaryConditions(pbb_grid)
        v_bcs = VVelocityBoundaryConditions(pbb_grid)
        w_bcs = WVelocityBoundaryConditions(pbb_grid)
        T_bcs = TracerBoundaryConditions(pbb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa ZFBC
        @test u_bcs.y.right isa ZFBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa NFBC
        @test v_bcs.y.right isa NFBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa ZFBC
        @test w_bcs.y.right isa ZFBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa ZFBC
        @test T_bcs.y.right isa ZFBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        # Triply bounded. Oceanographers call this a "Basin", engineers call it a "Box"
        bbb_topology = (Bounded, Bounded, Bounded)
        bbb_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=bbb_topology)

        u_bcs = UVelocityBoundaryConditions(bbb_grid)
        v_bcs = VVelocityBoundaryConditions(bbb_grid)
        w_bcs = WVelocityBoundaryConditions(bbb_grid)
        T_bcs = TracerBoundaryConditions(bbb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa NFBC
        @test u_bcs.x.right isa NFBC
        @test u_bcs.y.left  isa ZFBC
        @test u_bcs.y.right isa ZFBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa ZFBC
        @test v_bcs.x.right isa ZFBC
        @test v_bcs.y.left  isa NFBC
        @test v_bcs.y.right isa NFBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa ZFBC
        @test w_bcs.x.right isa ZFBC
        @test w_bcs.y.left  isa ZFBC
        @test w_bcs.y.right isa ZFBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa ZFBC
        @test T_bcs.x.right isa ZFBC
        @test T_bcs.y.left  isa ZFBC
        @test T_bcs.y.right isa ZFBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        grid = bbb_grid

        u_bcs = UVelocityBoundaryConditions(grid;
                                              east = BoundaryCondition(NormalFlow, simple_bc),
                                              west = BoundaryCondition(NormalFlow, simple_bc),
                                            bottom = BoundaryCondition(Value, simple_bc),
                                               top = BoundaryCondition(Value, simple_bc),
                                             north = BoundaryCondition(Value, simple_bc),
                                             south = BoundaryCondition(Value, simple_bc)
                                           )

        @test bc_location(u_bcs.x.east.condition) === (:x, Cell, Cell)
        @test bc_location(u_bcs.x.west.condition) === (:x, Cell, Cell)
        @test bc_location(u_bcs.y.north.condition) === (:y, Face, Cell)
        @test bc_location(u_bcs.y.south.condition) === (:y, Face, Cell)
        @test bc_location(u_bcs.z.top.condition) === (:z, Face, Cell)
        @test bc_location(u_bcs.z.bottom.condition) === (:z, Face, Cell)

        v_bcs = VVelocityBoundaryConditions(grid;
                                              east = BoundaryCondition(Value, simple_bc),
                                              west = BoundaryCondition(Value, simple_bc),
                                            bottom = BoundaryCondition(NormalFlow, simple_bc),
                                               top = BoundaryCondition(NormalFlow, simple_bc),
                                             north = BoundaryCondition(Value, simple_bc),
                                             south = BoundaryCondition(Value, simple_bc)
                                           )

        @test bc_location(v_bcs.x.east.condition) === (:x, Face, Cell)
        @test bc_location(v_bcs.x.west.condition) === (:x, Face, Cell)
        @test bc_location(v_bcs.y.north.condition) === (:y, Cell, Cell)
        @test bc_location(v_bcs.y.south.condition) === (:y, Cell, Cell)
        @test bc_location(v_bcs.z.top.condition) === (:z, Cell, Face)
        @test bc_location(v_bcs.z.bottom.condition) === (:z, Cell, Face)

        w_bcs = WVelocityBoundaryConditions(grid;
                                              east = BoundaryCondition(Value, simple_bc),
                                              west = BoundaryCondition(Value, simple_bc),
                                            bottom = BoundaryCondition(Value, simple_bc),
                                               top = BoundaryCondition(Value, simple_bc),
                                             north = BoundaryCondition(NormalFlow, simple_bc),
                                             south = BoundaryCondition(NormalFlow, simple_bc)
                                           )

        @test bc_location(w_bcs.x.east.condition) === (:x, Cell, Face)
        @test bc_location(w_bcs.x.west.condition) === (:x, Cell, Face)
        @test bc_location(w_bcs.y.north.condition) === (:y, Cell, Face)
        @test bc_location(w_bcs.y.south.condition) === (:y, Cell, Face)
        @test bc_location(w_bcs.z.top.condition) === (:z, Cell, Cell)
        @test bc_location(w_bcs.z.bottom.condition) === (:z, Cell, Cell)

        T_bcs = TracerBoundaryConditions(grid;
                                           east = BoundaryCondition(Value, simple_bc),
                                           west = BoundaryCondition(Value, simple_bc),
                                         bottom = BoundaryCondition(Value, simple_bc),
                                            top = BoundaryCondition(Value, simple_bc),
                                          north = BoundaryCondition(Value, simple_bc),
                                          south = BoundaryCondition(Value, simple_bc)
                                         )

        @test bc_location(T_bcs.x.east.condition) === (:x, Cell, Cell)
        @test bc_location(T_bcs.x.west.condition) === (:x, Cell, Cell)
        @test bc_location(T_bcs.y.north.condition) === (:y, Cell, Cell)
        @test bc_location(T_bcs.y.south.condition) === (:y, Cell, Cell)
        @test bc_location(T_bcs.z.top.condition) === (:z, Cell, Cell)
        @test bc_location(T_bcs.z.bottom.condition) === (:z, Cell, Cell)
    end
end


text:
@test isapprox(AutomotiveDrivingModels._mod2pi2(0.0), 0.0)
@test isapprox(AutomotiveDrivingModels._mod2pi2(0.5), 0.5)
@test isapprox(AutomotiveDrivingModels._mod2pi2(2pi + 1), 1.0)
@test isapprox(AutomotiveDrivingModels._mod2pi2(1 - 2pi), 1.0)


text:
import JuliaPackageWithRustDep
using Test

@testset "Rust FFI" begin

    # void
    JuliaPackageWithRustDep.rustdylib_printhello()

    # Int32
    @test JuliaPackageWithRustDep.rustdylib_abs_i32(Int32(10)) == 10
    @test JuliaPackageWithRustDep.rustdylib_abs_i32(Int32(-10)) == 10

    # Int64
    @test JuliaPackageWithRustDep.rustdylib_abs_i64(typemax(Int64)) == typemax(Int64)
    @test JuliaPackageWithRustDep.rustdylib_abs_i64(-9223372036854775807) == 9223372036854775807

    # Float32
    @test JuliaPackageWithRustDep.rustdylib_abs_f32(Float32(1.0)) == 1.0
    @test JuliaPackageWithRustDep.rustdylib_abs_f32(Float32(-1.0)) == 1.0

    # Float64
    @test JuliaPackageWithRustDep.rustdylib_abs_f64(Float64(1.0)) == 1.0
    @test JuliaPackageWithRustDep.rustdylib_abs_f64(Float64(-1.0)) == 1.0

    # Bool
    @test JuliaPackageWithRustDep.rustdylib_is_true_bool(true) == true
    @test JuliaPackageWithRustDep.rustdylib_is_true_bool(false) == false

    JuliaPackageWithRustDep.rustdylib_inspect_string("hey you")
    JuliaPackageWithRustDep.rustdylib_inspect_string("😻   é∈δ")

    println("Julia read: ", JuliaPackageWithRustDep.read_rust_owned_string())
end


text:
const modlobP = [
    -1.0,
    -0.81649658092772603273242802490196,
    -.44721359549995793928183473374626,
    0,
    0.44721359549995793928183473374626,
    0.81649658092772603273242802490196,
    1.0,
]

const modlobWei = [
    0.052380952380952380952380952380952,
    0.29387755102040816326530612244898,
    0.42517006802721088435374149659864,
    0.45714285714285714285714285714286,
]

const modlobW = [
    0.052380952380952380952380952380952,
    0.29387755102040816326530612244898,
    0.42517006802721088435374149659864,
    0.45714285714285714285714285714286,
    0.42517006802721088435374149659864,
    0.29387755102040816326530612244898,
    0.052380952380952380952380952380952,
]

const modlobNulw = [
    0.1148787002745345 -0.2954023721345174 0.4102810724090519 -0.4595148010981382
    0.2168775557284140 -0.4553480127399748 0.3463949695733805 0.0
    0.2957100706945793 -0.4080798975585194 -0.1182840282778318 0.4613077102835437
    0.3560582128767142 -0.1825453130808602 -0.4628897488019905 0.0
]
modlobNw = zeros(Float64, (4, 7))

for j = 1:4
    for k = 1:4
        modlobNw[j, k] = modlobNulw[j, k]
    end
    for k = 1:3
        modlobNw[j, 8-k] = modlobNulw[j, k]
    end
end
for k = 5:7
    modlobNw[2, k] = -modlobNw[2, k]
end
for k = 5:7
    modlobNw[4, k] = -modlobNw[4, k]
end

#Espelid 'modlob'
function modlob(f, a::T, b::T; tol = sqrt(eps(T))) where {T}
    #  %MODLOB Numerically evaluates an integral using adaptively
    #  % a Lobatto 7 point rule developed by Gander and Gautschi.
    #  %
    #  % Q=MODLOB('F',A,B) approximates the integral of F(X) from
    # % A to B to a relative precision TOL= 10^-3. 'F' is a
    #  % string containing the name of the function. The
    # % function F must return a vector of output values if
    #  % given a vector of input values.
    # % Q=MODLOB('F',A,B,TOL) integrates to a relative
    #  % error of TOL.
    #  %
    #  % Q=MODLOB('F',A,B,TOL,TRACE) displays the left
    #  % end point of the current interval, the interval
    #  % length, and the partial integral.
    #  %
    # % Q=MODLOB('F',A,B,TOL,TRACE,P1,P2,...) allows
    #  % coefficients P1, ... to be passed directly to the
    #  % function F: G=F(X,P1,P2,...). To use default values
    #  % for TOL or TRACE, one may pass the empty matrix ([]).
    #  %
    #  % Terje O. Espelid 31/01/02
    #  % Reference: Gander and Gautschi: Bit 2000
    #  % Berntsen and Espelid, TOMS 1991
    #  % Espelid, in Espelid and Genz, 1992.
    #  % Espelid, Tech report May 2002.

    if tol <= 0
        tol = eps(T)
    end
    # % 7 points Lobatto rule in the interval [-1,1]; we give only 4 points due to symmetry.
    # % 4 Orthogonal Nullrules of degree 5,4,3 and 2 are given. The odd numbered are symmetric.
    # % All even numbered are antisymmetric. We given first the five elements.
    # % Then we use the symmetry
    # % Thus: p is the points, w is the weights and
    # % nw is the nullrules associated with the 7 points Lobatto rule.
    h = (b - a) / 2
    center = (b + a) / 2
    x = Vector{T}(undef, 7)
    y = Vector{T}(undef, 7)
    is = zero(T)
    isabs = zero(T)
    for j = 1:length(x)
        x[j] = center + h * modlobP[j]
        y[j] = f(x[j])
        is += h * modlobW[j] * y[j]
        isabs += h * modlobW[j] * abs(y[j])
    end
    nEval = 32 #not more than 32 levels
    fa = y[1]
    fb = y[7]
    c = x[4]
    fc = y[4]
    if isnan(fa)
        throw(DomainError(string("function is NaN at ", a)))
    elseif isnan(fb)
        throw(DomainError(string("function is NaN at ", b)))
    elseif isnan(fc)
        throw(DomainError(string("function is NaN at ", c)))
    end
    # % Apply the seven point rule to compute an estimate of the integral  to be used in the
    # % error estimator.
    # % Define the noise level
    noise = 50 * eps(T) * isabs
    # % Check that the estimate is not too small.
    absis = abs(is)
    if absis <= noise
        is = b - a
        absis = is
    end
    # % If tol is chosen too small, redefine tol too be on the noise level.
    tol = max(tol, noise / absis)
    return lobattor(f, a, c, fa, fc, is, tol, nEval) +
           lobattor(f, c, b, fc, fb, is, tol, nEval)
end

function lobattor(f, a::T, b::T, fa::T, fb::T, is::T, tol::T, nEval::Int)::T where {T}
    # %LOBATTOR Recursive function used by Lobatto.
    # %
    # % Q = LOBATTOR('F',A,B,FA,FB,IS,TOL) tries to
    # % approximate the integral of F(X) from A to B to
    # % an appropriate relative error. The argument 'F' is
    # % a string containing the name of f. The remaining
    # % arguments are generated by MODLOB or by recursion.
    # %
    #
    h = (b - a) / 2
    center = (b + a) / 2
    x = (
        center + h * modlobP[1],
        center + h * modlobP[2],
        center + h * modlobP[3],
        center + h * modlobP[4],
        center + h * modlobP[5],
        center + h * modlobP[6],
        center + h * modlobP[7],
    )
    z = (x[2], x[3], x[4], x[5], x[6])
    y = (fa, f(z[1]), f(z[2]), f(z[3]), f(z[4]), f(z[5]), fb)
    c = x[4]
    fc = y[4]
    # i = h* (w' * y)
    # iabs = h* (w' * abs.(y))
    # E = abs(h * (modelobNw' * y))
    # Emin = min(E[2:4])
    i = zero(T)
    iabs = zero(T)

    E = zeros(T, 4)
    for j = 1:length(y)
    	i += h * modlobW[j] * y[j]
    	iabs += h * modlobW[j] * abs(y[j])
    	for k = 1:length(E)
    		E[k] += h * modlobNw[k,j] * y[j]
    	end
    end
    @. E = abs(E)
    noise = 50 * iabs * eps(T)
    Emin = minimum(E[2:4])
    rmax = 2 * one(T)
    # Compute the local error estimate
    if Emin != 0
        # r=E(1:3)./E(2:4);
        # if (sum(isinf(r))>0)
        # rmax=2;
        # else
        # rmax=max(r);
        r1, r2, r3 = E[1] / E[2], E[2] / E[3], E[3] / E[4]
        if isinf(r1) || isinf(r2) || isinf(r3) || isnan(r1) || isnan(r2) || isnan(r3)
            #rmax = 2
        else
            rmax = max(r1, r2, r3)
        end
    end
    err = zero(T)
    C = 10
    if rmax > 1
        err = C * maximum(E)
    elseif one(T) / 2 < rmax
        err = C * rmax * E[2]
    else
        err = C * (2* rmax)^4 * rmax * E[2]
    end
    # % If the best null-rules are on the noise level, then put the error to zero.
    if (E[1] < noise) && (E[2] < noise)
        err = zero(T)
    end
    # % Check if the error is small enough or the interval is too small.
    if (err <= max(tol * abs(is), noise)) || (x[2] <= a) || (b <= x[6])
        if ((x[2] <= a) || (b <= x[6]))
            #("Interval too small: required tolerance may not be met.")
        end
        return i
    end
    if nEval == 0
        return i
    end
    if isnan(fa) || isnan(fb) || isnan(fc)
        return T(NaN)
    end
    return lobattor(f, a, c, fa, fc, is, tol, nEval-1) +
           lobattor(f, c, b, fc, fb, is, tol, nEval-1)
end


text:
function lglnodes(N,a,b)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %
# % lglnodes.m
# %
# % Computes the Legendre-Gauss-Lobatto nodes, weights and the LGL Vandermonde
# % matrix. The LGL nodes are the zeros of (1-x^2)*P'_N(x). Useful for numerical
# % integration and spectral methods.
# %
# % Reference on LGL nodes and weights:
# %   C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods
# %   in Fluid Dynamics," Section 2.3. Springer-Verlag 1987
# %
# % Written by Greg von Winckel - 04/17/2004
# % Contact: gregvw@chtm.unm.edu
# %
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# % Truncation + 1
	N1=N+1;

	#% Use the Chebyshev-Gauss-Lobatto nodes as the first guess
	x=cos.(pi*collect(0:N)/N);

	# The Legendre Vandermonde Matrix
	P=zeros(N1,N1);

	# 	% Compute P_(N) using the recursion relation
	# 	% Compute its first and second derivatives and
	# 	% update x using the Newton-Raphson method.

	xold=2;

	while maximum(abs.(x .- xold))>eps()

	    xold=x;

	    P[:,1] .= 1;    P[:,2] = x;

	    for k=2:N
		P[:,k+1]=( (2*k-1)*x.*P[:,k]-(k-1)*P[:,k-1] )/k;
	    end

	    x=xold-( x.*P[:,N1]-P[:,N] )./( N1*P[:,N1] );

	end

	w=2 ./ (N*N1*P[:,N1].^2)/2 .* (b-a);
	x = (x .+ 1)*(b-a)/2.0 .+ a
	return x,w,P
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
using ThermalMemoryEraser
using Test
import ThermalMemoryEraser: potential, virtual_workrate, erased

include("spring_force.jl")

function multi_ensemble()

    nsample = 2
    T = 1.0
    dt = 1.e-2
    τ = 1.0
    f = SpringForce(1.0)

    h = Heun(step=dt, temperature=T)

    init = [10.0, 5.0]

    configurations = [Configuration(x, f, h, τ) for x in init]

    Ensemble(nsample, f, h, init, configurations)
end

function main()
    ens = multi_ensemble()

    er = erasure_rate(ens)
    @info "error rate" er

    w, q = workheatlist(ens)
    @info work=w, heat=q
    @test w !== nothing
    @test q !== nothing
end

@test potential(SpringForce(1), 1, 0.2) !== nothing

main()


text:
# ===== jacobian ===== #
import ForwardDiff: JACOBIAN_ERROR, reshape_jacobian, chunksize

function jacobian(f, z::AbstractArray{Complex{V}},
            cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, z)) where {T, V}
    #CHK && checktag(T, f, z)
    checktag(T, f, z)
    if chunksize(cfg) == length(z)
      vector_mode_jacobian(f, z, cfg)
    else
      chunk_mode_jacobian(f, z, cfg)
    end
end

function extract_jacobian!(::Type{T}, dz::AbstractArray, dzstar::AbstractArray,
                ydual::AbstractArray{<:ComplexDual{T,V,N}}, n) where {T,V,N}
    @assert size(dz) == size(dzstar) "Inconsistent result sizes"

    dz_reshaped = reshape(dz, length(ydual), n)
    dzstar_reshaped = reshape(dzstar, length(ydual), n)
    for col in 1:size(dz_reshaped, 2), row in 1:size(dz_reshaped, 1)
        dz_reshaped[row, col], dzstar_reshaped[row, col] = dz_partials(T, ydual[row], col)
    end
    return dz, dzstar
end

function vector_mode_jacobian(f::F, z, cfg::ComplexGradientConfig{T,V,N},
            evalfcn=ForwardDiff.vector_mode_dual_eval) where {F,T,V,N}
    ydual = evalfcn(f, z, cfg)
    dz = similar(ydual, Complex{valtype(eltype(ydual))}, length(ydual), N)
    dzstar = similar(ydual, Complex{valtype(eltype(ydual))}, length(ydual), N)
    extract_jacobian!(T, dz, dzstar, ydual, N)
    return dz, dzstar
end

function extract_jacobian_chunk!(::Type{T}, dz, dzstar, ydual, index, chunksize) where {T}
    ydual_reshaped = vec(ydual)
    offset = index - 1
    irange = 1:chunksize
    col = irange .+ offset

    # Use closure to avoid GPU broadcasting with Type
    dz_partials_wrap(ydual, nrange) = dz_partials(T, ydual, nrange)
    dztmp = collect(Iterators.flatten(dz_partials_wrap.(ydual_reshaped, transpose(irange))))
    len, ylen = length(dztmp), length(ydual_reshaped)

    dz[:,col] .= reshape(view(dztmp,1:2:len-1),ylen,chunksize)
    dzstar[:,col] .= reshape(view(dztmp,2:2:len),ylen,chunksize)

    #for icol in irange, row in 1:size(dz, 1)
    #    dz[row, icol+offset], dzstar[row, icol+offset] = dz_partials(T, ydual[row], icol)
    #end
    return dz, dzstar
end

@eval function chunk_mode_jacobian(f::F, z, cfg::ComplexGradientConfig{T,V,N}) where {F,T,V,N}
    $(chunk_mode_gradient_expr(quote
                                   zdual = cfg.duals
                                   seed!(zdual, z)
                               end,
                               quote
                                 dz = similar(z, Complex{valtype(ydual)},length(ydual), zlen)
                                 dzstar = similar(z, Complex{valtype(ydual)},length(ydual), zlen)
                                 dz_reshaped = reshape_jacobian(dz, ydual, zdual)
                                 dzstar_reshaped = reshape_jacobian(dzstar, ydual, zdual)
                               end,
                               :(dz, dzstar),
                               :(ydual = f(zdual)),
                               :(ydual isa AbstractArray || throw(JACOBIAN_ERROR)),
                               :(seed!(zdual, z, index, rseeds, iseeds, chunksize)),
                               :(seed!(zdual, z, index)),
                               :(extract_jacobian_chunk!(T, dz_reshaped, dzstar_reshaped, ydual, index, chunksize)),
                               :()))
end


text:
# ---
# title: 540. Single Element in a Sorted Array
# id: problem540
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/single-element-in-a-sorted-array/description/>
# hidden: true
# ---
#
# You are given a sorted array consisting of only integers where every element
# appears exactly twice, except for one element which appears exactly once. Find
# this single element that appears only once.
#
# **Follow up:** Your solution should run in O(log n) time and O(1) space.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,2,3,3,4,4,8,8]
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,3,7,7,10,11,11]
#     Output: 10
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `0 <= nums[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct Spline
    α::AbstractVector
    k::Integer
end

#=
\[f(x)=\sum_i\alpha_iB_{i,k,t}(x)\equiv B(x)\vec{\alpha} = \langle B|\vec{\alpha}\rangle.\]
=#

(S::Spline)(Bᵢ::AbstractMatrix) = Bᵢ*S.α
function (S::Spline)(B::Vector{M}) where {M<:AbstractMatrix}
    kB = length(B)
    S.k > kB && error("No support for splines of order $(S.k) in basis of order $(kB)")
    S(B[S.k])
end
(S::Spline)(B::Basis) =  S(B.B)

export Spline


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

import Base: _uv_hook_close, unsafe_convert,
    lock, trylock, unlock, islocked

export SpinLock, RecursiveSpinLock, Mutex

##########################################
# Atomic Locks
##########################################

"""
    AbstractLock

Abstract supertype describing types that
implement the thread-safe synchronization primitives:
`lock`, `trylock`, `unlock`, and `islocked`
"""
abstract type AbstractLock end

# Test-and-test-and-set spin locks are quickest up to about 30ish
# contending threads. If you have more contention than that, perhaps
# a lock is the wrong way to synchronize.
"""
    TatasLock()

See SpinLock.
"""
struct TatasLock <: AbstractLock
    handle::Atomic{Int}
    TatasLock() = new(Atomic{Int}(0))
end

"""
    SpinLock()

Creates a non-reentrant lock.
Recursive use will result in a deadlock.
Each `lock` must be matched with an `unlock`.

Test-and-test-and-set spin locks are quickest up to about 30ish
contending threads. If you have more contention than that, perhaps
a lock is the wrong way to synchronize.

See also RecursiveSpinLock for a version that permits recursion.

See also Mutex for a more efficient version on one core or if the lock may be held for a considerable length of time.
"""
const SpinLock = TatasLock

function lock(l::TatasLock)
    while true
        if l.handle[] == 0
            p = atomic_xchg!(l.handle, 1)
            if p == 0
                return
            end
        end
        ccall(:jl_cpu_pause, Void, ())
        # Temporary solution before we have gc transition support in codegen.
        ccall(:jl_gc_safepoint, Void, ())
    end
end

function trylock(l::TatasLock)
    if l.handle[] == 0
        return atomic_xchg!(l.handle, 1) == 0
    end
    return false
end

function unlock(l::TatasLock)
    l.handle[] = 0
    ccall(:jl_cpu_wake, Void, ())
    return
end

function islocked(l::TatasLock)
    return l.handle[] != 0
end

"""
    RecursiveTatasLock()

See RecursiveSpinLock.
"""
struct RecursiveTatasLock <: AbstractLock
    ownertid::Atomic{Int16}
    handle::Atomic{Int}
    RecursiveTatasLock() = new(Atomic{Int16}(0), Atomic{Int}(0))
end

"""
    RecursiveSpinLock()

Creates a reentrant lock.
The same thread can acquire the lock as many times as required.
Each `lock` must be matched with an `unlock`.

See also SpinLock for a slightly faster version.

See also Mutex for a more efficient version on one core or if the lock may be held for a considerable length of time.
"""
const RecursiveSpinLock = RecursiveTatasLock

function lock(l::RecursiveTatasLock)
    if l.ownertid[] == threadid()
        l.handle[] += 1
        return
    end
    while true
        if l.handle[] == 0
            if atomic_cas!(l.handle, 0, 1) == 0
                l.ownertid[] = threadid()
                return
            end
        end
        ccall(:jl_cpu_pause, Void, ())
        # Temporary solution before we have gc transition support in codegen.
        ccall(:jl_gc_safepoint, Void, ())
    end
end

function trylock(l::RecursiveTatasLock)
    if l.ownertid[] == threadid()
        l.handle[] += 1
        return true
    end
    if l.handle[] == 0
        if atomic_cas!(l.handle, 0, 1) == 0
            l.ownertid[] = threadid()
            return true
        end
        return false
    end
    return false
end

function unlock(l::RecursiveTatasLock)
    @assert(l.ownertid[] == threadid(), "unlock from wrong thread")
    @assert(l.handle[] != 0, "unlock count must match lock count")
    if l.handle[] == 1
        l.ownertid[] = 0
        l.handle[] = 0
        ccall(:jl_cpu_wake, Void, ())
    else
        l.handle[] -= 1
    end
    return
end

function islocked(l::RecursiveTatasLock)
    return l.handle[] != 0
end

##########################################
# System Mutexes
##########################################

# These are mutexes from libuv. We're doing some error checking (and
# paying for it in overhead), but regardless, in some situations,
# passing a bad parameter will cause an abort.

# TODO: how defensive to get, and how to turn it off?
# TODO: how to catch an abort?

const UV_MUTEX_SIZE = ccall(:jl_sizeof_uv_mutex, Cint, ())

"""
    Mutex()

These are standard system mutexes for locking critical sections of logic.

On Windows, this is a critical section object,
on pthreads, this is a `pthread_mutex_t`.

See also SpinLock for a lighter-weight lock.
"""
mutable struct Mutex <: AbstractLock
    ownertid::Int16
    handle::Ptr{Void}
    function Mutex()
        m = new(zero(Int16), Libc.malloc(UV_MUTEX_SIZE))
        ccall(:uv_mutex_init, Void, (Ptr{Void},), m.handle)
        finalizer(m, _uv_hook_close)
        return m
    end
end

unsafe_convert(::Type{Ptr{Void}}, m::Mutex) = m.handle

function _uv_hook_close(x::Mutex)
    h = x.handle
    if h != C_NULL
        x.handle = C_NULL
        ccall(:uv_mutex_destroy, Void, (Ptr{Void},), h)
        Libc.free(h)
        nothing
    end
end

function lock(m::Mutex)
    if m.ownertid == threadid()
        return
    end
    # Temporary solution before we have gc transition support in codegen.
    # This could mess up gc state when we add codegen support.
    gc_state = ccall(:jl_gc_safe_enter, Int8, ())
    ccall(:uv_mutex_lock, Void, (Ptr{Void},), m)
    ccall(:jl_gc_safe_leave, Void, (Int8,), gc_state)
    m.ownertid = threadid()
    return
end

function trylock(m::Mutex)
    if m.ownertid == threadid()
        return true
    end
    r = ccall(:uv_mutex_trylock, Cint, (Ptr{Void},), m)
    if r == 0
        m.ownertid = threadid()
    end
    return r == 0
end

function unlock(m::Mutex)
    @assert(m.ownertid == threadid(), "unlock from wrong thread")
    m.ownertid = 0
    ccall(:uv_mutex_unlock, Void, (Ptr{Void},), m)
    return
end

function islocked(m::Mutex)
    return m.ownertid != 0
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
sumabsdiff(z1::Vector{T}, z2::Vector{T}) where T <: Integer = sum(abs.(z1 - z2))

function meanabsdiff(z1::Vector{T}, z2::Vector{T}) where T <: Integer
  return sumabsdiff(z1, z2) / size(z1, 1)
end

function closest_column(z1::Vector{T}, Z2::Matrix{T};
                        z_diff=sumabsdiff) where T<: Integer

  # Compute distance between each column
  d = [z_diff(z1, Z2[:, k2]) for k2 in 1:size(Z2, 2)]

  k2 = argmin(d)
  return k2, d[k2]
end

"""
compute similarity of two feature matrices

Arguments
=========
- Z1: point estimate of Z from FAM (with K1 columns)
- W1: point estimate of W from FAM (with K1 columns)
- Z2: point estimate of Z from FAM (with K2 columns)
"""
function similarity_FM(Z1::Matrix{T}, W1::Matrix{F}, Z2::Matrix{T};
                       z_diff=sumabsdiff) where {T <: Integer,
                                                 F <: AbstractFloat}
  # cache dimensions
  I, K1 = size(W1)
  J = size(Z1, 1)
  K2 = size(Z2, 2)

  # Make sure dimensions match
  @assert size(Z1, 1) == size(Z2, 1) == J
  @assert size(Z1, 2) == size(W1, 2) == K1

  # Iterate over columns of Z1 -- the larger matrix
  @assert K1 >= K2

  function engine(z1::Matrix{T}, w1::Vector{F}, z2::Matrix{T}; s::Float64=0.0)
    if length(z1) == 0
      return s
    elseif length(z2) == 0
      return engine(z1[:, 2:end], w1[2:end], Z2, s=s)
    else
      k2, d = closest_column(z1[:, 1], z2, z_diff=z_diff)
      engine(z1[:, 2:end], w1[2:end], z2[:, 1:end .!= k2], s=s+d*w1[1])
    end
  end

  s = sum(begin
            ord_1 = sortperm(W1[i, :], rev=true)
            engine(Z1[:, ord_1], W1[i, ord_1], Z2)
          end for i in 1:I)

  return s
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
#!/usr/bin/env julia

# This is an example invocation of `debootstrap` to generate a Debian/Ubuntu-based rootfs
using Scratch, Pkg, Pkg.Artifacts, ghr_jll, SHA, Dates

# Utility functions
getuid() = ccall(:getuid, Cint, ())
getgid() = ccall(:getgid, Cint, ())

function debootstrap(name::String; release::String="buster", variant::String="minbase",
                     packages::Vector{String}=String[], force::Bool=false)
    if Sys.which("debootstrap") === nothing
        error("Must install `debootstrap`!")
    end

    tarball_path = joinpath(@get_scratch!("rootfs-images"), "$(name).tar.gz")
    if !force && isfile(tarball_path)
        @error("Refusing to overwrite tarball without `force` set", tarball_path)
        error()
    end

    artifact_hash = create_artifact() do rootfs
        packages_string = join(push!(packages, "locales"), ",")
        @info("Running debootstrap", release, variant, packages)
        run(`sudo debootstrap --variant=$(variant) --include=$(packages_string) $(release) "$(rootfs)"`)

        # Remove special `dev` files
        @info("Cleaning up `/dev`")
        for f in readdir(joinpath(rootfs, "dev"); join=true)
            # Keep the symlinks around (such as `/dev/fd`), as they're useful
            if !islink(f)
                run(`sudo rm -rf "$(f)"`)
            end
        end

        # take ownership of the entire rootfs
        @info("Chown'ing rootfs")
        run(`sudo chown $(getuid()):$(getgid()) -R "$(rootfs)"`)

        # Write out rootfs-info to contain a minimally-identifying string
        open(joinpath(rootfs, "etc", "rootfs-info"), write=true) do io
            write(io, """
            rootfs_type=debootstrap
            release=$(release)
            variant=$(variant)
            packages=$(packages_string)
            build_date=$(Dates.now())
            """)
        end

        # Write out a reasonable default resolv.conf
        open(joinpath(rootfs, "etc", "resolv.conf"), write=true) do io
            write(io, """
            nameserver 1.1.1.1
            nameserver 8.8.8.8
            """)
        end

        # Remove `_apt` user so that `apt` doesn't try to `setgroups()`
        @info("Removing `_apt` user")
        open(joinpath(rootfs, "etc", "passwd"), write=true, read=true) do io
            filtered_lines = filter(l -> !startswith(l, "_apt:"), readlines(io))
            truncate(io, 0)
            seek(io, 0)
            for l in filtered_lines
                println(io, l)
            end
        end

        # Set up the one true locale
        @info("Setting up UTF-8 locale")
        open(joinpath(rootfs, "etc", "locale.gen"), "a") do io
            println(io, "en_US.UTF-8 UTF-8")
        end
        run(`sudo chroot --userspec=$(getuid()):$(getgid()) $(rootfs) locale-gen`)
    end

    # Archive it into a `.tar.gz` file
    @info("Archiving", tarball_path, artifact_hash)
    archive_artifact(artifact_hash, tarball_path)

    return tarball_path
end

function upload_rootfs_image(tarball_path::String; github_repo::String="JuliaCI/rootfs-images")
    # Upload it to `github_repo`
    tag_name = "v1"
    tarball_url = "https://github.com/$(github_repo)/releases/download/$(tag_name)/$(basename(tarball_path))"
    @info("Uploading to $(github_repo)@$(tag_name)", tarball_url)
    run(`$(ghr_jll.ghr()) -u $(dirname(github_repo)) -r $(basename(github_repo)) -replace $(tag_name) $(tarball_path)`)
    return tarball_url
end

max_stars_repo_name: chengchingwen/julia

max_issues_repo_name: chengchingwen/julia

max_forks_repo_name: GregPlowman/julia


text:
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4908336/
struct SegmentedLSE{T, U} <: AggregationFunction
    p::T
    C::U
end

Flux.@functor SegmentedLSE

SegmentedLSE(d::Int) = SegmentedLSE(randn(Float32, d), zeros(Float32, d))

Base.show(io::IO, n::SegmentedLSE) = print(io, "SegmentedLSE($(length(n.p)))\n")
modelprint(io::IO, n::SegmentedLSE; pad=[]) = paddedprint(io, "SegmentedLSE($(length(n.p)))")

(m::SegmentedLSE)(x::MaybeMatrix, bags::AbstractBags, w=nothing) = segmented_lse_optim(x, m.C, m.p, bags)
function (m::SegmentedLSE)(x::AbstractMatrix, bags::AbstractBags, w::AggregationWeights, mask::AbstractVector)
    segmented_lse_optim(x .+ typemin(T) * mask', m.C, m.p, bags)
end

segmented_lse_optim(x::Missing, C::AbstractVector, p::AbstractVector, bags::AbstractBags) = segmented_lse_forw(x, C, bags)
function segmented_lse_optim(x::AbstractMatrix, C::AbstractVector, p::AbstractVector, bags::AbstractBags)
    a = p .* x
    m = maximum(a, dims=2)
    y = (m .+ segmented_lse_forw(a .- m, C, bags)) ./ p
end

segmented_lse_forw(::Missing, C::AbstractVector, bags::AbstractBags) = repeat(C, 1, length(bags))
function segmented_lse_forw(a::AbstractMatrix, C::AbstractVector, bags::AbstractBags)
    y = zeros(eltype(a), size(a, 1), length(bags))
    e = exp.(a)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                y[i, bi] = C[i]
            end
        else
            for j in b
                for i in 1:size(a, 1)
                    y[i, bi] += e[i, j]
                end
            end
            lb = log(max(1, length(b)))
            for i in 1:size(a, 1)
                y[i, bi] = log(y[i, bi]) - lb
            end
        end
    end
    y
end

function segmented_lse_back(Δ, y, a, C, bags)
    da = zero(a)
    s = zero(C)
    dC = zero(C)
    e = exp.(a)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                dC[i] += Δ[i, bi]
            end
        else
            s .= 0
            for j in b
                for i in 1:size(a, 1)
                    da[i, j] = Δ[i, bi] * e[i, j]
                    s[i] += e[i, j]
                end
            end
            da[:, b] ./= s
        end
    end
    da, dC, nothing, nothing
end

function segmented_lse_back(Δ, C, bags)
    dC = zero(C)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(C)
            dC[i] += Δ[i, bi]
        end
    end
    nothing, dC, nothing, nothing
end

@adjoint function segmented_lse_forw(a::AbstractMatrix, C::AbstractVector, bags::AbstractBags)
    y = segmented_lse_forw(a, C, bags)
    grad = Δ -> segmented_lse_back(Δ, y, a, C, bags)
    y, grad
end

@adjoint function segmented_lse_forw(a::Missing, C::AbstractVector, bags::AbstractBags)
    y = segmented_lse_forw(a, C, bags)
    grad = Δ -> segmented_lse_back(Δ, C, bags)
    y, grad
end


text:
@testset "item.jl" begin
    sup = Supplier(100, 10, leadtime = LeadTime(1, 3, 10))
    inv = Inventory(3, 50)
    item = Item(inv, sup)
    @test state(item) == [50,10,10,10]
    @test InventoryModels.inventory_position(item) == 80
    @test state_size(item) == 4
    @test action_size(item) == 2
    dest = []
    InventoryModels.pull!(item, 30, dest)
    item([30])
    InventoryModels.dispatch!(item)
    @test InventoryModels.reward!(item) == -(100 + 10*30 + 30*3 + 50)
    @test state(item) == [30,10,10,30]
    @test dest == [30, inv]
    InventoryModels.pull!(item, 40, dest)
    item([0])
    InventoryModels.dispatch!(item)
    @test InventoryModels.reward!(item) == -(40)
    @test state(item) == [0,10,30,0]
    @test dest == [30, inv, 40, inv]
    InventoryModels.pull!(item, 20, dest)
    item([0])
    InventoryModels.dispatch!(item)
    @test InventoryModels.reward!(item) == -(30)
    @test state(item) == [0,30,0,0]
    @test dest == [30, inv, 40, inv, 10, inv]
    reset!(item)
    @test state(item) == [50,10,10,10]
    @test InventoryModels.children(item) == ()

    item2 = deepcopy(item)
    inv3 = Inventory(7, 30)
    ass3 = Assembly(100, 1, item => 2, item2 => 4)
    item3 = Item(inv3, ass3)
    @test state(item3) == [30]
    @test InventoryModels.children(item3) == (item, item2)
    dest = []
    InventoryModels.pull!(item3, 20, dest)
    item3([10])
    item2([0])
    item([0])
    InventoryModels.dispatch!(item)
    InventoryModels.dispatch!(item2)
    InventoryModels.dispatch!(item3)
    InventoryModels.reward!(item)
    InventoryModels.reward!(item2)
    InventoryModels.reward!(item3)
    @test state(item) == [40,10,10,0]
    @test state(item2) == [20,10,10,0]
    @test state(item3) == [20]
    @test dest == [20, inv3]
end


text:
using DashBootstrapComponents, DashHtmlComponents

modal = html_div([
    dbc_button("Open modal", id="open", n_clicks=0),
    dbc_modal(
        [
            dbc_modalheader("Header"),
            dbc_modalbody("This is the content of the modal"),
            dbc_modalfooter(
                dbc_button("Close", id="close", className="ml-auto", n_clicks=0),
            ),
        ],
        id="modal",
        is_open=false,
    ),
]);

callback!(
    app,
    Output("modal", "is_open"),
    Input("open", "n_clicks"),
    Input("close", "n_clicks"),
    State("modal", "is_open"),
) do n1, n2, is_open
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using InvertedIndices, StatsBase

export imcmc_prop_sample_edit, imcmc_prop_sample_edit_informed
export imcmc_prop_sample_edit!, accept_reject_edit!, delete_insert!

function delete_insert!(
    x::Path,
    δ::Int, d::Int,
    ind_del::AbstractArray{Int},
    ind_add::AbstractArray{Int},
    vals::AbstractArray{Int})

    @views for (i, index) in enumerate(ind_del[1:d])
        deleteat!(x, index - i + 1)
    end
    @views for (index, val) in zip(ind_add[1:(δ-d)], vals[1:(δ-d)])
        # @show i, index, val
        insert!(x, index, val)
    end

end

# Alternative which just ennacts via the vectors given *no indexing)
function delete_insert!(
    x::Path,
    ind_del::AbstractArray{Int},
    ind_add::AbstractArray{Int},
    vals::AbstractArray{Int})

    @views for (i, index) in enumerate(ind_del)
        deleteat!(x, index - i + 1)
    end
    @views for (index, val) in zip(ind_add, vals)
        # @show i, index, val
        insert!(x, index, val)
    end

end

function imcmc_prop_sample_edit!(
    I_curr::Path{T}, I_prop::Path{T},
    δ::Int, # Number of edits
    d::Int, # Number of deletions (addititions = δ - d)
    vertex_set::Vector{T},
    ind_del::AbstractArray{Int}, ind_add::AbstractArray{Int}, vals::AbstractArray{T} # Storage
    ) where {T<:Union{Int, String}}

    n = length(I_curr)
    m = n + δ - 2*d

    # Sample indexing info and new entries (all in-place)
    @views StatsBase.seqsample_a!(1:n, ind_del[1:d])
    @views StatsBase.seqsample_a!(1:m, ind_add[1:(δ-d)])
    @views sample!(vertex_set, vals[1:(δ-d)])

    delete_insert!(I_prop, δ, d, ind_del, ind_add, vals)

    log_ratio = log(min(n-1, δ)+1) - log(min(m-1, δ)+1) + (δ - 2*d) * log(length(vertex_set))

    return log_ratio
end

ceil(Int, 4/2)

function imcmc_prop_sample_edit_informed(
    I_curr::Path{T},
    δ::Int, # Number of edits
    d::Int, # Number of deletions (addititions = δ - d)
    vertex_dist::DiscreteUnivariateDistribution
    ) where {T<:Union{Int, String}}

    n = length(I_curr)
    m = n + δ - 2*d
    ind_del = StatsBase.seqsample_a!(1:n, zeros(Int, d))
    ind_add = StatsBase.seqsample_a!(1:m, zeros(Int, δ - d))

    I_prop = Vector{T}(undef, m)
    I_prop[Not(ind_add)] = I_curr[Not(ind_del)]
    I_prop[ind_add] = rand(vertex_dist, δ-d)

    log_ratio = log(min(n-1, δ)+1) - log(min(m-1, δ)+1)

    for i in ind_del
        @inbounds log_ratio += logpdf(vertex_dist, I_curr[i])
    end
    for i in ind_add
        @inbounds log_ratio += -logpdf(vertex_dist, I_prop[i])
    end

    return I_prop, log_ratio
end

function accept_reject_edit!(
    I_curr::Path{T}, I_prop::Path{T},
    model::SPF{T},
    ν::Int,
    ind_del::AbstractArray{Int}, ind_add::AbstractArray{Int}, vals::AbstractArray{T}
    ) where {T<:Union{Int,String}}

    n = length(I_curr)
    δ = rand(1:ν) # Num. edits
    d = rand(0:min(n-1,δ)) # Num. deletions
    m = n + δ - 2*d

    # Set-up views
    ind_del = view(ind_del, 1:d)
    ind_add = view(ind_add, 1:(δ-d))
    vals = view(vals, 1:(δ-d))

    # Sample indexing info and new entries (all in-place)
    StatsBase.seqsample_a!(1:n, ind_del)
    StatsBase.seqsample_a!(1:m, ind_add)
    sample!(model.V, vals)

    delete_insert!(I_prop, δ, d, ind_del, ind_add, vals)

    log_ratio = log(min(n-1, δ)+1) - log(min(m-1, δ)+1) + (δ - 2*d) * log(length(model.V))

    log_α = -model.γ * (
        model.dist(I_prop, model.mode)
        - model.dist(I_curr, model.mode)
    ) + log_ratio

    if log(rand()) < log_α
        copy!(I_curr, I_prop)
        return 1
    else
        copy!(I_prop, I_curr)
        return 0
    end
end

function (mcmc::SpfInvolutiveMcmcEdit)(
    model::SPF{T};
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::Path{T}=model.mode
    ) where {T<:Union{Int,String}}

    # Find required length of chain
    req_samples = burn_in + 1 + (desired_samples - 1) * lag
    sample = fill(T[], req_samples)

    I_curr = mcmc.curr
    I_prop = mcmc.prop

    copy!(I_curr, init) # Initialisation
    copy!(I_prop, I_curr)

    # ind_del_mem = zeros(Int, mcmc.ν)
    # ind_add_mem = zeros(Int, mcmc.ν)
    # vals_mem = rand(model.V, mcmc.ν)
    # dist_store1 = zeros(Int, model.K+1)
    # dist_store2 = zeros(Int, model.K+1)

    lb(n::Int, δ::Int) = max(0, ceil(Int, (n + δ - model.K)/2))
    ub(n::Int, δ::Int) = min(n-1, δ)

    count=0

    for i in 1:req_samples

        n = length(I_curr)
        δ = rand(1:mcmc.ν) # Num. edits
        a, b = (lb(n, δ), ub(n, δ))
        # @show δ, a, b
        d = rand(a:b) # Num. deletions
        m = n + δ - 2*d

        # Set-up views
        ind_del = view(mcmc.ind_del, 1:d)
        ind_add = view(mcmc.ind_add, 1:(δ-d))
        vals = view(mcmc.vals, 1:(δ-d))

        # Sample indexing info and new entries (all in-place)
        StatsBase.seqsample_a!(1:n, ind_del)
        StatsBase.seqsample_a!(1:m, ind_add)
        sample!(model.V, vals)

        # Enact delete insertion of entries
        delete_insert!(I_prop, δ, d, ind_del, ind_add, vals)

        # Evaluate log ratio (iMCMC ratio of auxiliary distributions)
        log_ratio = (log(b - a + 1) - log(ub(m,δ) - lb(m,δ) + 1)
            + (δ - 2*d) * log(length(model.V))
            )
        # @show I_prop, log_ratio
        # Evaluate log acceptance probabilility
        log_α = -model.γ * (
            model.dist(I_prop, model.mode)
            - model.dist(I_curr, model.mode)
        ) + log_ratio

        if log(rand()) < log_α
            copy!(I_curr, I_prop)
            count += 1
        else
            copy!(I_prop, I_curr)
        end

        sample[i] = copy(I_curr)
    end
    output = SpfMcmcOutput(
        model,
        sample[(burn_in+1):lag:end],
        count/req_samples
        )
    return output
end

function StatsBase.sample!(
    mcmc::SpfInvolutiveMcmcEdit,
    model::SPF{T},
    output::Vector{Path{T}};
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::Path{T}=model.mode
    ) where {T<:Union{Int,String}}

    # Find required length of chain
    req_samples = burn_in + 1 + (length(output) - 1) * lag
    I_curr = mcmc.curr
    I_prop = mcmc.prop

    copy!(I_curr, init) # Initialisation
    copy!(I_prop, I_curr)

    count=1

    for i in 1:req_samples

        n = length(I_curr)
        δ = rand(1:mcmc.ν) # Num. edits
        d = rand(0:min(n-1,δ)) # Num. deletions
        m = n + δ - 2*d

        # Set-up views
        ind_del = view(mcmc.ind_del, 1:d)
        ind_add = view(mcmc.ind_add, 1:(δ-d))
        vals = view(mcmc.vals, 1:(δ-d))

        # Sample indexing info and new entries (all in-place)
        StatsBase.seqsample_a!(1:n, ind_del)
        StatsBase.seqsample_a!(1:m, ind_add)
        sample!(model.V, vals)

        # Enact delete insertion of entries
        delete_insert!(I_prop, δ, d, ind_del, ind_add, vals)

        # Evaluate log ratio (iMCMC ratio of auxiliary distributions)
        log_ratio = log(min(n-1, δ)+1) - log(min(m-1, δ)+1) + (δ - 2*d) * log(length(model.V))

        # Evaluate log acceptance probabilility
        log_α = -model.γ * (
            model.dist(I_prop, model.mode)
            - model.dist(I_curr, model.mode)
        ) + log_ratio

        if log(rand()) < log_α
            copy!(I_curr, I_prop)
        else
            copy!(I_prop, I_curr)
        end

        if (i > burn_in) & (((i-1) % lag)==0)
            output[count] = copy(I_curr)
            count += 1
        end
    end
end


text:
module Day05

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day05PuzzleInfo <: PuzzleInfo end

"""
Day 5 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day05PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 5 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day05PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Line segment joining two points in two-dimensional spatial domain.

$(FIELDS)
"""
struct LineSegment{T}
    """x-coordinate of first point."""
    x1::T
    """y-coordinate of first point."""
    y1::T
    """x-coordinate of second point."""
    x2::T
    """y-coordinate of second point."""
    y2::T
end

"""
    $(FUNCTIONNAME)(linesegmentstring)

Construct a `LineSegment` object from a string specification `"x1,y1 -> x2,y2"`.
"""
function parselinesegmentstring(linesegmentstring::AbstractString)
    point1string, point2string = split(linesegmentstring, " -> ")
    x1, y1 = [parse(Int, num) for num in split(point1string, ',')]
    x2, y2 = [parse(Int, num) for num in split(point2string, ',')]
    LineSegment{Int}(x1, y1, x2, y2)
end

"""
    $(FUNCTIONNAME)(linefrequencies, min_x, min_y, segment, processdiagonals)

Increment indices in `linefrequencies` representing line segment `segment`.

The value `linefrequencies[i, j]` represents the number of lines covering the point with
coordinates `(min_x + i - 1, min_y + i - 1)`. If `processdiagonals` is `false` (the
default) only horizontal and vertical line segments are recorded, if `true` diagonal
line segments are also recorded.
"""
function addlinesegment!(
    linefrequencies::AbstractArray{S, 2},
    min_x::T,
    min_y::T,
    segment::LineSegment{T},
    processdiagonals::Bool=false
) where {S<:Integer, T<:Integer}
    if segment.x1 == segment.x2 || segment.y1 == segment.y2
        linefrequencies[
            (min(segment.x1, segment.x2) - min_x + 1):(
                max(segment.x1, segment.x2) - min_x + 1),
            (min(segment.y1, segment.y2) - min_y + 1):(
                max(segment.y1, segment.y2) - min_y + 1)
        ] .+= 1
    elseif processdiagonals
        @assert abs(segment.x2 - segment.x1) == abs(segment.y2 - segment.y1)
        for i in 0:abs(segment.x2 - segment.x1)
            linefrequencies[
                segment.x1 + sign(segment.x2 - segment.x1) * i - min_x + 1,
                segment.y1 + sign(segment.y2 - segment.y1) * i - min_x + 1,
            ] += 1
        end
    end
end

function solve(info::Day05PuzzleInfo)
    linesegments = open(info.datapath, "r") do input
        [parselinesegmentstring(segmentstring) for segmentstring in eachline(input)]
    end
    min_x = minimum(l -> min(l.x1, l.x2), linesegments)
    max_x = maximum(l -> max(l.x1, l.x2), linesegments)
    min_y = minimum(l -> min(l.y1, l.y2), linesegments)
    max_y = maximum(l -> max(l.y1, l.y2), linesegments)
    linefrequencies = zeros(Int, max_x - min_x + 1, max_y - min_y + 1)
    for linesegment in linesegments
        addlinesegment!(
            linefrequencies, min_x, min_y, linesegment, isa(info, Part2PuzzleInfo)
        )
    end
    return sum(linefrequencies .>= 2)
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day05.txt"))) == 5
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day05.txt"))) == 12
end

end # module


text:
using Test

ENV["JULIA_PKG_SERVER"] = ""
using PkgServerClient
using PkgServerClient: registry, registry_response_time

function query_upstream(startup_path)
    if ispath(startup_path)
        startup_lines = readlines(startup_path)
    else
        startup_lines = String[]
    end

    pkg_matches = map(x->match(PkgServerClient.regex_PKG_SERVER, x), startup_lines)
    pkg_indices = findall(x->!isnothing(x), pkg_matches)

    if isempty(pkg_indices)
        return ""
    else
        # only last match counts
        return pkg_matches[last(pkg_indices)].captures[2]
    end
end

@testset "PkgServerClient.jl" begin
    @test !isempty(ENV["JULIA_PKG_SERVER"])

    @test all(keys(registry)) do k
        PkgServerClient.set_mirror(k)
        ENV["JULIA_PKG_SERVER"] == registry[k].url
    end

    println("Reigstry response time:")
    foreach(println, registry_response_time())

    # make sure the fallback "JuliaLang" is used when network issue happens
    @test registry_response_time(timeout=0.001)["JuliaLang"] == 0

    @test_throws ArgumentError PkgServerClient.set_mirror("UNKNOWN_SERVER_ABCDE")
    @test_throws ArgumentError PkgServerClient.generate_startup("UNKNOWN_SERVER_ABCDE")

    # check `generate_startup`
    mktempdir() do tmp
        pushfirst!(DEPOT_PATH, tmp)
        default_dot_julia = first(DEPOT_PATH)
        config_path = joinpath(default_dot_julia, "config")
        if !ispath(config_path)
            mkpath(config_path)
        end

        startup_path = joinpath(config_path, "startup.jl")
        @test !isfile(startup_path)
        PkgServerClient.generate_startup()
        @test !isempty(query_upstream(startup_path))
        PkgServerClient.generate_startup("BFSU")
        @test query_upstream(startup_path) == "https://mirrors.bfsu.edu.cn/julia"
    end
end


text:
# convert image to semilogy scaled image
function logimage(A,x)
    @assert size(A,1)==length(x)
    B = copy(A)
    y = logspace(log10(x[1]+0.1),log10(x[end]),length(x))
    for i=1:length(x)
        ind = searchsortedlast(x, y[i])
        ind < 1 && (ind=1)
        B[i,:] = A[ind,:]
    end
    return B
end

# bins of averages on log scale
function logbinavg(p::AbstractArray, f::AbstractArray, bins::Integer)
    @assert length(p) == length(f)
    n = length(p)
    fls = logspace(log10(1), log10(f[end]), bins)
    x = zeros(bins)
    j = 1
    for i = 1:bins
        sum = 0
        c = 0
        while( j <= n && f[j] <= fls[i] )
            sum += p[j]
            c += 1
            j += 1
        end
        if c != 0
            x[i] = sum/c
        else
            x[i] = -200dB
        end
    end

    return x, fls
end

function plotall(s::MCAudio; pgwin=tukey(nsamples(s),0.0003), pgnfft=nextfastfft(nsamples(s)), pgbins::Integer=600)

    fs = s.fs
    n = nsamples(s)
    x = s[:left]

    wl = 256*2
    wl > n && (wl=n>>2)
    gain = 10
    dbr = 80

    # audio
    sv = x
    st = [0:1/fs:(n-1)/fs]

    # periodogram
    P0 = 1
    P = periodogram(x; onesided=true, nfft=pgnfft, fs=fs, window=pgwin)
    Pp, Pf = logbinavg(power(P), freq(P), pgbins)
    Pp = 10*log10(Pp./P0)

    #specrogram
    SG = spectrogram(x, wl, wl>>1; onesided=true, fs=fs, window=hanning)
    SGp = 10*log10(power(SG))
    SGp[find(x->(x<-dbr),SGp)] = -dbr
    SGp[find(x->(x>-gain),SGp)] = -gain
    SGt = time(SG)
    SGf = freq(SG)

    # PLOT

    f, ax = plt.subplots(3, 1, sharex=false)
    ax[1][:plot](st, sv, "b")
    ax[1][:set_xlim]([0,st[end]])
    ax[1][:set_ylim]([-1.1,1.1])
    #ax[1][:set_xlabel]("seconds")

    ax[2][:imshow](flipud(logimage(SGp,SGf)),aspect="auto",interpolation="bilinear",extent=[SGt[1],SGt[end],SGf[1]+0.1,SGf[end]])
    #ax[2][:pcolormesh](repmat(vec(SGt),1,length(SGf)),repmat(vec(SGf),length(SGt),1),flipud(SGp))
    ax[2][:set_ylim]([60,SGf[end]])
    ax[2][:set_yscale]("log")

    ax[2][:set_ylabel]("Hz")
    #ax[2][:set_xlabel]("seconds")

    ax[3][:semilogx](Pf, Pp, "b")
    ax[3][:set_xlim]([10,Pf[end]])
    ax[3][:set_ylim]([-100,2])
    #ax[3][:set_xlabel]("Hz")
    ax[3][:set_ylabel]("dB")
    ax[3][:grid](b=true, axis="x", which="minor")
    ax[3][:grid](b=true, axis="y")
    #f[:set_size_inches](10,7)
    draw()
    return f, ax

end


text:
function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.tRb] = 5.000e+00
    p[C.Skp2] = 1.000e+00
    p[C.kSyE2f] = 3.000e-02
    p[C.kSyE2fE2f] = 4.000e-02
    p[C.jSyE2f] = 2.000e-01
    p[C.kAsRbE2f] = 5.000e+00
    p[C.kDsRbE2f] = 5.000e-03
    p[C.kDeE2f] = 5.000e-02
    p[C.kPhRbCd] = 2.000e-01
    p[C.kPhRbCe] = 3.000e-01
    p[C.kDpRb] = 5.000e-02
    p[C.kDeP21] = 2.500e-03
    p[C.kDeP21Cy] = 7.000e-03
    p[C.kSyCe] = 1.000e-02
    p[C.kAsCyP21] = 1.000e+00
    p[C.kDsCyP21] = 5.000e-02
    p[C.kDeCe] = 4.000e-03
    p[C.tP21] = 1.500e+00

    p[C.mitogen] = 1.000e+01
    p[C.DUSP] = 1.000e+00
    p[C.pAkt] = 1.000e+00
    p[C.kfRSK] = 1.735e+00
    p[C.krRSK] = 9.034e-01
    p[C.kSyCd] = 4.500e-3
    p[C.kDeCd] = 9.000e-3
    p[C.kon] = 1.000e-02
    p[C.koff] = 2.281e-03
    p[C.kf] = 4.650e-02
    p[C.kr] = 9.564e-02
    p[C.tERK] = 7.236e+01
    p[C.KmRSK] = 9.936e+01
    p[C.kfFOS] = 6.149e-03
    p[C.krFOS] = 1.475e-02
    p[C.KmFOS] = 1.123e-01

    return p
end


text:
include("dependencies_for_runtests.jl")

@testset "Utils" begin
    @info "Testing utils..."

    @testset "prettytime" begin
        @test prettytime(0) == "0 seconds"
        @test prettytime(35e-15) == "3.500e-14 seconds"

        @test prettytime(1e-9) == "1 ns"
        @test prettytime(1e-6) == "1 μs"
        @test prettytime(1e-3) == "1 ms"

        @test prettytime(second) == "1 second"
        @test prettytime(minute) == "1 minute"
        @test prettytime(hour) == "1 hour"
        @test prettytime(day) == "1 day"
        @test prettytime(year) == "1 year"

        @test prettytime(2second) == "2 seconds"
        @test prettytime(4minute) == "4 minutes"
        @test prettytime(8hour) == "8 hours"
        @test prettytime(16day) == "16 days"
        @test prettytime(32year) == "32 years"

        @test prettytime(13.7seconds) == "13.700 seconds"
        @test prettytime(6.666minutes) == "6.666 minutes"
        @test prettytime(1.234hour) == "1.234 hours"
        @test prettytime(40.5days) == "40.500 days"
        @test prettytime(5.0001years) == "5.000 years"
    end
end


text:
using Mosek

"""
    basicPortfolio(n::Int, x0 :: Vector{Float64})

Implements the simple Portfolio Selection example maximizing return subject to a limited risk.

    maximize
        mu' * x
    such that
        sum(x) = budget
        GT * x - s = 0
        s >= || t ||
        s <= gamma

Parameters:

- `w`     Initial uninvested wealth
- `x0`    Initial invested wealth for each asset
- `mu`    Expected returns for each asset
- `gamma` Risk limit as maximum standard deviation
- `GT`    Factorized covariance matrix
"""
function basicPortfolio(w :: Float64, x0 :: Vector{Float64}, mu :: Vector{Float64}, GT :: Array{Float64,2}, gamma :: Float64)
    n = length(mu)
    maketask() do t
        putstreamfunc(t,MSK_STREAM_LOG,msg -> print(msg))

        budget = w + sum(x0)

        # variables
        appendvars(t,2*n + 1)
        x_start  = 1
        x_end    = n
        s_start  = n+1
        s_end    = n+1
        t_start  = n+2
        t_end    = t_start+n-1

        putobjsense(t,MSK_OBJECTIVE_SENSE_MAXIMIZE)
        putclist(t,Int32[Int32(x_start):Int32(x_end) ...], mu)

        putvarboundslice(t,Int32(x_start), Int32(x_end+1),
                         fill(MSK_BK_LO,n),
                         zeros(Float64,n),
                         fill(Inf,n))

        putvarboundslice(t,Int32(t_start), Int32(t_end+1),
                         fill(MSK_BK_LO,n),
                         zeros(Float64,n),
                         fill(Inf,n))

        for j in 1:n
            putvarname(t,Int32(x_start+j-1), "x$j")
            putvarname(t,Int32(t_start+j-1), "t$j")
        end
        putvarname(t,Int32(s_start),"s")

        # constraints
        appendcons(t,Int32(1+n))

        # sum(x) = w+sum(x0)
        putconbound(t,Int32(1),MSK_BK_FX,w + sum(x0),w + sum(x0))
        putconname(t,Int32(1),"budget")
        putarow(t,Int32(1),Int32[Int32(x_start):Int32(x_end)...],ones(Float64,n))

        # GT * x - s = 0
        putconboundslice(t,2,n+2,
                         fill(MSK_BK_FX, n),
                         zeros(Float64,n),
                         zeros(Float64,n))

        for i in 1:n
            putconname(t,Int32(i+1),"GT[$i]")
            putarow(t,Int32(i+1),Int32[Int32(x_start):Int32(x_end)...],GT[i,:])
            putaij(t,Int32(i+1),s_start+i-1, -1.0)
        end

        # s > || t ||
        appendcone(t,MSK_CT_QUAD, 0.0, [ s_start, t_start:t_end... ])
        putconename(t,Int32(1),"stddev")

        # s = gamma
        putvarbound(t,Int32(s_start),MSK_BK_FX,gamma,gamma)

        # Turn all log output off.
        putintparam(t,MSK_IPAR_LOG, Int32(0))

        optimize(t,"mosek://solve.mosek.com:30080")

        # Display the solution summary for quick inspection of results.
        solutionsummary(t, MSK_STREAM_MSG)

        xxvector = getxx(t,MSK_SOL_ITR)

        #expectedReturn = sum(xxvector[x_start:x_end] .* mu)
        #standardDeviation = xxvector[s_start]

        xxvector[x_start:x_end],xxvector[s_start]
    end
end

# default data
gamma = 0.05
mu = [0.1073, 0.0737, 0.0627]
GT = [0.1667 0.0232  0.0013
      0.0000 0.1033 -0.0022
      0.0000 0.0000  0.0338]
x0 = [0.0, 0.0, 0.0]
w = 1.0

x,s = basicPortfolio(w,x0,mu,GT,gamma)

expectedReturn = sum(x .* mu)
println("------------------------------------------------------------")
println("Expected return $expectedReturn for gamma $s")
println("Optimal investment: $x")
println("------------------------------------------------------------")


text:
module Summary

import DivPopGa.ClusteredGa: Chromosome
import Statistics

function summary(population::Array{Chromosome, 1})::Dict{String, Any}
    costs = map(c -> c.cost, population)
    sorted_pop     = sort(population, by = ch -> ch.cost)
    mean_cost      = Statistics.mean(costs)
    median_cost    = Statistics.median(costs)
    min_cost       = minimum(costs)
    max_cost       = maximum(costs)
    q25_cost       = Statistics.quantile(costs, 0.25)
    q75_cost       = Statistics.quantile(costs, 0.75)
    best           = sorted_pop[1]

    result = Dict(
        "mean_cost"       => mean_cost,
        "min_cost"        => min_cost,
        "q25_cost"        => q25_cost,
        "median_cost"     => median_cost,
        "q75_cost"        => q75_cost,
        "max_cost"        => max_cost,
        "best_solution"   => best
    )

    return result
end

end # End of module

max_stars_repo_name: akadal/divpopga

max_issues_repo_name: akadal/divpopga

max_forks_repo_name: akadal/divpopga


text:
module Ingenuous

using LinearAlgebra

include("clustering.jl")
include("projection.jl")
include("quality.jl")
include("models.jl")

end # module

max_stars_repo_name: richardtjornhammar/ingenuous

max_issues_repo_name: richardtjornhammar/ingenuous

max_forks_repo_name: richardtjornhammar/ingenuous


text:
"""
    longdisplay(x,rows=1000,mimetype="text/plain")
    ld(args...)

Display `x` with `rows` overriding the number of rows of the display.
This can be used to display all the rows of a `DataFrame` instance,
nicely formatted, at the REPL. It does not work for arrays.
"""
longdisplay(x,rows=1000,mimetype="text/plain") = longdisplay(STDOUT,x,rows,mimetype)

"""
    longdisplay(io::IO, x,rows=1000,mimetype="text/plain")

Send output to `io` instead of `STDOUT`.
"""
function longdisplay(io::IO, x, rows=1000, mimetype="text/plain")
    dispsize = displaysize()
    newdispsize = (rows,dispsize[2])
    show(IOContext(io, displaysize=newdispsize), mimetype, x)
end

const ld = longdisplay

"""
    @withw f(filename,x,y,...)

Open file `filename` for writing and call `f(io,x,y,...)` where
`io` is the output stream. The stream is closed when the call to `f`
exits.

# Example
julia> @withw println("outfile.txt", "hello")
"""
macro withw(fcall)
    (fn,fcall.args[2]) = (fcall.args[2],:io)
    quote
        open($fn,"w") do io
            $(fcall)
        end
    end
end


text:
####################################################

#  RRT.jl

#  Implementing PRM and functions for plotting PRM
#  Parameters incude num samples, connection distance, a list of obstacles, and start and goal states
#  Graph search implemented is Astar. Edge cost is assumed to be equal to euclidean distance.
#  To use, call ` include("PRM.jl") ` from any Julia file in the same folder.
#  nouyang 2017

####################################################

using Plots
using DataStructures
using GeometryTypes
using Distributions

module algT
    using GeometryTypes
    # export GraphNode, Edge, Obstacle, Room

    struct GraphNode
        id::Int64
        state::Point{2, Float64}
    end

    struct Edge
        startID::Int64
        endID::Int64
    end

    struct Obstacle
        id::Int64
        rect::HyperRectangle{2, Vec}
    end

    struct Room
        width::Int64
        height::Int64
        walls::Vector{LineSegment}
        obstacles::Vector{HyperRectangle}
    end

    struct AlgParameters
        numSamples::Int64
        connectRadius::Int64
    end

    struct queueTmp
        node::algT.GraphNode
        statesList::Vector{algT.GraphNode}
        cost::Int64
    end

    struct Graph
        startstate::Point{2,Float64}
        goalstate::Point{2,Float64}
        nodeslist::Vector{GraphNode}
        edgeslist::Vector{Edge}
    end

    Base.isless(q1::queueTmp, q2::queueTmp) = q1.cost < q2.cost

end

module algfxn
    using GeometryTypes
    using Distributions
    using algT

    function sampleFree(roomW, roomH, obstacles)
        pt_rand = Void

        sampledFree = false
        while !sampledFree
            xrand = rand(Uniform(0, roomW))
            yrand = rand(Uniform(0, roomH))
            pt = Point(xrand, yrand)
            if !algfxn.isCollidingNode(pt, obstacles) # this does NOT check if point is on a line (e.g. a wall)
                pt_rand = Point(xrand, yrand)
                sampledFree = true
            end
        end

        return pt_rand
    end

    function steer(fromPt, toPt, connectRadius)
        dist = min_euclidean(Vec(fromPt), Vec(toPt))
        # If pt_rand is too far from its nearest node, "truncate" it to be closer
        if dist <= connectRadius
            return toPt
        else
            x1,y1 = fromPt
            x2,y2 = toPt

            theta = atan2((y2-y1) , (x2-x1))
            newX = x1 + connectRadius * cos(theta)
            newY = y1 + connectRadius * sin(theta)

            newPt = Point(newX, newY)
            return newPt
        end
    end

	function ccw(A,B,C)
        # determines direction of lines formed by three points
        return (C[2]-A[2]) * (B[1]-A[1]) > (B[2]-A[2]) * (C[1]-A[1])
    end

    function intersects(line1, line2) #no ":" at the end!
        A, B = line1[1], line1[2]
        C, D = line2[1], line2[2]
        return ( (ccw(A, C, D) != ccw(B, C, D)) && ccw(A, B, C) != ccw(A, B, D))
    end

    function decompRect(r::HyperRectangle) #GeometryTypes.HyperRectangle{2,Float64}
        corners = decompose(Point{2, Float64}, r)
        corners = [Point(pt) for pt in corners]
        lineBottom  = LineSegment(corners[1], corners[2])
        lineTop     = LineSegment(corners[3], corners[4])
        lineLeft    = LineSegment(corners[1], corners[3])
        lineRight   = LineSegment(corners[2], corners[4])
        lines = Vector{LineSegment}()
        push!(lines, lineTop, lineRight, lineBottom, lineLeft)
        return lines
    end

    #function isCollidingNode(node::algT.GraphNode, obsList::Vector{algT.Obstacle}
    function isCollidingNode(node::Point{2, Float64}, obsList::Vector{HyperRectangle})
        for obs in obsList
            if contains(obs, node)
                return true
            end
        end
        return false
    end

    function isCollidingObstacles(line::LineSegment, obsList::Vector{HyperRectangle})
        for obs in obsList
            rectLines = decompRect(obs)
            for rectline in rectLines
                if intersects(line, rectline)[1]
                    return true
                end
            end
        end
        return false
    end

    function isCollidingWalls(line::LineSegment, walls::Vector{LineSegment})
        for wall in walls
            if intersects(line, wall)
                return true
            end
        end
        return false
    end

    function findNearestNodes(nodestate, nodeslist, maxDist)
        # given maxDist, return all nodes within that distance of node
        # list is NOT sorted by distance
        nearestNodes = Vector{Tuple{algT.GraphNode, Float64}}()
        for n in nodeslist
            dist = min_euclidean(Vec(nodestate), Vec(n.state))
            if dist < maxDist
                push!(nearestNodes, (n, dist))
            end
        end
        return nearestNodes
    end

    function nearestN(pt::Point{2,Float64}, nodeslist)
        #
        # This loops through all nodes
        nearestDist = 99999999;
        nearestNode = Void;

        for n in nodeslist
            dist = min_euclidean(Vec(pt), Vec(n.state))
            if dist < nearestDist
                nearestDist = dist
                nearestNode = n
            end
        end
        return nearestNode
    end

    function inGoalRegion(pt, pt_goal)
        if pt == pt_goal
            return true
        else
            dist = min_euclidean(Vec(pt), Vec(pt_goal))
            if dist <= 3 #HARDCODED REGION #TODO
                return true
            end
        end

        return false
    end

    function costPath(solPath)
        pathcost = 0
        for i in 2:length(solPath)
            curN  = solPath[i].state
            prevN = solPath[i-1].state
            pathcost += min_euclidean(Vec(curN), Vec(prevN))
        end
        return pathcost
    end

    function findNode(nodeID, nodeslist)
        # Assumes we have not removed any nodes, hence we can sort by ID
        index = findfirst([node.id for node in nodeslist], nodeID)
        return nodeslist[index]
    end

end

module plotfxn

    using GeometryTypes
    using Plots
    using algfxn
    using algT

    ###  Plots.jl recipes

    @recipe function f(r::HyperRectangle)
        points = decompose(Point{2,Float64}, r)
        rectpoints = points[[1,2,4,3],:]
        xs = [pt[1] for pt in rectpoints];
        ys = [pt[2] for pt in rectpoints];
        seriestype := :shape
        color := :orange
        linecolor := :black
        #m = (:black, stroke(0))
        s = Shape(xs[:], ys[:])
    end

    @recipe function f(pt::Point)
        xs = [pt[1]]
        ys = [pt[2]]
        seriestype --> :scatter
        color := :orange
        markersize := 3
        xs, ys
    end

    @recipe function f(l::LineSegment)
        xs = [ l[1][1], l[2][1] ]
        ys = [ l[1][2], l[2][2] ]
        # seriestype = :line
        color := :darkblue
        lw := 3
        xs, ys
    end

    @recipe function f(rectList::Vector{<:HyperRectangle})
        for r in rectList
            @series begin
                r
            end
        end
    end

    @recipe function f(ptList::Vector{<:Point})
        for p in ptList
            @series begin
                p
            end
        end
    end

    @recipe function f(lineList::Vector{<:LineSegment})
        for l in lineList
            @series begin
                l
            end
        end

    end

    function plotRoom(room)
        aPlot = plot() #Todo! this assumes plotroom is first thing called()
        roomWidth, roomHeight, walls, obstacles = room.width, room.height, room.walls, room.obstacles
        plot!(aPlot, walls, color =:black)
        plot!(aPlot, obstacles, fillalpha=0.5)
        return aPlot
    end

    function plotRRT(roomPlot, graph, solPath, title)
        startstate, goalstate, nodeslist, edgeslist = graph.startstate, graph.goalstate, graph.nodeslist, graph.edgeslist

        x = [n.state[1] for n in nodeslist]
        y = [n.state[2] for n in nodeslist]
        scatter!(roomPlot, x,y, color=:black)

        edgeXs, edgeYs = [], []
        for e in edgeslist
            startN = algfxn.findNode(e.startID, nodeslist)
            endN = algfxn.findNode(e.endID, nodeslist)
            x1,y1 = startN.state[1], startN.state[2]
            x2,y2 = endN.state[1], endN.state[2]
            push!(edgeXs, x1, x2, NaN) #the NaNs, keep spaces between edges correctly unplotted
            push!(edgeYs, y1, y2, NaN)
        end

        plot!(roomPlot, edgeXs, edgeYs, color=:tan, linewidth=1)

        rrtPlot = plotSolPath(roomPlot, solPath)

        title!(rrtPlot, title)
        plot!(rrtPlot, legend=false, size=(600,600), xaxis=((-2,22), 0:1:20 ),
              yaxis=((-2,22), 0:1:20), foreground_color_grid= :black, titlefont=font("Arial", 10))

        return rrtPlot
    end

    function plotSolPath(aPlot, solPath)
        #print("\n ---Solution Path----- \n")
        if solPath != Void
            xPath = [n.state[1] for n in solPath]
            yPath = [n.state[2] for n in solPath]
            xstart, ystart = solPath[1].state

            xend, yend = solPath[end].state
            scatter!(aPlot, [xstart], [ystart],
                     markercolor= :red, markershape = :circle,  markersize = 6, markerstrokealpha = 0.5, markerstrokewidth=1)
            scatter!(aPlot, [xend], [yend],
                     markerstrokecolor = :green, markershape = :dtriangle,  markersize = 5, markerstrokealpha = 1, markerstrokewidth=5)
            plot!(aPlot, xPath, yPath, color = :orchid, linewidth=3, fillalpha = 0.3)
        else
           # #print("\n --- No solution path found ----- \n")
        end

        return aPlot
    end

end

# pts are pure Points, nodes have ID information
function rrtPlan(room, parameters, startstate, goalstate, obstaclesList)
    roomWidth, roomHeight, walls, obstacles = room.width, room.height, room.walls, room.obstacles
    numPts, connectRadius = parameters.numSamples, parameters.connectRadius

    nodeslist = Vector{algT.GraphNode}()
    edgeslist = Vector{algT.Edge}()

    startNode = algT.GraphNode(0, startstate)
    push!(nodeslist, startNode)

    currID = 1 #current node ID

    isPathFound = false

    n_new = algT.GraphNode(999999, goalstate)
    for i in 1:numPts
        pt_rand = algfxn.sampleFree(roomWidth, roomHeight, obstacles)

        n_nearest = algfxn.nearestN(pt_rand, nodeslist)
        pt_new = algfxn.steer(n_nearest.state, pt_rand, connectRadius) #no node ID yet
        candidateEdge = LineSegment(n_nearest.state, pt_new)

        if (    !algfxn.isCollidingObstacles(candidateEdge, obstacles) && !algfxn.isCollidingWalls(candidateEdge, walls) )
            n_new = algT.GraphNode(currID, pt_new)
            currID += 1
            newEdge = algT.Edge(n_nearest.id , n_new.id)
            push!(nodeslist, n_new)
            push!(edgeslist, newEdge)

            if algfxn.inGoalRegion(pt_new, goalstate)
                #print("\n\n - n_new $n_new , pt_new $pt_new -\n\n")
                isPathFound = true
                break
            end
        end
    end

    print("\n\n - pathfound $isPathFound - \n\n")

    pathNodes = Vector{algT.GraphNode}()

    goalNode = algT.GraphNode(999999, goalstate)
    push!(pathNodes, goalNode)

    if isPathFound
        n_last = n_new
        push!(pathNodes, n_last)

        currPathID = n_last.id

        while currPathID != 0
            foundParent = false

            node = Void
            for edge in edgeslist
                if edge.endID == currPathID
                    node = algfxn.findNode(edge.startID, nodeslist)
                    push!(pathNodes, node)
                    break
                end
            end
            currPathID = node.id
        end
    else
        # if no feasible path was found, just return the start and end nodes so we can plot them
    end

    push!(pathNodes, startNode)

    finalPathCost = algfxn.costPath(pathNodes)
    return (nodeslist, edgeslist, isPathFound, pathNodes, finalPathCost)

end