text:
using Distributed
addprocs()

using Gadfly, Compose, DataFrames, CSV, Colors, ColorSchemes
@everywhere using Distributions, StatsBase, Distances, SharedArrays

function gen_colors(n)
    cs = distinguishable_colors(n,
        [colorant"#66c2a5", colorant"#fc8d62", colorant"#8da0cb", colorant"#e78ac3",
            colorant"#a6d854", colorant"#ffd92f", colorant"#e5c494", colorant"#b3b3b3"],
        lchoices=Float64[58, 45, 72.5, 90],
        transform=c->deuteranopic(c, 0.1),
        cchoices=Float64[20,40],
        hchoices=[75,51,35,120,180,210,270,310]
    )
    convert(Vector{Color}, cs)
end

function bc_upd(ϵ::Float64, ev_start::Bool, n_agents::Int64, averaging)
    bc_ar = ev_start == true ? collect(range(1/n_agents, stop=1 - 1/n_agents, length=n_agents)) : rand(n_agents)
    v = zeros(Float64, n_agents)
    j = 1
    @views while true
        for i in 1:n_agents
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= ϵ + eps()) # we add eps() to avoid numerical issues
            v[i] = averaging(bc_ar[:, j][fa])
        end
        bc_ar = hcat(bc_ar, v)
        if bc_ar[:, j] == bc_ar[:, j + 1]
            break
        end
        j += 1
    end
    return bc_ar
end

ϵ = .1

res = bc_upd(ϵ, true, 50, mean)

df = DataFrame(res')
names!(df, [Symbol("$i") for i in 1:size(res, 1)])
df = stack(df)
df[:steps] = repeat(0:(size(res, 2) - 1), outer=size(res, 1))

plot(df, x=:steps, y=:value, color=:variable, Geom.point, Geom.line,
    Coord.cartesian(xmax=size(res, 2)),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    Guide.title("ϵ = $ϵ"),
    Theme(key_position=:none, panel_fill="white", point_size=1.25pt, line_width=1.5pt))

function bc_upd_rad(ϵ::Float64, ev_start::Bool, r_opinion::Float64, n_normals::Int64, n_radicals::Int64, averaging)
    st = ev_start == true ? collect(range(1/n_normals, stop=1 - 1/n_normals, length=n_normals)) : rand(n_normals)
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(rd, st)
    v = vcat(rd, st)
    j = 1
    @views while true
        for i in (n_radicals + 1):(n_radicals + n_normals)
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= ϵ + eps())
            v[i] = averaging(bc_ar[:, j][fa])
        end
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    return bc_ar
end

res_rad = bc_upd_rad(.3, true, 0.9, 50, 21, mean)

df = DataFrame(res_rad')
names!(df, [Symbol("$i") for i in 1:size(res_rad, 1)])
df = stack(df)
df[:steps] = repeat(0:size(res_rad, 2) - 1, outer=size(res_rad, 1));

plot(df, x=:steps, y=:value, color=:variable, Geom.point, Geom.line,
    Coord.cartesian(xmax=size(res_rad, 2)),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    Theme(key_position=:none, panel_fill="white", point_size=1.25pt, line_width=1.5pt))

@everywhere function radicalized_normals_count(ϵ::Float64,
                                               ev_start::Bool,
                                               r_opinion::Float64,
                                               n_normals::Int64,
                                               n_radicals::Int64,
                                               averaging)
    st = ev_start == true ? collect(range(1/n_normals, stop=1 - 1/n_normals, length=n_normals)) : rand(n_normals)
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(rd, st)
    v = vcat(rd, st)
    j = 1
    @views while true
        @inbounds for i in (n_radicals + 1):(n_radicals + n_normals)
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= ϵ + eps())
            v[i] = averaging(bc_ar[fa, j])
        end
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    return sum(@. abs(bc_ar[:, end] - r_opinion) <= 10^-3) - n_radicals
end

function radical_count_parallel(ev_start::Bool, r_opinion::Float64, n_normals::Int64, averaging)
    res = Array{Int64,2}(undef, 50, 50) # creates a 50 * 50 grid, each cell corresponding to one combination of ϵ and number of radicals
    for j in 1:50
        @inbounds res[:, j] = pmap(i -> radicalized_normals_count(i/100, ev_start, r_opinion, n_normals, j, averaging), 1:50)
    end
    return res
end

out = radical_count_parallel(true, 1.0, 50, mean);

function xlabelname(x)
    n = x/100
    return "$n"
end

function ylabelname(x)
    n = abs(x - 50)
    return "$n"
end

ticks = collect(0:10:50);

Gadfly.spy(rotl90(out),
    Guide.ColorKey(title="Converted\ntruth-seekers"),
    Guide.xlabel("ϵ"), Guide.ylabel("Campaigners"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p), minvalue=0, maxvalue=50),
    Theme(grid_color=colorant"white", panel_fill="white")
)

@everywhere function radicalized_normals_average(ϵ::Float64,
                                                 r_opinion::Float64,
                                                 n_normals::Int64,
                                                 n_radicals::Int64,
                                                 averaging)
    st = rand(n_normals)
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(rd, st)
    v = vcat(rd, st)
    j = 1
    @views while true
        for i in (n_radicals + 1):(n_radicals + n_normals)
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= ϵ)
            @inbounds v[i] = averaging(bc_ar[fa, j])
        end
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    return sum(@. abs(@view(bc_ar[:, end]) - r_opinion) <= 10^-3) - n_radicals
end

function radical_average_parallel(r_opinion::Float64, n_normals::Int64, n_radicals::Int64, averaging)
    res = Array{Int64,2}(undef, 50, n_radicals)
    for j in 1:n_radicals
        @inbounds res[:, j] = pmap(i -> radicalized_normals_average(i/100, r_opinion, n_normals, j, averaging), 1:50)
    end
    return res
end

function rand_sim(r_opinion::Float64, n_normals::Int64, n_radicals::Int64, n_sim::Int64, averaging)
    rand_res = Array{Int64,3}(undef, 50, n_radicals, n_sim)
    for i in 1:n_sim
        @inbounds rand_res[:, :, i] = radical_average_parallel(r_opinion, n_normals, n_radicals, averaging)
    end
    return rand_res
end

rand_res = rand_sim(1.0, 50, 50, 100, mean)

out_av = mean(rand_res, dims=3)

Gadfly.spy(rotl90(out_av[:, :, 1]),
    Guide.ColorKey(title="Converted\ntruth-seekers"),
    Guide.xlabel("ϵ"), Guide.ylabel("Campaigners"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p), minvalue=0, maxvalue=50),
    Theme(grid_color=colorant"white", panel_fill="white")
)

out_var = mapslices(variation, rand_res, dims=3)
mean(filter(x -> !isnan(x),out_var)), std(filter(x -> !isnan(x), out_var))

@everywhere function mixed_community_sim(ϵ::Float64,
                                         ev_start::Bool,
                                         r_opinion::Float64,
                                         n_truthseekers::Int64,
                                         n_indifs::Int64,
                                         n_radicals::Int64,
                                         α::Float64,
                                         τ::Float64,
                                         averaging_op,
                                         averaging_eps)
    n_normals = n_truthseekers + n_indifs
    st = ev_start == true ? collect(range(1/n_normals, stop=1 - 1/n_normals, length=n_normals)) : rand(n_normals)
    #δ = ev_start == true ? 10^-12 : 0.0
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(st, rd)
    v = vcat(st, rd)
    #eps = ϵ + δ
    j = 1
    @views while true
        for i in 1:n_normals
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= ϵ + eps())
            @inbounds v[i] = i <= n_truthseekers ? (1 - α)*averaging_op(bc_ar[fa, j]) + α*τ : averaging_op(bc_ar[fa, j])
        end
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    res = bc_ar[1:n_truthseekers + n_indifs, end]
    res_t = sum(@. abs(res[1:n_truthseekers] - τ) <= 10^-3)
    res_i = sum(@. abs(res[n_truthseekers + 1:n_normals] - τ) <= 10^-3)
    return res_t/n_truthseekers, res_i/n_indifs
end

@everywhere function sim_mixed(eps::Float64,
                               r_opinion::Float64,
                               n_indifs::Int64,
                               n_radicals::Int64,
                               α::Float64,
                               τ::Float64)
    mm = Array{Float64,2}(undef, 100, 2)
    for i in 1:100
        mm[i, :] = [mixed_community_sim(eps, false, r_opinion, 50, n_indifs, n_radicals, α, τ, mean, mean)...]
    end
    return mean(mm, dims=1)
end

τ = .65
ρ = .9
α = .5
ϵ = .1

truth_res = SharedArray{Float64,2}(50, 50)
indif_res = SharedArray{Float64,2}(50, 50)

@sync @distributed for j in 1:50
    for i in 1:50
        truth_res[i, j], indif_res[i, j] = sim_mixed(ϵ, ρ, i, j, α, τ)
    end
end

Gadfly.spy(truth_res[end:-1:1, :],
    Guide.ColorKey(title="Truth seekers\nbelieving τ (proportion)"),
    Guide.xlabel("Campaigners"), Guide.ylabel("Free riders"),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p), minvalue=0, maxvalue=maximum(truth_res)),
    Guide.title("τ = $τ, ρ = $ρ, ϵ = $ϵ, α = $α"),
    Theme(grid_color=colorant"white", panel_fill="white")
)

Gadfly.spy(indif_res[end:-1:1, :],
    Guide.ColorKey(title="Free riders\nbelieving τ (proportion)"),
    Guide.xlabel("Campaigners"), Guide.ylabel("Free riders"),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p), minvalue=0, maxvalue=maximum(indif_res)),
    Guide.title("τ = $τ, ρ = $ρ, ϵ = $ϵ, α = $α"),
    Theme(grid_color=colorant"white", panel_fill="white")
)

@everywhere function sim_mixed_rad(ϵ::Float64,
                                   r_opinion::Float64,
                                   n_indifs::Int64,
                                   n_radicals::Int64,
                                   α::Float64,
                                   τ::Float64) # does basically the same as the previous functions, except that we now average over 50 simulations
    mm = Array{Float64,2}(undef, 50, 2)
    for i in 1:50
        mm[i, :] = [mixed_community_sim(ϵ, false, r_opinion, 50, n_indifs, n_radicals, α, τ, mean, mean)...]
    end
    return mean(mm, dims=1)
end

function away_from_truth(τ::Float64, ρ::Float64, α::Float64, n_indifs::Int64)
    res = Array{Float64,2}(undef, 50, 50)
    for j in 1:50
        @inbounds res[:, j] = pmap(i -> sim_mixed_rad(i/100, ρ, n_indifs, j, α, τ)[1], 1:50)
    end
    return res
end

out3 = away_from_truth(.1, .3, .5, 0)
out5 = away_from_truth(.1, .5, .5, 0)
out7 = away_from_truth(.1, .7, .5, 0)

Gadfly.spy(rotl90(out3)*50,
    Guide.ColorKey(title="Truth-seekers believing\nthe truth"),
    Guide.xlabel("ϵ"), Guide.ylabel("Campaigners"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p), minvalue=0, maxvalue=50),
    Guide.title("τ = 0.1, ρ = 0.3, α = 0.5"),
    Theme(grid_color=colorant"white", panel_fill="white")
)

out3_50 = away_from_truth(.1, .3, .5, 50)
out5_50 = away_from_truth(.1, .5, .5, 50)
out7_50 = away_from_truth(.1, .7, .5, 50);

Gadfly.spy(rotl90(out7 - out7_50)*50,
    Guide.ColorKey(title="Truth-seekers believing\nthe truth"),
    Guide.xlabel("ϵ"), Guide.ylabel("Campaigners"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p), minvalue=0, maxvalue=50),
    Guide.title("τ = $τ, ρ = $ρ, α = $α, #FR = 50"),
    Theme(grid_color=colorant"white", panel_fill="white")
)

max_stars_repo_name: IgorDouven/Misinformation

max_issues_repo_name: IgorDouven/Misinformation

max_forks_repo_name: IgorDouven/Misinformation


text:
function attach_constraints_rows(cellvec,cellconstr,cellmask=Fill(true,length(cellconstr)))
  apply(ConstrainRowsKernel(),cellvec,cellconstr,cellmask)
end

function attach_constraints_cols(cellmat,cellconstr,cellmask=Fill(true,length(cellconstr)))
  cellconstr_t = apply(transpose,cellconstr)
  apply(ConstrainColsKernel(),cellmat,cellconstr_t,cellmask)
end

struct ConstrainRowsKernel <: Kernel end

function Arrays.kernel_cache(k::ConstrainRowsKernel,array::AbstractArray,constr,mask)
  kernel_cache(MulKernel(),constr,array)
end

@inline function Arrays.apply_kernel!(cache,k::ConstrainRowsKernel,array::AbstractArray,constr,mask)
  if mask
    apply_kernel!(cache,MulKernel(),constr,array)
  else
    array
  end
end

function Arrays.kernel_cache(k::ConstrainRowsKernel,matvec::Tuple,constr,mask)
  mat, vec = matvec
  cmat = kernel_cache(k,mat,constr,mask)
  cvec = kernel_cache(k,vec,constr,mask)
  (cmat,cvec)
end

@inline function Arrays.apply_kernel!(cache,k::ConstrainRowsKernel,matvec::Tuple,constr,mask)
  if mask
    cmat, cvec = cache
    mat, vec = matvec
    _mat = apply_kernel!(cmat,k,mat,constr,mask)
    _vec = apply_kernel!(cvec,k,vec,constr,mask)
    (_mat,_vec)
  else
    matvec
  end
end

struct ConstrainColsKernel <: Kernel end

function Arrays.kernel_cache(k::ConstrainColsKernel,array::AbstractArray,constr_t,mask)
  kernel_cache(MulKernel(),array,constr_t)
end

@inline function Arrays.apply_kernel!(cache,k::ConstrainColsKernel,array::AbstractArray,constr_t,mask)
  if mask
    apply_kernel!(cache,MulKernel(),array,constr_t)
  else
    array
  end
end

function Arrays.kernel_cache(k::ConstrainColsKernel,matvec::Tuple,constr_t,mask)
  mat, vec = matvec
  kernel_cache(k,mat,constr_t,mask)
end

@inline function Arrays.apply_kernel!(cache,k::ConstrainColsKernel,matvec::Tuple,constr_t,mask)
  if mask
    mat, vec = matvec
    _mat = apply_kernel!(cache,k,mat,constr_t,mask)
    (_mat,vec)
  else
    matvec
  end
end

function merge_cell_constraints_at_skeleton(cL,cR,axesL_rows,axesR_rows,axesL_cols,axesR_cols)
  blocks = (cL,cR)
  blockids = [(1,1),(2,2)]
  axs_rows = create_array_of_blocked_axes(axesL_rows,axesR_rows)
  axs_cols = create_array_of_blocked_axes(axesL_cols,axesR_cols)
  axs = apply((r,c) -> (r[1],c[1]),axs_rows,axs_cols)
  VectorOfBlockArrayCoo(blocks,blockids,axs)
end

function identity_constraints(cell_axes)
  apply(IdentityConstraintKernel(),cell_axes)
end

struct IdentityConstraintKernel <: Kernel end

function Arrays.kernel_cache(k::IdentityConstraintKernel,axs)
  n = length(axs[1])
  a = zeros(n,n)
  CachedArray(a)
end

function Arrays.apply_kernel!(cache,k::IdentityConstraintKernel,axs)
  n = length(axs[1])
  setsize!(cache,(n,n))
  a = cache.array
  fill!(a,zero(eltype(a)))
  o = one(eltype(a))
  @inbounds for i in 1:size(a,1)
    a[i,i] = o
  end
  a
end

ext: jl


text:
export endEncoding!

abstract type MtlCommandEncoder end

Base.:(==)(a::T, b::T) where {T <: MtlCommandEncoder} = a.handle == b.handle
Base.hash(q::MtlCommandEncoder, h::UInt) = hash(q.handle, h)

function unsafe_destroy!(cce::MtlCommandEncoder)
    mtRelease(cce.handle)
end

endEncoding!(ce::MtlCommandEncoder) = mtCommandEncoderEndEncoding(ce.handle)
Base.close(ce::MtlCommandEncoder) = endEncoding!(ce)


text:
# This file contains tools to iterate through the files that store the block
# chain on the hard drive.

# an iterator/stream chimera

# TODO: add fields for the current file and(? maybe not and) position in the
# file for faster comparison
mutable struct BCIterator
    io::IOStream
end

# This function does not make sense for i > 0, because we cannot be sure that in
# file i - 1 does not contain any blocks that come before the first block and
# that the first block of file i is actually the first block of the remaining
# chain.
function BCIterator(i::Integer)
    BCIterator(open(get_block_chain_file_path(i)))
end

BCIterator() = BCIterator(0)

# function Base.getindex(x::BCIterator, d::Symbol)
#     if     d == :file_name  x.io.name[end-12:end-1]
#     elseif d == :file_path  x.io.name[7:end-1]
#     elseif d == :file_index get_file_num(x)
#     else                    throw(KeyError(d))
#     end
# end

# # TODO: getproperty
# Base.getindex(x::BCIterator, ::Type{Val{:file_name}}) = x.io.name[end-12:end-1]
# Base.getindex(x::BCIterator, ::Type{Val{:file_path}}) = x.io.name[7:end-1]
# Base.getindex(x::BCIterator, ::Type{Val{:file_index}}) = x.io.name[7:end-1]

function is_last_file(x::BCIterator)
    get_num_block_chain_files() - 1 == get_file_num(x)
end

function open_next_file(x::BCIterator)
    if is_last_file(x)
        # TODO: don't use errors for programming logic
        throw(NoMoreFilesError())
    end
    new_file_idx = get_file_num(x) + 1
    # println(new_file_idx)
    close(x.io)
    x.io = open(get_block_chain_file_path(new_file_idx))
    return x
end

function get_file_num(x::BCIterator)
    parse(Int, x.io.name[end-9:end-5])
end

function get_file_pos(x::BCIterator)
    position(x.io)
end

function Base.:<(x::BCIterator, y::BCIterator)
    fn_x, fn_y = get_file_num(x), get_file_num(y)

    (fn_x < fn_y) ||
        ((fn_x == fn_y) &&
         get_file_pos(x) < get_file_pos(y))
end
Base.:>(x::BCIterator, y::BCIterator) = y < x
function Base.:<=(x::BCIterator, y::BCIterator)
    fn_x, fn_y = get_file_num(x), get_file_num(y)

    (fn_x < fn_y) ||
        ((fn_x == fn_y) &&
         get_file_pos(x) <= get_file_pos(y))
end
Base.:>=(x::BCIterator, y::BCIterator) = y <= x
function Base.:!=(x::BCIterator, y::BCIterator)
    (get_file_num(x) != get_file_num(y)) ||
        (get_file_pos(x) != get_file_pos(y))
end
function Base.:(==)(x::BCIterator, y::BCIterator)
    (get_file_num(x) == get_file_num(y)) &&
        (get_file_pos(x) == get_file_pos(y))
end

Base.close(x::BCIterator)      = close(x.io)
Base.read(x::BCIterator, T)    = read(x.io, T)
Base.read(x::BCIterator, T, n) = read(x.io, T, n)
Base.eof(x::BCIterator)        = eof(x.io)
Base.mark(x::BCIterator)       = mark(x.io)
Base.reset(x::BCIterator)      = reset(x.io)
Base.seekend(x::BCIterator)    = seekend(x.io)
Base.seekstart(x::BCIterator)  = seekstart(x.io)
Base.seek(x::BCIterator, pos)  = seek(x.io, pos)
Base.skip(x::BCIterator, off)  = skip(x.io, off)

# TODO: these functions do parsing and some of them are called all the time,
# the results should probably be cached
function get_block_chain_file_names()
    [x for x in readdir(DIR) if occursin(r"^blk", x)]
end
function get_block_chain_file_paths()
    [joinpath(DIR, x)
     for x in readdir(DIR) if occursin(r"^blk", x)]
end
function get_num_block_chain_files()
    length(get_block_chain_file_names())
end
function get_block_chain_file_name(i::Integer)
    "blk" * @sprintf("%05d", i) * ".dat"
end
function get_block_chain_file_path(i)
    joinpath(DIR, get_block_chain_file_name(i))
end

# TODO: should this function return something?
# function check_magic_bytes(io, the_magic::UInt32)
#     nextfourbytes = read(io, UInt32)
#     if nextfourbytes == zero(nextfourbytes)
#         nextfourbytes = read(io, UInt32)
#     end
#     if nextfourbytes != the_magic
#         throw(MagicBytesError(nextfourbytes))
#     end
#     nothing
# end
# check_magic_bytes(io) = check_magic_bytes(io, MAGIC)
function check_magic_bytes(io, the_magic::NTuple{4, UInt8})
    # at height 543617 there where some extra zero bytes after the block:
    while read(io, UInt8) == 0x00 end
    skip(io, -1)
    (read(io, UInt8) == the_magic[1] &&
     read(io, UInt8) == the_magic[2] &&
     read(io, UInt8) == the_magic[3] &&
     read(io, UInt8) == the_magic[4]) ||
     throw(MagicBytesError(0xffffffff))

    nothing
end
check_magic_bytes(io) = check_magic_bytes(io, to_byte_tuple(MAGIC))

"""
    seek_magic_bytes!(io, the_magic)

Advance `io` to next `the_magic`.
"""
function seek_magic_bytes!(io::IO, the_magic::NTuple{4, UInt8})
    @label a
    while read(io, UInt8) != the_magic[1] end
    if    read(io, UInt8) != the_magic[2] @goto a end
    if    read(io, UInt8) != the_magic[3] @goto a end
    if    read(io, UInt8) != the_magic[4] @goto a end
end

# I don't think this makes much sense here, they are about the same speed,
# @btime is not possible here, because it the function mutates io in place and we will reach EOF.
const CHAR_JUMP_TABLE = ntuple(i -> findfirst(isequal(UInt8(i - 1)), to_byte_tuple(MAGIC)) |>
                               x -> x == nothing ? UInt8(sizeof(MAGIC)) : UInt8(x - 1),
                               typemax(UInt8) + 1)

"""
    seek_magic_bytes_boyer_moore!(io, the_magic)

Advance `io` to next `the_magic`. Alternative implementation to `seek_magic_bytes!`
"""
function seek_magic_bytes_boyer_moore!(io::IO, the_magic::NTuple{4, UInt8})
    @label a
    # bad character rule
    current_byte = read(io, UInt8)
    if current_byte != the_magic[4]
        jump_length = CHAR_JUMP_TABLE[current_byte + 1]
        skip(io, jump_length)
        @goto a
    else
        skip(io, -2)
        current_byte = read(io, UInt8)
        if current_byte != the_magic[3]
            skip(io, 5)
            @goto a
        end

        skip(io, -2)
        current_byte = read(io, UInt8)
        if current_byte != the_magic[2]
            skip(io, 6)
            @goto a
        end

        skip(io, -2)
        current_byte = read(io, UInt8)
        if current_byte != the_magic[1]
            skip(io, 7)
            @goto a
        end

        skip(io, sizeof(the_magic) - 1)
    end
end


text:
using Documenter
using Imogen

DocMeta.setdocmeta!(Imogen, :DocTestSetup, :(using Imogen); recursive=true)
makedocs(
    sitename = "Imogen",
    format = Documenter.HTML(),
    modules = [Imogen],
    authors = "Douglas G. Moore",
    pages = Any[
        "Home" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/dglmoore/Imogen.jl.git"
)


text:
#MDDatasets: "cross" measurements
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#Constants used to filter out undesired crossings:
#-------------------------------------------------------------------------------
#RISE & FALL are mutually exclusive (but user can want both):
const XINGTYPE_RISE = UInt(0x1)
const XINGTYPE_FALL = UInt(0x2)
#const XINGTYPE_DIRMASK = UInt(0x3)
#FLAT & SING are mutually excluive (but user can want both):
const XINGTYPE_SING = UInt(0x4) #Does not cross @ single point
const XINGTYPE_FLAT = UInt(0x8) #Does not cross @ single point
#THRU & REV are mutually exclusive (but user can want both):
const XINGTYPE_THRU = UInt(0x10) #Goes through
const XINGTYPE_REV = UInt(0x20) #Reverses out.  Does not fully cross
const XINGTYPE_FIRSTLAST = UInt(0x40) #First/last point leaving/entering zero state

const XINGTYPE_ALL =
	XINGTYPE_RISE|XINGTYPE_FALL |
	XINGTYPE_SING|XINGTYPE_FLAT |
	XINGTYPE_THRU|XINGTYPE_REV|XINGTYPE_FIRSTLAST

#==Types
===============================================================================#

struct CrossType
	v::UInt
end

function CrossType(;rise=true, fall=true, sing=true, flat=true, thru=true, rev=false, firstlast=false)
	result = rise*XINGTYPE_RISE | fall*XINGTYPE_FALL |
	         sing*XINGTYPE_SING | flat*XINGTYPE_FLAT |
	         thru*XINGTYPE_THRU | rev*XINGTYPE_REV | firstlast*XINGTYPE_FIRSTLAST
	return CrossType(result)
end

function CrossType(id::Symbol)
	if :all == id
		return CrossType(rise=true, fall=true, sing=true, flat=true, thru=true, rev=true, firstlast=true)
	elseif :default == id
		return CrossType()
	elseif :rise == id
		return CrossType(rise=true, fall=false)
	elseif :fall == id
		return CrossType(rise=false, fall=true)
	elseif :risefall == id
		return CrossType(rise=true, fall=true)
	else
		throw("Unknown crossing-type preset: $id")
	end
end

#==Main algorithm
===============================================================================#

#Finds all zero crossing indices in a dataset, up to nmax.
#nmax = 0: find all crossings
#-------------------------------------------------------------------------------
function icross(d::DataF1, nmax::Integer, xstart::Real, allow::CrossType)
	#TODO: make into function if can force inline
	#(SGNTOXINGDIR >> (1+Int(sgncur)))&0x3 returns the type of xing:
	SGNTOXINGDIR = XINGTYPE_FALL|XINGTYPE_RISE<<2 #WANTCONST
	EMPTYRESULT = Limits1D{Int}[] #WANTCONST
	allow = allow.v #Get value

	validate(d); #Expensive, but might avoid headaches
	x = d.x; y = d.y #shortcuts
	ny = length(y)
	if ny < 1; return EMPTYRESULT; end
	xstart = max(x[1], xstart) #Simplify algorithm below
	if xstart > x[end]; return EMPTYRESULT; end
	nmax = nmax<1 ? ny : min(nmax, ny)
	idx = Vector{Limits1D{Int}}(undef, nmax) #resultant array of indices
	n = 0 #Index into idx[]
	i = 1 #Index into x/y[]
	while x[i] < xstart #Fast-forward to start point
		i+=1
	end
	#Here: x[i] >= xstart

	if x[i] > xstart
		i -= 1 #ok: xstart = max(x[1], xstart)
	end
	istart = i
	ystart = value(d, x=xstart) #TODO: direct interp, instead of using value()
	sgncur = sign(ystart)
	lastnzpos = 0; #Pretend like this is the start of the data
	if (sgncur != 0); lastnzpos = i; end

	#Move up to first non-zero value... if not already done:
	while lastnzpos < 1 && i < ny
		i+=1
		sgncur = sign(y[i])
		if (sgncur != 0); lastnzpos = i; end
	end
	if lastnzpos < 1; return EMPTYRESULT; end
	#Here: y[i] @ first nonzero value since xstart

	#Register first crossing if all initial values were zeros since xstart:
	if 0 == ystart
		xingdir = (SGNTOXINGDIR >> (1+Int(sgncur)))&0x3
		xingtype = XINGTYPE_FIRSTLAST
		xingtype |= xingdir

		if allow & xingtype == xingtype
			needsinterp = (1==lastnzpos-istart && y[istart] != 0)
			n+=1
			if needsinterp
				idx[n] = Limits1D(istart, istart+1)
			else
				idx[n] = Limits1D(lastnzpos-1, lastnzpos-1)
			end
		end
	end

	sgnenter=sgncur #Not really needed; mostly a declaration
	posfirstzero = 0 #No longer reading a zero sequence
	sgnprev = sgncur
	while i < ny
		i+=1
		sgncur = sign(y[i])
		if sgncur != sgnprev
			xingdir = (SGNTOXINGDIR >> (1+Int(sgncur)))&0x3
			if posfirstzero > 0
				xingtype = (posfirstzero == i-1 ? XINGTYPE_SING : XINGTYPE_FLAT)
				xingtype |= (0 == sgnenter+sgncur ? XINGTYPE_THRU : XINGTYPE_REV)
				xingtype |= xingdir
				if allow & xingtype == xingtype
					n+=1
					idx[n] = Limits1D(posfirstzero, i-1)
				end
				posfirstzero = 0
			elseif 0==sgncur
				sgnenter=sgnprev
				posfirstzero = i
			else
				xingtype = XINGTYPE_THRU|XINGTYPE_SING
				xingtype |= (SGNTOXINGDIR >> (1+Int(sgncur)))&0x3
				xingtype |= xingdir
				if allow & xingtype == xingtype
					n+=1
					idx[n] = Limits1D(i-1, i)
				end
			end
			if n >= nmax; break; end
			sgnprev = sgncur
		end
	end

	reachedend = (i >= ny)
	if reachedend && posfirstzero>0
		xingdir = (SGNTOXINGDIR >> (1-Int(sgnenter)))&0x3
		xingtype = XINGTYPE_FIRSTLAST
		xingtype |= xingdir
		if allow & xingtype == xingtype
			n+=1
			idx[n] = Limits1D(posfirstzero, posfirstzero)
		end
	end

	return resize!(idx, n)
end

#TODO: what about infinitiy? convert(Float32,NaN)
#-------------------------------------------------------------------------------
function xveccross(d::DataF1{TX,TY}, nmax::Integer,
	xstart::Real, allow::CrossType) where {TX<:Number, TY<:Number}
	idx = icross(d, nmax, xstart, allow)
	TR = typeof(one(promote_type(TX,TY))/2) #TODO: is there better way?
	result = Vector{TR}(undef, length(idx))
	for i in 1:length(idx)
		rng = idx[i]
		x1 = d.x[rng.min]; y1 = d.y[rng.min]
		if zero(TY) == y1
			result[i] = (d.x[rng.max] + x1)/2
		else
			Δy = d.y[rng.max] - y1
			Δx = d.x[rng.max] - x1
			result[i] = (0-y1)*(Δx/Δy)+x1
		end
	end
	return result
end

#-------------------------------------------------------------------------------
function xcross(d::DataF1; nmax::Integer=0, xstart::Real=-Inf,
	allow::CrossType=CrossType())
	x = xveccross(d, nmax, xstart, allow)
	return DataF1(x, x)
end
function xcross(::DS{:event}, d::DataF1, args...; kwargs...)
	d = xcross(d, args...;kwargs...)
	return DataF1(collect(1:length(d.x)), d.y)
end

#xcross1: return a single crossing point (new name for type stability)
#-------------------------------------------------------------------------------
function xcross1(d::DataF1; n::Integer=1, xstart::Real=-Inf,
	allow::CrossType=CrossType())
	n = max(n, 1)
	x = xveccross(d, n, xstart, allow)
	if length(x) < n
		return convert(eltype(x), NaN) #TODO: Will fail with int.  Use NA.
	else
		return x[n]
	end
end

#TODO: Make more efficient (don't use "value")
#-------------------------------------------------------------------------------
function ycross(d1::DataF1{TX,TY}, d2::T; nmax::Integer=0, xstart::Real=-Inf,
	allow::CrossType=CrossType()) where {TX<:Number, TY<:Number, T<:DF1_Num}
	x = xveccross(d1-d2, nmax, xstart, allow)
	TR = typeof(one(promote_type(TX,TY))/2) #TODO: is there better way?
	y = Vector{TR}(undef, length(x))
	for i in 1:length(x)
		y[i] = value(d1, x=x[i])
	end
	return DataF1(x, y)
end
function ycross(::DS{:event}, d1::DataF1, d2::T, args...; kwargs...) where T<:DF1_Num
	d = ycross(args...;kwargs...)
	return DataF1(collect(1:length(d.x)), d.y)
end

#ycross1: return a single crossing point (new name for type stability)
#-------------------------------------------------------------------------------
function ycross1(d1::DataF1{TX,TY}, d2::T; n::Integer=1,	xstart::Real=-Inf,
	allow::CrossType=CrossType()) where {TX<:Number, TY<:Number, T<:DF1_Num}
	n = max(n, 1)
	x = xveccross(d1-d2, n, xstart, allow)
	TR = typeof(one(promote_type(TX,TY))/2) #TODO: is there better way?
	if length(x) < n
		return convert(TR, NaN) #TODO: Will fail with int.  Use NA.
	else
		return value(d1, x=x[n])
	end
end

#Last line


text:
u1 = UncertainValue(Normal, 0, 0.1)
u2 = UncertainValue(Uniform, rand(100))
u3 = UncertainValue(rand(100))
u4 = UncertainValue(Normal, 5, 0.1)
u5 = UncertainValue(Normal, 10, 1.5)

vals = [u1, u2, u3, u4, u5, 6]
wts = [1, 3, 3, 1, 0.5, 0.1]

pop = UncertainScalarPopulation(vals, wts)
@test resample(pop) isa Real
@test resample(pop, 10) isa Vector{T} where {T <: Real}

x = UncertainScalarPopulation([u1, u2, u3, u4, u5], rand(5))
vals = [pop, u5, x, u4, u5, 6]
wts = rand(length(vals))

pop2 = UncertainScalarPopulation(vals, wts)
@test resample(pop2) isa Real
@test resample(pop2, 100000) isa Vector{T} where {T <: Real}

# Create a population consisting of a mixture of different types of uncertain values,
# both theoretical, fitted from empirical distributions, and discrete populations.
u1 = UncertainValue(Normal, 0, 0.1)
u2 = UncertainValue(Uniform, rand(100))
u3 = UncertainValue(rand(100))
u4 = UncertainValue(Normal, 5, 0.1)
u5 = UncertainValue(Normal, 10, 1.5)

wts = [1, 3, 3, 1, 0.5]
vals = [u1, u2, u3, u4, u5]
pop3 = UncertainScalarPopulation(vals, wts)

@test resample(pop3) isa Real
@test resample(pop3, 100) isa Vector{T} where {T <: Real}


text:
using OnlineStatsBase, LearnBase, OrderedCollections, StatsBase, Statistics, Dates, Test

import OnlineStatsBase: _fit!, _merge!, Mean, Variance

include("test_stats.jl")

#-----------------------------------------------------------------------# Random
@testset "Random Stuff" begin
string(Series(Mean(), Variance()))
string(Mean())
string(McclainWeight())
@test Mean() != Variance()
@test !(Mean() == Variance())
@test Mean() == merge(Mean(), Mean())
@test Mean() == merge(Mean(), fit!(Variance(), 1:5))
@test_throws Exception fit!(Mean(), "abc")
@test_throws Exception OnlineStatsBase._fit!(Mean(), "a")
@test OnlineStatsBase.input(Mean()) == Number
end

#-----------------------------------------------------------------------# Weight
@testset "Weight" begin
function test_weight(w::OnlineStatsBase.Weight, f::Function)
    println("  > $w")
    @test w == copy(w)
    for i in 1:20
        @test w(i) == f(i)
    end
end
test_weight(@inferred(EqualWeight()),               i -> 1 / i)
test_weight(@inferred(ExponentialWeight(.1)),       i -> ifelse(i==1, 1.0, .1))
test_weight(@inferred(LearningRate(.6)),            i -> 1 / i^.6)
test_weight(@inferred(LearningRate2(.5)),           i -> 1 / (1 + .5*(i-1)))
test_weight(@inferred(HarmonicWeight(4.)),          i -> 4 / (4 + i - 1))

@test ExponentialWeight(20) == ExponentialWeight(2 / 21)

@testset "McclainWeight" begin
    w = McclainWeight(.1)
    for i in 2:100
        @test .1 < w(i) < 1
    end
end
@testset "first weight is one" begin
    for w in [EqualWeight(), ExponentialWeight(), LearningRate(), LearningRate2(),
              HarmonicWeight(), McclainWeight()]
        @test w(1) == 1
    end
end
end #Weight

#-----------------------------------------------------------------------# Iteration
@testset "Iteration" begin
    x, y = randn(100,10), randn(100)
    @test fit!(Counter(AbstractVector), OnlineStatsBase.eachrow(x)).n == 100
    @test fit!(Counter(AbstractVector), OnlineStatsBase.eachcol(x)).n == 10
    @test fit!(Counter(Tuple),  zip(OnlineStatsBase.eachrow(x), y)).n == 100
    @test fit!(Counter(Tuple),  zip(OnlineStatsBase.eachcol(x), y)).n == 10

    for (j, xj) in enumerate(OnlineStatsBase.eachcol(x))
        @test xj == x[:, j]
    end
    for (i, xi) in enumerate(OnlineStatsBase.eachrow(x))
        @test xi == x[i, :]
    end
    @inferred OnlineStatsBase.eachrow(x)
    @inferred OnlineStatsBase.eachcol(x)
    @inferred zip(OnlineStatsBase.eachrow(x), y)
    @inferred zip(OnlineStatsBase.eachcol(x), y)

    @test length(OnlineStatsBase.eachcol(x)) == 10
    @test length(OnlineStatsBase.eachrow(x)) == 100
    @test length(zip(OnlineStatsBase.eachrow(x), y)) == 100
end


text:
"""
    DeclareMissings(; values)

DeclareMissings (or replace) various missing data values with `missing`.
This is useful for downstream imputation methods that assume missing data is represented by
a `missing`.

!!! In-place methods are only applicable for datasets which already `allowmissing`.

# Keyword Arguments
* `value::Tuple`: A tuple of values that should be considered `missing`

# Example
```jldoctest
julia> using Impute: DeclareMissings, apply

julia> M = [1.0 2.0 -9999.0 NaN 5.0; 1.1 2.2 3.3 0.0 5.5]
2×5 Matrix{Float64}:
 1.0  2.0  -9999.0  NaN    5.0
 1.1  2.2      3.3    0.0  5.5

julia> apply(M, DeclareMissings(; values=(NaN, -9999.0, 0.0)))
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5
```
"""
struct DeclareMissings{T<:Tuple}
    values::T
end

function DeclareMissings(; values)
    T = isa(values, Tuple) ? values : tuple(values)
    return DeclareMissings{typeof(T)}(T)
end

apply!(data::AbstractArray{Missing}, imp::DeclareMissings) = data

# Primary definition just calls `replace!`
function apply!(data::AbstractArray{Union{T, Missing}}, imp::DeclareMissings) where T
    # Reduce the possible set of values to those that could actually be found in the data
    # Useful, if we declare a `Replace` imputor that should be applied to multiple datasets.
    Base.replace!(data, (v => missing for v in imp.values if v isa T)...)
end

# Most of the time the in-place methods won't work because we need to change the
# eltype with allowmissing
apply(data::AbstractArray, imp::DeclareMissings) = apply!(allowmissing(data), imp)

# Custom implementation of a non-mutating impute for tables
function apply(table, imp::DeclareMissings)
    istable(table) || throw(MethodError(apply, (table, imp)))

    ctable = Tables.columns(table)

    cnames = Tuple(propertynames(ctable))
    cdata = Tuple(apply(getproperty(ctable, cname), imp) for cname in cnames)
    # Reconstruct as a ColumnTable
    result = NamedTuple{cnames}(cdata)

    # If our input was a ColumnTable just return the result. We can also do the same for
    if isa(table, Tables.ColumnTable)
        return result
    else
        return materializer(table)(result)
    end
end

# Specialcase for rowtable
function apply(data::T, imp::DeclareMissings) where T <: AbstractVector{<:NamedTuple}
    # We use columntable here so that we don't call `materialize` more often than needed.
    return materializer(data)(apply(Tables.columntable(data), imp))
end


text:
# This file was generated, do not modify it. # hide
folds(X::AbstractNode, nfolds) = node(XX->folds(XX, nfolds), X)
f = folds(X, 3)
f()


text:
#!/usr/bin/env julia
# Copyright 2017 Christian Henning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

workspace() # In case, this code block is copied to a notebook.

using MNIST
using PyPlot

include("network.jl")

# Specify dataset.
train_in, train_out = MNIST.traindata()
test_in, test_out = MNIST.testdata()
train_in /= 255
test_in /= 255
train_out = one_hot(Array{Int, 1}(train_out) + 1)
test_out = one_hot(Array{Int, 1}(test_out) + 1)

# Downsample data from 28x28 to 14x14.
if true
    train_in = reshape(train_in, 28, 28, size(train_in, 2))
    test_in = reshape(test_in, 28, 28, size(test_in, 2))

    train_in = train_in[2:2:end, 2:2:end, :]
    test_in = test_in[2:2:end, 2:2:end, :]

    train_in = reshape(train_in, 14*14, size(train_in, 3))
    test_in = reshape(test_in, 14*14, size(test_in, 3))
end

# Add weights for neural biases.
train_in = neural_bias_data_transform(train_in)
test_in = neural_bias_data_transform(test_in)

in_size = size(train_in)[1]
out_size = 10

# Define hyperparameters.
hparams = Hyperparameters()
hparams.λ = 1
hparams.γ = 0.99
hparams.αₛ = 2
hparams.αₚ = 2
hparams.rᵦ = 0.0
hparams.η = 0.9
hparams.num_epochs = 100
hparams.batch_size = 10
hparams.sample_size = 100
hparams.clip_gradients = true
hparams.clipping_value = 0.1
hparams.rate_decay = true
hparams.rate_decay_interval = 1000

# Weight initializer.
#initializer(dims...) = 1e-3 * (2*rand(dims) - 1) / prod(dims)
# We use Xavier initialization instead.

# Build network.
network = Network(hparams = hparams)

loss_params = LossParameters()
loss_params.sparsity = 0.5

layer0 = input_layer("in", in_size)
add_layer(network, layer0)

initializer = xavier_initializer(layer0.size, 300)
layer1 = fully_connected_layer("fc1", layer0, 300, loss_params=loss_params,
                               initializer=initializer)
add_regularizer(layer1, sparsity_regularizer)
add_regularizer(layer1, participation_regularizer)
add_layer(network, layer1)

initializer = xavier_initializer(layer1.size, 100)
layer2 = fully_connected_layer("fc2", layer1, 100, loss_params=loss_params,
                               initializer=initializer)
add_regularizer(layer2, sparsity_regularizer)
add_regularizer(layer2, participation_regularizer)
add_layer(network, layer2)

loss_params = LossParameters()
loss_params.sparsity = 0.9
initializer = xavier_initializer(layer2.size, out_size)
layer3 = fully_connected_layer("out", layer2, out_size, loss_params=loss_params,
                               initializer=initializer)
add_regularizer(layer3, sparsity_regularizer)
add_regularizer(layer3, participation_regularizer)
add_layer(network, layer3)

println("Predictions in the untrained network (to check initialization):")
println(predict(network, test_in[:,1:1]))
println(predict(network, test_in[:,2:2]))
println(predict(network, test_in[:,3:3]))
println()

# Train Network.
T, E = train_network(network, train_in, train_out = train_out,
                    test_in = test_in[:, 1:200],
                    test_out = test_out[:, 1:200],
                    eval_interval = 10)

#println("\nThe final participation trace in the output layer:")
#println(layer3.loss_memory["y_trace"])
#println()

# Plot training development.
title("Development of training loss")
xlabel("Epoch")
ylabel("Train Loss")
plot(1:length(T), T, label = "Train Loss")
legend()
show()

title("Development of test evaluation during training")
xlabel("Epoch")
ylabel("Test Eval")
plot(10:10:length(T), E, label="Test Accuracy")
axhline(y=0.0, color="k", linestyle="--")
axhline(y=1.0, color="k", linestyle="--")
legend()
ax = gca()
ax[:set_xlim]([0, length(T)])
show()

# Plot downsampling examples.
#subplot(221)
#axis("off")
#title("Original")
#imshow(reshape(test_in[1:end-1,1], 28, 28), vmin=-1.0, vmax=1.0)
#subplot(222)
#axis("off")
#title("Downsampled")
#imshow(reshape(test_in[1:end-1,1], 28, 28)[2:2:end, 2:2:end], vmin=-1.0,
#   vmax=1.0)
#
#subplot(223)
#axis("off")
#imshow(reshape(test_in[1:end-1,2], 28, 28), vmin=-1.0, vmax=1.0)
#subplot(224)
#axis("off")
#imshow(reshape(test_in[1:end-1,2], 28, 28)[2:2:end, 2:2:end], vmin=-1.0,
#   vmax=1.0)
#
#suptitle("Downsampling example")
#show()

for ii = 1:16
    pred = predict(network, test_in[:,ii:ii])

    subplot(4,8,(ii-1)*2 + 1)
    imshow(reshape(test_in[1:end-1,ii], 14, 14), vmin=-1.0, vmax=1.0)
    axis("off")

    subplot(4,8,(ii-1)*2 + 2)
    bar(1:10, squeeze(pred, 2), align="center")
    ax = gca()
    ax[:xaxis][:set_ticks](1:10)
end

title("Example of input (left) and output response (right)")
show()


text:
# Interacting with the Petsc Options Database
export PetscOptionsSetValue, PetscOptionsClearValue, PetscOptionsView, PetscSetOptions, PetscClearOptions

"""
  Typedef of PetscOptions
"""
const PetscOptions = Ptr{Void}

"""
  PetscOptionsSetValue

  **Inputs**

   * arg1: the key (string)
   * arg2: the value (string)
   * arg3: the PetscOptions object, defaults to the global options database
"""
function PetscOptionsSetValue(arg1::AbstractString,arg2::AbstractString, arg3::PetscOptions=C_NULL)
    ccall((:PetscOptionsSetValue,petsc),PetscErrorCode,(PetscOptions, Cstring, Cstring),arg3, arg1,arg2)
end

"""
  PetscOptionsView

  **Inputs**

   * arg1: a PetscViewer, defaults to Petsc's stdout
   * arg2: the PetscOptions object, defaults to the global options databse
"""
function PetscOptionsView(arg1::PetscViewer=C_NULL, arg2::PetscOptions=C_NULL)
    ccall((:PetscOptionsView,petsc),PetscErrorCode,(PetscOptions, PetscViewer,),arg2, arg1)
end

"""
  PetscOptionsClearValue

  **Inputs**

   * arg1: the key (string)
   * arg2: the PetscOptions object, defaults to the global options databse
"""
function PetscOptionsClearValue(arg1::AbstractString, arg3::PetscOptions=C_NULL)
    ccall((:PetscOptionsClearValue,petsc),PetscErrorCode,(PetscOptions, Cstring,),arg3, arg1)
end

"""
  Convenience wrapper for using a dictionary to set options
"""
function PetscSetOptions(opts::Dict)

  for (key, value) in opts
    PetscOptionsSetValue(key, value)
  end

end

"""
  Convenience wrapper for using a dictionary to clear options (only the keys
  are used).
"""
function PetscClearOptions(opts::Dict)

  for key in keys(opts)
    PetscOptionsClearValue(key)
  end

end


text:
module CBLS

using Constraints: sequential_tasks
using Constraints
using ConstraintDomains
using JuMP
using Lazy
using LocalSearchSolvers
using MathOptInterface

# Const
const LS = LocalSearchSolvers
const MOI = MathOptInterface
const MOIU = MOI.Utilities

# MOI functions
const SVF = MOI.SingleVariable
const VOV = MOI.VectorOfVariables
const OF = MOI.ObjectiveFunction

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

# MOI types
const VAR_TYPES = Union{MOI.ZeroOne, MOI.Integer}

# Export: domain
export DiscreteSet

# Export: Constraints
export AllDifferent
export AllEqual
export AllEqualParam
export AlwaysTrue
export DistDifferent
export Eq
export Error
export LessThanParam
export MinusEqualParam
export Ordered
export Predicate
export SequentialTasks
export SumEqualParam

#Export: Scalar objective function
export ScalarFunction

# Include
include("MOI_wrapper.jl")
include("attributes.jl")
include("variables.jl")
include("constraints.jl")
include("objectives.jl")
include("results.jl")

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -16.957145096415093 <= q <= 17.94995800016443)

                     add_NL_constraint(m, :(gelu(0.510608367415756 + 0.13576703913191768*gelu(-0.023366318069798275 + -0.6972645804796413*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.3585964327428499*gelu(-0.17574028037161238 + -0.933468976339102*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.2511872166134532*gelu(0.9643857894710246 + 0.9348568947078797*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) + gelu(-0.7493006135688143 + 0.2483293387413994*gelu(-0.023366318069798275 + -0.6972645804796413*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + -0.7456041150082044*gelu(-0.17574028037161238 + -0.933468976339102*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.49415154189847854*gelu(0.9643857894710246 + 0.9348568947078797*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) + gelu(0.3105901915386422 + 0.5384893864274902*gelu(-0.023366318069798275 + -0.6972645804796413*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.5405413918049242*gelu(-0.17574028037161238 + -0.933468976339102*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + -0.18067211353899282*gelu(0.9643857894710246 + 0.9348568947078797*gelu(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*gelu(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*gelu(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function summarize_regression_test(field_names, fields, correct_fields)
    for (field_name, φ, φ_c) in zip(field_names, fields, correct_fields)
        Δ = Array(φ) .- φ_c

        Δ_min      = minimum(Δ)
        Δ_max      = maximum(Δ)
        Δ_mean     = mean(Δ)
        Δ_abs_mean = mean(abs, Δ)
        Δ_std      = std(Δ)

        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, Δ_min, Δ_max, Δ_mean, Δ_abs_mean, Δ_std))
    end
end

include("regression_tests/thermal_bubble_regression_test.jl")
include("regression_tests/rayleigh_benard_regression_test.jl")
include("regression_tests/ocean_large_eddy_simulation_regression_test.jl")

@testset "Regression" begin
    println("Running regression tests...")

    for arch in archs
        @testset "Thermal bubble [$(typeof(arch))]" begin
            println("  Testing thermal bubble regression [$(typeof(arch))]")
            run_thermal_bubble_regression_test(arch)
        end

        @testset "Rayleigh–Bénard tracer [$(typeof(arch))]" begin
            println("  Testing Rayleigh–Bénard tracer regression [$(typeof(arch))]")
            run_rayleigh_benard_regression_test(arch)
        end

        @testset "Ocean large eddy simulation [$(typeof(arch))]" begin
            for closure in (AnisotropicMinimumDissipation(), ConstantSmagorinsky())
                closurename = string(typeof(closure).name.wrapper)
                println("  Testing oceanic large eddy simulation regression [$closurename, $(typeof(arch))]")
                run_ocean_large_eddy_simulation_regression_test(arch, closure)
            end
        end
    end
end


text:
using Flux
using Flux: @epochs
mutable struct NeuralSurrogate{X,Y,M,L,O,P,N,A,U} <: AbstractSurrogate
    x::X
    y::Y
    model::M
    loss::L
    opt::O
    ps::P
    n_echos::N
    lb::A
    ub::U
 end

function NeuralSurrogate(x,y,lb::Number,ub::Number,model,loss,opt,n_echos)
    X = reshape(x,length(x),1)
    data = Tuple{Array{eltype(X[1]),1},eltype(y[1])}[]
    for i in 1:length(x)
        push!(data, ([X[i]], y[i]))
    end
    ps = Flux.params(model)
    @epochs n_echos Flux.train!(loss, ps, data, opt)
    return NeuralSurrogate(x,y,model,loss,opt,ps,n_echos,lb,ub)
end

function (my_neural::NeuralSurrogate)(val::Number)
    return Tracker.data(first(my_neural.model([val])))
end

function add_point!(my_n::NeuralSurrogate,x_new,y_new)
    if length(my_n.lb) == 1
        #1D
        X = reshape(my_n.x,length(my_n.x),1)
        data = Tuple{Array{eltype(X[1]),1},eltype(my_n.y[1])}[]
        for i in 1:length(my_n.x)
            push!(data, ([X[i]], my_n.y[i]))
        end
        for j = 1:length(x_new)
            push!(data,([x_new[j]],y_new[j]))
        end
        my_n.x = vcat(my_n.x,x_new)
        my_n.y = vcat(my_n.y,y_new)
        @epochs my_n.n_echos Flux.train!(my_n.loss, my_n.ps, data, my_n.opt)
    else
        #ND

        n_previous = length(my_n.x)
        a = vcat(my_n.x,x_new)
        n_after = length(a)
        dim_new = n_after - n_previous
        n = length(my_n.x)
        d = length(my_n.x[1])
        tot_dim = n + dim_new
        X = Array{Float64,2}(undef,tot_dim,d)
        data = Tuple{Array{eltype(X[1]),1},eltype(my_n.y[1])}[]
        for j = 1:n
            X[j,:] = vec(collect(my_n.x[j]))
            push!(data, (X[j,:], my_n.y[j]))
        end
        if dim_new == 1
            X[n+1,:] = vec(collect(x_new))
            push!(data, (X[n+1,:], y_new))
        else
            i = 1
            for j = n+1:tot_dim
                X[j,:] = vec(collect(x_new[i]))
                push!(data,(X[j,:], y_new[i]))
                i = i + 1
            end
        end
        my_n.x = vcat(my_n.x,x_new)
        my_n.y = vcat(my_n.y,y_new)
        @epochs my_n.n_echos Flux.train!(my_n.loss, my_n.ps, data, my_n.opt)
    end
    nothing
end

"""
NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos)

- model: Flux layers
- loss: loss function
- opt: optimization function

"""
function NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos)
    X = Array{Float64,2}(undef,length(x),length(x[1]))
    for j = 1:length(x)
        X[j,:] = vec(collect(x[j]))
    end
    data = Tuple{Array{eltype(X[1]),1},eltype(y[1])}[]
    for i in 1:size(x,1)
        push!(data, (X[i,:], y[i]))
    end
    ps = Flux.params(model)
    @epochs n_echos Flux.train!(loss, ps, data, opt)
    return NeuralSurrogate(x,y,model,loss,opt,ps,n_echos,lb,ub)

end

function (my_neural::NeuralSurrogate)(val)
    v = collect(val)
    return Tracker.data(first(my_neural.model(reshape(v,length(v),1))))
end


text:
"""
    MatchingResult{U}

A type representing the result of a matching algorithm.

# Fields

- `weight::U`: total weight of the matching
- `mate::Vector{Int}`: pairwise assignment.

`mate[i] = j` if vertex `i` is matched to vertex `j`, and `mate[i] = -1` for unmatched vertices.
"""
struct MatchingResult{U<:Real}
    weight::U
    mate::Vector{Int}
end

## Maximum weight matching

"""
    maximum_weight_matching(g, w; optimizer)

Given a graph `g` and an matrix `w` of edge weights, return a matching  ([`MatchingResult`](@ref) object) with the maximum total weight.

If no weight matrix is given, all edges will be considered to have weight 1
(results in max cardinality matching).
Edges in `g` that are not present in `w` will not be considered for the matching.

A JuMP-compatible solver must be provided with the `optimizer` argument.

The efficiency of the algorithm depends on the input graph:
  - If the graph is bipartite, then the LP relaxation is integral.
  - If the graph is not bipartite, then it requires a MIP solver and the computation time may grow exponentially.
"""
function maximum_weight_matching(
    g::Graph, w::AbstractMatrix{U}=default_weights(g); optimizer
) where {U<:Real}
    model = Model(optimizer)
    n = nv(g)
    edge_list = collect(edges(g))

    # put the edge weights in w in the right order to be compatible with edge_list
    for j in 1:n
        for i in 1:n
            if i > j && w[i, j] > zero(U) && w[j, i] < w[i, j]
                w[j, i] = w[i, j]
            end
            if Edge(i, j) ∉ edge_list
                w[i, j] = zero(U)
            end
        end
    end

    if is_bipartite(g)
        @variable(model, x[edge_list] >= 0) # no need to enforce integrality
    else
        @variable(model, x[edge_list] >= 0, Int) # requires MIP solver
    end
    @objective(model, Max, sum(x[e] * w[src(e), dst(e)] for e in edge_list))

    @constraint(model, c1[i=1:n], sum(x[Edge(minmax(i, j))] for j in neighbors(g, i)) <= 1)
    optimize!(model)
    status = JuMP.termination_status(model)
    status != MOI.OPTIMAL && error("JuMP solver failed to find optimal solution.")
    solution = value.(x)
    cost = objective_value(model)

    mate = fill(-1, n)
    for e in edge_list
        if solution[e] >= 1 - 1e-5 # Some tolerance to numerical approximations by the solver.
            mate[src(e)] = dst(e)
            mate[dst(e)] = src(e)
        end
    end

    return MatchingResult(cost, mate)
end

## Maximum weight maximal matching

"""
    maximum_weight_maximal_matching(g, w[, cutoff]; optimizer)

Given a bipartite graph `g` and a matrix `w` of edge weights, return a matching ([`MatchingResult`](@ref) object) with the maximum total weight among the ones that contain the largest number of edges.

If no weight matrix is given, all edges will be considered to have weight 1
(results in max cardinality matching).
Edges in `g` that are not present in `w` will not be considered for the matching.

A JuMP-compatible solver must be provided with the `optimizer` argument.

A `cutoff` argument can be given to reduce computation time by excluding edges with weights lower than the specified value.

The algorithm relies on a linear relaxation on of the matching problem, which is
guaranteed to have integer solution on bipartite graphs.
"""
function maximum_weight_maximal_matching(
    g::Graph, w::AbstractMatrix{T}=default_weights(g); optimizer
) where {T<:Real}
    # TODO support for graphs with zero degree nodes
    # TODO apply separately on each connected component
    bpmap = bipartite_map(g)
    length(bpmap) != nv(g) && error("Graph is not bipartite")
    v1 = findall(isequal(1), bpmap)
    v2 = findall(isequal(2), bpmap)
    if length(v1) > length(v2)
        v1, v2 = v2, v1
    end

    nedg = 0
    edgemap = Dict{Edge,Int}()

    for j in 1:size(w, 2)
        for i in 1:size(w, 1)
            if w[i, j] > 0.0
                nedg += 1
                edgemap[Edge(i, j)] = nedg
                edgemap[Edge(j, i)] = nedg
            end
        end
    end

    model = Model(optimizer)
    @variable(model, x[1:length(w)] >= 0)

    for i in v1
        idx = Vector{Int}()
        for j in neighbors(g, i)
            if haskey(edgemap, Edge(i, j))
                push!(idx, edgemap[Edge(i, j)])
            end
        end
        if length(idx) > 0
            @constraint(model, sum(x[id] for id in idx) == 1)
        end
    end

    for j in v2
        idx = Vector{Int}()
        for i in neighbors(g, j)
            if haskey(edgemap, Edge(i, j))
                push!(idx, edgemap[Edge(i, j)])
            end
        end

        if length(idx) > 0
            @constraint(model, sum(x[id] for id in idx) <= 1)
        end
    end

    @objective(model, Max, sum(w[src(e), dst(e)] * x[edgemap[e]] for e in keys(edgemap)))

    optimize!(model)
    status = JuMP.termination_status(model)
    status != MOI.OPTIMAL && error("JuMP solver failed to find optimal solution.")
    sol = JuMP.value.(x)

    all(Bool[s == 1 || s == 0 for s in sol]) || error("Found non-integer solution.")

    cost = JuMP.objective_value(model)

    mate = fill(-1, nv(g))
    for e in edges(g)
        if w[src(e), dst(e)] > zero(T)
            inmatch = convert(Bool, sol[edgemap[e]])
            if inmatch
                mate[src(e)] = dst(e)
                mate[dst(e)] = src(e)
            end
        end
    end

    return MatchingResult(cost, mate)
end

"""
    cutoff_weights(w, cutoff)

Copy the weights matrix `w` with all elements below `cutoff` set to 0.
"""
function cutoff_weights(w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
    wnew = copy(w)
    for j in 1:size(w,2)
        for i in 1:size(w,1)
            if wnew[i,j] < cutoff
                wnew[i,j] = zero(T)
            end
        end
    end
    wnew
end

function maximum_weight_maximal_matching(
    g::Graph, w::AbstractMatrix{T}, cutoff::R; optimizer
) where {T<:Real,R<:Real}
    return maximum_weight_maximal_matching(g, cutoff_weights(w, cutoff); optimizer)
end


text:
using .Tracker: Tracker,
                TrackedReal,
                TrackedVector,
                TrackedMatrix,
                TrackedArray,
                TrackedVecOrMat,
                @grad,
                track,
                data,
                param

using Compat: eachcol
using LinearAlgebra

maporbroadcast(f, x::TrackedArray...) = f.(x...)
function maporbroadcast(
    f,
    x1::TrackedArray{T, N},
    x::AbstractArray{<:TrackedReal}...,
) where {T, N}
    return f.(convert(Array{TrackedReal{T}, N}, x1), x...)
end

_eps(::Type{<:TrackedReal{T}}) where {T} = _eps(T)
function Base.minimum(d::LocationScale{<:TrackedReal})
    m = minimum(d.ρ)
    if isfinite(m)
        return d.μ + d.σ * m
    else
        return m
    end
end
function Base.maximum(d::LocationScale{<:TrackedReal})
    m = maximum(d.ρ)
    if isfinite(m)
        return d.μ + d.σ * m
    else
        return m
    end
end

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:TrackerAD}, Inverse{<:ADBijector{<:TrackerAD}}},
    x::Real
)
    return data(Tracker.gradient(b, x)[1])
end
function jacobian(
    b::Union{<:ADBijector{<:TrackerAD}, Inverse{<:ADBijector{<:TrackerAD}}},
    x::AbstractVector{<:Real}
)
    # We extract `data` so that we don't return a `Tracked` type
    return data(Tracker.jacobian(b, x))
end

# implementations for Shift bijector
function _logabsdetjac_shift(a::TrackedReal, x::Real, ::Val{0})
    return tracker_shift_logabsdetjac(a, x, Val(0))
end
function _logabsdetjac_shift(a::TrackedReal, x::AbstractVector{<:Real}, ::Val{0})
    return tracker_shift_logabsdetjac(a, x, Val(0))
end
function _logabsdetjac_shift(
    a::Union{TrackedReal, TrackedVector{<:Real}},
    x::AbstractVector{<:Real},
    ::Val{1}
)
    return tracker_shift_logabsdetjac(a, x, Val(1))
end
function _logabsdetjac_shift(
    a::Union{TrackedReal, TrackedVector{<:Real}},
    x::AbstractMatrix{<:Real},
    ::Val{1}
)
    return tracker_shift_logabsdetjac(a, x, Val(1))
end
function tracker_shift_logabsdetjac(a, x, ::Val{N}) where {N}
    return param(_logabsdetjac_shift(data(a), data(x), Val(N)))
end

# Log bijector

@grad function logabsdetjac(b::Log{1}, x::AbstractVector)
    return -sum(log, data(x)), Δ -> (nothing, -Δ ./ data(x))
end
@grad function logabsdetjac(b::Log{1}, x::AbstractMatrix)
    return -vec(sum(log, data(x); dims = 1)), Δ -> (nothing, .- Δ' ./ data(x))
end
@grad function logabsdetjac(b::Log{2}, x::AbstractMatrix)
    return -sum(log, data(x)), Δ -> (nothing, -Δ ./ data(x))
end

# implementations for Scale bijector
# Adjoints for 0-dim and 1-dim `Scale` using `Real`
function _logabsdetjac_scale(a::TrackedReal, x::Real, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::Real, ::Val{0})
    return _logabsdetjac_scale(data(a), data(x), Val(0)), Δ -> (inv(data(a)) .* Δ, nothing, nothing)
end
# Need to treat `AbstractVector` and `AbstractMatrix` separately due to ambiguity errors
function _logabsdetjac_scale(a::TrackedReal, x::AbstractVector, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0})
    da = data(a)
    J = fill(inv.(da), length(x))
    return _logabsdetjac_scale(da, data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedReal, x::AbstractMatrix, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{0})
    da = data(a)
    J = fill(size(x, 1) / da, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
# adjoints for 1-dim and 2-dim `Scale` using `AbstractVector`
function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    return track(_logabsdetjac_scale, a, data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    da = data(a)
    J = inv.(da)
    return _logabsdetjac_scale(da, data(x), Val(1)), Δ -> (J .* Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    return track(_logabsdetjac_scale, a, data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    da = data(a)
    Jᵀ = repeat(inv.(da), 1, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end
# TODO: implement analytical gradient for scaling a vector using a matrix
# function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     track(_logabsdetjac_scale, a, data(x), Val{1})
# end
# @grad function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     throw
# end
# implementations for Stacked bijector
function logabsdetjac(b::Stacked, x::TrackedMatrix{<:Real})
    return map(eachcol(x)) do c
        logabsdetjac(b, c)
    end
end
# TODO: implement custom adjoint since we can exploit block-diagonal nature of `Stacked`
function (sb::Stacked)(x::TrackedMatrix{<:Real})
    return eachcolmaphcat(sb, x)
end
# Simplex adjoints
function _simplex_bijector(X::TrackedVecOrMat, b::SimplexBijector{1})
    return track(_simplex_bijector, X, b)
end
function _simplex_inv_bijector(Y::TrackedVecOrMat, b::SimplexBijector{1})
    return track(_simplex_inv_bijector, Y, b)
end
@grad function _simplex_bijector(X::AbstractVector, b::SimplexBijector{1})
    Xd = data(X)
    return _simplex_bijector(Xd, b), Δ -> (simplex_link_jacobian(Xd)' * Δ, nothing)
end
@grad function _simplex_inv_bijector(Y::AbstractVector, b::SimplexBijector{1})
    Yd = data(Y)
    return _simplex_inv_bijector(Yd, b), Δ -> (simplex_invlink_jacobian(Yd)' * Δ, nothing)
end

@grad function _simplex_bijector(X::AbstractMatrix, b::SimplexBijector{1})
    Xd = data(X)
    return _simplex_bijector(Xd, b), Δ -> begin
        maphcat(eachcol(Xd), eachcol(Δ)) do c1, c2
            simplex_link_jacobian(c1)' * c2
        end, nothing
    end
end
@grad function _simplex_inv_bijector(Y::AbstractMatrix, b::SimplexBijector{1})
    Yd = data(Y)
    return _simplex_inv_bijector(Yd, b), Δ -> begin
        maphcat(eachcol(Yd), eachcol(Δ)) do c1, c2
            simplex_invlink_jacobian(c1)' * c2
        end, nothing
    end
end

replace_diag(::typeof(log), X::TrackedMatrix) = track(replace_diag, log, X)
@grad function replace_diag(::typeof(log), X)
    Xd = data(X)
    f(i, j) = i == j ? log(Xd[i, j]) : Xd[i, j]
    out = f.(1:size(Xd, 1), (1:size(Xd, 2))')
    out, ∇ -> begin
        g(i, j) = i == j ? ∇[i, j]/Xd[i, j] : ∇[i, j]
        return (nothing, g.(1:size(Xd, 1), (1:size(Xd, 2))'))
    end
end

replace_diag(::typeof(exp), X::TrackedMatrix) = track(replace_diag, exp, X)
@grad function replace_diag(::typeof(exp), X)
    Xd = data(X)
    f(i, j) = ifelse(i == j, exp(Xd[i, j]), Xd[i, j])
    out = f.(1:size(Xd, 1), (1:size(Xd, 2))')
    out, ∇ -> begin
        g(i, j) = ifelse(i == j, ∇[i, j]*exp(Xd[i, j]), ∇[i, j])
        return (nothing, g.(1:size(Xd, 1), (1:size(Xd, 2))'))
    end
end

logabsdetjac(b::SimplexBijector{1}, x::TrackedVecOrMat) = track(logabsdetjac, b, x)
@grad function logabsdetjac(b::SimplexBijector{1}, x::AbstractVector)
    xd = data(x)
    return logabsdetjac(b, xd), Δ -> begin
        (nothing, simplex_logabsdetjac_gradient(xd) * Δ)
    end
end
@grad function logabsdetjac(b::SimplexBijector{1}, x::AbstractMatrix)
    xd = data(x)
    return logabsdetjac(b, xd), Δ -> begin
        (nothing, maphcat(eachcol(xd), Δ) do c, g
            simplex_logabsdetjac_gradient(c) * g
        end)
    end
end

for header in [
    (:(α_::TrackedReal), :β, :z_0, :(z::AbstractVector)),
    (:α_, :(β::TrackedReal), :z_0, :(z::AbstractVector)),
    (:α_, :β, :(z_0::TrackedVector), :(z::AbstractVector)),
    (:α_, :β, :z_0, :(z::TrackedVector)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::AbstractVector)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::AbstractVector)),
    (:(α_::TrackedReal), :β, :z_0, :(z::TrackedVecOrMat)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::AbstractVector)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::TrackedVector)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::TrackedVector)),
    (:α_, :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedVector)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedVector)),
]
    @eval begin
        function _radial_transform($(header...))
            α = softplus(α_)            # from A.2
            β_hat = -α + softplus(β)    # from A.2
            if β_hat isa TrackedReal
                TV = vectorof(typeof(β_hat))
                T = vectorof(typeof(β_hat))
            elseif z_0 isa TrackedVector
                TV = typeof(z_0)
                T = typeof(z_0)
            else
                T = TV = typeof(z)
            end
            Tr = promote_type(eltype(z), eltype(z_0))
            r::Tr = norm((z .- z_0)::TV)
            transformed::T = z .+ β_hat ./ (α .+ r') .* (z .- z_0)   # from eq(14)
            return (transformed = transformed, α = α, β_hat = β_hat, r = r)
        end
    end
end

for header in [
    (:(α_::TrackedReal), :β, :z_0, :(z::AbstractMatrix)),
    (:α_, :(β::TrackedReal), :z_0, :(z::AbstractMatrix)),
    (:α_, :β, :(z_0::TrackedVector), :(z::AbstractMatrix)),
    (:α_, :β, :z_0, :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::AbstractMatrix)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::AbstractMatrix)),
    (:(α_::TrackedReal), :β, :z_0, :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::AbstractMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::TrackedMatrix)),
    (:α_, :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedMatrix)),
]
    @eval begin
        function _radial_transform($(header...))
            α = softplus(α_)            # from A.2
            β_hat = -α + softplus(β)    # from A.2
            if β_hat isa TrackedReal
                TV = vectorof(typeof(β_hat))
                T = matrixof(TV)
            elseif z_0 isa TrackedVector
                TV = typeof(z_0)
                T = matrixof(TV)
            else
                T = typeof(z)
                TV = vectorof(T)
            end
            r::TV = eachcolnorm(z .- z_0)
            transformed::T = z .+ β_hat ./ (α .+ r') .* (z .- z_0)   # from eq(14)
            return (transformed = transformed, α = α, β_hat = β_hat, r = r)
        end
    end
end
eachcolnorm(X) = map(norm, eachcol(X))
eachcolnorm(X::TrackedMatrix) = track(eachcolnorm, X)
@grad function eachcolnorm(X)
    Xd = data(X)
    y = map(norm, eachcol(Xd))
    y, Δ -> begin
        (Xd .* (Δ ./ y)',)
    end
end

function matrixof(::Type{<:TrackedArray{T,1,Vector{T}}}) where {T<:Real}
    return TrackedArray{T,2,Matrix{T}}
end
function matrixof(::Type{TrackedReal{T}}) where {T<:Real}
    return TrackedArray{T,2,Matrix{T}}
end
function vectorof(::Type{<:TrackedArray{T,2,Matrix{T}}}) where {T<:Real}
    return TrackedArray{T,1,Vector{T}}
end
function vectorof(::Type{TrackedReal{T}}) where {T<:Real}
    return TrackedArray{T,1,Vector{T}}
end

(b::Exp{0})(x::TrackedVector) = exp.(x)::vectorof(float(eltype(x)))
(b::Exp{1})(x::TrackedVector) = exp.(x)::vectorof(float(eltype(x)))
(b::Exp{1})(x::TrackedMatrix) = exp.(x)::matrixof(float(eltype(x)))
(b::Exp{2})(x::TrackedMatrix) = exp.(x)::matrixof(float(eltype(x)))

(b::Log{0})(x::TrackedVector) = log.(x)::vectorof(float(eltype(x)))
(b::Log{1})(x::TrackedVector) = log.(x)::vectorof(float(eltype(x)))
(b::Log{1})(x::TrackedMatrix) = log.(x)::matrixof(float(eltype(x)))
(b::Log{2})(x::TrackedMatrix) = log.(x)::matrixof(float(eltype(x)))

logabsdetjac(b::Log{0}, x::TrackedVector) = .-log.(x)::vectorof(float(eltype(x)))
logabsdetjac(b::Log{1}, x::TrackedMatrix) = - vec(sum(log.(x); dims = 1))

getpd(X::TrackedMatrix) = track(getpd, X)
@grad function getpd(X::AbstractMatrix)
    Xd = data(X)
    return LowerTriangular(Xd) * LowerTriangular(Xd)', Δ -> begin
        Xl = LowerTriangular(Xd)
        return (LowerTriangular(Δ' * Xl + Δ * Xl),)
    end
end

lower(A::TrackedMatrix) = track(lower, A)
@grad function lower(A::AbstractMatrix)
    Ad = data(A)
    return lower(Ad), Δ -> (lower(Δ),)
end

_inv_link_chol_lkj(y::TrackedMatrix) = track(_inv_link_chol_lkj, y)
@grad function _inv_link_chol_lkj(y_tracked)
    y = data(y_tracked)

    K = LinearAlgebra.checksquare(y)

    w = similar(y)

    z_mat = similar(y) # cache for adjoint
    tmp_mat = similar(y)

    @inbounds for j in 1:K
        w[1, j] = 1
        for i in 2:j
            z = tanh(y[i-1, j])
            tmp = w[i-1, j]

            z_mat[i, j] = z
            tmp_mat[i, j] = tmp

            w[i-1, j] = z * tmp
            w[i, j] = tmp * sqrt(1 - z^2)
        end
        for i in (j+1):K
            w[i, j] = 0
        end
    end

    function pullback_inv_link_chol_lkj(Δw)
        LinearAlgebra.checksquare(Δw)

        Δy = zero(y)

        @inbounds for j in 1:K
            Δtmp = Δw[j,j]
            for i in j:-1:2
                Δz = Δw[i-1, j] * tmp_mat[i, j] - Δtmp * tmp_mat[i, j] / sqrt(1 - z_mat[i, j]^2) * z_mat[i, j]
                Δy[i-1, j] = Δz / cosh(y[i-1, j])^2
                Δtmp = Δw[i-1, j] * z_mat[i, j] + Δtmp * sqrt(1 - z_mat[i, j]^2)
            end
        end

        return (Δy,)
    end

    return w, pullback_inv_link_chol_lkj
end

_link_chol_lkj(w::TrackedMatrix) = track(_link_chol_lkj, w)
@grad function _link_chol_lkj(w_tracked)
    w = data(w_tracked)

    K = LinearAlgebra.checksquare(w)

    z = similar(w)

    @inbounds z[1, 1] = 0

    tmp_mat = similar(w) # cache for pullback.

    @inbounds for j=2:K
        z[1, j] = atanh(w[1, j])
        tmp = sqrt(1 - w[1, j]^2)
        tmp_mat[1, j] = tmp
        for i in 2:(j - 1)
            p = w[i, j] / tmp
            tmp *= sqrt(1 - p^2)
            tmp_mat[i, j] = tmp
            z[i, j] = atanh(p)
        end
        z[j, j] = 0
    end

    function pullback_link_chol_lkj(Δz)
        LinearAlgebra.checksquare(Δz)

        Δw = similar(w)

        @inbounds Δw[1,1] = zero(eltype(Δz))

        @inbounds for j=2:K
            Δw[j, j] = 0
            Δtmp = zero(eltype(Δz)) # Δtmp_mat[j-1,j]
            for i in (j-1):-1:2
                p = w[i, j] / tmp_mat[i-1, j]
                ftmp = sqrt(1 - p^2)
                d_ftmp_p = -p / ftmp
                d_p_tmp = -w[i,j] / tmp_mat[i-1, j]^2

                Δp = Δz[i,j] / (1-p^2) + Δtmp * tmp_mat[i-1, j] * d_ftmp_p
                Δw[i, j] = Δp / tmp_mat[i-1, j]
                Δtmp = Δp * d_p_tmp + Δtmp * ftmp # update to "previous" Δtmp
            end
            Δw[1, j] = Δz[1, j] / (1-w[1,j]^2) - Δtmp / sqrt(1 - w[1,j]^2) * w[1,j]
        end

        return (Δw,)
    end

    return z, pullback_link_chol_lkj
end