text:
include("SCM.jl")

# function to select the variables to intervene on.
"""
    inter_var_selection(scm::LinAddSCM, target::Int64, E::Int64, share::Float64; overlap = 0)

Select variables to intervene on in each interventional environment. The first environment being considered to be observational, it splits the variables in E-1 parts, corresponding to each interventional environment. Return array containing each of these E-1 groups of variables. It is possible to let the algo choose only `share`% of the variables to be intervened on in at least one environment. Also, `overlap` can be use to make two different environments to share common intervened variables.

* `scm`    : linear structural causal model
* `target` : index of target variable
* `E`      : nbr of environments (including observatinal one, hence E > 1)
* `share`  : share of the variables to get intervened on in at least one environment
* `overlap`: every two concecutive env have in common `overlap`% of the intervened variables
"""
function inter_var_selection(scm::SCM, target::Int64, E::Int64, share::Float64; overlap = 0.2)
    # init
    @assert E > 1 <= target <= scm.p && 0 <= share <= 1 >= overlap >= 0
    E₀ = E - 1 # nbr of interventional envs
    p::Int64 = floor(share*(scm.p - 1)) # nbr of vars to intervene on, -1 is for the target
    vars = deleteat!(Vector(1:scm.p), target)
    inter_variables::Array{Array{Int64,1},1} = []

    # generic fct, test if var is intervened on
    test_inter(x, l, u) = (l * p <= x-1 < u * p) || (l * p <= x-1 - p) || (x-1 + p < u * p)

    # random selection
    vars = vars[sortperm(rand(scm.p-1))[1:p]]
    cst = 1 / E₀
    for e in 1:E₀
        l = (e-1) * cst - overlap / 2
        u = e * cst + overlap / 2
        push!(inter_variables, vars[filter(x -> test_inter(x, l, u), 1:p)])
    end

    return inter_variables
end

# Choose which variables to invene on in each environment with a simple hack:
# With probability 1/3 intervene on 1 variable, 1/3 intervene on 3 variables,
# 1/3 intervene on 5 variables.
function simple_var_selection(scm::SCM, target::Int64, to_intervene = [1,3,5])
    @assert E > 1 <= target <= scm.p
    E₀ = E - 1 # nbr of interventional envs
    p = scm.p - 1 # nbr of vars to intervene on, -1 is for the target
    vars = deleteat!(Vector(1:scm.p), target)
    inter_variables::Array{Array{Int64,1},1} = []
    for e in 1:E₀
        # with probability 1/3 select either 1, 3, or 5 variables to intervene in env. e
        num_intervene = to_intervene[rand(Categorical([1/3,1/3,1/3]))]
        push!(inter_variables, sample(vars, num_intervene, replace=false))
    end
    return inter_variables
end

"""
    inter_distr_selection(inter_variables::Array{Array{Int64,1},1} ; lb=-2, ub=2, var_min=0.5, var_max=2)

Generate the noise distributions (~N) for the stochastic interventions. Each of these distr. is Normal(μ,σ), and we randomly select μ and σ. Return a list of tuples (one per interventional environments) of the same type as `inter_vect` in `stochastic_inter_sim`.

* `inter_variables`   : output of `inter_var_selection`, i.e. array of arrays containing the variables to intervene on in each env.
* `lb`                : lower bound for unif distr. generating a μ (actually we also randomly flip the sign of μ after)
* `ub`                : upper bound for unif distr. generating a μ
* `var_min`           : lower bound for unif distr. generating a σ^2
* `var_max`           : upper bound for unif distr. generating a σ^2
"""
function inter_distr_selection(inter_variables::Array{Array{Int64,1},1} ;
                                    lb=-2, ub=2, var_min=0.5^2, var_max=2^2)
    # init
    @assert 0 <= var_min <= var_max
    E = length(inter_variables)
    nSim = sum(length.(inter_variables))
    inter_distr::Array{Array{T,1},1} where T <: UnivariateDistribution = []

    # draw params
    means = (rand(nSim) .* (ub - lb) .+ lb) .* rand([-1,1], nSim)
    vars = rand(nSim) .* (var_max - var_min) .+ var_min

    # construct array of distr
    idx = 1
    for e in 1:E
        inter_dist_sub::Array{T,1} where T <:UnivariateDistribution = []
        for var in inter_variables[e]
            push!(inter_dist_sub, Normal(means[idx], sqrt(vars[idx])))
            idx = idx + 1
        end
        push!(inter_distr, inter_dist_sub)
    end

    return [tuple_ for tuple_ in zip(inter_variables, inter_distr)]
end

"""
    inter_noise_selection(inter_variables::Array{Array{Int64,1},1} ;
                                    additive = true, prob_cst = 1/3, lb = 0, ub = 1)

Generate the distributions of the shift/rescale variables (i.e. A or C) for mechanism changes. Each of these distr. is either uniform or constant. Return a list of tuples (one per interventional environments) of the same type as `noise_change` in `mech_changes_sim`.

* `inter_variables`   : output of `inter_var_selection`, i.e. array of arrays containing the variables to intervene on in each env.
* `additive`          : whether the noise is shifted (true) or rescaled (false)
* `prob_cst`          : whether the shift/rescale variables are constant or not
* `lb`                : lower bound of interval where the bounds for the uniform distr. are randomly chosen (uniformly)
* `ub`                : upper bound... Note: if additive, we randomly flip the interval around zero, else we flip around 1
* `flip`              : whether we randomly flip the intervals as mentionned above
"""
function inter_noise_selection(inter_variables::Array{Array{Int64,1},1} ;
                                    additive = true, prob_cst = 1/3, lb = 0, ub = 1, flip = true)
    # init
    @assert 0 <= prob_cst <= 1 && !(!additive && flip) || 0 < lb <= ub
    E = length(inter_variables)
    nSim = sum(length.(inter_variables))
    inter_distr::Array{Array{UnivariateDistribution,1},1} = []

    # draw params
    bds_1 = rand(nSim) .* (ub - lb) .+ lb
    bds_2 = rand(nSim) .* (ub - lb) .+ lb
    cst_probs = rand(nSim)
    switches = flip ? rand([0,1], nSim) : zeros(nSim)
    lower_bds = min.(bds_1, bds_2) .+ switches .* (.- min.(bds_1, bds_2) .+ (additive ? .- max.(bds_1, bds_2) : max.(bds_1, bds_2).^(-1)))
    upper_bds = max.(bds_1, bds_2) .+ switches .* (.- max.(bds_1, bds_2) .+ (additive ? .- min.(bds_1, bds_2) : min.(bds_1, bds_2).^(-1)))

    # construct array of distr
    idx = 1
    for e in 1:E
        inter_dist_sub::Array{UnivariateDistribution ,1} = []
        for var in inter_variables[e]
            noise_distr = Normal(0,1)
            if cst_probs[idx] <= prob_cst
                noise_distr = Normal((lower_bds[idx] + upper_bds[idx])/2,0)
            else
                noise_distr = Uniform(lower_bds[idx], upper_bds[idx])
            end
            push!(inter_dist_sub, noise_distr)
            idx = idx + 1
        end
        push!(inter_distr, inter_dist_sub)
    end

    return [tuple_ for tuple_ in zip(inter_variables, inter_distr)]
end

function simple_noise_selection(inter_variables::Array{Array{Int64,1},1}, d=TDist(4))
    meanshift = [0,0.1,0.2,0.5,1,2,5,10]
    strength = [0.1,0.2,0.5,1,2,5,10]
    E = length(inter_variables)
    inter_distr = Vector{Vector{<:LocationScale{Float64, <:UnivariateDistribution}}}(undef, E) #TDist{Float64}
    for e in 1:E
        num_int = length(inter_variables[e])
        tmp_dist::LocationScale{Float64, <:UnivariateDistribution} = rand(strength)*d + rand(meanshift)
        inter_dist_sub::Vector{<:LocationScale{Float64, <:UnivariateDistribution}} = [tmp_dist]#(undef, num_int) # TDist{Float64}
        for i in 2:num_int
#             inter_dist_sub[i] = rand(strength)*d + rand(meanshift)
            push!(inter_dist_sub, rand(strength)*d + rand(meanshift))
        end
        inter_distr[e] =  inter_dist_sub
    end
    return [tuple_ for tuple_ in zip(inter_variables, inter_distr)]
end

"""
    inter_coeff_selection(scm::LinAddSCM, inter_variables::Array{Array{Int64,1},1} ;
                                prob_unchanged = 2/3, redraw = false, lb=-2, ub=2, variance = 1)

Generate the new coeefficents from mechanism changes. It eithers generate completely new coefficent from scratch or just add some 'noise' to the existing coeff (this 'noise' is here constant when generating samples). With a certain probability it will not change the coeff. Return a list of tuples (one per interventional environments) of the same type as `coeff_change` in `mech_changes_sim`.

* `scm`               : the linear structural model
* `inter_variables`   : output of `inter_var_selection`, i.e. array of arrays containing the variables to intervene on in each env.
* `prob_unchanged`    : probability that a coeff is eventually unchanged
* `redraw`            : whether we redraw totally new coeff (true) or change the existing coeff by adding some noise (false)
* `lb`                : lower bound of uniform distr used to draw new coeff (when redraw = true). We also randomly flip signs after.
* `ub`                : upper bound of uniform distr used to draw new coeff.
* `variance`          : the noise term added to the coeff when redraw = false is drawn from N(0, `variance` )
"""
function inter_coeff_selection(scm::LinAddSCM, inter_variables::Array{Array{Int64,1},1} ;
                                prob_unchanged = 2/3, redraw = false, lb=-2, ub=2, variance = 1)
    @assert 0 <= prob_unchanged <= 1
    E = length(inter_variables)
    nSim = sum(map(x -> sum(scm.B[x,:] .!= 0), inter_variables))
    coeff_mat::Array{Array{Float64, 2}}= []

    # draw params
    unchange_probs = rand(nSim)
    coeff_changes = redraw ? (rand(nSim) .* (ub - lb) .+ lb) .* rand([-1,1],nSim) : rand(Normal(0,variance), nSim)

    # construct array of distr
    idx = 1
    for e in 1:E
        coeff_mat_sub::Array{Float64, 2} = zeros(0,scm.p)
        for var in inter_variables[e]
            vect = scm.B[var, :]
            for j in filter(i -> vect[i] != 0, 1:scm.p)
                if unchange_probs[idx] > prob_unchanged
                    vect[j] = redraw ? coeff_changes[idx] : vect[j] + coeff_changes[idx]
                end
                idx = idx +1
            end
            coeff_mat_sub = vcat(coeff_mat_sub, vect')
        end
        push!(coeff_mat, coeff_mat_sub)
    end

    return [tuple_ for tuple_ in zip(inter_variables, coeff_mat)]
end

# struct for linear structural models together with all its modifications under different interventional evironments and the information needed to generate samples from under both observational and interventional settings.
struct InterLinAddSCM
    scm::LinAddSCM
    is_mech_change::Bool
    inter_vect_list::Array{Tuple{Vector{Int64}, Vector{UnivariateDistribution }},1} # for stochastic_inter_sim
    noise_change_list::Array{Tuple{Vector{Int64}, Vector{UnivariateDistribution }},1}  # for mech_changes_sim
    coeff_change_list::Array{Tuple{Array{Int64,1},Array{Float64,2}},1} # for mech_changes_sim
    E # number of environments
    target # target variable
    additive # when the intervention is a mechanism change, if true -> shifted noise, otherwise rescaled by a random variable
    function InterLinAddSCM(scm::LinAddSCM, is_mech_change::Bool, target::Int64; additive = true,
                                inter_vect_list =[], #::Array{<:Tuple{Vector{Int64}, Vector{<:UnivariateDistribution}},1} = [],
                                noise_change_list = [],#::Array{<:Tuple{Vector{Int64}, Vector{<:UnivariateDistribution}},1} = [],
                                coeff_change_list = [])#::Array{<:Tuple{Array{Int64,1}, Array{Float64,2}},1} = [])
        @assert (!is_mech_change || length(noise_change_list) == length(coeff_change_list)) && 1 <= target <= scm.p
        E = is_mech_change ? length(noise_change_list) + 1 : length(inter_vect_list) + 1
        new(scm, is_mech_change, inter_vect_list, noise_change_list, coeff_change_list, E, target,additive)
    end
end

# struct for non-linear structural models together with all its modifications under different interventional evironments and the information needed to generate samples from under both observational and interventional settings.
struct InterNonLinAddSCM
    scm::NonLinAddSCM
    is_mech_change::Bool
    inter_vect_list#::Array{Tuple{Vector{Int64}, Vector{T}},1} where T <:LocationScale{Float64, TDist{Float64}}# for stochastic_inter_sim
    noise_change_list#::Array{Tuple{Vector{Int64}, Vector{T}},1} where T <:LocationScale{Float64, TDist{Float64}}# for mech_changes_sim
    E # number of environments
    target # target variable
    additive # when the intervention is a mechanism change, if true -> shifted noise, otherwise rescaled by a random variable
    function InterNonLinAddSCM(scm::NonLinAddSCM, is_mech_change::Bool, target::Int64; additive = true,
                                inter_vect_list =[], #::Array{<:Tuple{Vector{Int64}, Vector{<:UnivariateDistribution}},1} = [],
                                noise_change_list = [])
        @assert 1 <= target <= scm.p
        E = is_mech_change ? length(noise_change_list) + 1 : length(inter_vect_list) + 1
        new(scm, is_mech_change, inter_vect_list, noise_change_list, E, target,additive)
    end
end

"""
    generate_interLinAddSCM(is_mech_change::Bool, E::Int64, share::Float64; scm = nothing,
        high_dim_SCM_params = Dict("p_anterior" => 50, "p_posterior" => 50, "k" => 20, "parents_nbr" => 20, "lb" => -2, "ub" => 2, "var_min"    => 0.5^2, "var_max" => 2^2), overlap = 0.2, additive = true, inter_noise_selection_params = Dict("prob_cst" => 1/3, "lb" => 0, "ub" => 1),
        inter_coeff_selection_params = Dict("prob_unchanged" => 2/3, "redraw" => false, "lb" => -2, "ub" => 2, "variance" => 1),
        inter_distr_selection_params = Dict("lb" => -2, "ub" => 2, "var_min" => 0.5^2, "var_max" => 2^2))

Generate from scratch a new linear structural model and different interventional environment from it. Returns an InterLinAddSCM object. Interventional environments can only be either fully made of stochastic interventions or fully made of mechanism changes. The first environment is always considered to be purely observational.

* `is_mech_change`                : whether the interventions are mechanism changes (true) or not (false)
* `E`                             : number of environments, E >= 1
* `share`                         : between 0 and 1, share of the non-target variables that are intervened on at least once
* `scm`                           : can provide an SCM instead of generating one (generate one if `scm` = nothing)
* `target`                        : target variable (only when providing an SCM in `scm` field)
* `high_dim_SCM_params`           : parameters for the function `high_dim_SCM`
* `overlap`                       : between 0 and 1, overlap of inter. variables between two env, see `inter_var_selection`
* `additive`                      : for mech change, whether the noise change is additive (true) or multiplicative (false)
* `inter_noise_selection_params`  : parameters for the function `inter_noise_selection`
* `inter_coeff_selection_params`  : parameters for the function `inter_coeff_selection`
* `inter_distr_selection_params`  : parameters for the function `inter_distr_selection`
"""
function generate_interLinAddSCM(is_mech_change::Bool, E::Int64, share::Float64; scm = nothing, target = -1,
        high_dim_SCM_params = Dict("p_anterior" => 50, "p_posterior" => 50, "k" => 20, "parents_nbr" => 20, "lb" => -2, "ub" => 2,
        "var_min" => 0.5^2, "var_max" => 2^2), overlap = 0.2, additive = true,
        inter_noise_selection_params = Dict("prob_cst" => 1/3, "lb" => 0, "ub" => 1, "flip" => true),
        inter_coeff_selection_params = Dict("prob_unchanged" => 2/3, "redraw" => false, "lb" => -2, "ub" => 2, "variance" => 1),
        inter_distr_selection_params = Dict("lb" => -2, "ub" => 2, "var_min" => 0.5^2, "var_max" => 2^2))
    # check
    @assert E >= 1 && 0 <= share <= 1

    # create the linear SCM if it is not already given
    if scm == nothing
        hdsp = high_dim_SCM_params
        p_anterior, p_posterior, k, parents_nbr = hdsp["p_anterior"], hdsp["p_posterior"], hdsp["k"], hdsp["parents_nbr"]
        lb, ub, var_min, var_max = hdsp["lb"], hdsp["ub"], hdsp["var_min"], hdsp["var_max"]
        scm = high_dim_SCM(p_anterior,p_posterior,k,parents_nbr; lb=lb, ub=ub, var_min=var_min, var_max=var_max)
        target = p_anterior + 1
    else
        @assert typeof(scm) == LinAddSCM && 1 <= target <= scm.p
    end

    # create the interventional envs if needed
    if E > 1
        # generate variables to intervene on by (interventional) environments
        inter_variables = inter_var_selection(scm, target, E, share; overlap = overlap)

        # generate now the params for each interventional env, and return a InterLinAddSCM
        if is_mech_change
            # generate the noise scaling/shifting variables per intervention
            insp = inter_noise_selection_params
            prob_cst, lb, ub, flip = insp["prob_cst"], insp["lb"], insp["ub"], insp["flip"]
            noise_change_list = inter_noise_selection(inter_variables;
                                    additive = additive, prob_cst = prob_cst, lb = lb, ub = ub, flip = flip)

            # generate the potential coefficient changes per intervention
            icsp = inter_coeff_selection_params
            prob_unchanged, redraw, lb, ub, variance = icsp["prob_unchanged"], icsp["redraw"], icsp["lb"], icsp["ub"], icsp["variance"]
            coeff_change_list = inter_coeff_selection(scm, inter_variables;
                                    prob_unchanged = prob_unchanged, redraw = redraw, lb=lb, ub=lb, variance = variance)

            return InterLinAddSCM(scm,is_mech_change,target;
                                    additive = additive, noise_change_list = noise_change_list, coeff_change_list = coeff_change_list)
        else
            # generate the do-intervention noise distributions
            idsp = inter_distr_selection_params
            lb, ub, var_min, var_max = idsp["lb"], idsp["ub"], idsp["var_min"], idsp["var_max"]
            inter_vect_list = inter_distr_selection(inter_variables; lb=lb, ub=ub, var_min=var_min, var_max=var_max)

            return InterLinAddSCM(scm,is_mech_change,target; additive = additive, inter_vect_list = inter_vect_list)
        end
    else
        return InterLinAddSCM(scm,is_mech_change,target)
    end
end

"""
    simulate(inter_scm::InterLinAddSCM, nₑs, d::UnivariateDistribution)

Simulate observational and interventional data from an interventional linear SCM object, i.e. InterLinAddSCM.

* `inter_scm`   : the InterLinAddSCM object
* `nₑs`         : array of the number of samples from environment
* `ds`          : vector of (observational) distributions of the noises, assumed to be of mean = 0 and variance = 1, e.g. ds=[Gamma(), Exponential(), ...]
"""
function simulate(inter_scm::InterLinAddSCM, nₑs, ds::Vector{T}) where T <: UnivariateDistribution
    # check & init
    E = inter_scm.E
    @assert length(nₑs) == E
    p = inter_scm.scm.p
    target = inter_scm.target
    is_mech_change = inter_scm.is_mech_change
    scm = inter_scm.scm
    additive = inter_scm.additive

    # simulate data
    obs_noise_distr = rand(ds, p)

    # observational data
    sim_data = simulate(scm, obs_noise_distr, nₑs[1])

    # interventional data
    for e in 1:(E-1)
        if is_mech_change
            noise_change = inter_scm.noise_change_list[e]
            coeff_change = inter_scm.coeff_change_list[e]
            sim_data = vcat(sim_data, mech_changes_sim(scm, obs_noise_distr, noise_change, coeff_change, nₑs[e+1]; additive = additive))
        else
            inter_vect = inter_scm.inter_vect_list[e]
            sim_data = vcat(sim_data, stochastic_inter_sim(scm, obs_noise_distr, inter_vect, nₑs[e+1]))
        end
    end

    # separate target (y) from the rest (X)
    X, y = sim_data[:,1:p .!= target], sim_data[:,target]

    return (X,y)
end

"""
    simulate(inter_scm::InterNonLinAddSCM, nₑs, d::UnivariateDistribution)

Simulate observational and interventional data from an interventional non-linear SCM object, i.e. InterNonLinAddSCM.

* `inter_scm`   : the InterNonLinAddSCM object
* `nₑs`         : array of the number of samples from environment
* `ds`          : vector of (observational) distributions of the noises, assumed to be of mean = 0 and variance = 1, e.g. ds=[Gamma(), Exponential(), ...]
"""
function simulate(inter_scm::InterNonLinAddSCM, nₑs, ds::Vector{T}) where T <: UnivariateDistribution
    # check & init
    E = inter_scm.E
    @assert length(nₑs) == E
    p = inter_scm.scm.p
    target = inter_scm.target
    is_mech_change = inter_scm.is_mech_change
    scm = inter_scm.scm
    additive = inter_scm.additive

    # simulate data
    obs_noise_distr = rand(ds, p)

    # observational data
    sim_data = simulate(scm, obs_noise_distr, nₑs[1])

    # interventional data
    for e in 1:(E-1)
        if is_mech_change
            noise_change = inter_scm.noise_change_list[e]
            sim_data = vcat(sim_data, mech_changes_sim(scm, obs_noise_distr, noise_change, nₑs[e+1]; additive = additive))
        else
            inter_vect = inter_scm.inter_vect_list[e]
            sim_data = vcat(sim_data, stochastic_inter_sim(scm, obs_noise_distr, inter_vect, nₑs[e+1], scm.fs))
        end
    end

    # separate target (y) from the rest (X)
    X, y = sim_data[:,1:p .!= target], sim_data[:,target]

    return (X,y)
end


text:
using NNlib
import NNlib: softmax, ∇softmax, logsoftmax, ∇logsoftmax, conv, depthwiseconv, ∇conv_data, ∇depthwiseconv_data, maxpool, meanpool, σ, relu, batched_mul, batched_adjoint

drelu(x, Δ) = ifelse(x > 0, Δ, zero(x))

@adjoint function Base.Broadcast.broadcasted(::typeof(relu), x::Numeric)
  relu.(x), Δ -> (nothing, drelu.(x, Δ))
end

function dselu(x)
  λ = oftype(x/1, 1.0507009873554804934193349852946)
  α = oftype(x/1, 1.6732632423543772848170429916717)
  λ * ifelse(x > 0, 1, α * exp(x))
end

@adjoint function Base.Broadcast.broadcasted(::typeof(selu), x::Numeric)
  selu.(x), Δ -> (nothing, dselu.(x) .* Δ)
end

@adjoint function σ(x::Real)
    y = σ(x)
    return y, Δ -> (Δ * y * (1 - y),)
end

@adjoint softmax(xs; dims=1) = softmax(xs, dims=dims), Δ -> (∇softmax(Δ, xs, dims=dims),)

@adjoint logsoftmax(xs; dims=1) = logsoftmax(xs, dims=dims), Δ -> (∇logsoftmax(Δ, xs, dims=dims),)

@adjoint NNlib.DenseConvDims(args...; kwargs...) = NNlib.DenseConvDims(args...; kwargs...), _ -> nothing
@adjoint NNlib.DepthwiseConvDims(args...; kwargs...) = NNlib.DepthwiseConvDims(args...; kwargs...), _ -> nothing
@adjoint NNlib.PoolDims(args...; kwargs...) = NNlib.PoolDims(args...; kwargs...), _ -> nothing

colmajor(x) = colmajor(MemoryLayout(typeof(x)), x)
colmajor(_, x) = convert(Array, x)
colmajor(::AbstractColumnMajor, x) = x

@adjoint conv(x, w, cdims; kw...) =
  conv(x, w, cdims; kw...),
    Δ -> begin
       Δ = colmajor(Δ)
       return (
           NNlib.∇conv_data(Δ, w, cdims; kw...),
           NNlib.∇conv_filter(x, Δ, cdims; kw...),
           nothing,
       )
   end

@adjoint ∇conv_data(x, w, cdims; kw...) =
  ∇conv_data(x, w, cdims; kw...),
    Δ -> begin
       Δ = colmajor(Δ)
       return (
           NNlib.conv(Δ, w, cdims; kw...),
           NNlib.∇conv_filter(Δ, x, cdims; kw...),
           nothing,
       )
   end

@adjoint depthwiseconv(x, w, cdims; kw...) =
  depthwiseconv(x, w, cdims; kw...),
    Δ -> begin
       Δ = colmajor(Δ)
       return (
           NNlib.∇depthwiseconv_data(Δ, w, cdims; kw...),
           NNlib.∇depthwiseconv_filter(x, Δ, cdims; kw...),
           nothing,
       )
   end

@adjoint ∇depthwiseconv_data(x, w, cdims; kw...) =
  ∇depthwiseconv_data(x, w, cdims; kw...),
    Δ -> begin
       Δ = colmajor(Δ)
       return (
           NNlib.depthwiseconv(Δ, w, cdims; kw...),
           NNlib.∇depthwiseconv_filter(Δ, x, cdims; kw...),
           nothing,
       )
   end

@adjoint function maxpool(x, pdims::NNlib.PoolDims; kw...)
  y = maxpool(x, pdims; kw...)
  y, Δ -> (NNlib.∇maxpool(Δ, y, x, pdims; kw...), nothing)
end

@adjoint function meanpool(x, pdims::NNlib.PoolDims; kw...)
  y = meanpool(x, pdims; kw...)
  y, Δ -> (NNlib.∇meanpool(Δ, y, x, pdims; kw...), nothing)
end

@adjoint function batched_mul(A, B)
    C = batched_mul(A, B)
    C, Δ -> (batched_mul(Δ, batched_adjoint(B)), batched_mul(batched_adjoint(A), Δ))
end

ext: jl


text:
Cat = @theory begin
    1(A) => id(A).(A→A)
    id(A).(A→A) ⋅ f.(A→B) => f.(A→B)
    f.(A→B) ⋅ id(A).(A→A) => f.(A→B)
    f.(A→B)⋅id(B).(B→B) => f.(A→B)
    (f.(A→B)⋅g.(B→C))⋅h.(C→D) => f.(A→B)⋅(g.(B→C)⋅h.(C→D))
end

using Test

@testset "Cat" begin
@test @areequal Cat 1(B) id(B).(B→B)
@test @areequal Cat (id(A).(A→A)⋅f.(A→B)) f.(A→B)
@test @areequal Cat (1(A)⋅f.(A→B)) f.(A→B)
@test @areequal Cat (f.(A→B)⋅1(B)) f.(A→B)
@test @areequal Cat (f.(A→B)⋅g.(B→B))⋅g.(B→B) f.(A→B)⋅(g.(B→B)⋅g.(B→B))
@test @areequal Cat (f.(A→B)⋅g.(B→B))⋅h.(B→C) f.(A→B)⋅(g.(B→B)⋅h.(B→C))
# areequal(CatTyped, :(f.(A→B)⋅1(B)), :f; timeout=100)
end

MonCat = @theory begin
    1(A) => id(A).(A→A)
    id(A).(A→A) ⋅ f.(A→B) => f.(A→B)
    f.(A→B) ⋅ id(A).(A→A) => f.(A→B)
    f.(A→B)⋅id(B).(B→B) => f.(A→B)
    (f.(A→B)⋅g.(B→C))⋅h.(C→D) => f.(A→B)⋅(g.(B→C)⋅h.(C→D))

    munit() => munit()
    A⊗munit() => A
    munit()⊗A => A

    (f⊗g)⊗h => f⊗(g⊗h)  # associativity for both obs and homs
    f.(A→B)⊗g.(C→D) => (f⊗g).(A⊗C→B⊗D)
    id(A⊗B).(A⊗B→A⊗B) => id(A).(A→A)⊗id(B).(B→B)
end

@testset "MonCat" begin
    @test @areequal MonCat 1(A)⊗1(B) 1(A⊗B)
    @test @areequal MonCat (f.(A→B)⊗g.(B→C)) (f.(A→B)⊗g.(B→C))
    @test @areequal MonCat 1(A⊗B)⋅(f.(A→B)⊗g.(B→C)) (f.(A→B)⊗g.(B→C))
    @test @areequal MonCat (f.(A→B)⊗g.(B→C))⋅ 1(A⊗B) (f.(A→B)⊗g.(B→C))
    @test @areequal MonCat 1(A⊗B)⋅(f.(A→B)⊗g.(B→C)) (f⊗g).(A⊗B→B⊗C)
end


text:
#Plot basic Gaussian Process

using Gadfly, GaussianProcesses

#Training data
d, n = 2, 50         #Dimension and number of observations

x = 2π * rand(d, n)
y = vec(sin(x[1,:]).*sin(x[2,:])) + 0.05*rand(n)

mZero = MeanZero()                    #Zero mean function
#kern = Mat(5/2,[0.0,0.0],0.0) + SE(0.0,0.0)        #Matern 5/2 ARD kernel with parameters log(l₁) = 0 log(l₂) = 0 and log(σ) = 0 and SE Iso kernel with parameters log(ℓ) = 0 and log(σ) = 0

kern = SE(0.0, 0.0)

gp = GP(x,y,mZero,kern,-2.0)          # Fit the GP
optimize!(gp)                         # Optimize the hyperparameters

push!(gp, [[2.0; 3.0] [3.0, 4.0]], [2.0, 3.0])
push!(gp, [2.0; 3.0], 2.0)
plot(gp; clim=(-10.0, 10.0,-10.0,10.0)) #Plot the GP over range clim


text:
#RSDeltaSigmaPort display facilities
#-------------------------------------------------------------------------------

const PlotOrColl = Union{EasyPlot.Plot, EasyPlot.PlotCollection}

struct ImageRepr{T}
	obj::T #Object to display as image
	width::Float64 #Width of resultant image
	AR::Float64 #Aspect ratio
end

"""`ImageRepr{T}(obj::T; width=640, AR=1)`

Wrapper object requesting an image representation of something using a specific
width and aspect ratio.
"""
ImageRepr(obj; width=640, AR=1) = ImageRepr(obj, Float64(width), Float64(AR))

#==helper functions
===============================================================================#
#Convert EasyPlot.Plot -> EasyPlot.PlotCollection
_plt2coll(plot::EasyPlot.Plot) = push!(cons(:plot_collection), plot)

#=="show" interface
===============================================================================#
function _show(io::IO, mime::MIME, r::ImageRepr{EasyPlot.PlotCollection})
	w=round(Int, r.width); h = round(Int, w/r.AR)
	set = EasyPlot.set
	opt = EasyPlot.ShowOptions(dim=set(w=w, h=h))
	b = EasyPlot.getbuilder(:image, :InspectDR)
	nativeplot = EasyPlot.build(b, r.obj)
	EasyPlot._show(io, mime, opt, nativeplot)
	return nothing
end
function _show(io::IO, mime::MIME, r::ImageRepr{EasyPlot.Plot})
	pcoll = _plt2coll(r.obj)
	return _show(io, mime, ImageRepr(pcoll, width=r.width, AR=r.AR))
end

#Only publicly provide Base.show() for MIME"image/png":
Base.show(io::IO, mime::MIME"image/png", r::ImageRepr{T}) where T<:PlotOrColl = _show(io, mime, r)

function Base.show(io::IO, mime::MIME"text/plain", r::ImageRepr{EasyPlot.PlotCollection})
	#Do not warn. Jupyter calls function even if it doesn't use the results.
	#Might it be generating text for a hover-pop-up or something?
#	@warn("Trying to display inline plot on simple text display.")
	println(io, "[INLINE PLOT]: ", r.obj.title)
	return nothing
end
function Base.show(io::IO, mime::MIME"text/plain", r::ImageRepr{EasyPlot.Plot})
	println(io, "[INLINE PLOT]: ", r.obj.title)
	return nothing
end

#=="showable" interface
===============================================================================#
#Base.showable(mime::MIME"image/png", p::ImageRepr) = true

#==User-facing interface
===============================================================================#
"""`inlinedisp(plot; AR=1, scalew=1, maxw=900)`

Display plot as inline image.

# Arguments
 - `AR`: Aspect ratio (w/h) of generated image.
 - `maxw`: Reference width of image (typ. max width of display).
 - `scalew`: Amount by which to scale image before rendering it (`imagew=maxw*scalew`)
"""
inlinedisp(plot; AR=1, scalew=1, maxw=900) =
	display(ImageRepr(plot; AR=AR, width=maxw*scalew))

_write(mime::Symbol, filepath::String, nativeplot; opt_kwargs...) =
	_write(mime, filepath, ShowOptions(; opt_kwargs...), nativeplot)

"""`saveimage(mime::Symbol, filepath::String, plot; AR=1, width=900)`

Saves plot as an image.

# Arguments
 - `AR`: Aspect ratio (w/h) of generated image.
 - `maxw`: Reference width of image (typ. max width of display).
 - `scalew`: Amount by which to scale image before rendering it (`imagew=maxw*scalew`)
"""
function saveimage(mime::Symbol, filepath::String, plot::PlotOrColl; AR=1, width=900)
	ir = ImageRepr(plot; AR=AR, width=width)
	open(filepath, "w") do io
		_show(io, EasyPlot._getmime(mime), ir)
	end
end

function displaygui(plot::PlotOrColl)
	result = EasyPlot.displaygui(:InspectDR, plot)
	sleep(0.1) #Get GUI to refresh (Ideally: `doevents()`, or something)
	#NOTE: Appears sleep needs to be > 0 to refresh.
	return result
end

#Last line


text:
# sample a multinomial(n,p) r.v. using binomial r.v.s
@inline function sampleMultinomial!(n::Int64, p::Vector{Float64},
  result::Vector{Int64}, rng::RNG = GLOBAL_RNG) where RNG <: AbstractRNG
  nleft::Int64 = n
  tmp::Float64 = 0.0
  for i = 1:length(p)-1
    @inbounds result[i] = sampleBinomial(nleft, min(1.0, p[i] / (1.0 - tmp)), rng)
    @inbounds nleft -= result[i]
    @inbounds tmp += p[i]
  end
  @inbounds result[length(p)] = nleft;
  return
end

@inline function sampleMultinomial(n::Int64, p::Vector{Float64},
  rng::RNG = GLOBAL_RNG) where RNG <: AbstractRNG
  result::Vector{Int64} = Vector{Int64}(undef, length(p))
  sampleMultinomial!(n, p, result, rng)
  return result
end


text:
function get_body(obj::CelestialBodies, arg0::JString)
    return jcall(obj, "getBody", CelestialBody, (JString,), arg0)
end

function get_earth(obj::CelestialBodies)
    return jcall(obj, "getEarth", CelestialBody, ())
end

function get_earth_moon_barycenter(obj::CelestialBodies)
    return jcall(obj, "getEarthMoonBarycenter", CelestialBody, ())
end

function get_jupiter(obj::CelestialBodies)
    return jcall(obj, "getJupiter", CelestialBody, ())
end

function get_mars(obj::CelestialBodies)
    return jcall(obj, "getMars", CelestialBody, ())
end

function get_mercury(obj::CelestialBodies)
    return jcall(obj, "getMercury", CelestialBody, ())
end

function get_moon(obj::CelestialBodies)
    return jcall(obj, "getMoon", CelestialBody, ())
end

function get_neptune(obj::CelestialBodies)
    return jcall(obj, "getNeptune", CelestialBody, ())
end

function get_pluto(obj::CelestialBodies)
    return jcall(obj, "getPluto", CelestialBody, ())
end

function get_saturn(obj::CelestialBodies)
    return jcall(obj, "getSaturn", CelestialBody, ())
end

function get_solar_system_barycenter(obj::CelestialBodies)
    return jcall(obj, "getSolarSystemBarycenter", CelestialBody, ())
end

function get_sun(obj::CelestialBodies)
    return jcall(obj, "getSun", CelestialBody, ())
end

function get_uranus(obj::CelestialBodies)
    return jcall(obj, "getUranus", CelestialBody, ())
end

function get_venus(obj::CelestialBodies)
    return jcall(obj, "getVenus", CelestialBody, ())
end


text:
struct FiniteDifferenceBVPProblem{TF, TBCL, TBCR}
    a # (a,b) define the domain
    b
    Δx # mesh spacing
    n # number of points in the interior, excluding a and b
    x # interior mesh
    f::TF # value of f at the interior nodes
    A # the matrix
    rhs # right hand side of our system
    left_bc::TBCL
    right_bc::TBCR
end

struct DirichletBC
    u_bc
end

struct NeumannBC
    v_bc
end

struct PeriodicBC
end

struct SpectralBVPProblem{TF, TBCL, TBCR}
    a # (a,b) define the domain
    b
    Δx # mesh spacing
    n # number of points in the interior, excluding a and b
    N # number of wave numbers/number of points on [a,b)
    x # interior mesh
    f::TF # value of f at the interior nodes
    λ # eivenalues
    rhs # right hand side of our system
    left_bc::TBCL
    right_bc::TBCR
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
# # Sampling from grid approximate posterior
# using StatsBase
# using Distributions
# using Plots
# using StatsPlots

# %% 2.1
ways = [0, 3, 8, 9, 0]
ways = ways / sum(ways)

# %% 2.2
using Distributions
"""
With Distributions.jl working with distributions is a little different that with R. Instead of having `rbinom`, `dbinom` and `pbinom` we just got the `Binomial` distribution which we can work with.
d = Binomial(9, 0.5)
rand(d)     # singe random draw
rand(d, 10) # 10 random draws
pdf(d, 6)   # probability of getting 6
cdf(d, 6)   # cumulative probability of getting 6
Take a look at Distributions.jl docs to see what else you can do with distributions. https://juliastats.org/Distributions.jl/latest/index.html
"""

pdf(Binomial(9, 0.5), 6)

# %% 2.3
p_grid = range(0, 1, length = 20)
prior = ones(20)
likelihood = pdf.(Binomial.(9, p_grid), 6)
posterior = likelihood .* prior
posterior ./= sum(posterior)
"""
This is the next difference from Julia to R. You don't have automatic broadcasting (making functions work the same whether you run it on a single number or on an array of numbers). However, fixing this is very easy, you only have to annotate your function call with a dot: `f(single_number)` -> `f.(array_of_numbers)`.
That means the translation of the line `likelihood = pdf.(Binomial.(9, p_grid), 6)` goes something like: first, for every value in `p_grid` make a binomial distribution with that p value. For every distribution then take the pdf of 6. The same in the next line, this is elementwise multiplication.

Related advide: when you get an error of the kind "no method matching ...(::Float64, ::Array{Float64})" or similar, first check if all your dots are in order.
For more infos read https://julialang.org/blog/2017/01/moredots/
"""

# %% 2.4
using Plots
plot(p_grid, posterior, m = 3, legend = false,
    xlabel = "probability of water", ylabel = "posterior probability", title = "20 points")

# %% 2.5
prior = ifelse.(p_grid .< 0.5, 0, 1)
prior = @. exp(-5 * abs(p_grid - 0.5))  # @. means broadcast everything that follows

# %% 2.6
using Turing
include(srcdir("quap.jl"))
@model globethrowing(W, L) = begin
    p ~ Uniform(0, 1)
    W ~ Binomial(W + L, p)
end
m = globethrowing(6, 3)
r = quap(m)

d = MvNormal(r.coef, collect(reshape(r.vcov, 1, 1)))

s = rand(d, 10000)'

histogram!(collect(s), normalize = :probability)
r.vcov[1] |> sqrt

plot!(p_grid, posterior ./ 2.5, m = 3)

histogram!(rand(d, 10000), normalize = :probability)

https://hdfilme.cx/filme1/mara-und-der-feuerbringer-2015-hd-465-stream/deutsch


text:
using WordNet
using Test

function io_to_string(f::Function)
    buff = IOBuffer()
    f(buff)
    String(take!(buff))
end

include("test_pointer.jl")
include("test_db.jl")
include("test_lemma.jl")
include("test_synset.jl")
include("test_operations.jl")
include("test_sensekeys.jl")


text:
# Non-transducer reducing function combinators

# Some common definitions for TeeRF and ProductRF
abstract type AbstractMultiCastingRF{N} <: _Function end
# TOOD: Use `AbstractMultiCastingRF <: AbstractReduction` so that
# extra `BottomRF` wrapper can be removed.

start(rf::AbstractMultiCastingRF, init) = map(f -> start(f, init), rf.fs)

start(rf::AbstractMultiCastingRF{N}, init::NTuple{N,Any}) where {N} =
    map(start, rf.fs, init)
start(rf::AbstractMultiCastingRF, init::Tuple) = error(
    "`start($rf, init)` received an incompatible accumulator `init` of length ",
    length(init),
    ".\n",
    "It requires a $(length(rf.fs))-tuple instead.\n",
    "init = $init",
)

# Trying to detect bugs by narrowing method signatures when defining
# `next` for concrete reducing functions. Otherwise, it falls back to
# this method; e.g., foldl(TeeRF(min, max), [5, 2, 6, 8, 3])
next(rf::AbstractMultiCastingRF, accs, x) = error(
    "$rf received an incompatible accumulator of type $(typeof(accs))\n",
    "It requires a $(length(rf.fs))-tuple.\n",
    "Accumulator (first argument): $(accs)\n",
    "Input (second argument): $x",
)

@inline function complete(rf::AbstractMultiCastingRF, accs)
    results = map((f, a) -> complete(f, a), rf.fs, accs)
    if Baselet.Specialized.any(r -> r isa DefaultInitOf, results)  # is `any` too strict?
        return DefaultInitOf{typeof(rf)}()
    end
    return results
end

@inline completebasecase(rf::AbstractMultiCastingRF, accs) =
    map((f, a) -> completebasecase(f, a), rf.fs, accs)

@inline combine(rf::AbstractMultiCastingRF, lefts, rights) =
    map((f, l, r) -> combine(f, l, r), rf.fs, lefts, rights)

# # Something like these are needed for `AbstractMultiCastingRF <: AbstractReduction`:
# Transducer(::AbstractMultiCastingRF) = IdentityTransducer()
# has(rf::AbstractMultiCastingRF, ::Type{<:Transducer}) = false  # should be multicast?
# _realbottomrf(rf::AbstractMultiCastingRF) = rf
# _show_impl(io, m, rf::AbstractMultiCastingRF) = invoke(Base.show, Tuple{IO,Any}, io, rf)

# Adjoining init has to happen at "leafs":
_asmonoid(rf::AbstractMultiCastingRF) =
    default_constructorof(typeof(rf))(map(_asmonoid, rf.fs)::Tuple)

# See: [[./core.jl::TODOs for `Completing` and `skipcomplete`]]
Completing(rf::AbstractMultiCastingRF) =
    default_constructorof(typeof(rf))(map(Completing, rf.fs)::Tuple)

function Base.show(io::IO, rf::AbstractMultiCastingRF)
    @nospecialize
    T = default_constructorof(typeof(rf))
    if rf === T(rf.fs)
        print(io, T, rf.fs)
    else
        invoke(show, Tuple{IO,Any}, io, rf)
    end
end

"""
    TeeRF(reducing_functions::Tuple)
    TeeRF(reducing_functions...)

Combine multiple reducing functions into a new reducing function that
"multicast" the input to multiple reducing functions.

Roughly speaking, `TeeRF(op₁, op₂, ..., opₙ)` is equivalent to

```julia
((a₁, a₂, ..., aₙ), x) -> (op₁(a₁, x), op₂(a₂, x), ..., opₙ(aₙ, x))
```

For [`combine`](@ref), it behaves like [`ProductRF`](@ref).

!!! compat "Transducers.jl 0.4.32"

    New in version 0.4.32.

# Examples
```jldoctest
julia> using Transducers

julia> extrema′(xs, xf = Map(identity)) = foldl(TeeRF(min, max), xf, xs);

julia> extrema′([5, 2, 6, 8, 3])
(2, 8)
```

Note that the input is considered empty unless _all_ reducing
functions call their bottom reducing functions.  Specify `init` to
obain results even when the input collection is empty or all filtered
out.

```jldoctest; setup = :(using Transducers), filter = r"EmptyResultError: .*"
julia> filtering_max = Filter(isodd)'(max);

julia> foldl(TeeRF(min, filtering_max), Map(identity), [5, 2, 6, 8, 3])
(2, 5)

julia> foldl(TeeRF(min, filtering_max), Map(identity), 2:2:8)
ERROR: EmptyResultError: ...

julia> foldl(TeeRF(min, filtering_max), Map(identity), 2:2:8; init = Init)
(2, InitialValue(max))
```
"""
struct TeeRF{N,T<:NTuple{N,Any}} <: AbstractMultiCastingRF{N}
    fs::T
    TeeRF{N,T}(fs) where {N,T} = new{N,T}(fs)
end

TeeRF{N}(fs::NTuple{N,Any}) where {N} = TeeRF{N,typeof(fs)}(fs)
TeeRF(fs::NTuple{N,Any}) where {N} = TeeRF{N}(fs)
TeeRF(f, fs...) = TeeRF((f, fs...))

@inline next(rf::TeeRF{N}, accs::NTuple{N,Any}, x) where {N} =
    map((f, a) -> next(f, a, x), rf.fs, accs)
@inline (rf::TeeRF)(accs, x) = next(rf, accs, x)

"""
    ProductRF(reducing_functions::Tuple)
    ProductRF(reducing_functions...)

Combine `N` reducing functions into a new reducing function that work
on `N`-tuple.  The `i`-th reducing function recieves the `i`-th
element of the input tuple.

Roughly speaking, `ProductRF(op₁, op₂, ..., opₙ)` is equivalent to

```julia
((a₁, a₂, ..., aₙ), (b₁, b₂, ..., bₙ)) -> (op₁(a₁, b₁), op₂(a₂, b₂), ..., opₙ(aₙ, bₙ))
```

!!! compat "Transducers.jl 0.4.32"

    New in version 0.4.32.

# Examples

Like [`TeeRF`](@ref), `ProductRF` can be used to drive multiple reducing
functions.  `ProductRF` is more "low-level" in the sense that `TeeRF`
can be defined in terms of `ProductRF` (other direction is much harder):

```jldoctest
julia> using Transducers

julia> TeeRF′(fs...) = reducingfunction(
           Map(x -> ntuple(_ -> x, length(fs))),
           ProductRF(fs...),
       );

julia> foldl(TeeRF′(min, max), Map(identity), [5, 2, 6, 8, 3])
(2, 8)
```

`ProductRF` may be useful for handling pre-existing stream whose item
type is already a tuple:

```julia; setup = (using Transducers)
julia> foldl(ProductRF(&, +), Map(x -> (isodd(x), x)), [5, 2, 6, 8, 3])
(false, 24)

julia> foldl(TeeRF(reducingfunction(Map(isodd), &), +), Map(identity), [5, 2, 6, 8, 3])
(false, 24)
```
"""
struct ProductRF{N,T<:NTuple{N,Any}} <: AbstractMultiCastingRF{N}
    fs::T
    ProductRF{N,T}(fs) where {N,T} = new{N,T}(fs)
end

ProductRF{N}(fs::NTuple{N,Any}) where {N} = ProductRF{N,typeof(fs)}(fs)
ProductRF(fs::NTuple{N,Any}) where {N} = ProductRF{N}(fs)
ProductRF(f, fs...) = ProductRF((f, fs...))

@inline next(rf::ProductRF{N}, accs::NTuple{N,Any}, xs) where {N} =
    map(next, rf.fs, accs, xs)
@inline (rf::ProductRF)(accs, x) = next(rf, accs, x)

"""
    AdHocRF(next; oninit, start = identity, complete = identity, combine = next)

Define an ad-hoc reducing function `rf`.

!!! note

    Use [`wheninit`](@ref) etc. instead of this constructor.

# Arguments
- `next`: binary function

# Keyword Arguments
- `oninit`: nullary function that generates an initial value for `next`
- `start`: unary function that pre-process the initial value for `next`
- `complete`: unary function that post-process the accumulator
- `combine`: (approximately) associative binary function for combining
  multiple results of `next` (before post-processed by `complete`).

# Examples
```jldoctest
julia> using Transducers
       using Transducers: AdHocRF

julia> rf = AdHocRF(push!, combine = append!);

julia> foldxt(rf, Map(identity), 1:4; basesize = 1, init = OnInit(() -> []))
4-element Vector{Any}:
 1
 2
 3
 4
```
"""
struct AdHocRF{OnInit,Start,Next,CompleteBasecase,Complete,Combine} <: _Function
    oninit::OnInit
    start::Start
    next::Next
    completebasecase::CompleteBasecase
    complete::Complete
    combine::Combine

    function AdHocRF{OnInit,Start,Next,CompleteBasecase,Complete,Combine}(
        oninit,
        start,
        next,
        completebasecase,
        complete,
        combine,
    ) where {OnInit,Start,Next,CompleteBasecase,Complete,Combine}
        return new{OnInit,Start,Next,CompleteBasecase,Complete,Combine}(
            oninit,
            start,
            next,
            completebasecase,
            complete,
            combine,
        )
    end
end

# Capture T::Type as Type{T}
function AdHocRF(oninit, start, op, completebasecase, complete, combine)
    return AdHocRF{
        _typeof(oninit),
        _typeof(start),
        _typeof(op),
        _typeof(completebasecase),
        _typeof(complete),
        _typeof(combine),
    }(
        oninit,
        start,
        op,
        completebasecase,
        complete,
        combine,
    )
end

AdHocRF(
    op;
    oninit = nothing,
    start = identity,
    completebasecase = identity,
    complete = identity,
    combine = nothing,
) = AdHocRF(oninit, start, op, completebasecase, complete, combine)

AdHocRF(op::AdHocRF; kwargs...) = setproperties(op, values(kwargs))

@inline (rf::AdHocRF)(acc, x) = rf.next(acc, x)

@inline start(rf::AdHocRF, init::Union{InitOf,InitialValues.InitialValue}) =
    if rf.oninit === nothing
        rf.start(initialize(init, rf.next))
    else
        # No concrete `init` is given and `oninit`is defined.
        rf.start(rf.oninit())
    end
@inline start(rf::AdHocRF, init) = rf.start(initialize(init, rf.next))
@inline next(rf::AdHocRF, acc, x) = rf.next(acc, x)
@inline complete(rf::AdHocRF, acc) = rf.complete(acc)
@inline completebasecase(rf::AdHocRF, acc) = rf.completebasecase(acc)
@inline combine(rf::AdHocRF, a, b) = something(rf.combine, rf.next)(a, b)

_asmonoid(rf::AdHocRF) = @set rf.next = _asmonoid(rf.next)
Completing(rf::AdHocRF) = rf

wheninit(oninit, op) = AdHocRF(op; oninit = oninit)
whenstart(start, op) = AdHocRF(op; start = start)
whencompletebasecase(completebasecase, op) =
    AdHocRF(op; completebasecase = completebasecase)
whencomplete(complete, op) = AdHocRF(op; complete = complete)
whencombine(combine, op) = AdHocRF(op; combine = combine)

wheninit(oninit) = op -> wheninit(oninit, op)
whenstart(start) = op -> whenstart(start, op)
whencompletebasecase(completebasecase) = op -> whencompletebasecase(completebasecase, op)
whencomplete(complete) = op -> whencomplete(complete, op)
whencombine(combine) = op -> whencombine(combine, op)

"""
    wheninit(oninit, rf) -> rf′
    wheninit(oninit) -> rf -> rf′
    whenstart(start, rf) -> rf′
    whenstart(start) -> rf -> rf′
    whencomplete(complete, rf) -> rf′
    whencomplete(complete) -> rf -> rf′
    whencombine(combine, rf) -> rf′
    whencombine(combine) -> rf -> rf′

Add initialization/completion/merging phase to arbitrary reducing
function.

The functions passed to those combinators are used as follows in
`foldl`:

```julia
init′ = oninit()  # if oninit is given; otherwise standard `init`-preprocessing
acc = start(init′)
for x in collection
    acc += rf(acc, x)
end
result = acc
return complete(result)
```

In `foldxt`, a collection is split in multiple parts and then above
`foldl` except for `complete` is run on them, yielding multiple
`result`s which are combined by repeatedly calling `combine(result_1,
result_2)`.  Note that this allows non-associative function for `next`
while `combine` *must* be associative.

See also [`next`](@ref), [`start`](@ref), [`complete`](@ref), and
[`combine`](@ref).

# Arguments
- `rf`: reducing function
- `oninit`: nullary function that generates an initial value for `rf`
- `start`: unary function that pre-process the initial value for `rf`
- `complete`: unary function that post-process the accumulator
- `combine`: (approximately) associative binary function for combining
  multiple results of `rf` (before post-processed by `complete`).

# Extended help
## Examples

An example for using non-associative reducing function in `foldxt`:

```jldoctest wheninit
julia> using Transducers

julia> collector! = push! |> whencombine(append!) |> wheninit(() -> []);

julia> foldxt(collector!, Filter(isodd), 1:5; basesize = 1)
3-element Vector{Any}:
 1
 3
 5
```

More "tightly" typed vector can returned by using BangBang.jl interface:

```jldoctest wheninit
julia> collector!! = push!! |> whencombine(append!!) |> wheninit(Vector{Union{}});

julia> foldxt(collector!!, Filter(isodd), 1:5; basesize = 1)
3-element Vector{Int64}:
 1
 3
 5
```

Online averaging algorithm can be implemented, e.g., by:

```jldoctest wheninit
julia> averaging = function add_average((sum, count), x)
           (sum + x, count + 1)
       end |> wheninit() do
           (Init(+), 0)
       end |> whencombine() do (sum1, count1), (sum2, count2)
           (sum1 + sum2), (count1 + count2)
       end |> whencomplete() do (sum, count)
           sum / count
       end;

julia> foldl(averaging, Filter(isodd), 1:5)
3.0

julia> foldxt(averaging, Filter(isodd), 1:50; basesize = 1)
25.0
```

An alternative implementation is to use [`Map`](@ref) to construct a
singleton solution and then merge it into the accumulated solution:

```jldoctest wheninit
julia> averaging2 = function merge_average((sum1, count1), (sum2, count2))
           (sum1 + sum2, count1 + count2)
       end |> whencomplete() do (sum, count)
           sum / count
       end |> Map() do x
           (x, 1)
       end';  # `'` here is important

julia> foldl(averaging2, Filter(isodd), 1:5)
3.0

julia> foldxt(averaging2, Filter(isodd), 1:50; basesize = 1)
25.0
```
"""
(wheninit, whenstart, whencomplete, whencombine)

"""
    whencompletebasecase(completebasecase, rf) -> rf′
    whencompletebasecase(completebasecase) -> rf -> rf′

Add [`completebasecase`](@ref) protocol to arbitrary reducing function.

The function `completebasecase` is used as follows in the basecase
implementation of `reduce` as follows:

```julia
init′ = oninit()
acc = start(init′)
for x in collection
    acc += rf(acc, x)
end
result = completebasecase(acc)
return result
```

The `result₁` from basecase 1 and `result₂` from basecase 2 are combined
using [`combine`](@ref) protcol:

```julia
combine(result₁, result₂)
```

!!! note

    This function is an internal experimental interface for FoldsCUDA.
"""
whencompletebasecase


text:
import Flux

struct PointwiseFeedForward
    FirstFilter::Flux.Conv
    SecondFilter::Flux.Conv
end

PointwiseFeedForward(dims::Integer) = PointwiseFeedForward(
    Flux.Conv((1, 1), dims => 4 * dims, Flux.relu) |> Flux.gpu,
    Flux.Conv((1, 1), 4 * dims => dims) |> Flux.gpu
)
Flux.@treelike PointwiseFeedForward

function (m::PointwiseFeedForward)(x::AbstractArray{T, 3} where T)
    # x = permutedims(repeat(x, outer=[1, 1, 1, 1]), [1, 4, 2, 3])
    (T, D, N) = size(x)
    x = reshape(x, (T, 1, D, N))
    x = m.SecondFilter(m.FirstFilter(x))
    return dropdims(x, dims=2)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDA_Toolkit_jll
using Base
using Base: UUID, thisminor, thismajor
using LazyArtifacts
using Preferences
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDA_Toolkit")

#JLLWrappers.@generate_main_file("CUDA_Toolkit", UUID("6e453f8e-2bea-4beb-930a-5f3197ca40e4"))

"""
    is_available()
Return whether the artifact is available for the current platform.
"""
function is_available end

const PATH = Ref{String}("")
const LIBPATH = Ref{String}("")
# We put these inter-JLL-package API values here so that they are always defined, even if there
# is no underlying wrapper held within this JLL package.
const PATH_list = String[]
const LIBPATH_list = String[]

# We sub off to JLLWrappers' dev_jll, but avoid backedges
dev_jll() = Base.invokelatest(JLLWrappers.dev_jll, "CUDA_Toolkit")

const libcuda = Sys.iswindows() ? "nvcuda" : ( Sys.islinux() ? "libcuda.so.1" : "libcuda" )
function cuDriverGetVersion()
    version_ref = Ref{Cint}()
    status = ccall((:cuDriverGetVersion, libcuda), Cuint,
                   (Ptr{Cint},), version_ref)
    @assert status == 0
    major, ver = divrem(version_ref[], 1000)
    minor, patch = divrem(ver, 10)
    return VersionNumber(major, minor, patch)
end

function select_cuda_platform(download_info::Dict, platform::AbstractPlatform = HostPlatform())
    ps = collect(filter(p -> platforms_match(p, platform), keys(download_info)))

    if @has_preference("version")
        # select requested CUDA version
        preferred_version = parse(VersionNumber, @load_preference("version"))
        filter!(ps) do p
            @assert haskey(p.tags, "cuda")
            toolkit_version = parse(VersionNumber, p.tags["cuda"])
            thisminor(toolkit_version) == thisminor(preferred_version)
        end
    else
        # select compatible CUDA versions
        driver_version = cuDriverGetVersion()
        filter!(ps) do p
            @assert haskey(p.tags, "cuda")
            toolkit_version = parse(VersionNumber, p.tags["cuda"])
            # CUDA 11+ supports semantic versioning (dubbed Enhanced Compatibility)
            driver_version >= v"11" ? thismajor(toolkit_version) <= thismajor(driver_version) :
                                      thisminor(toolkit_version) <= thisminor(driver_version)
        end
    end

    if isempty(ps)
        return nothing
    end

    # chose the most recent CUDA version
    p = last(sort(ps, by = p -> parse(VersionNumber, p.tags["cuda"])))
    return download_info[p], parse(VersionNumber, p.tags["cuda"])
end

global best_wrapper, toolkit_version
# Load Artifacts.toml file and select best platform at compile-time, since this is
# running at toplevel, and therefore will be run completely at compile-time.  We use
# a `let` block here to avoid storing unnecessary data in our `.ji` files
best_wrapper, toolkit_version = let
    artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
    valid_wrappers = Dict{Platform,String}()
    artifacts = load_artifacts_toml(artifacts_toml; pkg_uuid=UUID("6e453f8e-2bea-4beb-930a-5f3197ca40e4"))["CUDA_Toolkit"]

    # Helper function to parse triplets for us
    parse_wrapper_platform(x) = parse(Platform, x)
    function make_wrapper_dict(dir, x)
        return Dict(
            parse_wrapper_platform(basename(f)[1:end-3]) =>
                joinpath(dir, "wrappers", f) for f in x
        )
    end

    # If it's a Dict, that means this is an AnyPlatform artifact, act accordingly:
    if isa(artifacts, Dict)
        joinpath(@__DIR__, "wrappers", "any.jl")
    else
        # Otherwise, it's a Vector, and we must select the best platform
        # First, find all wrappers on-disk, parse their platforms, and match:
        wrapper_files = filter(x -> endswith(x, ".jl"), readdir(joinpath(@__DIR__, "wrappers")))
        wrappers = make_wrapper_dict(@__DIR__, wrapper_files)
        for e in artifacts
            platform = unpack_platform(e, "CUDA_Toolkit_jll", artifacts_toml)

            # Filter platforms based on what wrappers we've generated on-disk.
            # Because the wrapper file naming strategy relies upon BB's assumptions of triplet
            # parsing, it can be somewhat fragile.  So we have loaded all the wrapper filenames
            # in, parsed them, and are now using `select_platform()` to match them, which is
            # much more robust.  This step avoids a disconnect between what is recorded in the
            # `Artifacts.toml` file and what wrappers are available on-disk.
            wrapper_file = select_platform(wrappers, platform)
            if wrapper_file !== nothing
                valid_wrappers[platform] = wrapper_file
            end
        end

        # From the available options, choose the best wrapper script
        select_cuda_platform(valid_wrappers)
    end
end

# Load in the wrapper, if it's not `nothing`!
if best_wrapper === nothing
    @debug(string("Unable to load ", "CUDA_Toolkit", "; unsupported platform ", triplet(HostPlatform())))
    is_available() = false
    version() = nothing
else
    Base.include(CUDA_Toolkit_jll, best_wrapper)
    is_available() = true
    version() = toolkit_version # XXX: set this in each wrapper, maybe also exposing the full version
end

end  # module CUDA_Toolkit_jll


text:
struct MappedDictionary{I, T, F, Maps <: Tuple{AbstractDictionary{<:I}, Vararg{AbstractDictionary{<:I}}}} <: AbstractDictionary{I, T}
    f::F
    dicts::Maps
end

Base.keys(d::MappedDictionary{I}) where {I} = keys(d.dicts[1])::AbstractIndices{I}

function Base.isassigned(d::MappedDictionary{I}, i::I) where {I}
    return all(Base.Fix2(isassigned, i), d.dicts)
end

function Base.isassigned(d::MappedDictionary{I, T, <:Any, <:Tuple{AbstractDictionary{<:I}}}, i::I) where {I, T}
    return isassigned(d.dicts[1], i)
end

@propagate_inbounds function Base.getindex(d::MappedDictionary{I, T, <:Any, <:Tuple{AbstractDictionary{<:I}}}, i::I) where {I, T}
    return d.f(d.dicts[1][i])::T
end

@inline function Base.getindex(d::MappedDictionary{I, T}, i::I) where {I, T}
    @boundscheck checkindex(d.dicts[1], i)
    return d.f(map(x -> @inbounds(x[i]), d.dicts)...)::T
end

# TODO FIXME what do about tokens when there is more than one mapped dictioanry? For now, we disable them...
istokenizable(d::MappedDictionary) = false
function istokenizable(d::MappedDictionary{I, T, <:Any, <:Tuple{AbstractDictionary{<:I}}}) where {I, T}
    return istokenizable(d.maps[1])
end

@propagate_inbounds function gettokenvalue(d::MappedDictionary{I, T, <:Any, <:Tuple{AbstractDictionary{<:I}}}, t) where {I, T}
    return d.f(gettokenvalue(d.dicts[1], t))
end

function istokenassigned(d::MappedDictionary{I, T, <:Any, <:Tuple{AbstractDictionary{<:I}}}, t) where {I, T}
    return istokenassigned(d.dicts[1], t)
end

Base.similar(dict::MappedDictionary, ::Type{T}, indices) where {T} = similar(parent(dict), T, indices)
Base.empty(dict::MappedDictionary, ::Type{I}, ::Type{T}) where {I, T} = similar(parent(dict), I, T)


text:
#module HPFData
#export read_data

import ..EchelleObservation, ..EchelleObservationSet, ..OrdersType
import ..AbstractInstrument, ..HPFType
#import ..hpf_flux_hdu, ..hpf_var_hdu, ..hpf_lambda_hdu, ..hpf_all_orders
using FITSIO

"""
read_data(filename; orders)
Example:
   data = read_data(filename_list[1])
"""
function read_data(filename::String; orders::OrdersType = hpf_all_orders )
   fits = FITS(filename)
   if typeof(orders) <: Integer
      orders = orders:orders
   end
   flux = read(fits[hpf_flux_hdu])[:,orders]
   var = read(fits[hpf_var_hdu])[:,orders]
   lambda = read(fits[hpf_lambda_hdu])[:,orders]
   metadata = Dict{Symbol,Any}(:Filename=>filename,:Instrument=>:HPF)
   return EchelleObservation{eltype(lambda),eltype(flux),eltype(var),HPFType}(flux=flux,var=var,lambda=lambda,metadata=metadata)
end

function read_data!(data::EchelleObservationSet{T1,T2,T3,HPFType}, i::Int, filename::String; orders::OrdersType = hpf_all_orders ) where {T1<:Number, T2<:Number, T3<:Number, T4<:AbstractInstrument}
   @assert 1 <= i <= size(data.flux,3)
   local fits
   try
      fits = FITS(filename)
   catch
      @error("# Can't read " * filename)
   end
   if typeof(orders) <: Integer
      orders = orders:orders
   end
   data.flux[:,:,i] = read(fits[hpf_flux_hdu])[:,orders]
   data.var[:,:,i] = read(fits[hpf_var_hdu])[:,orders]
   data.lambda[:,:,i] = read(fits[hpf_lambda_hdu])[:,orders]
   data.metadata_list[i] = Dict{Symbol,Any}(:Filename=>filename)
   return data
end

"""
read_data(filenames; orders)
Example:
   data = read_data(filename_list, orders=10:15)
"""
function read_data(filename_list::AbstractArray{String,1}; orders::OrdersType = hpf_all_orders )
   @assert length(filename_list) >= 1
   data_first = read_data(filename_list[1], orders=orders)
   data = EchelleObservationSet{Float64,Float32,Float32,HPFType}(size(data_first.flux), length(filename_list), metadata=Dict(:Instrument=>:HPF) )
   data.lambda[:,:,1] = data_first.lambda
   data.flux[:,:,1]   = data_first.flux
   data.var[:,:,1]    = data_first.var
   for (i,filename) in enumerate(filename_list)
          read_data!(data, i, filename, orders=orders)
   end
   return data
end

#end # module


text:
# ---
# title: 1104. Path In Zigzag Labelled Binary Tree
# id: problem1104
# author: Indigo
# date: 2022-03-08
# difficulty: Medium
# categories: Math, Tree
# link: <https://leetcode.com/problems/path-in-zigzag-labelled-binary-tree/description/>
# hidden: true
# ---
#
# In an infinite binary tree where every node has two children, the nodes are
# labelled in row order.
#
# In the odd numbered rows (ie., the first, third, fifth,...), the labelling is
# left to right, while in the even numbered rows (second, fourth, sixth,...),
# the labelling is right to left.
#
# ![](https://assets.leetcode.com/uploads/2019/06/24/tree.png)
#
# Given the `label` of a node in this tree, return the labels in the path from
# the root of the tree to the node with that `label`.
#
#
#
# **Example 1:**
#
#
#
#     Input: label = 14
#     Output: [1,3,4,14]
#
#
# **Example 2:**
#
#
#
#     Input: label = 26
#     Output: [1,2,6,10,26]
#
#
#
#
# **Constraints:**
#
#   * `1 <= label <= 10^6`
#
#
## @lc code=start
using LeetCode

function path_in_zig_zag_tree(n::Int)
    res = Int[]
    layer = floor(Int, log2(n)) + 1
    while layer != 0
        push!(res, n)
        b = (1 << layer) + (1 << (layer - 1)) - 1
        iseven(layer) && (n = b - n)
        n >>= 1
        layer -= 1
        iseven(layer) && (n = (b >> 1) - n)
    end
    reverse!(res)
    return res
end
## @lc code=end


text:
"""
    soille_watershed(img::Matrix, mask::BitArray{2}, connectity4::Bool = true)

Implements the watershed algorithm described in Soille, Pierre, and Luc M. Vincent.
"Determining watersheds in digital pictures via flooding simulations." Lausanne, D.L.
International Society for Optics and Photonics, 1990
"""
function soille_watershed(img::Matrix, mask::BitArray{2}, connectity4::Bool = true)
    @assert size(img)==size(mask)
    onimg(ci::CartesianIndex)::Bool =
        (ci[1] >= 1) && (ci[1] <= height(img)) && (ci[2] >= 1) && (ci[2] <= width(img))
    function neighbors4(ci::CartesianIndex)
        res = (
            CartesianIndex(ci[1] - 1, ci[2]),
            CartesianIndex(ci[1], ci[2] - 1),
            CartesianIndex(ci[1], ci[2] + 1),
            CartesianIndex(ci[1] + 1, ci[2]),
        )
        return (ci[1] > 1) && (ci[2] > 1) && (ci[1] < height(img)) && (ci[2] < width(img)) ? res : filter(onimg, res)
    end
    function neighbors8(ci::CartesianIndex)
        res = (
            CartesianIndex(ci[1] - 1, ci[2] - 1),
            CartesianIndex(ci[1] - 1, ci[2]),
            CartesianIndex(ci[1] - 1, ci[2] + 1),
            CartesianIndex(ci[1], ci[2] - 1),
            CartesianIndex(ci[1], ci[2] + 1),
            CartesianIndex(ci[1] + 1, ci[2] - 1),
            CartesianIndex(ci[1] + 1, ci[2]),
            CartesianIndex(ci[1] + 1, ci[2] + 1),
        )
        return (ci[1] > 1) && (ci[1] < height(img)) && (ci[2] > 1) && (ci[2] < width(img)) ? res : filter(onimg, res)
    end
    neighbors = connectity4 ? neighbors4 : neighbors8
    INQUEUE, MASK, INIT, WSHED = -3, -2, -1, 0
    tabLabels = fill!(zeros(Int16, size(img)), INIT)
    currentLabel, flag = 0, false
    pixelList = Tuple{CartesianIndex,typeof(img[1, 1])}[]
    for r = 1:height(img), c = 1:width(img)
        if mask[r, c]
            push!(pixelList, (CartesianIndex(r, c), img[r, c]))
        end
    end
    # Ascending order by intensity
    sort!(pixelList, lt = (p1, p2) -> p1[2] < p2[2])
    fifo = Deque{CartesianIndex}()
    currentIndex, heightIndex1, heightIndex2 = 1, 1, 1
    while currentIndex <= length(pixelList)
        h = pixelList[currentIndex][2]
        for pixelIndex = heightIndex1:length(pixelList)
            pij1, v1 = pixelList[pixelIndex]
            if v1 ≠ h
                heightIndex1 = pixelIndex
                break
            end
            tabLabels[pij1] = MASK
            for cuv1 in neighbors(pij1)
                # initialize queue with neighbors at level h of current basins or watersheds
                if (tabLabels[cuv1] >= WSHED) && mask[cuv1]
                    push!(fifo, pij1)
                    tabLabels[pij1] = INQUEUE
                    break
                end
            end
        end
        while !isempty(fifo)
            pij2 = popfirst!(fifo)
            @assert pij2 isa CartesianIndex
            for cuv2 in neighbors(pij2)
                # labeling current point by inspecting neighbors
                if mask[cuv2]
                    if (tabLabels[cuv2] > 0)  # i.e. the pixel belongs to an already labeled basin
                        if (tabLabels[pij2] == INQUEUE) || ((tabLabels[pij2] == WSHED) && flag)
                            tabLabels[pij2] = tabLabels[cuv2]
                        elseif (tabLabels[pij2] > 0) && (tabLabels[pij2] != tabLabels[cuv2])
                            tabLabels[pij2] = WSHED
                            flag = false
                        end
                    elseif (tabLabels[cuv2] == WSHED)
                        if tabLabels[pij2] == INQUEUE
                            tabLabels[pij2] = WSHED
                            flag = true
                        end
                    elseif tabLabels[cuv2] == MASK
                        tabLabels[cuv2] = INQUEUE
                        push!(fifo, cuv2)
                    end
                end
            end
        end
        # check for new minima at level h
        pixelIndex = heightIndex2
        while pixelIndex <= length(pixelList)
            pij3, v3 = pixelList[pixelIndex]
            if v3 ≠ h
                # this pixel is at level h+1
                heightIndex2 = pixelIndex
                break
            end
            if tabLabels[pij3] == MASK  # the pixel is inside a new minimum
                currentLabel += 1
                push!(fifo, pij3)
                tabLabels[pij3] = currentLabel
                while isempty(fifo) == false
                    p2 = popfirst!(fifo)
                    for cuv3 in neighbors(p2) # inspect neighbors of p2
                        if (tabLabels[cuv3] == MASK) && mask[cuv3]
                            push!(fifo, cuv3)
                            tabLabels[cuv3] = currentLabel
                        end
                    end
                end
            end
            currentIndex += 1
            pixelIndex += 1
        end
    end
    fp = zeros(UInt8, size(img))
    for c = 1:width(img), r = 1:height(img)
        fp[r, c] = tabLabels[r, c] == INIT ? 0 : tabLabels[r, c]
    end
    return fp
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
# https://github.com/nolta/Winston.jl

# credit goes to https://github.com/jverzani for contributing to the first draft of this backend implementation

supported_attrs(::WinstonBackend) = merge_with_base_supported([
    :annotations,
    :linecolor,
    :fillrange,
    :fillcolor,
    :label,
    :legend,
    :seriescolor, :seriesalpha,
    :linestyle,
    :linewidth,
    :markershape,
    :markercolor,
    :markersize,
    :bins,
    :title,
    :window_title,
    :guide, :lims, :scale,
  ])
supported_types(::WinstonBackend) = [:path, :scatter, :bar]
supported_styles(::WinstonBackend) = [:auto, :solid, :dash, :dot, :dashdot]
supported_markers(::WinstonBackend) = [:none, :auto, :rect, :circle, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5]
supported_scales(::WinstonBackend) = [:identity, :log10]
is_subplot_supported(::WinstonBackend) = false

# --------------------------------------------------------------------------------------

function _initialize_backend(::WinstonBackend; kw...)
  @eval begin
    # ENV["WINSTON_OUTPUT"] = "gtk"
    warn("Winston is no longer supported... many features will likely be broken.")
    import Winston, Gtk
    export Winston, Gtk
  end
end

# ---------------------------------------------------------------------------

## dictionaries for conversion of Plots.jl names to Winston ones.
const winston_linestyle =  KW(:solid=>"solid",
                                :dash=>"dash",
                                :dot=>"dotted",
                                :dashdot=>"dotdashed"
                               )

const winston_marker = KW(:none=>".",
                            :rect => "square",
                            :circle=>"circle",
                            :diamond=>"diamond",
                            :utriangle=>"triangle",
                            :dtriangle=>"down-triangle",
                            :cross => "plus",
                            :xcross => "cross",
                            :star5 => "asterisk"
                           )

function _before_update(plt::Plot{WinstonBackend})
  Winston.ghf(plt.o)
end

# ---------------------------------------------------------------------------

function _create_backend_figure(plt::Plot{WinstonBackend})
    Winston.FramedPlot(
        title = plt.attr[:title],
        xlabel = plt.attr[:xguide],
        ylabel = plt.attr[:yguide]
    )
end

copy_remove(d::KW, s::Symbol) = delete!(copy(d), s)

function addRegressionLineWinston(d::KW, wplt)
  xs, ys = regressionXY(d[:x], d[:y])
  Winston.add(wplt, Winston.Curve(xs, ys, kind="dotted"))
end

function getWinstonItems(plt::Plot)
  if isa(plt.o, Winston.FramedPlot)
    wplt = plt.o
    window, canvas = nothing, nothing
  else
    window, canvas, wplt = plt.o
  end
  window, canvas, wplt
end

function _series_added(plt::Plot{WinstonBackend}, series::Series)
    d = series.d
  window, canvas, wplt = getWinstonItems(plt)

  # until we call it normally, do the hack
  if d[:seriestype] == :bar
    d = barHack(;d...)
  end

  e = KW()
  e[:color] = getColor(d[:linecolor])
  e[:linewidth] = d[:linewidth]
  e[:kind] = winston_linestyle[d[:linestyle]]
  e[:symbolkind] = winston_marker[d[:markershape]]
  # markercolor     # same choices as `color`, or :match will set the color to be the same as `color`
  e[:symbolsize] = d[:markersize] / 5

  # pos             # (Int,Int), move the enclosing window to this position
  # screen          # Integer, move enclosing window to this screen number (for multiscreen desktops)

  ## lintype :path, :step, :stepinverted, :sticks, :dots, :none, :histogram2d, :hexbin, :histogram, :bar
  if d[:seriestype] == :none
    Winston.add(wplt, Winston.Points(d[:x], d[:y]; copy_remove(e, :kind)..., color=getColor(d[:markercolor])))

  elseif d[:seriestype] == :path
    x, y = d[:x], d[:y]
    Winston.add(wplt, Winston.Curve(x, y; e...))

    fillrange = d[:fillrange]
    if fillrange != nothing
      if isa(fillrange, AbstractVector)
        y2 = fillrange
      else
        y2 = Float64[fillrange for yi in y]
      end
      Winston.add(wplt, Winston.FillBetween(x, y, x, y2, fillcolor=getColor(d[:fillcolor])))
    end

  elseif d[:seriestype] == :scatter
    if d[:markershape] == :none
      d[:markershape] = :circle
    end

  # elseif d[:seriestype] == :step
  #     fn = Winston.XXX

  # elseif d[:seriestype] == :stepinverted
  #     fn = Winston.XXX

  elseif d[:seriestype] == :sticks
      Winston.add(wplt, Winston.Stems(d[:x], d[:y]; e...))

  # elseif d[:seriestype] == :dots
  #     fn = Winston.XXX

  # elseif d[:seriestype] == :histogram2d
  #     fn = Winston.XXX

  # elseif d[:seriestype] == :hexbin
  #     fn = Winston.XXX

  elseif d[:seriestype] == :histogram
      hst = hist(d[:y], d[:bins])
      Winston.add(wplt, Winston.Histogram(hst...; copy_remove(e, :bins)...))

  # elseif d[:seriestype] == :bar
  #     # fn = Winston.XXX

  else
    error("seriestype $(d[:seriestype]) not supported by Winston.")

  end

  # markershape
  if d[:markershape] != :none
    Winston.add(wplt, Winston.Points(d[:x], d[:y]; copy_remove(e, :kind)..., color=getColor(d[:markercolor])))
  end

  # optionally add a regression line
  d[:smooth] && d[:seriestype] != :histogram && addRegressionLineWinston(d, wplt)

  # push!(plt.seriesargs, d)
  # plt
end

# ----------------------------------------------------------------

const _winstonNames = KW(
    :xlims => :xrange,
    :ylims => :yrange,
    :xscale => :xlog,
    :yscale => :ylog,
  )

function _update_plot_object(plt::Plot{WinstonBackend}, d::KW)
  window, canvas, wplt = getWinstonItems(plt)
  for k in (:xguide, :yguide, :title, :xlims, :ylims)
    if haskey(d, k)
      Winston.setattr(wplt, string(get(_winstonNames, k, k)), d[k])
    end
  end

  for k in (:xscale, :yscale)
    if haskey(d, k)
      islogscale = d[k] == :log10
      Winston.setattr(wplt, (k == :xscale ? :xlog : :ylog), islogscale)
    end
  end

end

# ----------------------------------------------------------------

function createWinstonAnnotationObject(plt::Plot{WinstonBackend}, x, y, val::AbstractString)
  Winston.text(x, y, val)
end

function _add_annotations{X,Y,V}(plt::Plot{WinstonBackend}, anns::AVec{Tuple{X,Y,V}})
  for ann in anns
    createWinstonAnnotationObject(plt, ann...)
  end
end

# ----------------------------------------------------------------

# function _create_subplot(subplt::Subplot{WinstonBackend}, isbefore::Bool)
#   # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example
# end

# ----------------------------------------------------------------

function addWinstonLegend(plt::Plot, wplt)
  if plt.attr[:legend] != :none
    Winston.legend(wplt, [sd[:label] for sd in plt.seriesargs])
  end
end

function Base.show(io::IO, ::MIME"image/png", plt::AbstractPlot{WinstonBackend})
  window, canvas, wplt = getWinstonItems(plt)
  addWinstonLegend(plt, wplt)
  show(io, "image/png", wplt)
end

function Base.display(::PlotsDisplay, plt::Plot{WinstonBackend})

  window, canvas, wplt = getWinstonItems(plt)

  if window == nothing
    if Winston.output_surface != :gtk
      error("Gtk is the only supported display for Winston in Plots.  Set `output_surface = gtk` in src/Winston.ini")
    end
    # initialize window
    w,h = plt.attr[:size]
    canvas = Gtk.GtkCanvasLeaf()
    window = Gtk.GtkWindowLeaf(canvas, plt.attr[:window_title], w, h)
    plt.o = (window, canvas, wplt)
  end

  addWinstonLegend(plt, wplt)

  Winston.display(canvas, wplt)
  Gtk.showall(window)
end

# function Base.display(::PlotsDisplay, subplt::Subplot{WinstonBackend})
#   # TODO: display/show the Subplot object
# end


text:
### townsend.jl ---  Perform original Lomb–Scargle periodogram
#
# Copyright (C) 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Code:

struct LSPlan{T,A,B<:AbstractVector{T},C,D,E,F,G} <: PeriodogramPlan
    times::A
    signal::B
    freq::C
    sumw::D
    w::E
    X::B
    YY::T
    noise::F
    norm::Symbol
    P::G
end

# Original algorithm that doesn't take into account uncertainties and doesn't
# fit the mean of the signal.  This is implemented following the recipe by
# * Townsend, R. H. D. 2010, ApJS, 191, 247 (URL:
#   http://dx.doi.org/10.1088/0067-0049/191/2/247,
#   Bibcode: http://adsabs.harvard.edu/abs/2010ApJS..191..247T)
function _lombscargle_orig!(P::AbstractVector{T}, times::AbstractVector{<:Real},
                            X::AbstractVector{<:Real}, freqs::AbstractVector{<:Real},
                            XX::Real) where {T<:Complex}
    N = length(X)
    @inbounds Threads.@threads for n in eachindex(freqs)
        ω = freqs[n] * 2 * pi
        XC = XS = CC = CS = zero(T)
        @inbounds for j in eachindex(times)
            S, C = sincos(ω * times[j])
            XC  += X[j]*C
            XS  += X[j]*S
            CC  += C*C
            CS  += C*S
        end
        SS       = N - CC
        s_τ, c_τ = sincos(atan(CS, CC - N / 2) / 2)
        c_τ2     = c_τ*c_τ
        s_τ2     = s_τ*s_τ
        cs_τ_CS  = 2c_τ*s_τ*CS
        #=
        P[n] = sign(CC_τ)*(YC_τ)/sqrt(abs(2*CC_τ)) + 1.0im*sign(SS_τ)*(YS_τ)/sqrt(abs(2*SS_τ))
        P[n]     = (c_τ*XC + s_τ*XS)/sqrt(2*((c_τ2*CC + cs_τ_CS + s_τ2*SS))) +
                    1.0im*(c_τ*XS - s_τ*XC)/sqrt(2*(c_τ2*SS - cs_τ_CS + s_τ2*CC))
        =#
        P[n]     = (c_τ*XC + s_τ*XS)/sqrt(2*((c_τ2*CC + cs_τ_CS + s_τ2*SS))) +
                    1.0im*(c_τ*XS - s_τ*XC)/sqrt(2*(c_τ2*SS - cs_τ_CS + s_τ2*CC))
    end
    return P
end

_periodogram!(p::LSPlan) = _lombscargle_orig!(p.P, p.times, p.X, p.freq, p.YY)
_periodogram!(times, p::LSPlan) = _lombscargle_orig!(p.P, times, p.X, p.freq, p.YY)


text:
# Normal Arm Model

type Normal <: BanditArmBase
    armDist::Distributions.Normal

    function Normal( μ::Real, σ::Real )
        new(
            Distributions.Normal( μ, σ )
        )
    end
end

function pull!( arm::Normal )
    return rand( arm.armDist )
end

function tick!( arm::Normal )
    # Do nothing
    nothing
end

function reset!( arm::Normal )
    # Do nothing
    nothing
end


text:
using JetPack, Jets, Test

@testset "pad, 1D, end" begin
    A = JopPad(JetSpace(Float64, 5), 1:10)
    @test A*[1.0:5.0;] ≈ [1.0, 2.0, 3.0, 4.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    @test A' ∘ A * [1.0:5.0;] ≈ [1.0, 2.0, 3.0, 4.0, 5.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "pad, 1D, start" begin
    A = JopPad(JetSpace(Float64, 5), -4:5)
    @test A*[1.0:5.0;] ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
    @test A'∘A*[1.0:5.0;] ≈ [1.0, 2.0, 3.0, 4.0, 5.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "pad, 1D, start+end" begin
    A = JopPad(JetSpace(Float64, 5), -4:11)
    @test A*[1.0:5.0;] ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    @test A'∘A*[1.0:5.0;] ≈ [1.0, 2.0, 3.0, 4.0, 5.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "pad, 2D, all sides" begin
    A = JopPad(JetSpace(Float64, 3, 2), 0:5, -1:3)
    @test A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈
    [0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 1.0 2.0 0.0;
     0.0 0.0 3.0 4.0 0.0;
     0.0 0.0 5.0 6.0 0.0;
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0]
    @test A'∘A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈ [1.0 2.0;3.0 4.0;5.0 6.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "truncate, 2D, all sides" begin
    A = JopPad(JetSpace(Float64, 6, 5), 2:4, 3:4)
    @test A*
    [0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 1.0 2.0 0.0;
     0.0 0.0 3.0 4.0 0.0;
     0.0 0.0 5.0 6.0 0.0;
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0] ≈
    [1.0 2.0;
     3.0 4.0;
     5.0 6.0]
    @test A'∘A*
    [0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 1.0 2.0 0.0;
     0.0 0.0 3.0 4.0 0.0;
     0.0 0.0 5.0 6.0 0.0;
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0] ≈
    [0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 1.0 2.0 0.0;
     0.0 0.0 3.0 4.0 0.0;
     0.0 0.0 5.0 6.0 0.0;
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "pad, 2D, top and sides" begin
    A = JopPad(JetSpace(Float64, 3, 2), -2:3, -1:3)
    @test A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈
    [0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 1.0 2.0 0.0;
     0.0 0.0 3.0 4.0 0.0;
     0.0 0.0 5.0 6.0 0.0]
    @test A'∘A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈ [1.0 2.0;3.0 4.0;5.0 6.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "extend pad, 1D, end" begin
    A = JopPad(JetSpace(Float64, 5), 1:10, extend=true)
    @test A*[1.0:5.0;] ≈ [1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0]
    @test A'∘A*[1.0:5.0;] ≈ [1.0, 2.0, 3.0, 4.0, 30.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "extend pad, 1D, start" begin
    A = JopPad(JetSpace(Float64, 5), -4:5, extend=true)
    @test A*[1.0:5.0;] ≈ [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 3.0, 4.0, 5.0]
    @test A'∘A*[1.0:5.0;] ≈ [6.0, 2.0, 3.0, 4.0, 5.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "extend pad, 1D, end" begin
    A = JopPad(JetSpace(Float64, 5), -4:11, extend=true)
    @test A*collect(1.0:5.0) ≈ [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0]
    @test A'∘A*[1.0:5.0;] ≈ [6.0, 2.0, 3.0, 4.0, 35.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "extend pad, 2D" begin
    A = JopPad(JetSpace(Float64, 3, 2), 0:5, -1:3, extend=true)
    @test A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈
    [1.0 1.0 1.0 2.0 2.0;
     1.0 1.0 1.0 2.0 2.0;
     3.0 3.0 3.0 4.0 4.0;
     5.0 5.0 5.0 6.0 6.0;
     5.0 5.0 5.0 6.0 6.0;
     5.0 5.0 5.0 6.0 6.0]
    @test A'∘A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈ [6.0 8.0;9.0 8.0;45.0 36.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end

@testset "extend pad and truc, 2D" begin
    A = JopPad(JetSpace(Float64, 3, 2), -2:3, -1:3, extend=true)
    @test A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈
    [1.0 1.0 1.0 2.0 2.0;
     1.0 1.0 1.0 2.0 2.0;
     1.0 1.0 1.0 2.0 2.0;
     1.0 1.0 1.0 2.0 2.0;
     3.0 3.0 3.0 4.0 4.0;
     5.0 5.0 5.0 6.0 6.0]
    @test A'∘A*[1.0 2.0;3.0 4.0;5.0 6.0] ≈ [12.0 16.0;9.0 8.0;15.0 12.0]
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test lhs ≈ rhs
end