text:
"""
    barplot(text, heights; kw...)

# Description

Draws a horizontal barplot. It uses the first parameter (`text`)
to denote the names for the bars, and the second parameter
(`heights`) as their values. This means that the two vectors have
to have the same length. Alternatively, one can specify a barplot
using a dictionary `dict`. In that case, the keys will be used as
the names and the values, which have to be numeric, will be used
as the heights of the bars.

# Usage

    barplot(text, heights; $(keywords((border = :barplot, color = :green,), remove = (:xlim, :ylim, :xscale, :yscale, :height, :grid), add = (:symbols,))))

    barplot(dict; kw...)

# Arguments

$(arguments(
    (
        text = "the labels / captions of the bars",
        heights = "the values / heights of the bars",
        dict = "a dictonary in which the keys will be used as `text` and the values will be used as `heights`",
        xscale = "`Function` or `Symbol` to transform the bar length before plotting: this effectively scales the `x`-axis without influencing the captions of the individual bars (use `xscale = :log10` for logscale)",
    ); remove = (:xlim, :ylim, :xscale, :yscale, :height, :grid), add = (:symbols,),
))

# Returns

A plot object of type `Plot{BarplotGraphics}`.

# Author(s)

- Christof Stocker (github.com/Evizero)

# Examples

```julia-repl
julia> barplot(["Paris", "New York", "Moskau", "Madrid"],
               [2.244, 8.406, 11.92, 3.165],
               xlabel = "population [in mil]")
            ┌                                        ┐
      Paris ┤■■■■■■ 2.244
   New York ┤■■■■■■■■■■■■■■■■■■■■■■■ 8.406
     Moskau ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 11.92
     Madrid ┤■■■■■■■■■ 3.165
            └                                        ┘
                       population [in mil]
```

# See also

[`Plot`](@ref), [`histogram`](@ref), [`BarplotGraphics`](@ref)
"""
function barplot(
    text::AbstractVector{<:AbstractString},
    heights::AbstractVector{<:Number};
    border = :barplot,
    color = :green,
    out_stream::Union{Nothing,IO} = nothing,
    width::Int = out_stream_width(out_stream),
    symb = nothing,  # deprecated
    symbols = KEYWORDS.symbols,
    xscale = KEYWORDS.xscale,
    xlabel = transform_name(xscale),
    kw...,
)
    length(text) == length(heights) ||
        throw(DimensionMismatch("The given vectors must be of the same length"))
    minimum(heights) >= 0 || throw(
        ArgumentError("All values have to be positive. Negative bars are not supported."),
    )

    if any('\n' in t for t in text)
        _text = eltype(text)[]
        _heights = eltype(heights)[]
        for (t, h) in zip(text, heights)
            lines = split(t, '\n')
            if (n = length(lines)) > 1
                append!(_text, lines)
                for i in 1:n
                    push!(_heights, i == n ? h : -1)
                end
            else
                push!(_text, t)
                push!(_heights, h)
            end
        end
        text, heights = _text, _heights
    end

    area = BarplotGraphics(
        heights,
        width,
        xscale;
        color = color,
        symbols = _handle_deprecated_symb(symb, symbols),
    )
    plot = Plot(area; border = border, xlabel = xlabel, kw...)
    for i in 1:length(text)
        label!(plot, :l, i, text[i])
    end
    plot
end

barplot(dict::Dict{T,N}; kw...) where {T,N<:Number} =
    barplot(sorted_keys_values(dict)...; kw...)
barplot(label, height::Number; kw...) = barplot([string(label)], [height]; kw...)

function barplot(text::AbstractVector, heights::AbstractVector{<:Number}; kw...)
    text_str = map(string, text)
    barplot(text_str, heights; kw...)
end

"""
    barplot!(plot, text, heights)

Mutating variant of `barplot`, in which the first parameter
(`plot`) specifies the existing plot to draw additional bars on.

See [`barplot`](@ref) for more information.
"""
function barplot!(
    plot::Plot{<:BarplotGraphics},
    text::AbstractVector{<:AbstractString},
    heights::AbstractVector{<:Number},
)
    length(text) == length(heights) ||
        throw(DimensionMismatch("The given vectors must be of the same length"))
    isempty(text) && throw(ArgumentError("Can't append empty array to barplot"))
    curidx = nrows(plot.graphics)
    addrow!(plot.graphics, heights)
    for i in 1:length(heights)
        label!(plot, :l, curidx + i, text[i])
    end
    plot
end

barplot!(plot::Plot{<:BarplotGraphics}, dict::Dict{T,N}) where {T,N<:Number} =
    barplot!(plot, sorted_keys_values(dict)...)

function barplot!(
    plot::Plot{<:BarplotGraphics},
    text::AbstractVector,
    heights::AbstractVector{<:Number},
)
    text_str = map(string, text)
    barplot!(plot, text_str, heights)
end

function barplot!(plot::Plot{<:BarplotGraphics}, label, heights::Number)
    curidx = nrows(plot.graphics)
    addrow!(plot.graphics, heights)
    label!(plot, :l, curidx + 1, string(label))
    plot
end


text:
__precompile__()
module Remove_eq

export detect_eq_kurtosis, stalta, remove_eq, s_whiten, s_whiten!

using SeisIO, SeisNoise, FFTW, DSP, StatsBase, Plots, JLD2, SeisIO, Printf

"""
    detect_eq_kurtosis(data::SeisChannel,tw::Float64=60.0, threshold::Float64=3.0, overlap::Float64=30)

find earthquake by kurtosis threshold

# Input:
    - `data::SeisChannel`    : SeisData from SeisIO
    - `tw::Float64`  : time window to evaluate if earthquake is contained.
    - `threshold::Float64` : kurtosis threshold: if kurtosis > threshold, the time window contains earthquake
    - `overlap::Float64`           : overlap of time window to control the margin of earthquake removal. (large overlap assigns large margin before and after earthquake.)

    kurtosis evaluation following Baillard et al.(2013)
"""
function detect_eq_kurtosis(data::SeisChannel; tw::Float64=60.0, kurtosis_threshold::Float64=3.0, overlap::Float64=30)

    #convert window lengths from seconds to samples
    twsize = trunc(Int, tw * data.fs)
    overlapsize = trunc(Int, overlap * data.fs)

    #calculate how much to move beginning of window each iteration
    slide = twsize-overlapsize

    #kurtosis of timeseries
    ku1 = data.misc["kurtosis"][:]
    #reset long window counter and triggers for current channel
    i = 0

    #loop through current channel by sliding
    while i < length(ku1) - twsize

        #check if last window and change long window length if so
        if length(ku1) - i < twsize
            twsize = length(ku1)-i
        end

        #define chunk of data based on long window length and calculate long-term average
        twtrace = @views ku1[i+1:i+twsize]

        if !isnothing(findfirst(x -> x > kurtosis_threshold, twtrace))
            #this time window includes earthquake
            for tt= i+1:i+twsize
                data.misc["eqtimewindow"][tt] = false
            end
        end

        #advance long window
        i = i + slide

    end

    return data

end

"""
    detect_eq_stalta(data::SeisChannel,tw::Float64=60.0, threshold::Float64=3.0, overlap::Float64=30)

find earthquake and tremors by STA/LTA

# Input:
    - `data::SeisChannel`  : SeisChannel from SeisIO
    - `longwin::Float64`   : long time window
    - `shortwin::Float64`  : short time window
    - `threshold::Float64` : STA/LTA threshold: if STA/LTA > threshold, the time window contains earthquake
    - `overlap::Float64`   : overlap of time window to control the margin of earthquake removal. (large overlap assigns large margin before and after earthquake.)
    - `stalta_absoluteclip::Float64`   : clip if maximum absolute value exceeds this number (for the purpose of removing incoherent noise)

    original code written by Seth Olinger. For our purpose, overlap is applied for short time window
"""
function detect_eq_stalta(data::SeisChannel,longWinLength::Float64, shortWinLength::Float64, threshold::Float64, overlap::Float64, stalta_absoluteclip::Float64,
                            InputDict::Dict, tstamp::String, st1::String; datagap_eps::Float64=1e-8)

    #convert window lengths from seconds to samples
    longWin = trunc(Int,longWinLength * data.fs)
    shortWin = trunc(Int,shortWinLength * data.fs)
    overlapWin = trunc(Int,overlap * data.fs)
    trace = @view data.x[:]

    #buffer for sta/lta value
    staltatrace = zeros(length(trace))

    #calculate how much to move beginning of window each iteration
    slide = longWin

    #save weight
    eqweight = deepcopy(data.misc["eqtimewindow"])
    #print("before manipulation"); println(count(eqweight))

    # manipulate weight to avoid underestimation of lta due to data gap
    # threshold zero signal by 0.001% of lta within the timeseries, which encompasses small computational error such as 1e-21
    ltaall = StatsBase.mean(abs.(trace))
    for i = 1:length(trace)
        if isless(abs(trace[i]), ltaall*datagap_eps)
            # this is regarded as zero signal (data gap)
            eqweight[i] = false
        end
    end

    #print("after manipulation"); println(count(eqweight))

    #define long window counter and empty trigger vector
    i = 1

    #loop through current channel by sliding
    while i < length(trace)

        #check if last window and change long window length if so
        if length(trace) - i < longWin
            longWin = length(trace)-i
        end

        #define chunk of data based on long window length and calculate long-term average
        longTrace = trace[i:i+longWin]

        lta = StatsBase.mean(abs.(longTrace), weights(eqweight[i:i+longWin]))

        #reset short window counter
        n = 0
        breakflag = false
        #loop through long window in short windows
        while n <= longWin - shortWin

            #define chunk of data based on short window length and calculate short-term average
            shortTrace = @view trace[i+n:i+n+shortWin]
            shortTraceWeight = @view eqweight[i+n:i+n+shortWin]
            #sta = mean(abs.(shortTrace))
            sta = StatsBase.mean(abs.(shortTrace), weights(shortTraceWeight))
            #sta = StatsBase.mean(abs.(shortTrace))
            stamax = maximum(abs.(shortTrace))
            #calculate sta/lta ration
            staLta = sta/lta

            #record detection time if sta/lta ratio exceeds threshold
            if staLta > threshold || stamax > stalta_absoluteclip
                #triggers[i+n] = 1
                #this time window includes earthquake
                for tt= i+n:i+n+shortWin
                    data.misc["eqtimewindow"][tt] = false
                end
            end

            if breakflag
                break;
            end

            #advance short window
            n = n + shortWin - overlapWin

            #adjust for the last edge of long window
            if n >= longWin - shortWin
                n = longWin - shortWin
                breakflag = true
            end

            staltatrace[i+n+shortWin] = staLta

        end

        #advance long window
        i = i + slide

    end

    if InputDict["dumptraces"]
        #dump sta/lta trace
        mkpath(InputDict["dumppath"])
        tstamp_fname = replace(tstamp, ":" => "-")
        fname_out = join([tstamp_fname, st1, "stalta","dat"], '.')
        fopath = InputDict["dumppath"]*"/"*fname_out
        open(fopath, "w") do io
            for i in 1:length(staltatrace)
                write(io, @sprintf("%12.8f\n", staltatrace[i]))
            end
        end
    end
    return data

end

"""
    remove_eq(data::SeisChannel)

remove earthquake by kurtosis and STA/LTA threshold

# Input:
    - `data::SeisData`    : SeisData from SeisIO

"""
function remove_eq(data::SeisChannel, data_origin::SeisChannel, plot_kurtosis_α::Float64, max_taper_dur::Float64,
    plot_boxheight::Float64, plot_span::Int64, plot_fmt::String, fodir::String, tstamp::String, tvec::Array{Float64,1}, IsSaveFig::Bool)

    eqidlist = data.misc["eqtimewindow"][:]
    nx = length(data.x)

    i = 1

    t1 = []
    t2 = []
    y1 = []
    y2 = []

    tt1 = 0
    tt2 = 0
    tt3 = 0
    tt4 = 0
    tt5 = 0

    while i <= nx
        if !eqidlist[i]
            push!(t1, tvec[i])

            t1id = i

            #find next id
            tt1 += @elapsed nexttrueid = findfirst(x -> x == true, eqidlist[t1id:end])

            if isnothing(nexttrueid)
                # all data is removed
                t2id = nx
                iinc = nx
            else
                t2id = (t1id - 1) + nexttrueid - 1 # first false id to end false id
                iinc = t2id - i + 1 # first true as end of false + 1
            end

            push!(t2, tvec[t2id])

            max_wintaper_duration = Int(data.fs*max_taper_dur)

            # compute α given maximum taper length
            eq_length = t2id-t1id+1
            taper_length = 2*max_wintaper_duration
            tukey_length = eq_length + taper_length
            invert_tukey_α = taper_length/tukey_length

            # define inverted tukey window
            tt2 += @elapsed invtukeywin = -tukey(Int(tukey_length), invert_tukey_α) .+ 1

            # slice tukey window if it exceeds array bounds
            tt3 += @elapsed if t1id < max_wintaper_duration + 1
                left_overflow = (max_wintaper_duration-t1id)+1
                invtukeywin = @views invtukeywin[left_overflow+1:end]
                # leftmost t
                left = 1
            else
                # full taper length
                left = t1id - max_wintaper_duration
            end

            tt4 += @elapsed if (nx - t2id + 1) < max_wintaper_duration + 1
                right_overflow = (max_wintaper_duration-(nx - t2id + 1))+1
                invtukeywin = @views invtukeywin[1:end-right_overflow]
                # rightmost t
                right = nx
            else
                # full taper length
                right = t2id + max_wintaper_duration
            end

            # apply tukey window
            tt5 += @elapsed data.x[left:right] .*= invtukeywin

            #boxsize
            push!(y1, -plot_boxheight)
            push!(y2, plot_boxheight)

        else
            iinc = 1
        end

        i += iinc

    end

    #println([tt1, tt2, tt3, tt4, tt5])

    if IsSaveFig

        figdir = joinpath(fodir, "fig_removalEQ")
        if !ispath(figdir); mkpath(figdir); end

        Plots.plot(bg=:white)
        normalized_amp = 0.5 * maximum(data_origin.x[1:plot_span:end])

        Plots.plot!(tvec[1:plot_span:end], data_origin.x[1:plot_span:end]./ normalized_amp,
                    label="raw data", color="black")

        Plots.plot!(tvec[1:plot_span:end], data.x[1:plot_span:end]./ normalized_amp,
                    label="after removal", color="blue")

        #to plot kurtosis computed points
        kurtid = findall(x -> !iszero(data.misc["kurtosis"][x]), 1:length(data.misc["kurtosis"]))
        if isempty(filter(!isnan, data.misc["kurtosis"][kurtid]))
            kurtnormalize = 1.0
        else
            kurtnormalize = maximum(filter(!isnan, data.misc["kurtosis"][kurtid]))
        end

        Plots.plot!(tvec[kurtid], data.misc["kurtosis"][kurtid] ./ kurtnormalize .* plot_kurtosis_α,
                    color="red",
                    label="kurtosis")

        p = Plots.plot!(xlabel = "[hours]",
                    ylabel = "[m/s]",
                    title =  @sprintf("%s %s", data.id, tstamp),
                    size = (1200, 800))

        figname = @sprintf("%s/%s_%s.%s", figdir, data.id, tstamp, plot_fmt)
        Plots.savefig(p, figname)

        # save jld2 fig file
        # trace1 = scatter(;x=tvec[1:plot_span:end], y=data_origin.x[1:plot_span:end]./ normalized_amp,
        #  mode="lines", name="raw data", line_color="black")
        #
        # trace2 = scatter(;x=tvec[1:plot_span:end], y=data.x[1:plot_span:end]./ normalized_amp,
        #  mode="lines", name="after remove", line_color="blue")
        #
        # #to plot kurtosis computed points
        # kurtid = findall(x -> !iszero(data.misc["kurtosis"][x]), 1:length(data.misc["kurtosis"]))
        #
        # if isempty(filter(!isnan, data.misc["kurtosis"][kurtid]))
        #     kurtnormalize = 1.0
        # else
        #     kurtnormalize = maximum(filter(!isnan, data.misc["kurtosis"][kurtid]))
        # end
        #
        # trace3 = PlotlyJS.scatter(;x=tvec[kurtid], y=data.misc["kurtosis"][kurtid] ./ kurtnormalize .* plot_kurtosis_α,
        #  line_color="red", mode="lines", name="kurtosis")
        #
        # if !isempty(t1)
        #  shapes = PlotlyJS.rect(t1, t2, y1, y2; fillcolor="#ff99cc", opacity=0.3, line_width=0)
        #  layout = PlotlyJS.Layout(shapes=shapes, width=1200, height=600,
        #      xaxis=attr(title="Time [hour]"),
        #      yaxis=attr(title="Normalized velocity"),
        #      font =attr(size=13),
        #      showlegend=true,
        #      title = @sprintf("%s %s", data.id, tstamp))
        #
        #      p = PlotlyJS.plot([trace1; trace2; trace3],layout)
        #
        # else
        #  layout = PlotlyJS.Layout(width=1200, height=600,
        #      xaxis=attr(title="Time [hour]"),
        #      yaxis=attr(title="Normalized velocity"),
        #      font =attr(size=12),
        #      showlegend=true,
        #      title = @sprintf("%s %s", data.id, tstamp))
        #
        #  p = PlotlyJS.plot([trace1; trace2; trace3],layout)
        #
        # end
        #
        # figdir = joinpath(fodir, "fig")
        # mkpath(figdir)
        # figname = @sprintf("%s/%s_%s.%s", figdir, data.id, tstamp, plot_fmt)
        # ORCA.savefig(p, figname)
        # #display(p)
        # #println("press return for next plot...")
        # #readline()
    end

    return data

end

"""
    s_whiten!(data::SeisChannel)

Apply spectral whitening to seischannel trace
# Input:
    - `data::SeisChannel`    : SeisChannel from SeisIO
    - `freqmin::Real`: Pass band low corner frequency.
    - `freqmax::Real`: Pass band high corner frequency.
    - `pad::Int`: Number of tapering points outside whitening band.
"""
function s_whiten!(data::SeisChannel, freqmin::Float64, freqmax::Float64;pad::Int=50)

    # compute fft of time series
    FFT = rfft(data.x,1)

    # to use SeisNoise.whiten!() prepare (N, 2) array
    FFTtemp = Array{Complex{Float32}}(undef, length(FFT), 2)
    FFTtemp[:, 1] = FFT

    SeisNoise.whiten!(FFTtemp,freqmin,freqmax,data.fs, data.t[end, 1], pad=pad)

    data.x = irfft(FFTtemp[:,1],data.t[end, 1],1)
    return nothing

end
s_whiten(data::SeisChannel, freqmin::Float64, freqmax::Float64;pad::Int=50) = (U = deepcopy(data);
    s_whiten!(U,freqmin,freqmax,pad=pad);
    return U)

end

ext: jl

max_stars_repo_name: kura-okubo/SeisRemoveEQ

max_issues_repo_name: kura-okubo/SeisRemoveEQ

max_forks_repo_name: kura-okubo/SeisRemoveEQ


text:
using Test

include("./Tools.jl")

#evalpoly testing
coefs = [1.0,7.0,9.0]
x = 3.0
@test evalpoly(x,coefs) ≈ sum([x^(i-1)*coefs[i] for i = 1:length(coefs)])

#getIndexFirstGreater testing
a = [-4.0,1.0,3.0,7.0]
@test getBoundingIndsSorted(2.0,a) == (2,3)
@test getBoundingIndsSorted(-5.0,a) == (1,)
@test getBoundingIndsSorted(8.0,a) == (length(a),)

@test true


text:
# Autogenerated wrapper script for Gumbo_jll for aarch64-linux-musl
export libgumbo

JLLWrappers.@generate_wrapper_header("Gumbo")
JLLWrappers.@declare_library_product(libgumbo, "libgumbo.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgumbo,
        "lib/libgumbo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
a, b, c = Point(1, 2), Point(3, 7), Point(2, 4)
@show a b c
@show a + b
@show -a + b
@show a - b
@show a + b + c
@show a == b
@show a + a == c

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
for op = [:add, :sub, :max, :min]
    @eval function $(Symbol("scatter_", op, "!"))(ys::Matrix{T}, us::Array{T}, xs::Array{Int}) where T
        l = length(xs)
        s = size(ys, 1)

        Threads.@threads for num ∈ 1:(l*s)
            li = (num -1) ÷ s + 1
            i = (num - 1) % s + 1
            @inbounds ind = Tuple(CartesianIndices(xs)[li])
            @inbounds $(Symbol("atomic_", op, "!"))(
                pointer(ys,
                        Base._to_linear_index(ys,
                                              i,
                                              xs[li]
                                              )
                        ),
                us[i, ind...]
            )
        end

        return ys
    end

    @eval function $(Symbol("scatter_", op, "!"))(ys::Array{T}, us::Array{T}, xs::Array{<:Tuple}) where T
        l = length(xs)
        s = size(ys, 1)

        Threads.@threads for num = 1:(l*s)
            li = (num -1) ÷ s + 1
            i = (num - 1) % s + 1
            @inbounds ind = Tuple(CartesianIndices(xs)[li])
            @inbounds $(Symbol("atomic_", op, "!"))(
                pointer(ys,
                        Base._to_linear_index(ys,
                                              i,
                                              xs[li]...
                                              )
                        ),
                us[i, ind...]
            )
        end

        return ys
    end
end


text:
using Curiosity
using Curiosity: DPI, HordeDPI
using Curiosity: OpenWorld
using MinimalRLCore
const TDGWU = Curiosity.TwoDGridWorldUtils
include("../experiments/2d-gridworld.jl")
using Random
using StatsBase
using GVFHordes
using Statistics
using JLD2

using ProgressMeter

StatsBase.sample(p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(Weights([p.func(;state_t = s, action_t = a) for a in actions]))

StatsBase.sample(rng, p::GVFHordes.GVFParamFuncs.AbstractPolicy, s, actions) =
    sample(Weights([get(p;state_t = s, action_t = a) for a in actions]))

function gen_dpi(num_states=500, seed=1; start_type=:center, h_kwargs...)
    Random.seed!(seed)
    parsed = TwoDGridWorldExperiment.default_args()
    parsed["cumulant_schedule"] = "Constant"
    parsed["cumulant"] = 1.0
    parsed["demon_learner"] = "Q"

    dpis = DPI[]

    state_filter_func =  Curiosity.FeatureSubset(identity, 1:2)

    for gvf_i ∈ 1:4
        horde = Horde(
        [GVF(GVFParamFuncs.FeatureCumulant(gvf_i+2),
             TDGWU.GoalTermination(0.95),
             TDGWU.GoalPolicy(gvf_i, true))])

        cumulant_schedule = TDGWU.get_cumulant_schedule(parsed)

        env = OpenWorld(10, 10,
                        cumulant_schedule=cumulant_schedule,
                        start_type=start_type, normalized=true)

        ss = Vector{Float64}[]
        as = Int[]
        policy = TDGWU.NaiveGoalPolicy(gvf_i)

        while length(ss) < num_states
            s = MinimalRLCore.start!(env)
            push!(ss, state_filter_func(s))

            a = policy(s)
            term = false
            while term == false
            # for t ∈ 1:rand(1:15)
                s, r, term = MinimalRLCore.step!(env, a)
                a = policy(s)
                push!(ss, state_filter_func(s))
                if term
                    # println("termed")
                    break
                end
            end
        end

         push!(dpis, DPI(ss, policy; h_kwargs...))
    end

    return HordeDPI(dpis, state_filter_func)
end

function save_dpi(num_states, seed; start_type=:center, h_kwargs...)
    hdpi = gen_dpi(num_states, seed; start_type=start_type, h_kwargs...)
    @save "./src/data/dpi/2dOpenWorld_$(start_type).jld2" hdpi
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
mutable struct VectorWithEntryInserted{T,A} <: AbstractVector{T}
  a::A
  index::Int
  value::T
  function VectorWithEntryInserted(a::AbstractVector,index::Integer,value)
    A = typeof(a)
    T = eltype(a)
    @assert 1 <= index <= length(a)+1
    new{T,A}(a,index,value)
  end
end

Base.IndexStyle(::Type{<:VectorWithEntryInserted}) = IndexLinear()

function Base.getindex(v::VectorWithEntryInserted,i::Integer)
  i < v.index ? v.a[i] :
                (i==v.index ? v.value : v.a[i-1])
end

Base.size(v::VectorWithEntryInserted) = (length(v.a)+1,)

function array_cache(v::VectorWithEntryInserted)
  array_cache(v.a)
end

function getindex!(cache,v::VectorWithEntryInserted,i::Integer)
  i < v.index ? getindex!(cache,v.a,i) :
                (i==v.index ? v.value : getindex!(cache,v.a,i-1))
end

function Base.sum(a::VectorWithEntryInserted)
  sum(a.a) + a.value
end

function Base.setindex!(a::VectorWithEntryInserted,v,i::Integer)
   i < a.index ? (a.a[i] = v) :
                 (i==a.index ? a.value = v : a.a[i-1]=v)
end


text:
using BayesMultiFit
using Test
BAMF=BayesMultiFit
using ReversibleJumpMCMC
const RJMCMC = ReversibleJumpMCMC
using Random
using Distributions
using MicroscopePSFs
PSF=MicroscopePSFs

include("testhelpers.jl")

@testset "BayesMultiFit.jl" begin

    #TODO: update tests to use examples and test overall functionality
    datastate, psf, xystd,istd,split_std,bndpixels,prior_photons= gendatastate(Int32(1))

    # test genmodel! for DD type data for all BAMFData types
    DDdatatypelist=[BAMF.ArrayDD]
    expectDD=BAMF.ArrayDD(Int32(32))
    rjsDD = BAMF.RJStruct(32,psf,xystd,istd,split_std,expectDD,bndpixels,prior_photons)
    BAMF.genmodel!(datastate, rjsDD, expectDD)
    for datatype in DDdatatypelist
        data=genBAMFDD(datatype)
        BAMF.genmodel!(datastate, rjsDD, data)
        sz, dim, DD=DDinfo(data)
        @test sz == 32
        if length(dim)==3
            @test dim==(32, 32, 1)
            for i in 1:32, j in 1:32
                @test data.data[i, j, 1] ≈ expectDD.data[i, j] atol=1f-4
            end
        else
            @test dim==(32,32)
            #@test isapprox(data.data , expected.DD.data; atol=1f-4)
            #for i in 1:32, j in 1:32
            #    @test data.data[i, j] ≈ expectDD.data[i, j] atol=1f-4
            #end
        end
        @test isa(DD, datatype)

    end
end


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

export LAPACKException,
       ARPACKException,
       SingularException,
       PosDefException,
       RankDeficientException

struct LAPACKException <: Exception
    info::BlasInt
end

struct ARPACKException <: Exception
    info::String
end

function ARPACKException(i::Integer)
    if i == -8
        return ARPACKException("error return from calculation of a real Schur form.")
    elseif i == -9
        return ARPACKException("error return from calculation of eigenvectors.")
    elseif i == -14
        return ARPACKException("did not find any eigenvalues to sufficient accuracy. Try with a different starting vector or more Lanczos vectors by increasing the value of ncv.")
    end
    return ARPACKException("unspecified ARPACK error: $i")
end

struct SingularException <: Exception
    info::BlasInt
end

struct PosDefException <: Exception
    info::BlasInt
end
function Base.showerror(io::IO, ex::PosDefException)
    print(io, "PosDefException: matrix is not ")
    if ex.info == -1
        print(io, "Hermitian")
    else
        print(io, "positive definite")
    end
    print(io, "; Cholesky factorization failed.")
end

struct RankDeficientException <: Exception
    info::BlasInt
end

max_stars_repo_name: nextjournal/julia-lang

max_issues_repo_name: nextjournal/julia-lang

max_forks_repo_name: nextjournal/julia-lang


text:
struct Domain{FT<:Real}
    α :: FT
    β :: FT
    g :: FT
    H :: FT
end

function DomainProperties(
	FT = Float64;
	α :: Real = 0.2,
	β :: Real = 0.5,
	g :: Real = 1,
	H :: Real = 1
)
	return Domain{FT}(α,β,g,H)

end

function show(io::IO, D::Domain{FT}) where FT <: Real
    print(
		io,
		"The Domain{$FT} properties are as follows:\n",
		"     Rayleigh Damping (α) : ", D.α, '\n',
		"    Coriolis Gradient (β) : ", D.β, '\n',
		"              Gravity (g) : ", D.g, '\n',
		"       Surface Height (H) : ", D.H, '\n'
	)
end


text:
module TestLastitems

using DataTools
using Test
using Transducers

include("utils.jl")

@testset "lastitem" begin
    @test lastitem(3:7) === 7
    @test 3:7 |> Map(x -> x + 1) |> Filter(isodd) |> lastitem == 7
end

@testset "lastitems" begin
    @test lastitems(3:7, 2) ==ₜ view(3:7, 4:5)
    @test 3:7 |> Map(x -> x + 1) |> Filter(isodd) |> lastitems(2) == [5, 7]
end

end  # module


text:
__precompile__()
module CopyMethods

export copy_col!, copy_row!, unit

"""
    function copy_col!{T}(A::AbstractArray{T,2},col1::Int,col2::Int,row_min=0,row_max=Inf)

In matrix `A` copy column `col1` over `col2`. Optionally constrain the range from `row_min` to `row_max`.
__Warning:__ Bound checks eliminated for maximizing performance.
"""
function copy_col!{T}(A::AbstractArray{T,2},col1::Int,col2::Int,row_min=0,row_max=Inf)
  @inbounds @simd for i in max(1, row_min):min(size(A,1), row_max)
    A[i,col2] = A[i,col1]
  end
end

"""
    function copy_row!{T}(A::AbstractArray{T,2},row1::Int,row2::Int,col_min=0,col_max=Inf)

In matrix `A` copy row `row1` over `row2`. Optionally constrain the range from `col_min` to `col_max`.
__Warning:__ Bound checks eliminated for maximizing performance.
"""
function copy_row!{T}(A::AbstractArray{T,2},row1::Int,row2::Int,col_min=0,col_max=Inf)
  for i in max(1, col_min):min(size(A,2), col_max)
    @inbounds A[row2,i] = A[row1,i]
  end
end

function copy_row!{T}(A::AbstractArray{T,1},N::Int, row1::Int,row2::Int)
  if length(A) == 0 || max(row1, row2) > N
      return
  end
  # @printf("%s %s %s %s %s",row1, row2, N, length(A), length(A)/N)
  @inbounds @simd for i in 1:Int(length(A)/N)
    A[ row2+N*(i-1) ] = A[ row1+N*(i-1) ]
  end
end

unit(N,j) = map(x->(if x==j return 1 else 0 end), 1:N)

end


text: Dict{String,Any}("formfeed" => Dict{String,Any}("value" => "This string has a \f form feed character.","type" => "string"),"notunicode2" => Dict{String,Any}("value" => "This string does not have a unicode \\u escape.","type" => "string"),"backslash" => Dict{String,Any}("value" => "This string has a \\ backslash character.","type" => "string"),"notunicode3" => Dict{String,Any}("value" => "This string does not have a unicode \\u0075 escape.","type" => "string"),"notunicode4" => Dict{String,Any}("value" => "This string does not have a unicode \\u escape.","type" => "string"),"tab" => Dict{String,Any}("value" => "This string has a \t tab character.","type" => "string"),"carriage" => Dict{String,Any}("value" => "This string has a \r carriage return character.","type" => "string"),"quote" => Dict{String,Any}("value" => "This string has a \" quote character.","type" => "string"),"newline" => Dict{String,Any}("value" => "This string has a \n new line character.","type" => "string"),"notunicode1" => Dict{String,Any}("value" => "This string does not have a unicode \\u escape.","type" => "string"),"backspace" => Dict{String,Any}("value" => "This string has a \b backspace character.","type" => "string"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using BinaryBuilder

name = "LibSpatialIndex"
version = v"1.8.5"

# Collection of sources required to build LibSpatialIndex
sources = [
    ArchiveSource("http://download.osgeo.org/libspatialindex/spatialindex-src-1.8.5.tar.bz2",
        "31ec0a9305c3bd6b4ad60a5261cba5402366dd7d1969a8846099717778e9a50a"),
    DirectorySource("./patches"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir

cd spatialindex-src-*

patch < ${WORKSPACE}/srcdir/makefile.patch
rm Makefile.am.orig

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  patch < ${WORKSPACE}/srcdir/header-check.patch
fi

aclocal
autoconf
automake --add-missing --foreign

# Show options in the log
./configure --help

./configure --prefix=${prefix} --host=$target --build=${MACHTYPE} --enable-static=no
make
make install
install_license COPYING
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libspatialindex_c", :libspatialindex_c),
    LibraryProduct("libspatialindex", :libspatialindex),
]

# Dependencies that must be installed before this package can be built
dependencies = []

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
export Scalingtrain1st!

using Flux
using Flux: onehotbatch
using Base.Iterators
using StatsBase
using Zygote
using Zygote: @nograd
using LinearAlgebra
using CUDA
using Plots

function collectscales(model)
    for (i,cell) in enumerate(m.cells)
        for (j,mixedop) in enumerate(cell.mixedops)
            for (k,op) in enumerate(mixedop.ops)
                for (l,layer) in enumerate(op.op)
                    if typeof(layer) <: Flux.BatchNorm
                        @show (i,j,k,l,layer.λ)
                    end
                end
            end
        end
    end
end

connectstrings = vcat([[string(j,"-",i) for j = 1:i-1] for i = 3:6]...)
axlab = vcat([string("n",e) for e in connectstrings],[string("r",e) for e in connectstrings])

edgedict = Dict{String, Int64}()
for (i,strconn) in enumerate(connectstrings)
    edgedict[strconn] = i
end
primdict = Dict{String, Int64}()
for (i,prim) in enumerate(PRIMITIVES)
    primdict[prim] = i
end

function scalingupdate(model)
    reducecells = Array{Int64, 1}()
    normalcells = Array{Int64, 1}()
    scales = Array{Float32, 3}(undef, length(model.cells), length(model.cells[1].mixedops), length(model.cells[1].mixedops[1].ops)) |> gpu#cell id x mixop id x op id
    for (i,cell) in enumerate(model.cells)
        if cell.reduction
            push!(reducecells, i)
        else
            push!(normalcells, i)
        end
        for (j,mixedop) in enumerate(cell.mixedops)
            for (k,op) in enumerate(mixedop.ops)
                scales[i,j,k] = 0.0
                for (l,layer) in enumerate(op.op)
                    if typeof(layer) <: Flux.BatchNorm
                        #@show(k, mean(layer.γ .- 1.0))
                        scales[i,j,k] = mean(layer.γ .- 1.0)  #automatically gets only last BatchNorm in op
                    end
                end
            end
        end
    end
    reduce = scales[reducecells,:,:]
    normal = scales[normalcells,:,:]
    display(vcat(dropdims(mean(normal,dims=1),dims=1),dropdims(mean(reduce,dims=1),dims=1)))
    vcat(dropdims(mean(normal,dims=1),dims=1),dropdims(mean(reduce,dims=1),dims=1))'
end

function scalingreg(model)
    reducecells = Array{Int64, 1}()
    normalcells = Array{Int64, 1}()
    collectscales = 0.0
    for (i,cell) in enumerate(model.cells)
        for (j,mixedop) in enumerate(cell.mixedops)
            for (k,op) in enumerate(mixedop.ops)
                for (l,layer) in enumerate(op.op)
                    if typeof(layer) <: Flux.BatchNorm
                        collectscales += abs(mean(layer.γ))
                    end
                end
            end
        end
    end
    collectscales
end

i = 1
function Scalingtrain1st!(loss, model, train, val, opt_α, opt_w, λ=0.0001, losses=[0.0,0.0]; cbepoch = () -> (), cbbatch = () -> ())
    local train_loss
    local val_loss
    local gsα
    w = all_ws(model)
    α = all_αs(model)
    fakeg = all_αs(model).*0
    fakea = all_αs(model).*0
    global i
    for (train_batch, val_batch) in zip(TrainCuIterator(train), TrainCuIterator(val))
        gsw = gradient(w) do
            train_loss = loss(model, train_batch...) + λ*scalingreg(model)#sum([sum(abs, ss) for ss in scalingupdate(model)])
            return train_loss
        end
        losses[1] = train_loss
        foreach(CUDA.unsafe_free!, train_batch)
        Flux.Optimise.update!(opt_w, w, gsw)
        CUDA.reclaim()

        gsα = gradient(α) do
            val_loss = loss(model, val_batch...)
            return val_loss
        end
        losses[2] = val_loss
        foreach(CUDA.unsafe_free!, val_batch)
        CUDA.reclaim()
        GC.gc()
        gs = [gsα[x] for x in α]
        for (f,g) in zip(fakeg, gs)
            f .-= g
        end
        CUDA.reclaim()
        GC.gc()
        cbbatch()
    end
    p1 = heatmap(axlab, PRIMITIVES, scalingupdate(model)|>cpu, xrotation = 90, title = "Mean Scaling Delta")
    #ac2 = activationupdatema(model)
    #p2 = heatmap(axlab, PRIMITIVES, hcat((ac2 |> cpu)...), xrotation = 90, title = "Activation Mean of Absolute Value")
    #gs = [gsα[x] for x in α]
    p3 = heatmap(axlab, PRIMITIVES, hcat((fakeg |> cpu)...), xrotation = 90 , title = "Alpha Gradient Update")
    plot(p1,p3,layout = (2,1),size = (1200, 800))
    savefig(string("test/models/fig_scale", i, ".png"));

    i += 1
    cbepoch()
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
module IOCTL

if Sys.islinux()
    # Most generic definitions for ioctl for Linux
    # https://github.com/torvalds/linux/blob/b6839ef26e549de68c10359d45163b0cfb031183/tools/include/uapi/asm-generic/ioctls.h

    const TCGETS = 0x5401
    const TCSETS = 0x5402
    const TCSETSW = 0x5403
    const TCSETSF = 0x5404
    const TCGETA = 0x5405
    const TCSETA = 0x5406
    const TCSETAW = 0x5407
    const TCSETAF = 0x5408
    const TCSBRK = 0x5409
    const TCXONC = 0x540A
    const TCFLSH = 0x540B
    const TIOCEXCL = 0x540C
    const TIOCNXCL = 0x540D
    const TIOCSCTTY = 0x540E
    const TIOCGPGRP = 0x540F
    const TIOCSPGRP = 0x5410
    const TIOCOUTQ = 0x5411
    const TIOCSTI = 0x5412
    const TIOCGWINSZ = 0x5413
    const TIOCSWINSZ = 0x5414
    const TIOCMGET = 0x5415
    const TIOCMBIS = 0x5416
    const TIOCMBIC = 0x5417
    const TIOCMSET = 0x5418
    const TIOCGSOFTCAR = 0x5419
    const TIOCSSOFTCAR = 0x541A
    const FIONREAD = 0x541B
    const TIOCINQ = FIONREAD
    const TIOCLINUX = 0x541C
    const TIOCCONS = 0x541D
    const TIOCGSERIAL = 0x541E
    const TIOCSSERIAL = 0x541F
    const TIOCPKT = 0x5420
    const FIONBIO = 0x5421
    const TIOCNOTTY = 0x5422
    const TIOCSETD = 0x5423
    const TIOCGETD = 0x5424
    const TCSBRKP = 0x5425 # /* Needed for POSIX tcsendbreak() */
    const TIOCSBRK = 0x5427   # /* BSD compatibility */
    const TIOCCBRK = 0x5428   # /* BSD compatibility */
    const TIOCGSID = 0x5429   # /* Return the session ID of FD */
    const TIOCGRS485 = 0x542E
    const TIOCSRS485 = 0x542F
    const TIOCGPTN = 0 # TODO  # /* Get Pty Number (of pty-mux device) */
    const TIOCSPTLCK = 0 # TODO   # /* Lock/unlock Pty */
    const TIOCGDEV = 0 # TODO  # /* Get primary device node of /dev/console */
    const TCGETX = 0x5432  # /* SYS5 TCGETX compatibility */
    const TCSETX = 0x5433
    const TCSETXF = 0x5434
    const TCSETXW = 0x5435
    const TIOCSIG = 0 # TODO   # /* pty: generate signal */
    const TIOCVHANGUP = 0x5437
    const TIOCGPKT = 0 # TODO  # /* Get packet mode state */
    const TIOCGPTLCK = 0 # TODO  # /* Get Pty lock state */
    const TIOCGEXCL = 0 # TODO  # /* Get exclusive mode state */
    const TIOCGPTPEER = 0 # TODO  # /* Safely open the slave */
    const TIOCGISO7816 = 0 # TODO
    const TIOCSISO7816 = 0 # TODO

    const FIONCLEX = 0x5450
    const FIOCLEX = 0x5451
    const FIOASYNC = 0x5452
    const TIOCSERCONFIG = 0x5453
    const TIOCSERGWILD = 0x5454
    const TIOCSERSWILD = 0x5455
    const TIOCGLCKTRMIOS = 0x5456
    const TIOCSLCKTRMIOS = 0x5457
    const TIOCSERGSTRUCT = 0x5458  # /* For debugging only */
    const TIOCSERGETLSR =  0x5459  # /* Get line status register */
    const TIOCSERGETMULTI = 0x545A  # /* Get multiport config  */
    const TIOCSERSETMULTI = 0x545B  # /* Set multiport config */

    const TIOCMIWAIT = 0x545C # /* wait for a change on serial input line(s) */
    const TIOCGICOUNT = 0x545D # /* read serial port inline interrupt counts */
    const TIOCPKT_DATA =  0
    const TIOCPKT_FLUSHREAD =  1
    const TIOCPKT_FLUSHWRITE =  2
    const TIOCPKT_STOP =  4
    const TIOCPKT_START =  8
    const TIOCPKT_NOSTOP = 16
    const TIOCPKT_DOSTOP = 32
    const TIOCPKT_IOCTL = 64

    const TIOCSER_TEMT = 0x01 # /* Transmitter physically empty */

else

    # MACOS
    const TIOCGETP = Sys.islinux() ? 0 : 1074164744 # Get parameters -- V6/V7 gtty()
    const TIOCSETP = Sys.islinux() ? 1 : -2147060727 # Set parameters -- V6/V7 stty()

    const TIOCSETN = 2 # V7:   as above, but no flushtty
    const TIOCEXCL = 3 # V7: set exclusive use of tty
    const TIOCNXCL = 4 # V7: reset excl. use of tty
    const TIOCHPCL = 5 # V7: hang up on last close
    const TIOCFLUSH = 6 # V7: flush buffers

    const TIOCSTI = 7 # simulate terminal input
    const TIOCSBRK = 8 # set   break bit
    const TIOCCBRK = 9 # clear break bit
    const TIOCSDTR = 10 # set   data terminal ready
    const TIOCCDTR = 11 # clear data terminal ready
    const TIOCGPGRP = 12 # get pgrp of tty
    const TIOCSPGRP = 13 # set pgrp of tty

    const TCGETP = 1076130901
    const TCGETA = 1075082331
    const TCSETAW = -2146143143

    const FIONREAD = 22 # get # bytes to read

    const TIOCGETD = 23 # Get line discipline
    const TIOCSETD = 24 # Set line discipline

    const TIOCGWINSZ = 1074295912 # Get window size info
    const TIOCSWINSZ = 26 # Set window size info (maybe gen SIGWINCH)

end

# TODO: check that these values make sense

const TIOCGETC = Sys.islinux() ? 14 : 1074164754 # get special characters
const TIOCSETC = Sys.islinux() ? 15 : -2147060719 # set special characters
const TIOCLBIS = Sys.islinux() ? 16 : 0x8004747f # set   bits in local mode word
const TIOCLBIC = Sys.islinux() ? 17 : 0x8004747e # clear bits in local mode word
const TIOCLGET = Sys.islinux() ? 18 : 0x8004747c # get local mode mask
const TIOCLSET = Sys.islinux() ? 19 : 0x8004747d # set local mode mask
const TIOCSLTC = Sys.islinux() ? 20 : -2147060619 # set local special chars
const TIOCGLTC = Sys.islinux() ? 21 : 1074164852 # get local special chars

const ALLDELAY =0177400 # Delay algorithm selection
const BSDELAY = 0100000 # Select backspace delays
const BS0 = 0
const BS1 = 0100000
const VTDELAY = 040000 # for,-feed/v-tab delay
const FF0 = 0
const FF1 = 040000
const CRDELAY = 030000 # carriage-return delay
const CR0 = 0
const CR1 = 010000
const CR2 = 020000
const CR3 = 030000
const TBDELAY = 06000 # tab delays
const TAB0 = 0
const TAB1 = 01000
const TAB2 = 04000
const XTABS = 06000
const NLDELAY = 01400 # new-line delays
const NL0 = 0
const NL1 = 0400
const NL2 = 01000
const NL3 = 01400
const EVENP = 0200 # even parity allowed on input
const ODDP = 0100 # odd parity allowed on input
const RAW = 040 # wake on all chars, 8-bit input
const CRMOD = 020 # map CR->LF; echo LF or CR as CRLF
const ECHO = 010 # echo (full duplex)
const LCASE = 04 # map upper case to lower case
const CBREAK = 02 # return each char as soon as typed
const TANDEM = 01 # automatic flow control

const LCRTBS = 01 # Backspace on erase rather than echoing erase
const LPRTERA = 02 # Printing terminal erase mode
const LCRTERA = 04 # Erase char echoes as BS-SP-BS
const LTILDE = 010 # Convert ~ to ` on output (for Hazeltines)
const LMDMBUF = 020 # Stop/start output when carrier drops
const LLITOUT = 040 # Suppress output translations
const LTOSTOP = 0100 # Send SIGTTOU for background output
const LFLUSHO = 0200 # Output is being flushed
const LNOHANG = 0400 # Don't send hangup when carrier drops
const LETXACK = 01000 # Diablo style buffer hacking (??)
const LCRTKIL = 02000 # Use BS-SP-BS to erase entire line on line kill
const LINTRUP = 04000 # Generate SIGTINT when input ready to read
const LCTLECH = 010000 # Echo input control chars as ^X (DEL as ^?)
const LPENDIN = 020000 # Retype pending input at next read or input char
const LDECCTQ = 040000 # Only ^Q restarts after ^S, like DEC systems

const OTTYDISC = 0 # Old V7-style discipline (must be zero)
const NTTYDISC = 1 # New BSD-style discipline
const NETLDISC = 2 # high-speed "net" discipline (not supported)

struct winsize
    ws_row::UInt16
    ws_col::UInt16
    ws_xpixel::UInt16
    ws_ypixel::UInt16
end

struct tchars
    t_intrc::UInt8 # interrupt
    t_quitc::UInt8 # quit
    t_startc::UInt8 # start output
    t_stopc::UInt8 # stop output
    t_eofc::UInt8 # end-of-file
    t_brkc::UInt8 # input delimiter (like nl)
end

struct ltchars
    t_suspc::UInt8 # /* stop process signal */
    t_dsuspc::UInt8 # /* delayed stop process signal */
    t_rprntc::UInt8 # /* reprint line */
    t_flushc::UInt8 # /* flush output (toggles) */
    t_werasc::UInt8 # /* word erase */
    t_lnextc::UInt8 # /* literal next character */
end

_file_handle(s::Base.LibuvStream) = ccall(:jl_uv_file_handle, Base.OS_HANDLE, (Ptr{Cvoid},), s.handle)

function ioctl(fd::RawFD, parameter::Int)
    ws = Ref{winsize}()
    r = ccall(:ioctl, Cint, (Cint, Cint, Ptr{Cvoid}), fd, parameter, ws)
    r == -1 ? error("ioctl failed: $(Base.Libc.strerror())") : nothing
    return ws[]
end
ioctl(s::Base.LibuvStream, parameter) = ioctl(_file_handle(s), Int(parameter))
ioctl(f::Int, parameter) = ioctl(RawFD(f), parameter)

end


text:
# single qubit Cliffords

export localclifford

const C1 = Vector{Clifford{1}}(undef, 24)

# identity
C1[1]  = RI

# pi/2 and pi rotations about X,Y,Z
C1[2]  = exp(-im*pi/4*X)
C1[3]  = im * X
C1[4]  = exp(im*pi/4*X)
C1[5]  = exp(-im*pi/4*Y)
C1[6]  = im * Y
C1[7]  = exp(im*pi/4*Y)
C1[8]  = exp(-im*pi/4*Z)
C1[9]  = im * Z
C1[10] = exp(im*pi/4*Z)

# Hadamard class
C1[11] = exp(-im*pi/2/sqrt(2) * (X+Y))
C1[12] = exp(-im*pi/2/sqrt(2) * (X-Y))
C1[13] = exp(-im*pi/2/sqrt(2) * (X+Z)) # standard Hadamard
C1[14] = exp(-im*pi/2/sqrt(2) * (X-Z))
C1[15] = exp(-im*pi/2/sqrt(2) * (Y+Z))
C1[16] = exp(-im*pi/2/sqrt(2) * (Y-Z))

# Axis exchange class
C1[17] = exp(-1im*pi/3/sqrt(3) * (X+Y+Z))
C1[18] = exp(-2im*pi/3/sqrt(3) * (X+Y+Z))
C1[19] = exp(-1im*pi/3/sqrt(3) * (X-Y+Z))
C1[20] = exp(-2im*pi/3/sqrt(3) * (X-Y+Z))
C1[21] = exp(-1im*pi/3/sqrt(3) * (X+Y-Z))
C1[22] = exp(-2im*pi/3/sqrt(3) * (X+Y-Z))
C1[23] = exp(-1im*pi/3/sqrt(3) * (-X+Y+Z))
C1[24] = exp(-2im*pi/3/sqrt(3) * (-X+Y+Z))

const rC1 = Dict{Clifford{1},UInt}()

for (idx,v) in enumerate(C1)
    rC1[v] = idx
end

localclifford(i::Int) = C1[i]
localclifford(v::Vector) = kron(C1[v]...)

localcliffordindex(c::Clifford) = rC1[c]


text:
module RocketDelayOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: delay()" begin

    println("Testing: operator delay()")

    run_proxyshowcheck("Delay", delay(0), args = (check_subscription = true, ))

    run_testset([
        (
            source = of(2) |> delay(50),
            values = @ts(50 ~ [ 2 ] ~ c)
        ),
        (
            source = completed(Int) |> delay(50),
            values = @ts(50 ~ c)
        ),
        (
            source = faulted("e") |> delay(50),
            values = @ts(50 ~ e("e"))
        ),
        (
            source = never() |> delay(50),
            values = @ts()
        )
    ])

end

end


text:
"""
    IDNotFoundInBackbone

This exception is thrown when the ID is not found in the current version of the backbone.
"""
struct IDNotFoundInBackbone <: Exception
    id::Integer
end

function Base.showerror(io::IO, e::IDNotFoundInBackbone)
    message = "The TAX ID $(e.id) is not in the current version of the backbone."
    print(io, message)
end

"""
    NameHasNoDirectMatch

This exception is thrown when the name passed as an argument does not have a
direct match, in which case using `strict=false` to switch to fuzzy matching
may be advised.
"""
struct NameHasNoDirectMatch <: Exception
    name::String
end

function Base.showerror(io::IO, e::NameHasNoDirectMatch)
    message = "The name $(e.name) has no direct match\n"
    message *= "Possible fixes are\n"
    message *= "→ do a lowercase search with `lowercase=true`\n"
    message *= "→ do a fuzzy search with `strict=false`\n"
    message *= "→ look at the output of `similarnames(\"$(e.name)\")`"
    print(io, message)
end

"""
    NameHasMultipleMatches

This exception is thrown when the name is an "in-part" name, which is not a
valid node but an aggregation of multiple nodes. It is also thrown when the
name is valid for several nodes. The error message will return the taxa that
could be used instead. "Reptilia" is an example of a node that will throw
this exception (in-part name); "Mus" will throw this example as it is valid
subgenus of itself.
"""
struct NameHasMultipleMatches <: Exception
    name::String
    taxa::Vector{NCBITaxonomy.NCBITaxon}
end

function Base.showerror(io::IO, e::NameHasMultipleMatches)
    message = "The name $(e.name) is ambiguous.\n"
    message *= "In the NCBI taxonomy, it refers to\n"
    for taxa in e.taxa
        message *= "→ the $(rank(taxa)) $(taxa)\n"
    end
    message *= "Please pick the correct taxa\n"
    message *= "You can get a list of correct taxa using `alternativetaxa(\"$(e.name)\")`\n"
    print(io, message)
end


text:
@inline function tostr(val::Transition)
    hval = Cell(val >> BSHIFT_CELL)
    lval = Cell(val & BMASK_C)
    string(tostr(hval),tostr(lval))
end

@inline function tostr(val::Cell)
    string((Char((val >> BSHIFT_NORTH) & BMASK_P)),
           (Char((val >> BSHIFT_EAST) & BMASK_P)),
           (Char((val >> BSHIFT_SOUTH) & BMASK_P)),
           (Char((val >> BSHIFT_WEST) & BMASK_P)))
end

@inline function toval64(str::AbstractString)
    @assert(length(str)==sizeof(UInt64), length(str))
    hval = toval32(str[1:4])
    lval = toval32(str[5:8])
    Transition(hval) << BSHIFT_CELL|Transition(lval)
end

@inline function toval32(str::AbstractString)
    @assert(length(str)==sizeof(UInt32))
    Cell(str[1]) << BSHIFT_NORTH|
    Cell(str[2]) << BSHIFT_EAST|
    Cell(str[3]) << BSHIFT_SOUTH|
    Cell(str[4]) << BSHIFT_WEST
end

const loadfuncs = Dict(:CellSpace=>load_cell, :Rule=>load_rule)
function load(fname::String, typ::Symbol)
    if haskey(loadfuncs, typ)
        loadfuncs[typ](fname)
    else
        throw("no function for $typ")
    end
end

const updatefuncs = Dict(:Checkerboard=>update_checkerboard!, :Random=>update_random!)
function update!(cs::CellSpace, rule::Rule, typ::Symbol)
    if haskey(updatefuncs, typ)
        updatefuncs[typ](cs, rule)
    else
        throw("no function for $typ")
    end
end


text:
# The halfpsace functions are private methods, and all take a coordinate and a surface to determine which halfspace the coordinate is in.

# Plane halfspace determination
function halfspace(c::Coord, plane::Plane)
    d::Float64 = -dot(plane.normal, plane.point)
    half::Float64 = dot(plane.normal, c) + d
    if half <= 0
        return -1
    else
        return 1
    end
end

# Cone halfspace determination
function halfspace(c::Coord, cone::Cone)
	p_minus_c::Coord = c - cone.tip
    half::Float64 = dot(p_minus_c,cone.axis)^2 - dot(p_minus_c,p_minus_c) * (cos(cone.theta))^2
    if half <= 0
        return -1
    else
        return 1
    end
end

# Sphere halfspace determination
function halfspace(c::Coord, sphere::Sphere)
    half::Float64 = (c.x - sphere.center.x)^2 + (c.y - sphere.center.y)^2 + (c.z - sphere.center.z)^2 - sphere.radius^2
    if half <= 0
        return -1
    else
        return 1
    end
end

# Infinite cylinder halfspace
function halfspace(c::Coord, cyl::InfCylinder)
    tmp::Coord = cross((c-cyl.center), cyl.normal)
    half::Float64 = dot(tmp, tmp) - cyl.radius^2
    if half <= 0
        return -1
    else
        return 1
    end
end

function ^(a::Region, b::Region)
    if halfspace(_p, a.surface) == a.halfspace
        if halfspace(_p, b.surface) == b.halfspace
            return true
        end
    end
    return false
end

function ^(a::Region, b::Bool)
    if halfspace(_p, a.surface) == a.halfspace
        if b == true
            return true
        end
    end
    return false
end

function ^(b::Bool, a::Region)
    if halfspace(_p, a.surface) == a.halfspace
        if b == true
            return true
        end
    end
    return false
end

function |(a::Region, b::Region)
    if halfspace(_p, a.surface) == a.halfspace
        return true
    end
    if halfspace(_p, b.surface) == b.halfspace
        return true
    end
    return false
end

function |(a::Region, b::Bool)
    if halfspace(_p, a.surface) == a.halfspace
        return true
    end
    if b == true
        return true
    end
    return false
end

function |(b::Bool, a::Region)
    if halfspace(_p, a.surface) == a.halfspace
        return true
    end
    if b == true
        return true
    end
    return false
end

function ~(a::Region)
    b::Region = Region(a.surface, a.halfspace)
    if a.halfspace == -1
        b.halfspace = 1
    else
        b.halfspace = -1
    end
    return b
end


text:
using Luxor
using Lynx

include("utils.jl")

mutable struct SpriteSheet
    image::Image
    sprites::Dict{NTuple{4, Int}, Image}
    SpriteSheet() = new()
end

function SpriteSheet(image::Image)
    this = SpriteSheet()
    this.sprites = Dict{NTuple{4, Int}, Image}()
    this.image = image
    return this
end

SpriteSheet(path::AbstractString) = SpriteSheet(Image(path))

function Base.getindex(this::SpriteSheet, x::Int, y::Int, w::Int, h::Int)
    return get!(this.sprites, (x, y, w, h), this.image[x, y, w, h])
end

Lynx.init("SpriteSheet", 800, 600)

sprites = SpriteSheet(asset"fruits.png")
w, h = 128, 128
tiles = Tiler(800, 600, 3, 3)

function update(dt)
    background("#111")
    origin()
    for (pos, n) in tiles
        i, j = 16 .* (tiles.currentrow - 1, tiles.currentcol - 1)
        drawimage(sprites[i, j, 16, 16], pos.x, pos.y, w, h, centered = true)
    end
end

run!(update)


text:
using Distributions, StatsBase, Distributed

export draw_sample, draw_sample!, rand_multinomial_init

# Gibbs Move
# ----------

function imcmc_gibbs_update!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    i::Int,
    model::SIS,
    mcmc::SisMcmcInsertDeleteGibbs
    )
    # S_prop = copy(S_curr)

    # @inbounds m = rand_reflect(length(curr[i]), ν, 1, model.K_inner)

    @inbounds n = length(S_curr[i])
    δ = rand(1:mcmc.ν_gibbs)
    # a,b = (lb(n, δ, model), ub(n, δ))
    # @show n, δ, a, b
    d = rand(0:min(n,δ))
    m = n + δ - 2*d

    # Catch invalid proposal (m > K_inner). Here we imediately reject, making no changes.
    if (m > model.K_inner) | (m < 1)
        return 0
    end

    # @show m
    # 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)

    delete_insert!(S_prop[i], ind_del, ind_add, vals)

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

    # @show curr_dist, prop_dist
    @inbounds log_α = (
        -model.γ * (
            model.dist(model.mode, S_prop)-model.dist(model.mode, S_curr)
            ) + log_ratio
        )
    if log(rand()) < log_α
        # println("accepted")
        copy!(S_curr[i], S_prop[i])
        # @inbounds S_curr[i] = copy(S_prop[i])
        # println("$(i) value proposal $(tmp_prop) was accepted")
        return 1
    else
        copy!(S_prop[i], S_curr[i])
        # @inbounds S_prop[i] = copy(S_curr[i])
        # println("$(i) value proposal $(tmp_prop) was rejected")
        return 0
    end
    # @show S_curr
end

function imcmc_gibbs_scan!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    model::SIS,
    mcmc::SisMcmcInsertDeleteGibbs
    )
    count = 0
    N = length(S_curr)
    for i = 1:N
        count += imcmc_gibbs_update!(S_curr, S_prop, i, model, mcmc)
    end
    return count
end

# Sampler Functions
# -----------------

function draw_sample!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::SisMcmcInsertDeleteGibbs,
    model::SIS;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::InteractionSequence{Int}=get_init(model, mcmc.init)
    )

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], init[i])
        copy!(S_prop[i], init[i])
    end
    # @show S_curr, S_prop, init
    # S_curr = copy(init)
    # S_prop = copy(init)

    ind = 0
    sample_count = 1
    i = 0
    gibbs_scan_count = 0
    gibbs_tot_count = 0
    gibbs_acc_count = 0
    tr_dim_count = 0
    tr_dim_acc_count = 0

    # Bounds for uniform sampling
    lb(x::Vector{Path{Int}}) = max(1, length(x) - 1 )
    ub(x::Vector{Path{Int}}) = min(model.K_outer, length(x) + 1)

    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out[sample_count] = deepcopy(S_curr)
            sample_count += 1
        end
        # Gibbs scan
        if rand() < mcmc.β
            # println("Gibbs")
            gibbs_scan_count += 1
            gibbs_tot_count += length(S_curr)
            gibbs_acc_count += imcmc_gibbs_scan!(S_curr, S_prop, model, mcmc) # This enacts the scan, changing curr, and outputs number of accepted moves.
        # Else do insert or delete
        else
            # println("Transdim")
            tr_dim_count += 1
            tr_dim_acc_count += imcmc_trans_dim_accept_reject!(
                S_curr, S_prop,
                model, mcmc
            )
        end
    end
    # Send storage back
    # @show S_curr, S_prop
    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end
    return (
        gibbs_tot_count, gibbs_scan_count, gibbs_acc_count,
        tr_dim_count, tr_dim_acc_count
    )

end

function draw_sample(
    mcmc::SisMcmcInsertDeleteGibbs,
    model::SIS;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::Vector{Path{Int}}=get_init(model, mcmc.init)
    )
    sample_out = Vector{Vector{Path{Int}}}(undef, desired_samples)
    # @show sample_out
    draw_sample!(sample_out, mcmc, model, burn_in=burn_in, lag=lag, init=init)
    return sample_out

end

function (mcmc::SisMcmcInsertDeleteGibbs)(
    model::SIS;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::Vector{Path{Int}}=get_init(model, mcmc.init)
    )

    sample_out = Vector{Vector{Path{Int}}}(undef, desired_samples)
    # @show sample_out
    (
        gibbs_tot_count,
        gibbs_scan_count,
        gibbs_acc_count,
        tr_dim_count,
        tr_dim_acc_count
        ) = draw_sample!(sample_out, mcmc, model, burn_in=burn_in, lag=lag, init=init)

    p_measures = Dict(
            "Proportion Gibbs moves" => gibbs_scan_count/(tr_dim_count + gibbs_scan_count),
            "Trans-dimensional move acceptance probability" => tr_dim_acc_count/tr_dim_count,
            "Gibbs move acceptance probability" => gibbs_acc_count/gibbs_tot_count
        )
    output = SisMcmcOutput(
            model,
            sample_out,
            p_measures
            )

    return output

end


text:
using Mimi
using Test
using DataFrames
using VegaLite
using Electron

import Mimi:
    dataframe_or_scalar, _spec_for_item, menu_item_list, getdataframe, dim_names, time_labels

# Helper function returns true if VegaLite is verison 3 or above, and false otherwise
function _is_VegaLite_v3()
    return isdefined(VegaLite, :vlplot) ? true : false
end

include("../examples/compositecomp-model.jl") # constructs and runs model
# m's structure is as follows:
#
#          top
#        /    \
#       A       B
#     /  \     /  \
#    1    2   3    4

# 1.  dataframe helper functions
@test typeof(dataframe_or_scalar(m, :top, :fooA1)) == Float64 # same for fooA2, foo3, foo4
@test typeof(dataframe_or_scalar(m, :top, :var_3_1)) == DataFrame
@test typeof(dataframe_or_scalar(m, :top, :par_1_1)) == DataFrame

#2.  Specs and menu
items = [:fooA1, :fooA2, :foo3, :foo4, :var_3_1, :par_1_1]
for item in items
    static_spec = _spec_for_item(m, :top, item; interactive = false)
    interactive_spec = _spec_for_item(m, :top, item)
    if length(dim_names(m, :top, item)) == 0
        name =  string(:top, " : ", item, " = ", m[:top, item])
    else
        name = string(:top, " : ", item)
    end
    @test static_spec["name"] == interactive_spec["name"] == name
end

s = menu_item_list(m)
@test collect(keys(s)) == ["pars", "vars"]
@test length(collect(keys(s["pars"]))) == 5
@test length(collect(keys(s["vars"]))) == 1

#3.  explore(m::Model)
w = explore(m)
@test typeof(w) == Electron.Window
close(w)

#4.  Mim.plot(m::Model, comp_name::Symbol, datum_name::Symbol;
#       dim_name::Union{Nothing, Symbol} = nothing)
items = [:fooA1, :fooA2, :foo3, :foo4, :var_3_1, :par_1_1]
for item in items
    p_type = _is_VegaLite_v3() ? VegaLite.VLSpec : VegaLite.VLSpec{:plot}
    @test typeof(Mimi.plot(m, :top, item)) == p_type
end


text:
using Test, SafeTestsets

println("importing CVChannel")
@time using CVChannel

println("running ./test/CVChannel.jl")
@time begin # timing block

@testset "running ./test/runtests.jl" begin

    println("testing ./src/operations.jl")
    @time @safetestset "./test/operations.jl" begin include("operations.jl") end

    println("testing ./src/channels.jl")
    @time @safetestset "./test/channels.jl" begin include("channels.jl") end

    println("testing ./states.jl")
    @time @safetestset "./test/states.jl" begin include("states.jl") end

    println("testing ./src/optimizations.jl")
    @time @safetestset "./test/optimizations.jl" begin include("optimizations.jl") end

    println("testing ./src/optimizer_interface.jl")
    @time @safetestset "./test/optimizer_interface.jl" begin include("optimizer_interface.jl") end

    println("testing ./src/see-saw_optimization.jl")
    @time @safetestset "./test/see-saw_optimization.jl" begin include("see-saw_optimization.jl") end
end

println("\ntotal elapsed time :")
end # timing block

max_forks_repo_name: ChitambarLab/cv-channel


text:
function calc_order_fee(client, trades, config)
    com = 0
    for t in trades
        commission = typeof(t["commission"]) == Float64 ? t["commission"] : parse(Float64, t["commission"])
        if commission != 0
            if t["commissionAsset"] == config["farm_coin"]
                commission *= parse(Float64, t["price"])
            elseif t["commissionAsset"] != config["base_coin"]
                translator = string(t["commissionAsset"], config["base_coin"])
                weightedAvgPrice = try client.get_ticker(symbol=translator)["weightedAvgPrice"]
                catch
                    string(1/parse(Float64, client.get_ticker(symbol=string(config["farm_coin"], t["commissionAsset"]))["weightedAvgPrice"]))
                end
                if typeof(weightedAvgPrice) == String
                    weightedAvgPrice = parse(Float64, weightedAvgPrice)
                end
                commission *= weightedAvgPrice
            end
        end
        com += commission
    end
    return com
end
export calc_order_fee


text:
#RSDeltaSigma: Base functions to assist with filtering operations
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
"""`CoeffCSD`

A structure containing filter coefficents and their canonical-signed digit (CSD)
representation.
"""
struct CoeffCSD
	val
	csd
end

#==Accessors
===============================================================================#
csdsize2(f::Vector{CoeffCSD}) = [size(fi.csd,2) for fi in f] #VERIFYME

Base.values(f::Vector{CoeffCSD}) = [fi.val for fi in f]

#==Functions
===============================================================================#
"""`bquantize(x, nsd=3, abstol=eps, reltol=10*eps)`

Bidirectionally quantize a n by 1 vector x to nsd signed digits.

Terminate early if the error is less than the specified tolerances.

# Output
y is a CoeffCSD[] with the following fields:
 - `y[i].val` is the quantized value in floating-point form
 - `y[i].csd` is a 2-by-nsd (or less) matrix containing
   the powers of two (first row) and their signs (second row).

# See also:
[`bunquantize`](@ref)
"""
function bquantize(x, nsd::Int=3, abstol::Float64=eps(1.0), reltol::Float64=10*eps(1.0))
	n = length(x)
	q = zeros(2*n,nsd)
	y = Array{CoeffCSD}(undef, n)
	offset = -log2(0.75)

	for i in 1:n
		xp = x[i]
		val = 0
		csd = zeros(2,nsd)
		for j in 1:nsd
			error = abs(val-x[i])
			if error <= abstol || error <= abs(x[i])*reltol
				break
			end
			p = floor( log2(abs(xp)) + offset )
			p2 = 2^p
			sx = sign(xp)
			xp = xp - sx*p2
			val = val + sx*p2
			csd[1:2, j] = [p; sx]
		end
		y[i] = CoeffCSD(val, csd)
	end
	return y
end

"""`bunquantize(q)`

Calculate the value corresponding to a bidirectionally quantized quantity.
q is a 2n by m matrix containing the powers of 2 and their signs for each
quantized value.

# See also:
[`bquantize`](@ref)
"""
function bunquantize(q)
	q = 1.0*q #Ensure we have floating point values
	n = size(q,1)/2
	y = zeros(array_round(n),1)
	signs = 2:2:array_round(2*n)
	powers = signs .- 1
	for i in 1:size(q,2)
		y = y + 2 .^ q[powers,i] .* q[signs,i]
	end
	return y
end

function convertForm(f,q)
	n = length(f)
	F = CoeffCSD[]
	for i in 1:n
		val = f[i]
		csdrep = q[2*i-1:2*i,:]
		csd = csdrep[:,findall(csdrep[1,:] .!=0)]
		push!(F, CoeffCSD(val, csd))
	end
	return F
end

#Last line


text:
const SLEN = 3

autobet() = randbool(SLEN)
function autobet(ob::BitArray{1})
    rlen = length(ob)
    2 < rlen || return ~ob
    3 < rlen || return [~ob[2], ob[1:2]]
    opt = falses(rlen)
    opt[1] = true
    opt[end-1:end] = true
    ob != opt || return ~opt
    return opt
end
autobet(ob::Array{Bool,1}) = autobet(convert(BitArray{1}, ob))

function pgencode{T<:String}(a::T)
    b = uppercase(a)
    0 < length(b) || return trues(0)
    !ismatch(r"[^HT]+", b) || error(@sprintf "%s is not a HT sequence" a)
    b = split(b, "")
    b .== "H"
end
pgdecode(a::BitArray{1}) = join([i ? "H" : "T" for i in a], "")

function humanbet()
    b = ""
    while length(b) != SLEN  || ismatch(r"[^HT]+", b)
        print("Your bet? ")
        b = uppercase(chomp(readline()))
    end
    return b
end


text:
module Perceptron

using MetadataTools,  DocStringExtensions

export PerceptronClassifier, fit!, predict

type PerceptronClassifier{T}
    W::AbstractMatrix{T}
    b::AbstractVector{T}
    n_classes::Int
    n_features::Int
end

function Base.show{T}(io::IO, p::PerceptronClassifier{T})
    n_classes  = p.n_classes
    n_features = p.n_features
    print(io, "PerceptronClassifier{$T}(n_classes=$n_classes, n_features=$n_features)")
end

PerceptronClassifier(T::Type, n_classes::Int, n_features::Int) = PerceptronClassifier{T}(rand(T, n_features, n_classes),
                                                                                zeros(T, n_classes),
                                                                                n_classes,
                                                                                n_features)

"""
$(SIGNATURES)

Compute the accuracy betwwen `y` and `y_hat`.
"""
function accuracy(y::AbstractVector, y_hat::AbstractVector)
    acc = 0.
    @fastmath for k = 1:length(y)
            @inbounds  acc += y[k] == y_hat[k]
    end
    return acc/length(y_hat)
end

"""
$(SIGNATURES)

Predicts the class for a given input in a `PerceptronClassifier`.
The placeholder is used to avoid allocating memory for each matrix-vector multiplication.

- Returns the predicted class.
"""
function predict(h::PerceptronClassifier, x::AbstractVector, class_placeholder::AbstractVector)
    @fastmath class_placeholder .= At_mul_B!(class_placeholder, h.W, x) .+ h.b
    return indmax(class_placeholder)
end

"""
$(SIGNATURES)

Function to predict the class for a given example.

- Returns the predicted class.
"""
function predict(h::PerceptronClassifier, x::AbstractVector)
    score = h.W' * x .+ h.b
    return @fastmath indmax(score), score
end

"""
$(SIGNATURES)

Function to predict the class for a given input batch.
- Returns the predicted class.
"""
function predict(h::PerceptronClassifier, X::AbstractMatrix)
    predictions = zeros(Int32, size(X, 2))
    class_placeholder = zeros(eltype(h.W), h.n_classes)

    @inbounds for m in 1:length(predictions)
        predictions[m] = predict(h, view(X,:,m), class_placeholder)
    end
    return predictions
end

"""
$(SIGNATURES)

>    fit!(h::PerceptronClassifier,
>         X::Array,
>         y::Array;
>         n_epochs=50,
>         learning_rate=0.1,
>         print_flag=false,
>         compute_accuracy=true,
>         seed=srand(1234),
>         pocket=false,
>         shuffle_data=false)

##### Arguments

- **`h`**, (PerceptronClassifier{T} type), Multiclass perceptron.
- **`X`**, (Array{T,2} type), data contained in the columns of X.
- **`y`**, (Vector{T} type), class labels (as integers from 1 to n_classes).

##### Keyword arguments

- **`n_epochs`**, (Integer type), number of passes (epochs) through the data.
- **`learning_rate`**, (Float type), learning rate (The standard perceptron is with learning_rate=1.)
- **`compute_accuracy`**, (Bool type), if `true` the accuracy is computed at the end of every epoch.
- **`print_flag`**, (Bool type), if `true` the accuracy is printed at the end of every epoch.
- **`seed`**, (MersenneTwister type), seed for the permutation of the datapoints in case there the data is shuffled.
- **`pocket`** , (Bool type), if `true` the best weights are saved (in the pocket) during learning.
- **`shuffle_data`**, (Bool type),  if `true` the data is shuffled at every epoch (in reality we only shuffle indicies for performance).

"""
function fit!(h::PerceptronClassifier, X::AbstractArray, y::AbstractVector, scores::Array;
              n_epochs=50, learning_rate=1., print_flag=false,
              compute_accuracy=false, seed=srand(1234), pocket=false,
              shuffle_data=false)

    n_features, n_samples = size(X)
    @assert length(y) == n_samples

    T = eltype(X)
    learning_rate = T(learning_rate)

    class_placeholder = zeros(T, h.n_classes)
    y_preds = zeros(Int64, n_samples)
    data_indices = Array(1:n_samples)
    max_acc = zero(T)

    if pocket
        W_hat = zeros(T, h.n_features, h.n_classes)
        b_hat = zeros(T, h.n_classes)
    end

    @fastmath for epoch in 1:n_epochs

        n_mistakes = 0
        if shuffle_data
            shuffle!(seed, data_indices)
        end

        @inbounds for m in data_indices
            x = view(X, :, m)

            y_hat = predict(h, x, class_placeholder)

            if y[m] != y_hat
                n_mistakes += 1
                ####  wij ← wij − η (yj −tj) · xi
                h.W[:, y[m]]  .= h.W[:, y[m]]  .+ learning_rate .* x
                h.b[y[m]]     .= h.b[y[m]]     .+ learning_rate
                h.W[:, y_hat] .= h.W[:, y_hat] .- learning_rate .* x
                h.b[y_hat]    .= h.b[y_hat]    .- learning_rate
            end
        end

        if compute_accuracy
             @inbounds for m in  data_indices
                 y_preds[m] = predict(h, view(X, :, m), class_placeholder)
            end
            acc = accuracy(y, y_preds)
            push!(scores, acc)
        else
            acc = (n_samples - n_mistakes)/n_samples
            push!(scores, acc)
        end

        if pocket
            if acc > max_acc
                max_acc = acc
                copy!(W_hat, h.W)
                copy!(b_hat, h.b)
            end
        end

        if print_flag
            print("Epoch: $(epoch) \t Accuracy: $(round(acc,3))\r")
            flush(STDOUT)
        end
    end

end

end # module


text:
struct SplitHLAData <: AbstractHLAData
    name::String
    records::Vector{FASTX.FASTA.Record}
    hla_types::Vector{HLAType}
    tree::Union{Missing, PhylogeneticTree}
    stan_input::Dict{String, Any}
    allele_depth::Int
    replacement_mincount::Int
    original_data::HLAData
    split::UnitRange{Int64}
    split_idx::Int

    function SplitHLAData(
        name::String,
        records::Vector{FASTX.FASTA.Record},
        hla_types::Vector{HLAType},
        tree::Union{Missing, PhylogeneticTree},
        stan_input::Union{Missing, Dict{String, Any}},
        allele_depth::Int,
        replacement_mincount::Int,
        original_data::HLAData,
        split::UnitRange{Int64},
        split_idx::Int
    )

        if length(hla_types) != length(records)
            msg = "Vector of HLATypes has size $(length(hla_types)), " *
                "expected $(length(records))"
            error(msg)
        end

        if !ismissing(tree)
            m = Escape.matching(tree, records)
            m isa Exception && throw(m)
        end

        if !ismissing(stan_input)
            m = ErrorException("stan_input is not in a valid format.")
            is_valid(stan_input) || throw(m)
        end

        new(name, records, hla_types, tree, stan_input, allele_depth, replacement_mincount,
            original_data, split, split_idx)
    end
end

function split_hla_data(data::HLAData, n::Int = 10)
    split_data = SplitHLAData[]
    replacements = Escape.replacements(data)
    stan_input = Escape.stan_input(data)
    R = stan_input["R"]
    splits = Iterators.partition(1:R, ceil(Int, R // n))

    for (i, split) in collect(enumerate(splits))
        stan_input_split = copy(stan_input)
        in_split = findall(x -> x in split, stan_input_split["rs"])

        stan_input_split["Z"] = stan_input_split["Z"][split, :]
        stan_input_split["rs"] = map(x -> findfirst(y -> y == x, split),
        stan_input_split["rs"][in_split])
        stan_input_split["idx"] = stan_input_split["idx"][in_split]
        stan_input_split["y"] = stan_input_split["y"][in_split]
        stan_input_split["N"] = length(in_split)
        stan_input_split["R"] = length(split)
        stan_input_split["phy"] = stan_input_split["phy"][split, :]

        split_data_ = SplitHLAData(data.name, data.records, data.hla_types, data.tree,
            stan_input_split, allele_depth(data), replacement_mincount(data), data, split, i)

        @assert length(Escape.replacements(split_data_)) == stan_input_split["R"]

        push!(split_data, split_data_)
    end

    return split_data
end


text:
include("data/player.jl")

player = Player("test-player")

Player_AppendStatus!(player, Status("test", -1))

for i in player.statusEffects
    print(i)
end

max_stars_repo_name: Flobii/discord-bot

max_issues_repo_name: Flobii/discord-bot

max_forks_repo_name: Flobii/discord-bot


text:
N = Float64

c = zeros(N, 2)
E1, E2 = diagm(N[-1, 0.5]), N[1 1; 0.5 0.3]
E = [E1, E2];
F2 = N[-0.5 1]';
F = [F2];
G = diagm(N[0.3, 0.3]);
p = PolynomialZonotope(c, E, F, G)

@test dim(p) == 2
@test order(p) == 7//2
@test polynomial_order(p) == 2

# type-specific concrete methods
scale(N(0.5), p)
linear_map(N[1.0 2.0; 2.0 5.0], p)
z = Zonotope(N[1., 2.], N(1.) * eye(N, 2))
minkowski_sum(p, z)
minkowski_sum(z, p)


text:
#
# Correctness Tests
#
using LsqFit, Test, LinearAlgebra, Random
using OptimBase, Calculus
import NLSolversBase: OnceDifferentiable

my_tests = ["curve_fit.jl", "levenberg_marquardt.jl", "curve_fit_inplace.jl", "geodesic.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
## Expectation
## ===========
function sample_mean{T<:Real}(X::RandVar{T}; nsamples = 1000)
  mean(rand(X, nsamples))
end

# probability found using samples
function prob_sampled(X::RandVar; nsamples = 1000)
  ntrue = 0
  for i=1:nsamples if rand(X) ntrue += 1 end end
  ntrue/nsamples
end


text:
function dJv_seriesdk(k2::T,v::Int64) where {T <: Real}
# Use series expansion to compute J_v:
nmax = 100
n = 1; error = Inf; if k2 < 1; tol = eps(k2); else; tol = eps(inv(k2)); end
# Computing leading coefficient (n=0):
#coeff = 3pi/(2^(2+v)*factorial(v+2))
if k2 < 1
  coeff = .75*pi/exp(lfact(v+2))
#  println("coefficient: ",coeff)
# multiply by (2v-1)!!
  for i=1:v
    coeff *= (2i-1)/2
  end
# Add leading term to J_v:
  Jv = one(k2)*coeff
  dJvdk = one(k2)*coeff*(2v+1)
# Now, compute higher order terms until desired precision is reached:
  while n < nmax && abs(error) > tol
    coeff *= (2n-1)*(2(n+v)-1)/(2n*(2n+2v+4))*k2
    Jv += coeff
    dJvdk += coeff*(2n+2v+1)
#    error = coeff/Jv
    error = coeff
    n += 1
  end
  dJvdk *= k2^v
  Jv *= k2^v*sqrt(k2)
#  println("Jv: ",Jv," dJv/dk: ",dJvdk)
  return Jv,dJvdk
else # k^2 >= 1
  coeff = convert(typeof(k2),pi)
  # Compute (2v-1)!!/(2^v v!):
  for i=1:v
#    coeff *= 1.-.5/i
    coeff *= (2i-1)/(2i)
  end
  Jv = one(k2)*coeff
  dJvdk = zero(k2)
  k2inv = inv(k2); n=2
  while n < nmax && abs(error) > tol
#    coeff *= (1.-2.5/n)*(1.-.5/(n+v))*k2inv
    coeff *= (n-5)*(n+2v-1)/(n*(n+2v))
    coeff *= k2inv
    Jv += coeff
    dJvdk -= n*coeff
#    error = coeff/Jv
    error = coeff
    n += 2
  end
  dJvdk /= sqrt(k2)
  return Jv,dJvdk
end
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
"""
Store a diagonally-scaled Toeplitz∘Hankel matrix:

    DL(T∘H)DR

where the Hankel matrix `H` is non-negative definite. This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹.
"""
immutable ToeplitzHankelPlan{S}
    T::TriangularToeplitz{S}
    C::Vector{Vector{S}}
    DL::Vector{S}
    DR::Vector{S}
    ToeplitzHankelPlan(T,C,DL,DR)=new(T,C,DL,DR)
end

function ToeplitzHankelPlan(T::TriangularToeplitz,C::Vector,DL::AbstractVector,DR::AbstractVector)
    S=promote_type(eltype(T),eltype(C[1]),eltype(DL),eltype(DR))
    ToeplitzHankelPlan{S}(T,C,collect(S,DL),collect(S,DR))
end
ToeplitzHankelPlan(T::TriangularToeplitz,C::Matrix) =
    ToeplitzHankelPlan(T,C,ones(size(T,1)),ones(size(T,2)))
ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,DL::AbstractVector,DR::AbstractVector) =
    ToeplitzHankelPlan(T,partialchol(H),DL,DR)
ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,D::AbstractVector,DL::AbstractVector,DR::AbstractVector) =
    ToeplitzHankelPlan(T,partialchol(H,D),DL,DR)

*(P::ToeplitzHankelPlan,v::AbstractVector)=P.DL.*toeplitzcholmult(P.T,P.C,P.DR.*v)

function partialchol(H::Hankel)
    # Assumes positive definite
    σ=eltype(H)[]
    n=size(H,1)
    C=Vector{eltype(H)}[]
    v=[H[:,1];vec(H[end,2:end])]
    d=diag(H)
    @assert length(v) ≥ 2n-1
    reltol=maxabs(d)*eps(eltype(H))*log(n)
    for k=1:n
        mx,idx=findmax(d)
        if mx ≤ reltol break end
        push!(σ,inv(mx))
        push!(C,v[idx:n+idx-1])
        for j=1:k-1
            nCjidxσj = -C[j][idx]*σ[j]
            Base.axpy!(nCjidxσj, C[j], C[k])
        end
        @simd for p=1:n
            @inbounds d[p]-=C[k][p]^2/mx
        end
    end
    for k=1:length(σ) scale!(C[k],sqrt(σ[k])) end
    C
end

function partialchol(H::Hankel,D::AbstractVector)
    # Assumes positive definite
    T = promote_type(eltype(H),eltype(D))
    σ=T[]
    n=size(H,1)
    C=Vector{T}[]
    v=[H[:,1];vec(H[end,2:end])]
    d=diag(H).*D.^2
    @assert length(v) ≥ 2n-1
    reltol=maxabs(d)*eps(T)*log(n)
    for k=1:n
        mx,idx=findmax(d)
        if mx ≤ reltol break end
        push!(σ,inv(mx))
        push!(C,v[idx:n+idx-1].*D.*D[idx])
        for j=1:k-1
            nCjidxσj = -C[j][idx]*σ[j]
            Base.axpy!(nCjidxσj, C[j], C[k])
        end
        @simd for p=1:n
            @inbounds d[p]-=C[k][p]^2/mx
        end
    end
    for k=1:length(σ) scale!(C[k],sqrt(σ[k])) end
    C
end

function toeplitzcholmult(T,C,v)
    n,K = length(v),length(C)
    ret,temp1,temp2 = zero(v),zero(v),zero(v)
    un,ze = one(eltype(v)),zero(eltype(v))
    broadcast!(*, temp1, C[K], v)
    A_mul_B!(un, T, temp1, ze, temp2)
    broadcast!(*, ret, C[K], temp2)
    for k=K-1:-1:1
        broadcast!(*, temp1, C[k], v)
        A_mul_B!(un, T, temp1, ze, temp2)
        broadcast!(*, temp1, C[k], temp2)
        broadcast!(+, ret, ret, temp1)
    end
    ret
end

# Diagonally-scaled Toeplitz∘Hankel polynomial transforms

function leg2chebTH{S}(::Type{S},n)
    λ = Λ(0:half(S):n-1)
    t = zeros(S,n)
    t[1:2:end] = λ[1:2:n]
    T = TriangularToeplitz(2t/π,:U)
    H = Hankel(λ[1:n],λ[n:end])
    DL = ones(S,n)
    DL[1] /= 2
    T,H,DL
end

function cheb2legTH{S}(::Type{S},n)
    t = zeros(S,n-1)
    t[1:2:end] = Λ(0:one(S):div(n-2,2),-half(S),one(S))
    T = TriangularToeplitz(t,:U)
    h = Λ(1:half(S):n-1,zero(S),3half(S))
    H = Hankel(h[1:n-1],h[n-1:end])
    D = 1:one(S):n-1
    DL = (3half(S):n-half(S))./D
    DR = -(one(S):n-one(S))./4D
    T,H,D,DL,DR
end

function leg2chebuTH{S}(::Type{S},n)
    λ = Λ(0:half(S):n-1)
    t = zeros(S,n)
    t[1:2:end] = λ[1:2:n]./(((1:2:n)-2))
    T = TriangularToeplitz(-2t/π,:U)
    H = Hankel(λ[1:n]./((1:n)+1),λ[n:end]./((n:2n-1)+1))
    T,H
end

function ultra2ultraTH{S}(::Type{S},n,λ₁,λ₂)
    @assert abs(λ₁-λ₂) < 1
    DL = (zero(S):n-one(S))+λ₂
    jk = 0:half(S):n-1
    t = zeros(S,n)
    t[1:2:n] = Λ(jk,λ₁-λ₂,one(S))[1:2:n]
    T = TriangularToeplitz(scale!(inv(gamma(λ₁-λ₂)),t),:U)
    h = Λ(jk,λ₁,λ₂+one(S))
    scale!(gamma(λ₂)/gamma(λ₁),h)
    H = Hankel(h[1:n],h[n:end])
    DR = ones(S,n)
    T,H,DL,DR
end

function jac2jacTH{S}(::Type{S},n,α,β,γ,δ)
    @assert β == δ
    @assert abs(α-γ) < 1
    @assert α+β > -1
    jk = zero(S):n-one(S)
    DL = (2jk+γ+β+one(S)).*Λ(jk,γ+β+one(S),β+one(S))
    T = TriangularToeplitz(Λ(jk,α-γ,one(S)),:U)
    H = Hankel(Λ(jk,α+β+one(S),γ+β+two(S)),Λ(jk+n-one(S),α+β+one(S),γ+β+two(S)))
    DR = Λ(jk,β+one(S),α+β+one(S))/gamma(α-γ)
    T,H,DL,DR
end

immutable ChebyshevToLegendrePlan{TH}
    toeplitzhankel::TH
end

ChebyshevToLegendrePlan{S}(::Type{S},n) = ChebyshevToLegendrePlan(th_cheb2legplan(S,n))

function *(P::ChebyshevToLegendrePlan,v::AbstractVector)
    w = zero(v)
    S,n = eltype(v),length(v)
    w[1:2:end] = -one(S)./(one(S):two(S):n)./(-one(S):two(S):n-two(S))
    [dot(w,v);P.toeplitzhankel*view(v,2:n)]
end

th_leg2chebplan{S}(::Type{S},n)=ToeplitzHankelPlan(leg2chebTH(S,n)...,ones(S,n))
th_cheb2legplan{S}(::Type{S},n)=ChebyshevToLegendrePlan(ToeplitzHankelPlan(cheb2legTH(S,n)...))
th_leg2chebuplan{S}(::Type{S},n)=ToeplitzHankelPlan(leg2chebuTH(S,n)...,1:n,ones(S,n))
th_ultra2ultraplan{S}(::Type{S},n,λ₁,λ₂)=ToeplitzHankelPlan(ultra2ultraTH(S,n,λ₁,λ₂)...)
th_jac2jacplan{S}(::Type{S},n,α,β,γ,δ)=ToeplitzHankelPlan(jac2jacTH(S,n,α,β,γ,δ)...)

th_leg2cheb(v)=th_leg2chebplan(eltype(v),length(v))*v
th_cheb2leg(v) = th_cheb2legplan(eltype(v),length(v))*v
th_leg2chebu(v)=th_leg2chebuplan(eltype(v),length(v))*v
th_ultra2ultra(v,λ₁,λ₂)=th_ultra2ultraplan(eltype(v),length(v),λ₁,λ₂)*v
th_jac2jac(v,α,β,γ,δ)=th_jac2jacplan(eltype(v),length(v),α,β,γ,δ)*v


text:
"""
```julia
GibbsSampling(;ϵ::T=1e-5,nBurnin::Int=100,samplefrequency::Int=1)
```

Draw samples from the true posterior via Gibbs Sampling.

**Keywords arguments**
    - `ϵ::T` : convergence criteria
    - `nBurnin::Int` : Number of samples discarded before starting to save samples
    - `samplefrequency::Int` : Frequency of sampling
"""
mutable struct GibbsSampling{T<:Real,N} <: SamplingInference{T}
    nBurnin::Integer # Number of burnin samples
    samplefrequency::Integer # Frequency at which samples are saved
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 # Number of samples
    nMinibatch::Int64
    ρ::Float64
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    opt::NTuple{N,SOptimizer}
    sample_store::AbstractArray{T,3}
    xview::SubArray{T,2,Matrix{T}}#,Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true}
    yview::SubArray
    function GibbsSampling{T}(nBurnin::Int,samplefrequency::Int,ϵ::Real) where {T}
        @assert nBurnin >= 0 "nBurnin should be a positive integer"
        @assert samplefrequency >= 0 "samplefrequency should be a positive integer"
        return new{T,1}(nBurnin,samplefrequency,ϵ)
    end
    function GibbsSampling{T,1}(nBurnin::Int,samplefrequency::Int,ϵ::Real,nFeatures::Int,nSamples::Int,nMinibatch::Int,nLatent::Int) where {T}
        opts = ntuple(_->SOptimizer{T}(VanillaGradDescent(η=1.0)),nLatent)
        new{T,nLatent}(nBurnin,samplefrequency,ϵ,0,false,nSamples,nMinibatch,nSamples/nMinibatch,true,opts)
    end
end

function GibbsSampling(;ϵ::T=1e-5,nBurnin::Int=100,samplefrequency::Int=1) where {T<:Real}
    GibbsSampling{Float64}(nBurnin,samplefrequency,ϵ)
end

function Base.show(io::IO,inference::GibbsSampling{T}) where {T<:Real}
    print(io,"Gibbs Sampler")
end

function tuple_inference(i::TSamp,nLatent::Integer,nFeatures::Integer,nSamples::Integer) where {TSamp <: GibbsSampling}
    return TSamp(i.nBurnin,i.samplefrequency,i.ϵ,nFeatures,nSamples,nSamples,nLatent)
end

function init_sampler(inference::GibbsSampling{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer,cat_samples::Bool) where {T<:Real}
    if inference.nIter == 0 || !cat_samples
        inference.sample_store = zeros(T,nSamples,nFeatures,nLatent)
    else
        inference.sample_store = cat(inference.sample_store,zeros(T,nSamples,nFeatures,nLatent),dims=1)
    end
    return inference
end

function sample_parameters(model::MCGP{T,L,<:GibbsSampling},nSamples::Int,callback::Union{Nothing,Function},cat_samples::Bool) where {T,L}
    init_sampler(model.inference,model.nLatent,model.nFeatures,nSamples,cat_samples)
    computeMatrices!(model)
    for i in 1:(model.inference.nBurnin+nSamples*model.inference.samplefrequency)
        model.inference.nIter += 1
        sample_local!(model.likelihood,get_y(model),get_f(model))
        sample_global!.(
        ∇E_μ(model.likelihood,model.inference.opt[1],get_y(model)),
        ∇E_Σ(model.likelihood,model.inference.opt[1],get_y(model)),
        model.f)
        if model.inference.nIter > model.inference.nBurnin && (model.inference.nIter-model.inference.nBurnin)%model.inference.samplefrequency==0
            store_variables!(model.inference,get_f(model))
        end
    end
    symbols = ["f_"*string(i) for i in 1:model.nFeatures]
    chains = [Chains(reshape(model.inference.sample_store[:,:,i],:,model.nFeatures,1),symbols) for i in 1:model.nLatent]
end

sample_local!(l::Likelihood,y,f::Tuple{<:AbstractVector{T}}) where {T} =sample_local!(l,y,first(f))
set_ω!(l::Likelihood,ω) = l.θ .= ω
get_ω(l::Likelihood) = l.θ

function logpdf(model::AbstractGP{T,<:Likelihood,<:GibbsSampling}) where {T}
    return 0.0
end

function sample_global!(∇E_μ::AbstractVector,∇E_Σ::AbstractVector,gp::_MCGP{T}) where {T}
    # @info ∇E_Σ
    # @info eigvals(Matrix(gp.K))
    global Σ = inv(Symmetric(2.0*Diagonal(∇E_Σ)+inv(gp.K)))
    gp.f .= rand(MvNormal(Σ*(∇E_μ+gp.K\gp.μ₀),Σ))
    return nothing
end

function store_variables!(i::SamplingInference{T},fs) where {T}
    i.sample_store[(i.nIter-i.nBurnin)÷i.samplefrequency,:,:] .= hcat(fs...)
end

function post_process!(model::AbstractGP{T,<:Likelihood,<:GibbsSampling}) where {T}
    for k in 1:model.nLatent
        model.μ[k] = vec(mean(hcat(model.inference.sample_store[k]...),dims=2))
        model.Σ[k] = Symmetric(cov(hcat(model.inference.sample_store[k]...),dims=2))
    end
    nothing
end


text:
export colon_ideal

order(a::NfOrdGenFracIdlSet) = a.order

order(a::NfOrdGenFracIdl) = order(a.parent)

function colon_ideal(a::NfOrdGenIdl)
  B = basis(a)
  O = order(a)
  @vprint :NfOrd 1 "$(representation_mat(B[1]))\n"
  @vprint :NfOrd 1 FakeFmpqMat(representation_mat(B[1]),ZZ(1))*basis_mat_inv(a)
  m = to_fmpz_mat(FakeFmpqMat(representation_mat(B[1]),ZZ(1))*basis_mat_inv(a))
  for i in 2:degree(O)
    m = hcat(to_fmpz_mat(FakeFmpqMat(representation_mat(B[i]),ZZ(1))*basis_mat_inv(a)),m)
  end
  n = hnf(transpose(m))
  # n is upper right HNF
  n = transpose(sub(n,1:degree(O),1:degree(O)))
  b,d = pseudo_inv(n)
  return NfOrdGenFracIdl(O,hnf(FakeFmpqMat(b,d)))
end

function NfOrdGen(a::NfOrdGenFracIdl)
  z = NfOrdGen(nf(order(a)), a.basis_mat*order(a).basis_mat)
  return z
end


text:
module SmoothLivePlot

using Observables, WebIO

include("plotFunctions.jl")
export modifyPlotObject!

include("plotMacros.jl")
export @makeLivePlot

end


text:
using PkgSRA

"""
Steps 1 and 3: fit a distribution of models

calc_distribution_of_models(dat,
                             numerical_grad,
                             sindy_library;
                             U=nothing,
                             selection_dist=Normal(0,10),
                             use_clustering_minimization=false,
                             chain_opt=(iterations=200,
                                        num_training_pts=100,
                                        start_ind=101),
                                        val_list = Iterators.product(1:3,1:3))

Returns both the Turing chain and the sindycModel used as the
        priors for the coefficients.
"""
function calc_distribution_of_models(dat,
                                     numerical_grad,
                                     sindy_library;
                                     U=nothing,
                                     selection_dist=Normal(0,10),
                                     use_clustering_minimization=false,
                                     chain_opt=(iterations=200,
                                                num_training_pts=100,
                                                start_ind=101),
                                                val_list = Iterators.product(1:3,1:3)) # DEBUG
   # Use AIC to find the best model
   (best_sindy_model,best_criterion,all_criteria,all_models, best_index) =
       sindyc_ensemble(dat, numerical_grad, sindy_library, val_list,
                       U=U,
                       selection_criterion=my_aicc,
                       sparsification_mode="num_terms",
                       selection_dist=selection_dist,
                       use_clustering_minimization=use_clustering_minimization)

   # Use Turing to fit a distribution of models
   turing_model = convert_sindy_to_turing_enforce_zeros(best_sindy_model;
                                   dat_noise_prior=Normal(0.0, 20.0),
                                   coef_noise_std=1.0)
   chain = generate_chain(dat, numerical_grad,
                          turing_model;
                          chain_opt...)[1]

   return (chain=chain, best_sindy_model=best_sindy_model)
end

"""
Step 2: Get the distribution of residuals

calc_distribution_of_residuals(dat, numerical_grad, chain,
                                sample_ind,
                                sindy_template)

Returns:
        residual
        sample_gradients
        sample_noise
        dat_grad
"""
function calc_distribution_of_residuals(dat, numerical_grad, chain,
                                        sample_ind,
                                        sindy_template)
    sample_gradients3d, sample_noise2d = sample_sindy_posterior_grad(chain,
                dat, sample_ind, sindy_template)

    # Calculate the residuals
    sample_gradients = transpose(drop_all_1dims(mean(sample_gradients3d, dims=1)))
        # sample_gradients_std = transpose(drop_all_1dims(std(sample_gradients3d, dims=1)))
        sample_noise = mean(sample_noise2d)
    dat_grad = numerical_grad[:, sample_ind]

    residual = dat_grad .- sample_gradients

    return (residual=residual,
            sample_gradients=sample_gradients,
            sample_noise=sample_noise,
            dat_grad=dat_grad)
end

"""
Step 3: Subsample the data using the residuals

function subsample_using_residual(residual, noise;
                                  noise_factor=1.0,
                                  min_length=1,
                                  shorten_length=1)
return accepted_ind
"""
function subsample_using_residual(residual, noise;
                                  noise_factor=1.0,
                                  min_length=1,
                                  shorten_length=1)
    partial_accepted_ind = abs.(residual) .< (noise_factor*noise)
    partial_accepted_ind = vec(prod(Int.(partial_accepted_ind), dims=1))
    _, _, all_starts, all_ends = calc_contiguous_blocks(
            partial_accepted_ind,
            minimum_length=min_length + 2*shorten_length)

    # Get the "good" indices
    accepted_ind = []
    for (s, e) in zip(all_starts, all_ends)
        accepted_ind = vcat(accepted_ind,
                (s+shorten_length):(e-shorten_length))
    end
    return accepted_ind
end

"""
Optional step 0: Randomly subsample the initial data (for the naive model)

calc_best_random_subsample(model_template, dat2, numerical_grad2, sindy_library;
                                    num_pts=400,
                                    num_subsamples=10,
                                    val_list = Iterators.product(1:3,1:3),
                                    sindyc_ensemble_params=get_sindyc_ensemble_parameters(),
                                    use_control=true,
                                    max_error=nothing)
return (best_initial_subsample=best_initial_subsample,
        best_model=all_final_models2[tmp_ind],
        all_initial_subsamples=initial_subsamples,
        all_errL2=all_errL2)
"""
function calc_best_random_subsample(model_template, dat2,
                        numerical_grad2;
                        num_pts=400,
                        num_subsamples=20,
                        val_list = Iterators.product(1:3,1:3),
                        sindyc_ensemble_params=get_sindyc_ensemble_parameters(),
                        use_control=true,
                        max_error=0)
    # Initially, randomly subsample the data
    initial_subsamples = []
    for i in 1:num_subsamples
        push!(initial_subsamples, randperm(size(dat2,2))[1:num_pts])
    end

    println("Calculating best random subsample; control=$use_control)")

    # Just do SINDY here, not Turing yet
    sz = size(initial_subsamples)
    all_final_models2 = Vector{DynamicalSystemModel}(undef,sz)
    all_errL2 = zeros(sz)
    current_error = max_error+1
    max_iter = 2
    this_iter = 1
    if use_control; f = sindyc_ensemble; else; f = sindy_ensemble end
    while current_error > max_error
        if this_iter > 1
            println("Random subsampling gave a large error ($current_error), retrying")
        end
        for (i, inds) in enumerate(initial_subsamples)
            i%10==0 && println("Testing random subsample $i/$(length(initial_subsamples))")
            (all_final_models2[i], all_errL2[i], _, _, _) =
                 f(model_template, dat2[:,inds], numerical_grad2[:, inds],
                    val_list; sindyc_ensemble_params...)
            # if sindyc_ensemble_params[:selection_criterion] == my_aicc
            #     # Calculate real L2 error
            #     this_dat = all_final_models2[i](dat2, 0)
            #     all_errL2[i] = sum(abs, numerical_grad2.-this_dat)
            # else
            #     # Should be cross validation; already L2 error
            # end
        end
        current_error = minimum(all_errL2)
        this_iter += 1
        if max_error == 0 || this_iter >= max_iter
            break
        end
    end

    # Choose the best of the above models and subsamples
    min_err, tmp_ind = findmin(all_errL2)
    best_initial_subsample = initial_subsamples[tmp_ind]

    return (best_initial_subsample=best_initial_subsample,
            best_model=all_final_models2[tmp_ind],
            all_initial_subsamples=initial_subsamples,
            all_errL2=all_errL2)
end

#####
#####
#####
export calc_distribution_of_models,
        calc_distribution_of_residuals,
        subsample_using_residual,
        calc_best_random_subsample