text:
using ArgParse
include("clustering.jl")
include("datasets.jl")

settings = ArgParseSettings()
@add_arg_table settings begin
    "--problem"
    arg_type = String
    default = "iris"
    "--fname"
    arg_type = String
    default = "lif"
    "--seed"
    arg_type = Int
    default = 0
    "--logfile"
    arg_type = String
    default = "stdp.log"
    "--train_epochs"
    arg_type = Int
    default = 1
    "--weight_mean"
    arg_type = Float64
    default = 0.5
    "--weight_std"
    arg_type = Float64
    default = 0.1
    "--t_train"
    arg_type = Int
    default = 350
    "--t_blank"
    arg_type = Int
    default = 150
    "--fr"
    arg_type = Float64
    default = 65.0
    "--pre_target"
    arg_type = Float64
    default = 0.4
    "--stdp_lr"
    arg_type = Float64
    default = 0.0001
    "--stdp_mu"
    arg_type = Float64
    default = 2.0
    "--inhib_weight"
    arg_type = Float64
    default = 0.1
    "--ika"
    arg_type = Float64
    default = 0.02
    "--ikb"
    arg_type = Float64
    default = 0.2
    "--ikd"
    arg_type = Float64
    default = 2.0
end
settings

args = parse_args(settings)
Logging.configure(filename=args["logfile"], level=INFO)

X, Y = get_data(args["problem"])
n_cluster = length(unique(Y))

stdp_labels = stdp_cluster(
    X, Y, n_cluster, x->x; seed=args["seed"], logfile=args["logfile"],
    problem=args["problem"], fname=args["fname"],
    train_epochs=args["train_epochs"], weight_mean=args["weight_mean"],
    weight_std=args["weight_std"], t_train=args["t_train"],
    t_blank=args["t_blank"], fr=args["fr"], pre_target=args["pre_target"],
    stdp_lr=args["stdp_lr"], stdp_mu=args["stdp_mu"],
    inhib_weight=args["inhib_weight"], ika=args["ika"],
    ikb=args["ikb"], ikd=args["ikd"])

acc = randindex(stdp_labels, Y)

Logging.info(@sprintf("E%0.6f", -acc[1]))


text:
@defstruct MemoryDataLayer Layer (
  name :: AbstractString = "memory-data",
  (tops :: Vector{Symbol} = Symbol[:data,:label], length(tops) > 0),
  (batch_size :: Int = 0, batch_size > 0),
  (data :: Vector{Array} = Array[], length(data) == length(tops)),
  shuffle :: Bool = false,
  transformers :: Vector = [],
)
@characterize_layer(MemoryDataLayer,
  is_source => true
)

type MemoryDataLayerState <: LayerState
  layer :: MemoryDataLayer
  blobs :: Vector{Blob}
  epoch :: Int
  trans :: Vector{Vector{DataTransformerState}}

  curr_idx :: Int
  shuffle_idx :: Vector{Int}

  MemoryDataLayerState(backend::Backend, layer::MemoryDataLayer) = begin
    blobs = Array{Blob}(length(layer.tops))
    trans = Array{Vector{DataTransformerState}}(length(layer.tops))
    transformers = convert(Vector{@compat(Tuple{Symbol, DataTransformerType})}, layer.transformers)
    for i = 1:length(blobs)
      dims = tuple(size(layer.data[i])[1:end-1]..., layer.batch_size)
      idxs = map(x -> 1:x, dims)

      blobs[i] = make_blob(backend, eltype(layer.data[i]), dims...)
      trans[i] = [setup(backend, convert(DataTransformerType, t), blobs[i])
          for (k,t) in filter(kt -> kt[1] == layer.tops[i], transformers)]
    end

    if layer.shuffle
      shuffle_idx = randperm(size(layer.data[1])[end])
    else
      shuffle_idx = Int[]
    end

    new(layer, blobs, 0, trans, 1, shuffle_idx)
  end
end

function setup(backend::Backend, layer::MemoryDataLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  @assert length(inputs) == 0
  for i = 2:length(layer.data)
    @assert eltype(layer.data[i]) == eltype(layer.data[1])
  end

  state = MemoryDataLayerState(backend, layer)
  return state
end
function shutdown(backend::Backend, state::MemoryDataLayerState)
  map(destroy, state.blobs)
  map(ts -> map(t -> shutdown(backend, t), ts), state.trans)
end

function forward(backend::Backend, state::MemoryDataLayerState, inputs::Vector{Blob})
  n_done = 0
  while n_done < state.layer.batch_size
    n_remain = size(state.layer.data[1])[end] - state.curr_idx + 1
    if n_remain == 0
      state.curr_idx = 1
      n_remain = size(state.layer.data[1])[end]

      if state.layer.shuffle
        state.shuffle_idx = randperm(size(state.layer.data[1])[end])
      end
    end

    n1 = min(state.layer.batch_size - n_done, n_remain)
    for i = 1:length(state.blobs)
      dset = state.layer.data[i]
      idx = map(x -> 1:x, size(state.blobs[i])[1:end-1])
      if state.layer.shuffle
        the_data = dset[idx..., state.shuffle_idx[state.curr_idx:state.curr_idx+n1-1]]
      else
        the_data = dset[idx..., state.curr_idx:state.curr_idx+n1-1]
      end
      set_blob_data(the_data, state.blobs[i], n_done+1)
    end
    state.curr_idx += n1
    n_done += n1

    if state.curr_idx > size(state.layer.data[1])[end]
      state.epoch += 1
    end
  end

  for i = 1:length(state.blobs)
    for j = 1:length(state.trans[i])
      forward(backend, state.trans[i][j], state.blobs[i])
    end
  end
end

function backward(backend::Backend, state::MemoryDataLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end

ext: jl


text:
# Autogenerated wrapper script for OpenLSTO_jll for x86_64-w64-mingw32-cxx11
export m2do_3d_lsm, m2do_fea, m2do_lsm

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenLSTO")
JLLWrappers.@declare_library_product(m2do_3d_lsm, "m2do_3d_lsm.dll")
JLLWrappers.@declare_library_product(m2do_fea, "m2do_fea.dll")
JLLWrappers.@declare_library_product(m2do_lsm, "m2do_lsm.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        m2do_3d_lsm,
        "bin\\m2do_3d_lsm.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        m2do_fea,
        "bin\\m2do_fea.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        m2do_lsm,
        "bin\\m2do_lsm.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module BoundaryGridsTests

using Test
using Gridap

model = CartesianDiscreteModel(partition=(3,4,2))

tags = [23,24,25]
grid = BoundaryGrid(model,tags)
test_grid(grid,60,20)

tags = "boundary"
grid = BoundaryGrid(model,tags)

tags = ["physical_tag_23","physical_tag_24","physical_tag_25"]
grid = BoundaryGrid(model,tags)
test_grid(grid,60,20)

#writevtk(model,"model")

#writevtk(grid,"grid",celldata=["cell"=>grid.descriptor.facet_to_cell,"lfacet"=>grid.descriptor.facet_to_lfacet])

trian = Triangulation(grid)
@test isa(trian,BoundaryTriangulation)

end # module


text:
export read, write, existsfile, mkdir
export filenames, filepaths, dirnames, dirpaths
export readmat, writemat
# using MAT

existsfile(filename::AbstractString) = (s = stat(filename); s.inode!=0)

#import Base.mkdir
#function amkdir(a)
#    if a[1]=='/'
#        start = '/';
#    else
#        start = "";
#    end
#
#    if !existsfile(a)
#        parts = split(a,'/');
#        n = len(parts);
#        for i = 1:n
#      @show riffle(parts[1:i],'/')
#        d = [start join(riffle(parts[1:i],'/'))];
#            if !(existsfile(d))
#                Base.mkdir(d)
#            end
#        end
#    end
#end

import Base.read
read(filename::AbstractString) = open(read, filename)

import Base.write
write(data::AbstractString, filename::AbstractString) = write(data, filename, "w")
function write(data::AbstractString, filename::AbstractString, mode)
    io = open(filename, mode)
    finalizer(io,close)
    write(io,data)
    close(io)
end

function filedirnames(path = pwd(); selector = isdir, hidden = false, withpath = false,
    recursive = false, suffix = "", regex = Void)
    # @show path selector hidden withpath
    files = readdir(path)
    r = filter(x->selector(joinpath(path,x)) && (hidden || x[1]!='.'), files)
    # @show r
    r = sort(r)
    r = withpath ? map(r, x->joinpath(path,x)) : r
    f = (VERSION < v"0.5-") ? utf8 : AbstractString
    r = map(r, f)
    if recursive
        f(x) = filedirnames(x;
            selector = selector, hidden = hidden, withpath = withpath, recursive = true)
        r = @p dirpaths path | map f | flatten | concat r _
    end
    if !isempty(suffix)
        suffix = isa(suffix, AbstractString) ? [suffix] : suffix
        r = @p mapvec suffix (suf->@p filter r x->@p takelast x len(suf) | isequal suf) | flatten
    end
    if regex != Void
        r = @p filter r x->@p ismatch regex x
    end
    r
end
dirnames(path = pwd(); kargs...) = filedirnames(path; selector = isdir, kargs...)
dirpaths(path = pwd(); kargs...) = dirnames(path; withpath = true, kargs...)
filenames(path = pwd(); kargs...) = filedirnames(path; selector = x->!isdir(x), kargs...)
filepaths(path = pwd(); kargs...) = filenames(path; withpath = true, kargs...)

function readmat(filename)
    return matread(filename)
end

readmat(filename,variables) = readmat(filename, variables...)

function readmat(filename,variables...)
    r = Any[]
    mat = matopen(filename)
    try
        for v in variables
            r[v] = read(mat, v)
        end
    finally
        close(mat)
    end
end


text:
using IL_simulations, PDBTools

top_dir     =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/VSIL/ITP"
pdb_dir     =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/VSIL/PDB"
input_dir   =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/ubiquitin_files"
top_ion_dir =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/ion_files"

atoms   = readPDB("$input_dir/ubq.pdb")
MMP     = mass(atoms)

cation1 = "BMIM"
anion1  = "NO3"
anion2  = "DCA"

function molarmass(anion, cation, pdb_dir)
  atcation = readPDB("$pdb_dir/$(cation)_VSIL.pdb")
  atanion  = readPDB("$pdb_dir/$(anion)_VSIL.pdb")
  MMIL = mass(atcation) + mass(atanion)
  return MMIL
end

MMIL  = molarmass(anion1, cation1, pdb_dir)
MMIL2 = molarmass(anion2, cation1, pdb_dir)
c       =  0.10

## test with one IL/salt and a protein/polymer
#cd("PONEIL")
#newil = EMIMCL

cation = "EMIM"
anion  = "Br"
MMIL = molarmass(anion, cation, pdb_dir)
dir = pwd();
for i in ["0.50", "1.00", "1.50", "2.00", "2.50", "3.00"]
  for j in ["00", "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19","20"]
    cd("$dir/EMIMBr/$i/$j")
    conc = parse(Float64,i)
    data = IL_simulations.Data_il(protein="$input_dir/ubq.pdb", MMP = MMP, cation = cation, anion = anion,MM = MMIL,c=conc)
    nions, nwater, sides = prot_elec(data)
    IL_simulations.top(dict, top_dir, input_dir, data, nions, nwater)
    IL_simulations.pack_input(data,pdb_dir,nions,nwater,sides)
    IL_simulations.mdp_files_solute()
  end
end

#println("PONEIL")
## test with one IL/salt
#cd("../ONEIL")
#data = IL_simulations.Data_elec(cation = cation1, anion = anion1, MM = MMIL,c=c)
#nions, nwater, sides = sol_elec(data)
#IL_simulations.topelec(dict, top_dir, top_ion_dir, data, nions, nwater)
#IL_simulations.pack_input_sol(data,pdb_dir,nions,nwater,sides)
#IL_simulations.mdp_files()

#println("ONEIL")

# test with two ILs/salts and a protein/polymer -same cation
#cd("../PTWOIL")
#data = IL_simulations.Data_2il(protein="$input_dir/ubq.pdb", MMP = MMP, cation = cation1, anion1 = anion1, anion2 = anion2, MM1 = MMIL, MM2 = MMIL2,c=c)
#ncat, nan, nwater, sides = prot_melec(data)
#IL_simulations.top(dict, top_dir, input_dir, data, ncat, nan, nwater)
#IL_simulations.pack_input(data, pdb_dir, ncat, nan, nwater, sides)
#IL_simulations.mdp_files()

#println("PTWOIL")

### test with two ILs/salts - same cation
#cd("../TWOIL")
#data = IL_simulations.Data_twoelec(cation = cation1, anion1 = anion1, anion2 = anion2 , MM1 = MMIL, MM2 = MMIL2, c=c)
#ncat, nan, nwater, sides = sol_melec(data)
#IL_simulations.topelec(dict,top_dir, top_ion_dir, data, ncat, nan, nwater)
#IL_simulations.pack_input_elec(data, pdb_dir, ncat, nan, nwater, sides)
#IL_simulations.mdp_files()
#
#
#println("TWOIL")


text:
mutable struct Swarm{M,T}

    Particles::Vector{Particle} # Vector of Particles

    b::T                        # Global Best Objective Function Value
    d::SVector{M,T}             # Location of Global Best Objective Function

    N::Int                      # Neighborhood Size
    W::AbstractFloat            # Inertia
    c::Int                      # Adaptive Inertia Counter

    y1::AbstractFloat           # SelfAdjustmentWeight
    y2::AbstractFloat           # SocialAdjustmentWeight

end

##########################
#   Swarm Constructors   #
##########################
function Swarm(NDims,NParticals)
    temp = Vector{Particle}(undef, NParticals)
    for i in 1:NParticals
        temp[i] = Particle(NDims)
    end
    Swarm(temp, Inf, temp[1].p, 2, Inf, 0, 0.0, 0.0)
end

#####################
#   Swarm Methods   #
#####################

Base.length(S::Swarm) = length(S.Particles)

function Base.getindex(S::Swarm, i::Int)
    1 <= i <= length(S.Particles) || throw(BoundsError(S, i))
    return S.Particles[i]
end

function Base.setindex!(S::Swarm, v, i::Int)
    1 <= i <= length(S.Particles) || throw(BoundsError(S, i))
    S.Particles[i] = v
end

Base.lastindex(S::Swarm) = S[length(S)]

function fvalcheck(S::Swarm)
    n = length(S.Particles)
    @inbounds for i in 1:n
        if ifinf(S[i].fx) || isnan(S[i].fx)
            throw(ErrorException("Objective function returned Inf or NaN!"))
        end
    end
end

function setgbest!(S::Swarm)
    n = length(S)
    updated = false
    @inbounds for i in 1:n
        if S[i].fp < S.b
            S.b = S[i].fp
            S.d = deepcopy(S[i].p)
            updated = true
        end
    end
    return updated
end

function printstatus(S::Swarm, Time::AbstractFloat, Iteration::Int, StallCount::Int)
    # Time Elapsed, Iteration Number, Function Evaluations, Stall Iterations, Global Best
    fspec1 = FormatExpr("Time Elapsed: {1:f} sec, Iteration Number: {2:d}, Function Evaluations: {3:d}")
    fspec2 = FormatExpr("Stall Iterations: {1:d}, Global Best: {2:e}")
    printfmtln(fspec1, Time, Iteration, (Iteration + 1)*length(S))
    printfmtln(fspec2, StallCount, S.b)
    println(" ")
end

function updatevelocities!(S::Swarm)
    n = length(S.Particles)
    m = length(S.d)
    @inbounds for i in 1:n
        # Chose random subset of N Particles not including i
        ws = Weights([float(j != i) for j in 1:n])
        samp = sample(1:n, ws, S.N; replace=false, ordered=false)

        # Determine fbest(S)
        fbest = Inf
        best  = 0
        @inbounds for j in samp
            if S[j].fp < fbest
                fbest = S[j].fp
                best = j
            end
        end

        # Update i's velocity
        u1 = SVector{m}(rand(m))
        u2 = SVector{m}(rand(m))
        S[i].v = @. S.W*S[i].v + S.y1*u1*(S[i].p - S[i].x) + S.y2*u2*(S[best].p - S[i].x)
    end
end

function step!(S::Swarm)
    n = length(S.Particles)
    @inbounds for i in 1:n
        S[i].x = S[i].x .+ S[i].v
    end
end

function enforcebnds!(S::Swarm, Opts::PSOOptions)
    n = length(S)
    tempx = MVector{Opts.NDims}(zeros(Opts.NDims))
    tempv = MVector{Opts.NDims}(zeros(Opts.NDims))
    @inbounds for i in 1:n
        tempx .* 0.0
        tempv .* 0.0
        @inbounds for j in 1:Opts.NDims
            if S[i].x[j] > Opts.UpperBounds[j]
                tempx[j] = Opts.UpperBounds[j]
                tempv[j] = 0.0
            elseif S[i].x[j] < Opts.LowerBounds[j]
                tempx[j] = Opts.LowerBounds[j]
                tempv[j] = 0.0
            else
                tempx[j] = S[i].x[j]
                tempv[j] = S[i].v[j]
            end
        end
        S[i].x = SVector{Opts.NDims}(tempx)
        S[i].v = SVector{Opts.NDims}(tempv)
    end
end

function feval!(f, S::Swarm, Opts::PSOOptions)

    # Evaluate Objective Function for all Particles
    if Opts.UseParallel
        Threads.@threads for i in 1:Opts.SwarmSize
            S[i].fx = f(S[i].x)
        end
    else
        @inbounds for i in 1:Opts.SwarmSize
            S[i].fx = f(S[i].x)
        end
    end

    # Check Objective Function Values if Desired
    if Opts.FunValCheck == "on"
        fvalcheck(S)
    end

    # Update Each Particles Best Objective Function Value and its Location
    @inbounds for i in 1:Opts.SwarmSize
        if S[i].fx < S[i].fp
            S[i].p = deepcopy(S[i].x)
            S[i].fp = S[i].fx
        end
    end
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
### General utilities

"""
    top, success = chaintop(box)

Return a Box `top` that lies at the top of a chain containing
`box`. If we let `n = ndims(box)` and `p = degree(box)`, a chain has
the property that if one follows `ceil(Int, n/p)` splits downward from
`top`, at each stage choosing the last child, the splits will cover
each dimension once. (Depending on whether `p` divides `n` evenly,
some of the dimensions may be fictive.)  All nodes generated by splits
between the top and bottom of the chain are considered to be part of
the chain. For every node in the chain, its position satifies `x[i] ∈
{b[i], y[i]}` for a base point `b = position(top)` and "opposite
point" `y` which is the position of the final node in the chain.

`success` is true if `box` is a member of a chain. For a general CSp
tree, success is not guaranteed. However, if you only grow the tree
via [`addpoint!`](@ref), all nodes are members of one or more chains
and success is guaranteed.
"""
function chaintop(box::Box{p}) where p
    function process_split!(covered, nremaining, split)
        n = length(covered)
        nnovel = ntarget = p  # how many novel dimensions we're looking for
        for d in split.dims
            nnovel -= d > n || covered[d]
            ntarget -= d > n  # don't count fictive dimensions against the total
        end
        nnovel == ntarget || nnovel == nremaining || return false, nremaining
        # All of the dimensions in this split are novel, so we accept
        # the split and signal that we should keep going
        for d in split.dims
            d <= n || continue
            covered[d] = true
        end
        nremaining -= nnovel
        return true, nremaining
    end

    n = ndims(box)
    covered = falses(n)   # covered[d] signals whether dimension d has appeared
    nremaining = n        # number of dimensions that we haven't yet seen
    bottom = top = box
    # First, go down
    split = box.parent.split
    while nremaining > 0
        didsplit, nremaining = process_split!(covered, nremaining, split)
        didsplit || break
        bottom = split.others.children[end]
        isleaf(bottom) && break
        # If the tree was grown via `addpoint!`, ends of chains may be
        # marked by fictive dimensions. So look for those.
        hasfictive = false
        for d in split.dims
            if d > n
                hasfictive = true
            end
        end
        hasfictive && break
        split = bottom.split
    end
    top = box.parent
    nremaining == 0 && return top, true  # box was already at the top
    isroot(box) && return box, false  # didn't find enough, but we can't go up
    # Go up until all are covered or we encounter a conflict
    split = top.parent.split
    while nremaining > 0
        didsplit, nremaining = process_split!(covered, nremaining, split)
        didsplit || nremaining > 0 || break
        isroot(top) && return top, false
        top = top.parent
        split = top.parent.split
    end
    while nremaining > 0
        # We encountered a conflict before finishing all the
        # dimensions.  For a tree grown with `addpoint!`, this
        # indicates nodes currently between bottom and top are split
        # between two chains. (Since we look for fictive dimensions
        # and those are always at the end of a chain, this can only
        # happen in cases where p divides n evenly.)

        # To solve the problem, backtrack one level at a time until we
        # can complete the chain at the top.
        split = bottom.parent.split
        for d in split.dims
            covered[d] = false
        end
        nremaining += p
        split = top.parent.split
        didsplit, nremaining = process_split!(covered, nremaining, split)
        didsplit || nremaining > 0 || break
        isroot(top) && return top, false
        bottom = bottom.parent
        top = top.parent
    end
    return top, nremaining == 0
end

"""
    Cp = coefficients_p([f=value], box; skip=0)

When fitting a polynomial of degree `p = degree(box)` to the data
around `box`, compute those highest-order coefficients that can be
inferred directly via the `O(1)`-per-split algorithm. This computes
elements of `Cp` that correspond to splits present in `box`'s tree.
Any coefficients that can't be determined via the `O(1)` algorithm are
filled with NaN.

For a CS2 tree, `Cp` will store the off-diagonal elements of the
Hessian. The diagonals (and any undetermined off-diagonals) will be
NaN.

The optional argument `skip` allows you to omit early splits. This can
be useful if you know that some calculation will be adversely affected
by early splits. An example would be in estimating a residual (model
error), where the model's degrees of freedom allow it to fit early
splits exactly but for which there is no reason to believe that this
is anything other than overfitting.

Finally, there are many other variants:

    Cp = coefficients_p(f, box, callback)   # callback(box, val) is called for each box visited
    Cp, state = coefficients_p(f, box, iter, state, callback)   # for custom iteration over splits
    coefficients_p!(f, Cp, box)    # for pre-allocated outputs, can exploit sparsity
"""
coefficients_p(box::Box{p,T}) where {p,T} = coefficients_p(value, box) # supply the "readout" function

coefficients_p(f::Function, box::Box; skip::Int=0) =
    coefficients_p(f, box, (box,val)->nothing; skip=skip) # add the callback

function coefficients_p(f::Function, box::Box{p,T}, callback::Function;
                        skip::Int=0) where {p,T}
    # Create the iterator
    iter = splits(box)
    state = skipsplits(iter, skip)
    Cp, state = coefficients_p(f, box, iter, state, callback)
    return Cp
end

function coefficients_p(f::Function,
                        box::Box{p,T},
                        iter,
                        state,
                        callback::Function) where {p,T}
    # Allocate the output
    n = ndims(box)
    Cp = allocate_coefficients_p(n, box)
    return coefficients_p!(f, Cp, box, iter, state, callback)
end

coefficients_p!(Cp::AbstractArray, box::Box{p,T}; skip::Int=0) where {p,T} =
    coefficients_p!(value, Cp, box; skip=skip)

coefficients_p!(f::Function, Cp::AbstractArray, box::Box; skip::Int=0) =
    coefficients_p!(f, Cp, box, (box,val)->nothing; skip=skip)

coefficients_p!(Cp::AbstractArray, box::Box{p,T}, callback::Function;
    skip::Int=0) where {p,T} = coefficients_p!(value, Cp, box, callback; skip=skip)

function coefficients_p!(f::Function, Cp::AbstractArray, box::Box{p,T}, callback::Function;
                         skip::Int=0) where {p,T}
    iter = splits(box)
    state = skipsplits(iter, skip)
    Cp, state = coefficients_p!(f, Cp, box, iter, state, callback)
    return Cp
end

function coefficients_p!(Cp::AbstractArray{S,p},
                         box::Box{p,T},
                         iter,
                         state,
                         callback::Function) where {p,S,T}
    return coefficients_p!(value, Cp, box, iter, state, callback)
end

# This implementation makes it possible for callers to continue
# iterating where this left off. A practical application is a CS2 tree
# being used to construct a quadratic model in 2 dimensions: 1 split
# is enough to determine the off-diagonal coefficient, but a 2nd split
# is needed to determine the gradient and diagonal entries.
function coefficients_p!(f::Function,
                         Cp::Union{SymmetricArray{S,p},SymmetricMatrix{S}},
                         box::Box{p,T},
                         iter,
                         state,
                         callback::Function) where {p,S,T}
    function bitsum(i)
        s = Int(i & 0x01)
        while i != 0
            i = i >> 1
            s += i & 0x01
        end
        return s
    end
    bitsign(i) = isodd(bitsum(i)) ? -1 : 1

    # Calculate the number of coefficients that could be set (n choose p)
    n = ndims(box)
    ndims(Cp) == p || throw(DimensionMismatch("output should be $p-dimensional, got $(ndims(Cp))"))
    all(isequal(n), size(Cp)) || throw(DimensionMismatch("output should have size $n along each axis, got $(size(Cp))"))
    nremaining = nnz_offdiag_sym(Cp)
    fillnz!(Cp, NaN)
    # Iterate until we fill all coefficients or exhaust the tree
    result = iterate(iter, state)
    while result !== nothing
        split, state = result
        result = iterate(iter, state)
        dims = split.dims
        maximum(dims) > n && continue
        if isnan(Cp[dims...])
            denom = oneunit(T)
            for (d, x) in zip(split.dims, split.xs)
                denom *= x - position(split.self, d)
            end
            s = isodd(p) ? -1 : 1
            num = zero(T)
            for ichild = 0:maxchildren(box)-1
                child = getchild(split, ichild)
                fc = f(child)
                num += (s * bitsign(ichild)) * fc
                callback(child, fc)
            end
            Cp[dims...] = num/denom
            nremaining -= 1
        end
        nremaining <= 0 && break
    end
    return Cp, state
end

allocate_coefficients_p(n, box::Box{p,T}) where {p,T} =
    SymmetricArray(fill(T(NaN), ntuple(d->n, Val(p))))

function nnz_offdiag_sym(A)
    n = size(A, 1)
    nremaining = n
    for i = 1:ndims(A)-1
        nremaining *= (n-i)/(i+1)
    end
    return round(Int, nremaining)
end

function nnz_offdiag_sym(A::SymmetricArray{T,2,<:SparseMatrixCSC}) where T
    nremaining = 0
    S = A.data
    n = size(S, 1)
    for j = 1:n
        for idx = S.colptr[j]:S.colptr[j+1]-1
            i = S.rowval[idx]
            nremaining += i < j
        end
    end
    return nremaining
end

function nnz_offdiag_sym(A::SymTridiagonal)
    return length(A.ev)
end

fillnz!(A::AbstractArray, val) = fill!(A, val)
fillnz!(A::SymmetricArray, val) = fillnz!(A.data, val)
function fillnz!(A::SparseMatrixCSC, val)
    fill!(A.nzval, val)
    return A
end
function fillnz!(A::SymTridiagonal, val)
    fill!(A.dv, val)
    fill!(A.ev, val)
    return A
end

# Move to stdlib/Base:
function Base.setindex!(A::SymTridiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    if i == j
        @inbounds A.dv[i] = x
    elseif i == j+1
        @inbounds A.ev[j] = x
    elseif i+1 == j
        @inbounds A.ev[i] = x
    else
        throw(ArgumentError("cannot set off-diagonal entry ($i, $j)"))
    end
    return x
end


text:
struct WeightPointValue{W,P,V}
  weight::W
  point::P
  value::V
end

function quadpoints(chart::ReferenceSimplex{1}, rule)
    u, w = legendre(rule, 0.0, 1.0)
    [(neighborhood(chart, u[:,i]), w[i]) for i in eachindex(w)]
end

function quadpoints(chart::ReferenceSimplex{2}, rule)
    u, w = trgauss(rule)
    [(neighborhood(chart, u[:,i]), w[i]) for i in eachindex(w)]
end

"""
    pw = quadpoints(chart, rule)

Returns a collection of (point, weight) tuples corresponding to the numerical quadrature `rule`
defined on the domain of `chart`. The weights returned take into account the Jacobian determinant
resulting from mapping from the reference domain to the configuration space.

Functions can be integrated like:

```julia
PW = quadpoints(chart, rule)
I = sum(pw[2]*f(pw[1]) for pw in PW)
```
"""
function quadpoints(chart::Union{Simplex,subd_chart}, rule)
    PV = quadpoints(domain(chart), rule)
    map(PV) do pv
        q = neighborhood(chart, pv[1])
        w = jacobian(q)*pv[2]
        (q,w)
    end
end

# function quadpoints(chart::subd_chart, rule)
#     u, w = trgauss(rule)
#     PV = [(neighborhood(chart, u[:,i]), w[i]) for i in eachindex(w)]
#     map(PV) do pv
#         q = pv[1]
#         w = jacobian(q)*pv[2]
#         (q,w)
#     end
# end

"""
    quadpoints(refspace, charts, rules)

Computed a matrix of vectors containing (weight, point, value) triples that can
be used in numerical integration over the elements described by the charts. Internally,
this method used `quadpoints(chart, rule)` to retrieve the points and weights for
a certain quadrature rule over `chart`.
"""
function quadpoints(f, charts, rules)

    pw = quadpoints(charts[1], rules[1])
    P = typeof(pw[1][1])
    W = typeof(pw[1][2])
    V = typeof(f(pw[1][1]))

    WPV = WeightPointValue{W,P,V}

    qd = Array{Vector{WPV}}(undef,length(rules), length(charts))
    for j in eachindex(charts)
        for i in eachindex(rules)
            pw = quadpoints(charts[j], rules[i])
            qd[i,j] = Vector{WPV}(undef,length(pw))
            for k in eachindex(pw)
                qd[i,j][k] = WeightPointValue(pw[k][2],pw[k][1],f(pw[k][1]))
            end
        end
    end

  qd
end


text:
#Julia Vega Plotting

#Pkg.add("Vega")
#Pkg.add("Compat")
#Pkg.add("Patchwork")
using Vega
using Compat
using Patchwork
Patchwork.load_js_runtime()
stock = ["chairs", "tables", "desks", "rugs", "lamps"];
quantity = [15, 10, 10, 5, 20];
piechart(x = stock, y = quantity)

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
"""
    debug_mode() -> Bool

Determines if ChainRulesCore is in `debug_mode`.
Defaults to `false`, but if the user redefines it to return `true` then extra
information will be shown when errors occur.

Enable via:
```
ChainRulesCore.debug_mode() = true
```
"""
debug_mode() = false


text:
import Base: *, A_mul_B!, Ac_mul_B!, At_mul_B!
import Base: getindex, setindex!

type SparseMatrixCSR{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Int                  # Number of rows
    n::Int                  # Number of columns
    rowptr::Vector{Ti}      # Row i is in rowptr[i]:(rowptr[i+1]-1)
    colval::Vector{Ti}      # Column values of nonzeros
    nzval::Vector{Tv}       # Nonzero values
end

# constructor in the form of `sparse(I, J, V, m, n, combine) -> CSC` in Base
function SparseMatrixCSR{Tv,Ti<:Integer}(I::AbstractVector{Ti},
                                         J::AbstractVector{Ti},
                                         V::AbstractVector{Tv},
                                         nrow::Integer, ncol::Integer,
                                         combine::Union{Function,Base.Func}=Base.AddFun())
    # super simple solution for now. Just construct the CSC with the rows/cols
    # and I/J flipped then extract the fields
    csc = sparse(J, I, V, ncol, nrow, combine)
    rowptr = csc.colptr
    colval = csc.rowval
    nzval = csc.nzval
    SparseMatrixCSR(nrow, ncol, rowptr, colval, nzval)
end

copy(S::SparseMatrixCSR) =
    SparseMatrixCSR(S.m, S.n, copy(S.rowptr), copy(S.colval), copy(S.nzval))

spzeros_csr(m::Integer, n::Integer) = spzeros_csr(Float64, m, n)
spzeros_csr(Tv::Type, m::Integer, n::Integer) = spzeros_csr(Tv, Int, m, n)
function spzeros_csr(Tv::Type, Ti::Type, m::Integer, n::Integer)
    ((m < 0) || (n < 0)) && throw(ArgumentError("Invalid Array dimensions"))
    SparseMatrixCSR(m, n, ones(Ti, m+1), Array(Ti, 0), Array(Tv, 0))
end

# NOTE:  copied/pasted from base. Modified slightly
function Base.showarray(io::IO, S::SparseMatrixCSR;
                        header::Bool=true, limit::Bool=Base._limit_output,
                        rows = Base.tty_size()[1], repr=false)

    if header
        print(io, S.m, "x", S.n, " CSR sparse matrix with ", nnz(S), " ",
              eltype(S), " entries:")
    end

    if limit
        half_screen_rows = div(rows - 8, 2)
    else
        half_screen_rows = typemax(Int)
    end
    pad = ndigits(max(S.m, S.n))
    k = 0
    sep = "\n\t"
    for row = 1:S.m, k = S.rowptr[row] : (S.rowptr[row+1]-1)
        if k < half_screen_rows || k > nnz(S)-half_screen_rows
            print(io, sep, '[', rpad(row, pad), ", ", lpad(S.colval[k], pad), "]  =  ")
            showcompact(io, S.nzval[k])
        elseif k == half_screen_rows
            print(io, sep, '\u22ee')
        end
        k += 1
    end
end

# algorithm based on scipy/sparse/sparsetools/csr.h
function csr2csc{Tv,Ti}(indptr::Vector{Ti}, indval::Vector{Ti},
                        nzval::Vector{Tv}, m::Int, n::Int)
    # allocate
    Bptr = zeros(Ti, n+1)
    Bind = similar(indval)
    Bval = similar(nzval)

    numnz = indptr[m+1] - 1

    # get colptr by accumulating hits on colval and doing cumsum
    @inbounds Bptr[1] = 1
    @inbounds for i=1:numnz Bptr[indval[i]+1] += 1; end
    Bptr = cumsum(Bptr)

    @inbounds for row=1:m
        for j=indptr[row]:(indptr[row+1]-1)
            col = indval[j]
            dest = Bptr[col]

            Bind[dest] = row
            Bval[dest] = nzval[j]
            Bptr[col] += 1
        end
    end

    # fix up Bptr
    last = 1
    @inbounds for col=1:n
        temp = Bptr[col]
        Bptr[col] = last
        last = temp
    end
    Bptr, Bind, Bval
end

function Base.convert{Tv,Ti<:Integer}(::Type{SparseMatrixCSC},
                                      A::SparseMatrixCSR{Tv,Ti})
    Bptr, Bind, Bval = csr2csc(A.rowptr, A.colval, A.nzval, A.m, A.n)
    SparseMatrixCSC(A.m, A.n, Bptr, Bind, Bval)
end

function Base.convert{Tv,Ti<:Integer}(::Type{SparseMatrixCSR},
                                      B::SparseMatrixCSC{Tv,Ti})
    Aptr, Aind, Aval = csr2csc(B.colptr, B.rowval, B.nzval, B.n, B.m)
    SparseMatrixCSR(B.m, B.n, Aptr, Aind, Aval)
end

## getindex methods. NOTE: copied and modified from Base
getindex(A::SparseMatrixCSR, I::Tuple{Integer,Integer}) =
    getindex(A, I[1], I[2])

function getindex{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, i0::Integer, i1::Integer)
    !(1 <= i0 <= A.m && 1 <= i1 <= A.n) && throw(BoundsError())
    r1 = Int(A.rowptr[i0])
    r2 = Int(A.rowptr[i0+1]-1)
    (r1 > r2) && return zero(Tv)
    r1 = searchsortedfirst(A.colval, i1, r1, r2, Base.Order.Forward)
    ((r1 > r2) || (A.colval[r1] != i1)) ? zero(Tv) : A.nzval[r1]
end

getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{T}, j::Integer) =
    getindex(A, I, [j])
getindex{T<:Integer}(A::SparseMatrixCSR, i::Integer, J::AbstractVector{T}) =
    getindex(A, [i], J)

# Colon translation (this could be done more efficiently)
getindex(A::SparseMatrixCSR, ::Colon) = getindex(A, 1:length(A))
getindex(A::SparseMatrixCSR, ::Colon, ::Colon) = getindex(A, 1:size(A, 1), 1:size(A, 2))
getindex(A::SparseMatrixCSR, ::Colon, j) = getindex(A, 1:size(A, 1), j)
getindex(A::SparseMatrixCSR, i, ::Colon) = getindex(A, i, 1:size(A, 2))

function getindex_rows{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector)
    m, n = size(A)
    nI = length(I)

    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval

    rowptrB = Array(Ti, nI+1)
    rowptrB[1] = 1
    nnzB = 0

    @inbounds for i = 1:nI
        row = I[i]
        1 <= row <= m || throw(BoundsError())
        nnzB += rowptrA[row+1] - rowptrA[row]
        rowptrB[i+1] = nnzB + 1
    end

    colvalB = Array(Ti, nnzB)
    nzvalB  = Array(Tv, nnzB)
    ptrB = 0

    @inbounds for i = 1:nI
        row = I[i]
        for k = rowptrA[row]:rowptrA[row+1]-1
            ptrB += 1
            colvalB[ptrB] = colvalA[k]
            nzvalB[ptrB] = nzvalA[k]
        end
    end
    return SparseMatrixCSR(nI, n, rowptrB, colvalB, nzvalB)
end

function getindex{Tv,Ti<:Integer}(A::SparseMatrixCSR{Tv,Ti},
                                  I::AbstractVector, J::Range)
    # Ranges for indexing rows
    (m, n) = size(A)

    # whole rows:
    if J == 1:n
        return getindex_rows(A, I)
    end

    nI = length(I)
    nJ = length(J)
    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval
    rowptrS = Array(Ti, nI+1)
    rowptrS[1] = 1
    nnzS = 0

    # Form the structure of the result and compute space
    @inbounds for i = 1:nI
        row = I[i]
        @simd for k in rowptrA[row]:rowptrA[row+1]-1
            nnzS += colvalA[k] in J
        end
        rowptrS[i+1] = nnzS+1
    end

    # Populate the values in the result
    colvalS = Array(Ti, nnzS)
    nzvalS  = Array(Tv, nnzS)
    ptrS    = 1

    @inbounds for i = 1:nI
        row = I[i]
        for k = rowptrA[row]:rowptrA[row+1]-1
            colA = colvalA[k]
            j = Base.SparseMatrix.rangesearch(J, colA)
            if j > 0
                colvalS[ptrS] = j
                nzvalS[ptrS] = nzvalA[k]
                ptrS += 1
            end
        end
    end

    return SparseMatrixCSR(nI, nJ, rowptrS, colvalS, nzvalS)
end

function getindex_J_sorted{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector,
                                  J::AbstractVector)
    # sorted vectors for indexing columns
    # similar to getindex_general, but without the transpose trick
    m, n = size(A)

    nJ = length(J)
    avgN = div(nnz(A), m)

    # heuristics based on experiments
    alg = ((nJ - avgN) > 2^8)  ? 1 :
          ((avgN - nJ) > 2^10) ? 0 : 2

    return (alg == 0) ? getindex_J_sorted_bsearch_A(A, I, J) :
           (alg == 1) ? getindex_J_sorted_bsearch_J(A, I, J) :
                        getindex_J_sorted_linear(A, I, J)
end

function getindex_J_sorted_bsearch_A{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti},
                                            I::AbstractVector, J::AbstractVector)
    const nI = length(I)
    const nJ = length(J)

    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval
    rowptrS = Array(Ti, nI+1)
    rowptrS[1] = 1

    ptrS = 1
    # determine result size
    @inbounds for i = 1:nI
        row = I[i]
        ptrJ::Int = 1 # runs through J
        ptrA::Int = rowptrA[row]
        stopA::Int = rowptrA[row+1]-1
        if ptrA <= stopA
            while ptrJ <= nJ
                colJ = J[ptrJ]
                ptrJ += 1
                (colvalA[ptrA] > colJ) && continue
                ptrA = searchsortedfirst(colvalA, colJ, ptrA, stopA, Base.Order.Forward)
                (ptrA <= stopA) || break
                if colvalA[ptrA] == colJ
                    ptrS += 1
                end
            end
        end
        rowptrS[i+1] = ptrS
    end

    colvalS = Array(Ti, ptrS-1)
    nzvalS  = Array(Tv, ptrS-1)

    # fill the values
    ptrS = 1
    @inbounds for i = 1:nI
        row = I[i]
        ptrJ::Int = 1 # runs through J
        ptrA::Int = rowptrA[row]
        stopA::Int = rowptrA[row+1]-1
        if ptrA <= stopA
            while ptrJ <= nJ
                colJ = J[ptrJ]
                if colvalA[ptrA] <= colJ
                    ptrA = searchsortedfirst(colvalA, colJ, ptrA, stopA, Base.Order.Forward)
                    (ptrA <= stopA) || break
                    if colvalA[ptrA] == colJ
                        colvalS[ptrS] = ptrJ
                        nzvalS[ptrS] = nzvalA[ptrA]
                        ptrS += 1
                    end
                end
                ptrJ += 1
            end
        end
    end
    return SparseMatrixCSR(nI, nJ, rowptrS, colvalS, nzvalS)
end

function getindex_J_sorted_linear{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti},
                                         I::AbstractVector, J::AbstractVector)
    const nI = length(I)
    const nJ = length(J)

    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval
    rowptrS = Array(Ti, nI+1)
    rowptrS[1] = 1
    cacheJ = zeros(Int, A.n)

    ptrS   = 1
    # build the cache and determine result size
    @inbounds for i = 1:nI
        row = I[i]
        ptrJ::Int = 1 # runs through I
        ptrA::Int = rowptrA[row]
        stopA::Int = rowptrA[row+1]
        while ptrJ <= nJ && ptrA < stopA
            colA = colvalA[ptrA]
            colJ = J[ptrJ]

            if colJ > colA
                ptrA += 1
            elseif colJ < colA
                ptrJ += 1
            else
                (cacheJ[colA] == 0) && (cacheJ[colA] = ptrJ)
                ptrS += 1
                ptrJ += 1
            end
        end
        rowptrS[i+1] = ptrS
    end

    colvalS = Array(Ti, ptrS-1)
    nzvalS  = Array(Tv, ptrS-1)

    # fill the values
    ptrS = 1
    @inbounds for i = 1:nI
        row = I[i]
        ptrA::Int = rowptrA[row]
        stopA::Int = rowptrA[row+1]
        while ptrA < stopA
            colA = colvalA[ptrA]
            ptrJ = cacheJ[colA]
            if ptrJ > 0
                while ptrJ <= nJ && J[ptrJ] == colA
                    colvalS[ptrS] = ptrJ
                    nzvalS[ptrS] = nzvalA[ptrA]
                    ptrS += 1
                    ptrJ += 1
                end
            end
            ptrA += 1
        end
    end
    return SparseMatrixCSR(nI, nJ, rowptrS, colvalS, nzvalS)
end

function getindex_J_sorted_bsearch_J{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti},
                                            I::AbstractVector, J::AbstractVector)
    const nI = length(I)
    const nJ = length(J)

    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval
    rowptrS = Array(Ti, nI+1)
    rowptrS[1] = 1

    n = A.n

    # cacheJ is used first to store num occurrences of each col in rows of
    # interest and later to store position of first occurrence of each col in J
    cacheJ = zeros(Int, n)

    # count rows
    @inbounds for i = 1:nI
        row = I[i]
        for ptrA in rowptrA[row]:(rowptrA[row+1]-1)
            cacheJ[colvalA[ptrA]] += 1
        end
    end

    # fill cache and count nnz
    ptrS::Int = 0
    ptrJ::Int = 1
    @inbounds for i = 1:n
        rval = cacheJ[i]
        (rval == 0) && continue
        ptrJ = searchsortedfirst(J, i, ptrJ, nJ, Base.Order.Forward)
        cacheJ[i] = ptrJ
        while ptrJ <= nJ && J[ptrJ] == i
            ptrS += rval
            ptrJ += 1
        end
        if ptrJ > nJ
            @simd for j=(i+1):n; @inbounds cacheJ[i]=ptrJ; end
            break
        end
    end

    colvalS = Array(Ti, ptrS)
    nzvalS  = Array(Tv, ptrS)
    rowptrS[nI+1] = ptrS+1

    # fill the values
    ptrS = 1
    @inbounds for i = 1:nI
        row = I[i]
        ptrA::Int = rowptrA[row]
        stopA::Int = rowptrA[row+1]
        while ptrA < stopA
            colA = colvalA[ptrA]
            ptrJ = cacheJ[colA]
            (ptrJ > nJ) && break
            if ptrJ > 0
                while J[ptrJ] == colA
                    colvalS[ptrS] = ptrJ
                    nzvalS[ptrS] = nzvalA[ptrA]
                    ptrS += 1
                    ptrJ += 1
                    (ptrJ > nJ) && break
                end
            end
            ptrA += 1
        end
        rowptrS[i+1] = ptrS
    end
    return SparseMatrixCSR(nI, nI, rowptrS, colvalS, nzvalS)
end

function permute_cols!{Tv,Ti}(S::SparseMatrixCSR{Tv,Ti}, pJ::Vector{Int})
    (m, n) = size(S)
    rowptrS = S.rowptr; colvalS = S.colval; nzvalS = S.nzval
    # preallocate temporary sort space
    nc = min(nnz(S), n)
    colperm = Array(Int, nc)
    colvalTemp = Array(Ti, nc)
    nzvalTemp = Array(Tv, nc)

    for i in 1:m
        colrange = rowptrS[i]:(rowptrS[i+1]-1)
        nc = length(colrange)
        (nc > 0) || continue
        k = 1
        for j in colrange
            colA = colvalS[j]
            colvalTemp[k] = pJ[colA]
            nzvalTemp[k] = nzvalS[j]
            k += 1
        end
        sortperm!(pointer_to_array(pointer(colperm), nc), pointer_to_array(pointer(colvalTemp), nc))
        k = 1
        for j in colrange
            kperm = colperm[k]
            colvalS[j] = colvalTemp[kperm]
            nzvalS[j] = nzvalTemp[kperm]
            k += 1
        end
    end
    S
end

function getindex_general{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector,
                                 J::AbstractVector)
    pJ = sortperm(J)
    @inbounds J = J[pJ]
    permute_cols!(getindex_J_sorted(A, I, J), pJ)
end

# the general case:
function getindex{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector,
                         J::AbstractVector)
    (m, n) = size(A)

    if !isempty(J)
        minj, maxj = extrema(J)
        ((minj < 1) || (maxj > n)) && throw(BoundsError())
    end

    if !isempty(I)
        mini, maxi = extrema(I)
        ((mini < 1) || (maxi > m)) && throw(BoundsError())
    end

    if isempty(I) || isempty(J) || (0 == nnz(A))
        return spzeros(Tv, Ti, length(I), length(J))
    end

    if issorted(I)
        return getindex_J_sorted(A, I, J)
    else
        return getindex_general(A, I, J)
    end
end

# logical getindex - the methods immediately below are to avoid ambiguity errors
getindex(::SparseMatrixCSR, ::AbstractVector{Bool}, ::Range{Bool}) = error("Cannot index with Range{Bool}")
getindex{T<:Integer}(::SparseMatrixCSR, ::AbstractVector{T}, ::Range{Bool}) = error("Cannot index with Range{Bool}")
getindex(::SparseMatrixCSR, ::Range{Bool}, ::Range{Bool}) = error("cannot index with Range{Bool}")
getindex{T<:Integer}(::SparseMatrixCSR, ::Range{T}, ::Range{Bool}) = error("Cannot index with Range{Bool}")

getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::Range{T}) = A[find(I), J]
getindex{T<:Integer}(A::SparseMatrixCSR, I::Range{T}, J::AbstractVector{Bool}) = A[I,find(J)]
getindex(A::SparseMatrixCSR, I::Integer, J::AbstractVector{Bool}) = A[I,find(J)]
getindex(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::Integer) = A[find(I),J]
getindex(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = A[find(I),find(J)]
getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{T}, J::AbstractVector{Bool}) = A[I,find(J)]
getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::AbstractVector{T}) = A[find(I),J]

function getindex{Tv}(A::SparseMatrixCSR{Tv}, I::AbstractArray{Bool})
    checkbounds(A, I)
    n = sum(I)

    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval
    rowptrB = Int[1,n+1]
    colvalB = Array(Int, n)
    nzvalB = Array(Tv, n)
    c = 1
    colB = 1

    @inbounds for row in 1:A.m
        c1 = rowptrA[row]
        c2 = rowptrA[row+1]-1

        for col in 1:A.n
            if I[row, col]
                while (c1 <= c2) && (colvalA[c1] < col)
                    c1 += 1
                end
                if (c1 <= c2) && (colvalA[c1] == col)
                    nzvalB[c] = nzvalA[c1]
                    colvalB[c] = colB
                    c += 1
                end
                colB += 1
                (colB > n) && break
            end
        end
        (colB > n) && break
    end
    rowptrB[end] = c
    n = length(nzvalB)
    if n > (c-1)
        deleteat!(nzvalB, c:n)
        deleteat!(colvalB, c:n)
    end
    SparseMatrixCSR(1, n, rowptrB, colvalB, nzvalB)
end

function getindex{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractArray)
    szA = size(A);
    nA = szA[1] * szA[2]
    rowptrA = A.rowptr
    colvalA = A.colval
    nzvalA = A.nzval

    n = length(I)
    outm = size(I, 1)
    outn = size(I, 2)
    szB = (outm, outn)
    rowptrB = zeros(Int, outm+1)
    colvalB = Array(Int, n)
    nzvalB = Array(Tv, n)

    rowB = colB = 1
    rowptrB[rowB] = 1
    idxB = 1

    @inbounds for i in 1:n
        ((I[i] < 1) | (I[i] > nA)) && throw(BoundsError())
        rowA, colA = ind2sub(szA, I[i])
        rowB, colB = ind2sub(szB, i)

        for c in rowptrA[rowA]:(rowptrA[rowA+1]-1)
            if colvalA[c] == colA
                rowptrB[rowB+1] += 1
                colvalB[idxB] = colB
                nzvalB[idxB] = nzvalA[c]
                idxB += 1
                break
            end
        end
    end

    rowptrB = cumsum(rowptrB)

    if n > (idxB-1)
        deleteat!(nzvalB, idxB:n)
        deleteat!(colvalB, idxB:n)
    end

    SparseMatrixCSR(outm, outn, rowptrB, colvalB, nzvalB)
end

Base.size(S::SparseMatrixCSR) = (S.m, S.n)
Base.nnz(S::SparseMatrixCSR) = Int(S.rowptr[end] - 1)

function A_mul_B!{T}(α::Number, A::SparseMatrixCSR, B::StridedVecOrMat,
                     β::Number, C::StridedVecOrMat{T})
    A.m == size(C, 1) || throw(DimensionMismatch("size(C, 1) != size(A, 1)"))
    A.n == size(B, 1) || throw(DimensionMismatch("size(B, 1) != size(A, 2)"))
    size(C, 2) == size(B, 2) || throw(DimensionMismatch("size(C, 2) != size(B, 2)"))

    @inbounds begin
        for k=1:size(C, 2)
            for row=1:A.m
                temp = zero(T)
                @simd for j=A.rowptr[row]:A.rowptr[row+1]-1
                    temp += B[A.colval[j], k] * A.nzval[j]
                end
                C[row, k] = α*temp + β*C[row, k]
            end
        end
    end  # inbounds
    C
end

# API to be consistent with Base dense operations
A_mul_B!{T}(C::StridedVecOrMat{T}, A::SparseMatrixCSR, B::StridedVecOrMat) =
    A_mul_B!(one(T), A, B, zero(T), C)

function (*){T,S}(A::SparseMatrixCSR{T}, x::StridedVector{S})
    TS = promote_type(Base.LinAlg.arithtype(T), Base.LinAlg.arithtype(S))
    A_mul_B!(similar(x, TS, size(A, 1)), A, convert(AbstractVector{TS}, x))
end

function (*){T,S}(A::SparseMatrixCSR{T}, B::StridedMatrix{S})
    TS = promote_type(Base.LinAlg.arithtype(T), Base.LinAlg.arithtype(S))
    A_mul_B!(similar(B, TS, size(A, 1), size(B, 2)), A,
             convert(AbstractMatrix{TS}, B))
end

function row_kron{TvA,TiA,TvB,TiB}(A::SparseMatrixCSR{TvA,TiA},
                                   B::SparseMatrixCSR{TvB,TiB})
    TvS = promote_type(Base.LinAlg.arithtype(TvA), Base.LinAlg.arithtype(TvB))
    TiS = promote_type(TiA, TiB)

    mA, nA = size(A)
    mB, nB = size(B)

    mA == mB || error(DimensionMismatch("A and B must have same number of rows"))

    rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval
    rowptrB = B.rowptr; colvalB = B.colval; nzvalB = B.nzval

    rowptrS = Array(TiS, mA+1)
    colvalS = Array(TiS, nA*nB)  # overallocate then deleteat! to avoid cache misses
    nzvalS = Array(TvS, nA*nB)

    rowptrS[1] = 1
    ptrS = 1
    @inbounds for row = 1:mA
        cindA = rowptrA[row]:rowptrA[row+1]-1
        cindB = rowptrB[row]:rowptrB[row+1]-1

        for ca in cindA
            for cb in cindB
                colvalS[ptrS] = (colvalA[ca]-1)*nB + colvalB[cb]
                nzvalS[ptrS] = nzvalB[cb]*nzvalA[ca]
                ptrS += 1
            end
        end

        rowptrS[row+1] = ptrS
    end

    if nA*nB > (ptrS-1)
        deleteat!(nzvalS, ptrS:nA*nB)
        deleteat!(colvalS, ptrS:nA*nB)
    end

    SparseMatrixCSR(mA, nA*nB, rowptrS, colvalS, nzvalS)
end


text:
module GenericLinearAlgebra

import LinearAlgebra: mul!, ldiv!

include("juliaBLAS.jl")
include("lapack.jl")
include("cholesky.jl")
include("householder.jl")
include("qr.jl")
include("eigenSelfAdjoint.jl")
include("eigenGeneral.jl")
include("tridiag.jl")
include("svd.jl")
include("rectfullpacked.jl")
end


text:
module ZigZagBoomerang
using Random
using Requires
using Statistics
using ProgressMeter
using RandomNumbers.Xorshifts
using RandomNumbers: gen_seed
const Rng = Xoroshiro128Plus
Seed() = gen_seed(UInt64, 2)

#using AbstractMCMC

# ZigZag1d and Boomerang1d reference implementation
include("types.jl")
include("common.jl")
include("dynamics.jl")
export ZigZag1d, Boomerang1d, ZigZag, FactBoomerang
const LocalZigZag = ZigZag
export LocalZigZag, BouncyParticle, Boomerang
#include("laplace.jl")

include("poissontime.jl")
export poisson_time
include("jointflow.jl")
include("fact_samplers.jl")
include("not_fact_samplers.jl")

include("priorityqueue.jl")
#const SPriorityQueue = PriorityQueue
include("sfact.jl")
include("local.jl")

include("parallel.jl")
include("sfactiter.jl")
include("notfactiter.jl")

include("zigzagboom1d.jl")
export pdmp, spdmp, eventtime, eventposition

include("ss_fact.jl")
export sspdmp

include("ss_not_fact.jl")
export sticky_pdmp

include("trace.jl")
include("condition.jl")
include("discretise.jl")
const discretise = discretize
export discretise, discretize, sdiscretize, subtrace

function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" include("staticarrays.jl")
end

end


text:
module DescriptorSystems

using LinearAlgebra
using MatrixEquations
using MatrixPencils
using Polynomials
using Random

import LinearAlgebra: BlasFloat, BlasReal, BlasComplex, copy_oftype, transpose, adjoint, opnorm, normalize, rdiv!
import Base: +, -, *, /, \, (==), (!=), ^, isapprox, iszero, convert, promote_op, size, length, ndims,
             hcat, vcat, hvcat, inv, show, lastindex, require_one_based_indexing, print, show, one, zero, eltype
import MatrixPencils: isregular, rmeval
import Polynomials: AbstractRationalFunction, AbstractPolynomial, poles, isconstant, variable, degree, pqs

export DescriptorStateSpace, AbstractDescriptorStateSpace, dss, dssdata, rdss, rss, iszero, order
export RationalTransferFunction, rtf
export gminreal, gir, gir_lrtran, gbalmr, gsvselect, gss2ss, dss2ss
export gbilin, c2d, dss2rm, dss2pm, timeresp
export confmap, simplify, normalize, poles, gain, zpk, rtfbilin, numpoly, denpoly, isconstant, sampling_time
export blockdiag, eye, rcond
export gdual, ctranspose, inv, ldiv, rdiv, ginv
export append, series, parallel, horzcat, vertcat
export order, evalfr, dcgain, opnorm, freqresp, chess
export gpole, gzero, gpoleinfo, gzeroinfo, gnrank, isregular, isproper, isstable,
       glinfnorm, ghinfnorm, gl2norm, gh2norm, ghanorm, gnugap
export gsdec, grnull, glnull, grange, gcrange, grsol, glsol, grmcover1, grmcover2, glmcover1, glmcover2
export grcf, glcf, grcfid, glcfid, gnrcf, gnlcf, giofac, goifac, grsfg, glsfg
export gnehari, glinfldp, glasol, grasol
# export PencilStateSpace, pss, pssdata

abstract type AbstractDynamicalSystem end
abstract type AbstractLTISystem <: AbstractDynamicalSystem end
abstract type AbstractDescriptorStateSpace <: AbstractLTISystem end

include("types/DescriptorStateSpace.jl")
include("types/RationalTransferFunction.jl")
#include("types/PencilStateSpace.jl")
include("ginv.jl")
include("dss.jl")
include("connections.jl")
include("operations.jl")
include("conversions.jl")
include("order_reduction.jl")
include("analysis.jl")
include("decompositions.jl")
include("factorizations.jl")
include("covers.jl")
include("linsol.jl")
include("nullrange.jl")
include("model_matching.jl")
include("dstools.jl")
include("dsutils.jl")
end


text:
using Jacobi, Plots
gr(size=(700,700))

ProjDir = dirname(@__FILE__)
cd(ProjDir) #do

  if !isdefined(Main, :lgt)
    function lgt(Q, nx, lb, ub)
      z = zglj(Q)
      x = collect(linspace(lb, ub, nx))
      y = zeros(nx, Q)
      for k = 1:Q, i=1:nx
        y[i,k]=lagrange(k,x[i],z)
      end
      (x, y)
    end
  end

  Q = 5
  nx = 201
  lb = -1
  ub = 1
  (x, y) = lgt(Q, nx, lb, ub)

  p = plot()
  for k = 1:Q
    plot!(x, y[:, k], label="Q = $k")
  end

  plot(p)
  savefig("lagrange_n201.png")
  gui()

  Q = 5
  nx = 3
  lb = -1
  ub = 1
  (x, y) = lgt(Q, nx, lb, ub)
  x |> display
  println()
  y |> display
  println()

  p = plot()
  for k = 1:Q
    plot!(x, y[:, k], label="Q = $k")
  end

  plot(p)
  savefig("lagrange_n3.png")
  gui()

#end


text:
# This file was generated, do not modify it. # hide
hpn  = xx.Horsepower
Xnew = DataFrame(hp1=hpn, hp2=hpn.^2, hp3=hpn.^3)

yy1 = predict(lr1, Xnew)
yy2 = predict(lr2, Xnew)
yy3 = predict(lr3, Xnew)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy1, lw=3, label="Order 1")
plot(xx.Horsepower, yy2, lw=3, label="Order 2")
plot(xx.Horsepower, yy3, lw=3, label="Order 3")

legend(fontsize=14)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g3.svg")) # hide


text:
using Combinatorics

Tile = Matrix{Bool}

function extracttiles(lines::Vector{<:AbstractString})
    sep = findall(x -> ':' ∈ x, lines)
    # 6:end - 1 because the line is "Tile <id>:"
    ids = parse.(Int, [l[6:end - 1] for l in lines[sep]])
    # -2 because there is a blank line between the tiles
    ends = sep[2:end] .- 2
    push!(ends, length(lines))
    d = Dict{Int,Tile}()
    # the tiles are 10 by 10, hence the hardcoded 10 below
    for (id, start, stop) in zip(ids, sep .+ 1, ends)
        tile = lines[start:stop]
        d[id] = [t[i] == '#' for t in tile, i = 1:10]
    end
    return d
end

function pprint(tile::Tile)
    for l in eachrow(tile)
        line = join([i == false ? "░░" : "██" for i in l])
        println(line)
    end
end

function makeallpositions(tile::Tile)
    m = Tile[tile]
    # Rotations
    for i in 1:3
        push!(m, rotl90(tile, i))
    end
    # Flips
    reversed = reverse(tile, dims=1)
    push!(m, reversed)
    push!(m, rotl90(reversed, 1))
    push!(m, rotl90(reversed, 3))
    push!(m, reverse(tile, dims=2))
    return m
end

rightside(t::Tile) = view(t, :, 10)
leftside(t::Tile) = view(t, :, 1)
upside(t::Tile) = view(t, 1, :)
downside(t::Tile) = view(t, 10, :)

function findmatch(t1::Tile, t2::Tile)
    facing = [(leftside, rightside), (rightside, leftside), (upside, downside), (downside, upside)]
    for (op1, op2) in facing
        if op1(t1) == op2(t2)
            return true
        end
    end
    return false
end

function allpositions(tiles::Dict{Int,Tile})
    alltiles = Dict{Int,Vector{Tile}}()
    for (id, tile) in tiles
        alltiles[id] = makeallpositions(tile)
    end
    return alltiles
end

function findallmatches(all_tiles::Dict{Int,Vector{Tile}})
    # Init dictionary of matching tiles
    matches = Dict{Int,Vector{Int}}()
    for id in keys(all_tiles)
        matches[id] = Int[]
    end
    # The hard way, check all combinations of all possible positions
    for (id1, id2) in combinations(collect(keys(all_tiles)), 2)
        for t1 in all_tiles[id1], t2 in all_tiles[id2]
            if findmatch(t1, t2)
                push!(matches[id1], id2)
                push!(matches[id2], id1)
                break
            end
        end
    end
    return matches
end

function getcorners(matches::Dict{Int,Vector{Int}})
    return [k for k in keys(matches) if length(matches[k]) == 2]
end

const data = readlines("day20.txt")
const tiles = extracttiles(data)
const alltiles = allpositions(tiles)
const matches = findallmatches(alltiles)
const corners = getcorners(matches)

println("Part 1: $(prod(corners))")

# We managed to make part 1 without making the grid, but we need to go there for part 2.
# one grid for the actual tiles in the correct position, and one with the id, so we can
# use our already computed matches dictionary.
const grid = Matrix{Tile}(undef, 12, 12)
const gridid = Matrix{Int}(undef, 12, 12)

function removeid!(matches, id)
    # remove from matches when the tile is set on the grid, so we do not need to look it
    # up anymore.
    for (i, tileids) in matches
        matches[i] = [t for t in tileids if t != id]
    end
    return matches
end

function settileingrid!(pos, id, tile, grid, gridid, matches, alltiles)
    grid[pos...] = tile
    gridid[pos...] = id
    alltiles[id] = Tile[tile]
    removeid!(matches, id)
end

function startingpoint!(grid, gridid, alltiles, matches, corners)
    # select arbitrarily one corner
    id = corners[1]
    gridid[1, 1] = id
    # it is a corner so we know it can only be connected to 2 tiles
    id1, id2 = matches[id]
    for t in alltiles[id], t1 in alltiles[id1], t2 in alltiles[id2]
        if rightside(t) == leftside(t1) && downside(t) == upside(t2)
            settileingrid!((1, 1), id, t, grid, gridid, matches, alltiles)
            settileingrid!((1, 2), id1, t1, grid, gridid, matches, alltiles)
            settileingrid!((2, 1), id2, t2, grid, gridid, matches, alltiles)
        elseif rightside(t) == leftside(t2) && downside(t) == upside(t1)
            settileingrid!((1, 1), id, t, grid, gridid, matches, alltiles)
            settileingrid!((1, 2), id2, t2, grid, gridid, matches, alltiles)
            settileingrid!((2, 1), id1, t1, grid, gridid, matches, alltiles)
        end
    end
end

function setnextright!(pos, grid, gridid, alltiles, matches)
    # set the next id and tile in the grids for a given position
    i, j = pos
    lefttile = j > 1 ?   grid[i, j - 1] : nothing
    leftid   = j > 1 ? gridid[i, j - 1] : nothing
    uptile   = i > 1 ?   grid[i - 1, j] : nothing
    upid     = i > 1 ? gridid[i - 1, j] : nothing
    possiblenexts = leftid !== nothing ? matches[leftid] : matches[upid]
    for id1 in possiblenexts
        # shortcut here if the related id is not on the same line
        if length(matches[id1]) > 3
            continue
        end
        for t1 in alltiles[id1]
            if lefttile === nothing
                check = downside(uptile) == upside(t1)
            elseif uptile === nothing
                check = rightside(lefttile) == leftside(t1)
            else
                check = downside(uptile) == upside(t1) && rightside(lefttile) == leftside(t1)
            end
            if check
                settileingrid!(pos, id1, t1, grid, gridid, matches, alltiles)
                return nothing
            end
        end
    end
    return nothing
end

function fillgrids!(grid, gridid, matches, alltiles, corners)
    x, y = size(grid)
    startingpoint!(grid, gridid, alltiles, matches, corners)
    # The first line, we start at 3
    for i = 3:y
        setnextright!((1, i), grid, gridid, alltiles, matches)
    end
    # The 2nd line, we start at 2
    for i = 2:y
        setnextright!((2, i), grid, gridid, alltiles, matches)
    end
    # For all the other lines, we start at 1
    for j in 3:x, i in 1:y
        setnextright!((j, i), grid, gridid, alltiles, matches)
    end
    return nothing
end

function makeimage(grid)
    x, y = size(grid)
    grid2 = Tile(undef, x * 8, y * 8)
    for i in 1:x, j in 1:y
        grid2[(8 * (i - 1) + 1):(8 * (i - 1) + 8), (8 * (j - 1) + 1):(8 * (j - 1) + 8)] = grid[i, j][2:9, 2:9]
    end
    return grid2
end

function createseamonster()
    seamonster = split("""
                      #
    #    ##    ##    ###
     #  #  #  #  #  #   """, "\n")
    l = length(seamonster[1])
    seamonstergrid = Tile(undef, 3, l)
    for i in 1:3, j in 1:l
        seamonstergrid[i, j] = seamonster[i][j] == '#'
    end
    return seamonstergrid
end

function scanimage(image::Tile, seamonster::Tile)
    imx, imy = size(image)
    mx, my = size(seamonster)
    numones = count(isone, seamonster)
    scanned = Tile(undef, imx - mx + 1, imy - my + 1)

    for x in 1:(imx - mx + 1), y in 1:(imy - my + 1)
        scanned[x, y] = count(isone, image[x:x + mx - 1, y:y + my - 1] .* seamonster) == numones
    end
    return scanned
end

function part2(grid, gridid, matches, alltiles, corners)
    fillgrids!(grid, gridid, matches, alltiles, corners)
    image = makeimage(grid)
    seamonster = createseamonster()
    allmonsters = makeallpositions(seamonster)
    for monster in allmonsters
        scan = scanimage(image, monster)
        c = count(isone, scan)
        if c > 0
            return count(isone, image) - c * count(isone, monster)
        end
    end
end

println("Part 2: $(part2(grid, gridid, matches, alltiles, corners))")


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EstimationProblem(spatialdata, domain, targetvars)

A spatial estimation problem on a given `domain` in which the
variables to be estimated are listed in `targetvars`. The
data of the problem is stored in `spatialdata`.

## Examples

Create an estimation problem for rainfall precipitation measurements:

```julia
julia> EstimationProblem(spatialdata, domain, :precipitation)
```

Create an estimation problem for precipitation and CO₂:

```julia
julia> EstimationProblem(spatialdata, domain, (:precipitation,:CO₂))
```
"""
struct EstimationProblem{S<:AbstractData,
                         D<:AbstractDomain,
                         M<:AbstractMapper} <: AbstractProblem
  # input fields
  spatialdata::S
  domain::D
  mapper::M
  targetvars::Dict{Symbol,DataType}

  # state fields
  mappings::Dict{Symbol,Dict{Int,Int}}

  function EstimationProblem{S,D,M}(spatialdata, domain, targetvars,
                                    mapper) where {S<:AbstractData,
                                                   D<:AbstractDomain,
                                                   M<:AbstractMapper}
    probvnames = Tuple(keys(targetvars))
    datavnames = Tuple(keys(variables(spatialdata)))
    datacnames = coordnames(spatialdata)

    @assert !isempty(probvnames) && probvnames ⊆ datavnames "target variables must be present in spatial data"
    @assert isempty(probvnames ∩ datacnames) "target variables can't be coordinates"
    @assert ndims(domain) == length(datacnames) "data and domain must have the same number of dimensions"
    @assert coordtype(spatialdata) == coordtype(domain) "data and domain must have the same coordinate type"

    mappings = map(spatialdata, domain, probvnames, mapper)

    new(spatialdata, domain, mapper, targetvars, mappings)
  end
end

function EstimationProblem(spatialdata::S, domain::D, targetvarnames::NTuple{N,Symbol};
                           mapper::M=NearestMapper()) where {S<:AbstractData,
                                                            D<:AbstractDomain,
                                                            M<:AbstractMapper,
                                                            N}
  # build dictionary of target variables
  datavars = variables(spatialdata)
  targetvars = Dict(var => Base.nonmissingtype(T) for (var,T) in datavars if var ∈ targetvarnames)

  EstimationProblem{S,D,M}(spatialdata, domain, targetvars, mapper)
end

function EstimationProblem(spatialdata::S, domain::D, targetvarname::Symbol;
                           mapper::M=NearestMapper()) where {S<:AbstractData,
                                                            D<:AbstractDomain,
                                                            M<:AbstractMapper}
  EstimationProblem(spatialdata, domain, (targetvarname,); mapper=mapper)
end

"""
    data(problem)

Return the spatial data of the estimation `problem`.
"""
data(problem::EstimationProblem) = problem.spatialdata

"""
    domain(problem)

Return the spatial domain of the estimation `problem`.
"""
domain(problem::EstimationProblem) = problem.domain

"""
    mapper(problem)

Return the mapper of the estimation `problem`.
"""
mapper(problem::EstimationProblem) = problem.mapper

"""
    variables(problem)

Return the variable names of the estimation `problem` and their types.
"""
variables(problem::EstimationProblem) = problem.targetvars

"""
    datamap(problem, targetvar)

Return the mapping from domain locations to data locations for the
`targetvar` of the `problem`.
"""
datamap(problem::EstimationProblem, var::Symbol) = problem.mappings[var]

"""
    datamap(problem)

Return the mappings from domain locations to data locations for all
the variables of the `problem`.
"""
datamap(problem::EstimationProblem) = problem.mappings

# ------------
# IO methods
# ------------
function Base.show(io::IO, problem::EstimationProblem)
  dim = ndims(problem.domain)
  print(io, "$(dim)D EstimationProblem")
end

function Base.show(io::IO, ::MIME"text/plain", problem::EstimationProblem)
  vars = ["$var ($T)" for (var,T) in problem.targetvars]
  println(io, problem)
  println(io, "  data:      ", problem.spatialdata)
  println(io, "  domain:    ", problem.domain)
  print(  io, "  variables: ", join(vars, ", ", " and "))
end


text:
module Discord

using OpenTrick
using HTTP
using JSON

const URL = "https://discord.com/api"
const API_VERSION = 8 # Gateway & API version

include("helpers.jl")
include("websocket.jl")
include("heartbeat.jl")
include("event_handler.jl")
include("client.jl")
include("op_handler.jl")

end


text:
# models.jl
# define ODEs
# 5 models: resource dependent, resource independent
# no TL coupling, with TL coupling, with TX and TL coupling
# July 2020, Nadanai Laohakunakorn (nadanai.laohakunakorn@ed.ac.uk)

###############################################
######### 1. Resource dependent model #########
###############################################

function grads(du,u,params,t)

    R = u[1]
    dT = u[2]
    dG = u[3]
    mT = u[4]
    mG = u[5]
    pT = u[6]
    pG = u[7]

    alpha = params[1]
    beta = params[2]
    K = params[3]

    du[1] = -alpha*R/(R+K)*dT*pT-alpha*R/(R+K)*dG*pT-beta*R/(R+K)*mT-beta*R/(R+K)*mG
    du[2] = 0
    du[3] = 0
    du[4] = alpha*R/(R+K)*dT*pT
    du[5] = alpha*R/(R+K)*dG*pT
    du[6] = beta*R/(R+K)*mT
    du[7] = beta*R/(R+K)*mG

end

###############################################
######## 2. Resource independent model ########
###############################################

function grads_RI(du,u,params,t)

    R = u[1]
    dT = u[2]
    dG = u[3]
    mT = u[4]
    mG = u[5]
    pT = u[6]
    pG = u[7]
    TX = u[8]
    TL = u[9]

    alpha = params[1]
    beta = params[2]
    lamb1 = params[3]
    lamb2 = params[4]
    KTL = params[5]

    du[1] = 0
    du[2] = 0
    du[3] = 0
    du[4] = alpha*TX*dT*pT
    du[5] = alpha*TX*dG*pT
    du[6] = beta*TL*mT/(mT+mG+KTL)
    du[7] = beta*TL*mG/(mT+mG+KTL)
    du[8] = -lamb1*TX
    du[9] = -lamb2*TL

end

###############################################
## 3. Resource independent model  no loading ##
###############################################

function grads_RInl(du,u,params,t)

    R = u[1]
    dT = u[2]
    dG = u[3]
    mT = u[4]
    mG = u[5]
    pT = u[6]
    pG = u[7]
    TX = u[8]
    TL = u[9]

    alpha = params[1]
    beta = params[2]
    lamb1 = params[3]
    lamb2 = params[4]
    KTL = params[5]

    du[1] = 0
    du[2] = 0
    du[3] = 0
    du[4] = alpha*TX*dT*pT
    du[5] = alpha*TX*dG*pT
    du[6] = beta*TL*mT/(mT+KTL)
    du[7] = beta*TL*mG/(mG+KTL)
    du[8] = -lamb1*TX
    du[9] = -lamb2*TL

end

###############################################
######## 4. Coupled model TL saturation #######
###############################################

function grads_1RI(du,u,params,t)

    R = u[1]
    dT = u[2]
    dG = u[3]
    mT = u[4]
    mG = u[5]
    pT = u[6]
    pG = u[7]
    TX = u[8]
    TL = u[9]

    alpha = params[1]
    beta = params[2]
    lamb1 = params[3]
    lamb2 = params[4]
    KTL = params[5]
    K = params[6]

    du[1] = -alpha*TX*dT*pT*R/(R+K)-alpha*TX*dG*pT*R/(R+K)-beta*TL*mT/(mT+mG+KTL)*R/(R+K)-beta*TL*mG/(mT+mG+KTL)*R/(R+K)
    du[2] = 0
    du[3] = 0
    du[4] = alpha*TX*dT*pT*R/(R+K)
    du[5] = alpha*TX*dG*pT*R/(R+K)
    du[6] = beta*TL*mT/(mT+mG+KTL)*R/(R+K)
    du[7] = beta*TL*mG/(mT+mG+KTL)*R/(R+K)
    du[8] = -lamb1*TX
    du[9] = -lamb2*TL

end

###############################################
####### 5. Coupled model TXTL saturation ######
###############################################

function grads_1RIf(du,u,params,t)

    R = u[1]
    dT = u[2]
    dG = u[3]
    mT = u[4]
    mG = u[5]
    pT = u[6]
    pG = u[7]
    TX = u[8]
    TL = u[9]

    alpha = params[1]
    beta = params[2]
    lamb1 = params[3]
    lamb2 = params[4]
    KTL = params[5]
    K = params[6]
    KTX = params[7]

    du[1] = -alpha*TX*dT/(dT+dG+KTX)*pT*R/(R+K)-alpha*TX*dG/(dT+dG+KTX)*pT*R/(R+K)-beta*TL*mT/(mT+mG+KTL)*R/(R+K)-beta*TL*mG/(mT+mG+KTL)*R/(R+K)
    du[2] = 0
    du[3] = 0
    du[4] = alpha*TX*dT/(dT+dG+KTX)*pT*R/(R+K)
    du[5] = alpha*TX*dG/(dT+dG+KTX)*pT*R/(R+K)
    du[6] = beta*TL*mT/(mT+mG+KTL)*R/(R+K)
    du[7] = beta*TL*mG/(mT+mG+KTL)*R/(R+K)
    du[8] = -lamb1*TX
    du[9] = -lamb2*TL

end


text:
# TODO: respect ordering
module AtomixCUDA

using Atomix: AtomicRef
using CUDA: CUDA, CuArray, CuDeviceArray
using UnsafeAtomics: UnsafeAtomics, Ordering

const GenericCuArray{T,N} = Union{CuArray{T,N},CuDeviceArray{T,N}}
const CuAtomicRef{T,Ptr,Data<:GenericCuArray} = AtomicRef{T,Ptr,Data}

function UnsafeAtomics.load(ref::CuAtomicRef{T}, ord::Ordering) where {T}
    error("not implemented")
end

function UnsafeAtomics.store!(ref::CuAtomicRef{T}, v, ord::Ordering) where {T}
    error("not implemented")
end

@inline function UnsafeAtomics.cas!(
    ref::CuAtomicRef{T},
    expected,
    desired,
    success_ordering::Ordering,
    failure_ordering::Ordering,
) where {T}
    ptr = ref.ptr
    expected = convert(T, expected)
    desired = convert(T, desired)
    data = ref.data
    GC.@preserve data begin
        old = CUDA.atomic_cas!(ptr, expected, desired)
    end
    return (; old = old, success = old === expected)
end

@inline function UnsafeAtomics.modify!(ref::CuAtomicRef{T}, op::OP, x, ord) where {T,OP}
    x = convert(T, x)
    ptr = ref.ptr
    data = ref.data
    GC.@preserve data begin
        old = if op === (+)
            CUDA.atomic_add!(ptr, x)
        elseif op === (-)
            CUDA.atomic_sub!(ptr, x)
        elseif op === (&)
            CUDA.atomic_and!(ptr, x)
        elseif op === (|)
            CUDA.atomic_or!(ptr, x)
        elseif op === xor
            CUDA.atomic_xor!(ptr, x)
        elseif op === min
            CUDA.atomic_min!(ptr, x)
        elseif op === max
            CUDA.atomic_max!(ptr, x)
        else
            error("not implemented")
        end
    end
    return (old, op(old, x))
end

end  # module AtomixCUDA


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

using Markdown
using InteractiveUtils

# ╔═╡ 794a5356-6e05-44b4-ad9d-375a1fa3613b
begin
	using PlutoUI
	using DataFrames
	using CSV

	using Arrow
	using Serialization
	using JLSO
	using JSONTables
	using CodecZlib
	using ZipFile
	using JDF

	using StatsPlots # for charts
	using Mmap       # for compression
	using Test
end

# ╔═╡ b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
md"""
## Introduction DataFrames / 4 - Load and Save

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/04_loadsave.ipynb)

We do not cover all features of the packages. Please refer to their documentation to learn them.

Here we will load `CSV.jl` to read and write CSV files and `Arrow.jl`, `JLSO.jl`, and `serialization`, which allow us to work with a binary format and `JSONTables.jl` for JSON interaction. Finally we consider a custom `JDF.jl` format.

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ d6993a5b-8158-4ce3-accb-8e8c514229b3
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ de97ded0-3952-4188-9a43-a8f383e1ae65
md"""
### Load and save DataFrames

Let us create a simple DataFrame for testing purpose :
"""

# ╔═╡ 887ea560-b824-44bd-96ba-7756271993e7
x₀ = DataFrame(A=[true, false, true], B=[1, 2, missing],
	C=[missing, "b", "c"], D=['a', missing, 'c'])

# ╔═╡ 327fae3b-f00a-4bb4-ba95-a003e0740c58
md"""
and use `eltypes` to look at the columnwise types.
"""

# ╔═╡ b18bb84d-99f2-454e-b8bd-7c5e52b94141
eltype.(eachcol(x₀))

# ╔═╡ ca9927e5-3ff2-4c73-8206-4d7495668312
md"""
### CSV.jl

Let's use `CSV.jl` to save x₀ to disk; make sure `x0.csv` does not conflict with some file in your working directory.
"""

# ╔═╡ 696d306e-c000-47af-a904-eff748b317f3
CSV.write("./Data/x0.csv", x₀)

# ╔═╡ 157caf96-53ad-4ae6-aa1f-5d3fef0b178f
md"""
Now we can see how it was saved by reading `x0.csv`.
"""

# ╔═╡ c04fb113-c95f-4a23-8884-147cdc62ac5c
with_terminal() do
	print(read("./Data/x0.csv", String))
end

# ╔═╡ 06ee8941-57ef-4dfa-8385-861146755d53
md"""
We can also load it back, as usual running:
"""

# ╔═╡ ae8ac551-ea7f-4a9f-bfb4-e1b86677a6a8
y₀ = CSV.read("./Data/x0.csv", DataFrame)

# ╔═╡ 86253acd-78a2-44c2-931e-f2322ed67dc2
eltype.(eachcol(y₀))

# ╔═╡ 069ea341-4097-42f5-83ff-18548bd34de2
md"""
*Remark:* Note that when loading in a DataFrame from a CSV the column type for column :D has changed!
"""

# ╔═╡ a4615e9b-d3d9-443a-9146-e411f9bbd749
md"""
### Serialization, JDF.jl, and JLSO.j

#### Serialization

Now we use serialization to save x₀.

There are two ways to perform serialization.
  - The first way is to use the Serialization.serialize as below:
    **Note** *that in general, this process will not work if the reading and writing are done by different versions of Julia, or an instance of Julia with a different system image*.
"""

# ╔═╡ 8cc8ae4d-44df-495b-9469-6b467abd6add
open("./Data/x₀.bin", "w") do io
    serialize(io, x₀)
end

# ╔═╡ 8c7e6303-0d09-4993-943f-d9bd3a9c316b
md"""
Now we load back the saved file to `y₁` variable. Again `y₁` is identical to `x₀`. However, please beware that if you session does not have `DataFrames.jl` loaded, then it may not recognise the content as `DataFrames.jl`.
"""

# ╔═╡ aa8ffd59-e791-4dca-bce7-05a18dfda955
y₁ = open(deserialize, "./Data/x₀.bin")

# ╔═╡ b67f1614-90e2-4e4f-9890-9fd66d24cf5b
eltype.(eachcol(y₁))

# ╔═╡ 3d64997f-bfed-4c48-b962-d0bd6c7aa1ef
md"""
#### JDF.jl

[`JDF.jl`](https://github.com/xiaodaigh/JDF) is a relatively new package designed to serialize DataFrames. You can save a DataFrame with the `savejdf` function.
"""

# ╔═╡ 9be83f87-de31-4403-b567-4f488f88645e
savejdf("./Data/x₀.jdf", x₀);

# ╔═╡ 2717467d-c2c5-4006-bc6f-f5c5de3c6304
md"""
To load the saved JDF file, we use the `loadjdf` function
"""

# ╔═╡ 1f9c7d56-3471-4679-aa17-c5ba3b7c31c6
x₀_loaded = loadjdf("./Data/x₀.jdf") |> DataFrame

# ╔═╡ 332e66fe-5a54-4741-a09c-465afef318b7
@assert isequal(x₀_loaded, x₀)

# ╔═╡ 494e54f6-f29c-457a-9aba-0b7a3cd143b9
md"""
`JDF.jl` offers the ability to load only certain columns from disk to help with working with large files
"""

# ╔═╡ 2909e4a4-5355-4ca9-aca2-efc43d938b29
# set up a JDFFile which is a on disk representation of `x₀` backed by JDF.jl
x₀_ondisk = jdf"./Data/x₀.jdf"

# ╔═╡ 8beffe4b-051a-41e9-af93-10329564355f
md"""
We can see all the names of `x` without loading it into memory
"""

# ╔═╡ 43b9a909-dd99-4bf0-bf41-8a21c5127e24
names(x₀_ondisk)

# ╔═╡ 424dc5a0-f64d-4563-a8be-99ffbb1b211d
md"""
The below is an example of how to load only columns `:A` and `:D`
"""

# ╔═╡ 6109fe95-9cd4-47b7-a79d-387e7aac78db
xd = sloadjdf(x₀_ondisk; cols = ["A", "D"]) |> DataFrame

# ╔═╡ 7211a7cc-3ab1-4c4e-a858-60ee8d910054
eltype.(eachcol(xd))

# ╔═╡ c7ea6a2b-8719-4a1a-b2ce-34237bbbc625
md"""
#### JDF.jl vs others

`JDF.jl` is specialized to `DataFrames.jl` and only supports a restricted list of columns, so it cannot save dataframes with arbitrary column types. However, this also means that `JDF.jl` has specialised algorithms to serailize the type it supports to optimize speed, minimize disk usage, and reduce the chance of errors

The list support columns for `JDF` include

```Julia
WeakRefStrings.StringVector
Vector{T}, Vector{Union{Mising, T}}, Vector{Union{Nothing, T}}
CategoricalArrays.CategoricalVetors{T}
```

where `T` can be `String`, `Bool`, `Symbol`, `Char`, `TimeZones.ZonedDateTime` (experimental) and `isbits` types i.e. `UInt*`, `Int*`, `Float*`, and `Date*` types etc.
"""

# ╔═╡ 2516d88a-9141-4a25-90b3-fc212e6299f9
md"""
#### JLSO.jl

JLSO.jl

Another way to perform serialization is by using the [JLSO.jl](https://github.com/invenia/JLSO.jl) library:
"""

# ╔═╡ 381916dc-3610-4df4-9911-c822d28d3956
JLSO.save("./Data/x₀.jlso", :data => x₀)

# ╔═╡ db07eb7e-2fa7-44e7-adb6-c0254a83e0e3
md"""
Now we can load back the file to y₂.
"""

# ╔═╡ 4f63f693-8392-4e8e-a04e-00f45ba016a3
y₂ = JLSO.load("./Data/x₀.jlso")[:data]

# ╔═╡ a3203a29-33f5-457a-9a90-80c5e14422a6
eltype.(eachcol(y₂))

# ╔═╡ f12dc3e2-167d-4cad-a999-53c67d55af07
md"""
#### JSONTables.jl

Often you might need to read and write data stored in JSON format. `JSONTables.jl` provides a way to process them in row-oriented or column-oriented layout. We present both options below.
"""

# ╔═╡ 23b96b56-709f-4f5b-8bb9-961764f61faf
open(io -> arraytable(io, x₀), "./Data/x₁.json", "w")

# ╔═╡ 6223627b-7541-42f7-892a-80c86e3a95ec
open(io -> objecttable(io, x₀), "./Data/x₂.json", "w")

# ╔═╡ f9281efb-e8aa-4267-9d63-7094e6010315
with_terminal() do
	print(read("./Data/x₁.json", String))
end

# ╔═╡ f215042e-5dec-4e2f-940f-02467c0b3d94
with_terminal() do
	print(read("./Data/x₂.json", String))
end

# ╔═╡ ef097417-b572-4110-b943-7b2500c81838
yj₁ = open(jsontable, "./Data/x₁.json") |> DataFrame

# ╔═╡ 1aca803a-197a-47ac-8eb3-73ab3c40de47
eltype.(eachcol(yj₁))

# ╔═╡ 0e193005-cd0b-4c1c-aa1e-42cee639238c
yj₂ = open(jsontable, "./Data/x₁.json") |> DataFrame

# ╔═╡ 6dd53b4a-a936-40c3-873d-42abe4421189
eltype.(eachcol(yj₂))

# ╔═╡ dc97e056-2d8d-442e-90d2-b75931cecb33
md"""
#### Arrow.jl

Finally we use `Apache Arrow` format that allows, in particular, for data interchange with `R` or `Python`.
"""

# ╔═╡ d1a0d653-ac96-4712-bc26-430ba4436544
Arrow.write("./Data/x₀.arrow", x₀)

# ╔═╡ 2244a140-ec46-4ec5-84c6-fc40511953ad
yₐ = Arrow.Table("./Data/x₀.arrow") |> DataFrame

# ╔═╡ e392701a-75a4-465b-9547-a8d6a3bdc059
eltype.(eachcol(yₐ))

# ╔═╡ 2bf8f141-c809-4599-ace3-3fabd1bbc1d6
yₐ.A[1] = false   ## can be changed - used to be immutable

# ╔═╡ f00b7857-6393-4501-b0a8-bfd7d50667c4
yₐ.A

# ╔═╡ 613569c5-20f4-4463-a039-12fecab721eb
yₐ.B

# ╔═╡ edd10614-b548-4148-acf5-33e16ab0a8be
md"""
### Basic benchmarking

Next, we'll create some files, so be careful that you don't already have these files in your working directory!

In particular, we'll time how long it takes us to write a DataFrame with 10^3 rows and 10^5 columns.

"""

# ╔═╡ 3d478ed4-cf5a-44a1-9441-6da73f7986d2
begin
	bigdf = DataFrame(rand(Bool, 10^5, 500), :auto)
	bigdf[!, 1] = Int.(bigdf[!, 1])
	bigdf[!, 2] = bigdf[!, 2] .+ 0.5
	bigdf[!, 3] = string.(bigdf[!, 3], ", as string")

	bigdf = DataFrame(rand(Bool, 10^5, 500), :auto)
	bigdf[!, 1] = Int.(bigdf[!, 1])
	bigdf[!, 2] = bigdf[!, 2] .+ 0.5
	bigdf[!, 3] = string.(bigdf[!, 3], ", as string")

	with_terminal() do
		println("First run")
		println("CSV.jl")
		global csvwrite1 = @elapsed @time CSV.write("./Data/bigdf1.csv", bigdf)

		println("Serialization")
		global serializewrite1 = @elapsed @time open(io -> serialize(io, bigdf),
			"./Data/bigdf.bin", "w")

		println("JDF.jl")
		global jdfwrite1 = @elapsed @time savejdf("./Data/bigdf.jdf", bigdf)

		println("JLSO.jl")
		global jlsowrite1 = @elapsed @time JLSO.save("./Data/bigdf.jlso", :data => bigdf)

		println("Arrow.jl")
		global arrowwrite1 = @elapsed @time Arrow.write("./Data/bigdf.arrow", bigdf)

		println("JSONTables.jl arraytable")
		global jsontablesawrite1 = @elapsed @time open(io -> arraytable(io, bigdf),
			"./Data/bigdf1.json", "w")

		println("JSONTables.jl objecttable")
		global jsontablesowrite1 = @elapsed @time open(io -> objecttable(io, bigdf),
			"./Data/bigdf2.json", "w")

		println("Second run")
		println("CSV.jl")
		global csvwrite2 = @elapsed @time CSV.write("./Data/bigdf1.csv", bigdf)

		println("Serialization")
		global serializewrite2 = @elapsed @time open(io -> serialize(io, bigdf),
				"./Data/bigdf.bin", "w")

		println("JDF.jl")
		global jdfwrite2 = @elapsed @time savejdf("./Data/bigdf.jdf", bigdf)

		println("JLSO.jl")
		global jlsowrite2 = @elapsed @time JLSO.save("./Data/bigdf.jlso", :data => bigdf)

		println("Arrow.jl")
		global arrowwrite2 = @elapsed @time Arrow.write("./Data/bigdf.arrow", bigdf)

		println("JSONTables.jl arraytable")
		global jsontablesawrite2 = @elapsed @time open(io -> arraytable(io, bigdf),
			"./Data/bigdf1.json", "w")

		println("JSONTables.jl objecttable")
		global jsontablesowrite2 = @elapsed @time open(io -> objecttable(io, bigdf),
			"./Data/bigdf2.json", "w")
	end
end

# ╔═╡ 782f560a-9fd5-498b-9ff8-0215736b240e
groupedbar(
    # Exclude JSONTables.jl arraytable due to timing
	repeat(["CSV.jl", "Serialization", "JDF.jl", "JLSO.jl",
			"Arrow.jl", "JSONTables.jl\nobjecttable"], inner = 2),

	[csvwrite1, csvwrite2, serializewrite1, serializewrite1, jdfwrite1, jdfwrite2,
     jlsowrite1, jlsowrite2, arrowwrite1, arrowwrite2, jsontablesowrite2, jsontablesowrite2],

	group = repeat(["1st", "2nd"], outer = 6),
	ylab = "Second",
	title = "Write Performance\nDataFrame: bigdf\nSize: $(size(bigdf))"
)

# ╔═╡ 5ab7852b-ea4e-42de-829f-a101e58a0b93
begin
	data_files = ["./Data/bigdf1.csv", "./Data/bigdf.bin", "./Data/bigdf.arrow",
		"./Data/bigdf1.json", "./Data/bigdf2.json"]
	df = DataFrame(file = data_files, size = getfield.(stat.(data_files), :size))
	append!(df,
			DataFrame(file = "./Data/bigdf.jdf",
			size=reduce((x,y) -> x + y.size, stat.(joinpath.("./Data/bigdf.jdf", readdir("./Data/bigdf.jdf"))),
				init=0)))
	sort!(df, :size)
end

# ╔═╡ 175391bb-3cdc-409f-8963-8aa3d613b601
@df df plot(:file, :size/1024^2, seriestype=:bar, title = "Format File Size (MB)", label="Size", ylab="MB", legend=:topleft, xrotation=45)

# ╔═╡ 5e475f09-3177-45b5-a9a1-f25c1cfeef37
begin

	with_terminal() do
		println("First run")
		println("CSV.jl")
		global csvread1 = @elapsed @time CSV.read("./Data/bigdf1.csv", DataFrame)

		println("Serialization")
		global serializeread1 = @elapsed @time open(deserialize, "./Data/bigdf.bin")

		println("JDF.jl")
		global jdfread1 = @elapsed @time loadjdf("./Data/bigdf.jdf")

		println("JLSO.jl")
		global jlsoread1 = @elapsed @time JLSO.load("./Data/bigdf.jlso")

		println("Arrow.jl")
		global arrowread1 = @elapsed @time df_tmp = Arrow.Table("./Data/bigdf.arrow") |> DataFrame
		global arrowread1copy = @elapsed @time copy(df_tmp)

		println("JSONTables.jl arraytable")
		global jsontablesaread1 = @elapsed @time open(jsontable, "./Data/bigdf1.json")

		println("JSONTables.jl objecttable")
		global jsontablesoread1 = @elapsed @time open(jsontable, "./Data/bigdf2.json")

		println("Second run")
		global csvread2 = @elapsed @time CSV.read("./Data/bigdf1.csv", DataFrame)

		println("Serialization")
		global serializeread2 = @elapsed @time open(deserialize, "./Data/bigdf.bin")

		println("JDF.jl")
		global jdfread2 = @elapsed @time loadjdf("./Data/bigdf.jdf")

		println("JLSO.jl")
		global jlsoread2 = @elapsed @time JLSO.load("./Data/bigdf.jlso")

		println("Arrow.jl")
		global arrowread2 = @elapsed @time df_tmp = Arrow.Table("./Data/bigdf.arrow") |> DataFrame
		global arrowread2copy = @elapsed @time copy(df_tmp)

		println("JSONTables.jl arraytable")
		global jsontablesaread2 = @elapsed @time open(jsontable, "./Data/bigdf1.json")

		println("JSONTables.jl objecttable")
		global jsontablesoread2 = @elapsed @time open(jsontable, "./Data/bigdf2.json");
	end
end

# ╔═╡ 2fac7472-ddf7-477e-9142-0640bda99890
groupedbar(
    repeat(["CSV.jl", ".Serialization", "JDF.jl", "JLSO.jl", "Arrow.jl",
			"Arrow.jl\ncopy", "JSON\narraytable", "JSON\nobjecttable"], inner = 2),

    [csvread1, csvread2, serializeread1, serializeread2, jdfread1, jdfread2, jlsoread1, jlsoread2,
     arrowread1, arrowread2, arrowread1+arrowread1copy, arrowread2+arrowread2copy,
     jsontablesaread1, jsontablesaread2, jsontablesoread1, jsontablesoread2],

	group = repeat(["1st", "2nd"], outer = 8),
    ylab = "Second",
    title = "Read Performance\nDataFrame: bigdf\nSize: $(size(bigdf))",
	legend=:topleft
)

# ╔═╡ d2df242d-259c-4935-b71b-e1d444ac6f07
md"""
#### Using gzip compression

A common user requirement is to be able to load and save CSV that are compressed using `gzip`. Below we show how this can be accomplished using `CodecZlib.jl`. The same pattern is applicable to `JSONTables.jl` compression/decompression.

Again make sure that you do not have file named df_compress_test.csv.gz in your working directory

We first generate a random data frame

"""

# ╔═╡ b7b657d6-0723-4209-9ecb-efcc72be402f
dfₓ = DataFrame(rand(1:10, 10, 1000), :auto)

# ╔═╡ 8357bef1-d31c-4fa5-992f-f2a6eb6d9573
# GzipCompressorStream comes from CodecZlib

open("./Data/df_compress_test.csv.gz", "w") do io
    stream = GzipCompressorStream(io)
    CSV.write(stream, dfₓ)
    close(stream)
end

# ╔═╡ c33160d1-b12e-42ab-a3c3-b8d6bc58bd9d
df₂ = CSV.File(transcode(GzipDecompressor,
		Mmap.mmap("./Data/df_compress_test.csv.gz"))) |> DataFrame

# ╔═╡ 544fbeae-ca73-487d-92b5-8d619d9b8db7
@assert dfₓ == df₂

# ╔═╡ bcee6550-63e7-4ff2-a324-37bef45a28b9
md"""
#### Using zip files

Sometimes you may have files compressed inside a zip file.

In such a situation you may use `ZipFile.jl` in conjunction an an appropriate reader to read the files.

Here we first create a ZIP file and then read back its contents into a DataFrame.
"""

# ╔═╡ 50eee498-9548-4ee7-8375-90c17b376c0a
zdf₁ = DataFrame(rand(1:10, 3, 4), :auto)

# ╔═╡ 90eb146e-9dbf-45ef-a093-437fb938b10d
zdf₂ = DataFrame(rand(1:10, 3, 4), :auto)

# ╔═╡ f5cb08ba-814e-4a7b-b676-5967a352d2a7
md"""
And we show yet another way to write a DataFrame into a CSV file
"""

# ╔═╡ 223e5415-0ab4-4016-b506-658adc75c57b
begin
	# write a CSV file into the zip file
	w = ZipFile.Writer("./Data/x.zip")

	f₁ = ZipFile.addfile(w, "x₁.csv")
	write(f₁, sprint(show, "text/csv", zdf₁))

	# write a second CSV file into zip file
	f₂ = ZipFile.addfile(w, "x₂.csv", method=ZipFile.Deflate)
	write(f₂, sprint(show, "text/csv", zdf₂))

	close(w)
end

# ╔═╡ dd26958e-55f1-4b2e-b557-1bd24eb43fa0
md"""
Now we read the CSV we have written.
"""

# ╔═╡ 8b87edb5-d8a7-4234-ac8d-f9e0c8f5f8a0
z = ZipFile.Reader("./Data/x.zip");

# ╔═╡ f257c2c9-6cdf-4df8-8b53-7a8d3794b7c3
begin
	# find the index index of file called x₁.csv
	ix_xcsv₁ = findfirst(x -> x.name == "x₁.csv", z.files);

	# to read the x₁.csv file in the zip file
	rzdf₁ = read(z.files[ix_xcsv₁]) |> f -> CSV.read(f, DataFrame);
end

# ╔═╡ ec462215-2451-4e8e-b440-58503378faae
@assert rzdf₁ == zdf₁

# ╔═╡ 4fb9645f-ee4d-410a-aec3-06f13fdec761
begin
	# find the index index of file called x₂.csv
	ix_xcsv₂ = findfirst(x->x.name == "x₂.csv", z.files)

	# to read the x2.csv file in the zip file
	rzdf₂ = read(z.files[ix_xcsv₂]) |> f -> CSV.read(f, DataFrame);
end

# ╔═╡ cc833913-383f-46ac-8ff2-5b6c8c8cbff4
@assert rzdf₂ == zdf₂

# ╔═╡ c3f1fae0-4422-4023-b43f-de1651a89983
md"""
Note that once you read a given file from z object its stream is all used-up (it is at its end). Therefore to read it again you need to close z and open it again.

Also do not forget to close the zip file once done.
"""

# ╔═╡ 5a1bd5be-9329-44e8-848e-efb7bac1e5f4
close(z)

# ╔═╡ 1d8e3b27-8af6-4541-923a-a7744b271464
md"""
Finally, let's clean up. Do not run the next cell unless you are sure that it will not erase your important files.
"""

# ╔═╡ d775d1a6-906c-4013-8ab7-5951cb5fe071
foreach(f -> rm(string("./Data/", f), force=true),
	["x0.csv", "x₀.bin", "x₀.jlso", "x₁.json", "x₂.json", "bigdf1.csv", "bigdf.bin",
		"bigdf.jlso", "bigdf1.json", "bigdf2.json", "x.zip"])

# ╔═╡ 2ad38587-5e55-4e24-9c23-a3180dd1bf99
begin
	rm("bigdf.jdf", recursive=true, force=true)
	rm("x₀.jdf", recursive=true, force=true)
end

# ╔═╡ e0fa8b35-19fd-4155-959d-51c2817199e6
md"""
Note that we did not remove `x₀.arrow`, `bigdf.arrow` and `df_compress_test.csv.gz` - you have to do it manually, as these files are memory mapped.
"""

# ╔═╡ c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e
html"""

"""

# ╔═╡ Cell order:
# ╟─b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
# ╠═794a5356-6e05-44b4-ad9d-375a1fa3613b
# ╟─d6993a5b-8158-4ce3-accb-8e8c514229b3
# ╟─de97ded0-3952-4188-9a43-a8f383e1ae65
# ╠═887ea560-b824-44bd-96ba-7756271993e7
# ╟─327fae3b-f00a-4bb4-ba95-a003e0740c58
# ╠═b18bb84d-99f2-454e-b8bd-7c5e52b94141
# ╟─ca9927e5-3ff2-4c73-8206-4d7495668312
# ╠═696d306e-c000-47af-a904-eff748b317f3
# ╟─157caf96-53ad-4ae6-aa1f-5d3fef0b178f
# ╠═c04fb113-c95f-4a23-8884-147cdc62ac5c
# ╟─06ee8941-57ef-4dfa-8385-861146755d53
# ╠═ae8ac551-ea7f-4a9f-bfb4-e1b86677a6a8
# ╠═86253acd-78a2-44c2-931e-f2322ed67dc2
# ╟─069ea341-4097-42f5-83ff-18548bd34de2
# ╟─a4615e9b-d3d9-443a-9146-e411f9bbd749
# ╠═8cc8ae4d-44df-495b-9469-6b467abd6add
# ╟─8c7e6303-0d09-4993-943f-d9bd3a9c316b
# ╠═aa8ffd59-e791-4dca-bce7-05a18dfda955
# ╠═b67f1614-90e2-4e4f-9890-9fd66d24cf5b
# ╟─3d64997f-bfed-4c48-b962-d0bd6c7aa1ef
# ╠═9be83f87-de31-4403-b567-4f488f88645e
# ╟─2717467d-c2c5-4006-bc6f-f5c5de3c6304
# ╠═1f9c7d56-3471-4679-aa17-c5ba3b7c31c6
# ╠═332e66fe-5a54-4741-a09c-465afef318b7
# ╟─494e54f6-f29c-457a-9aba-0b7a3cd143b9
# ╠═2909e4a4-5355-4ca9-aca2-efc43d938b29
# ╟─8beffe4b-051a-41e9-af93-10329564355f
# ╠═43b9a909-dd99-4bf0-bf41-8a21c5127e24
# ╟─424dc5a0-f64d-4563-a8be-99ffbb1b211d
# ╠═6109fe95-9cd4-47b7-a79d-387e7aac78db
# ╠═7211a7cc-3ab1-4c4e-a858-60ee8d910054
# ╟─c7ea6a2b-8719-4a1a-b2ce-34237bbbc625
# ╟─2516d88a-9141-4a25-90b3-fc212e6299f9
# ╠═381916dc-3610-4df4-9911-c822d28d3956
# ╟─db07eb7e-2fa7-44e7-adb6-c0254a83e0e3
# ╠═4f63f693-8392-4e8e-a04e-00f45ba016a3
# ╠═a3203a29-33f5-457a-9a90-80c5e14422a6
# ╟─f12dc3e2-167d-4cad-a999-53c67d55af07
# ╠═23b96b56-709f-4f5b-8bb9-961764f61faf
# ╠═6223627b-7541-42f7-892a-80c86e3a95ec
# ╠═f9281efb-e8aa-4267-9d63-7094e6010315
# ╠═f215042e-5dec-4e2f-940f-02467c0b3d94
# ╠═ef097417-b572-4110-b943-7b2500c81838
# ╠═1aca803a-197a-47ac-8eb3-73ab3c40de47
# ╠═0e193005-cd0b-4c1c-aa1e-42cee639238c
# ╠═6dd53b4a-a936-40c3-873d-42abe4421189
# ╟─dc97e056-2d8d-442e-90d2-b75931cecb33
# ╠═d1a0d653-ac96-4712-bc26-430ba4436544
# ╠═2244a140-ec46-4ec5-84c6-fc40511953ad
# ╠═e392701a-75a4-465b-9547-a8d6a3bdc059
# ╠═2bf8f141-c809-4599-ace3-3fabd1bbc1d6
# ╠═f00b7857-6393-4501-b0a8-bfd7d50667c4
# ╠═613569c5-20f4-4463-a039-12fecab721eb
# ╟─edd10614-b548-4148-acf5-33e16ab0a8be
# ╠═3d478ed4-cf5a-44a1-9441-6da73f7986d2
# ╠═782f560a-9fd5-498b-9ff8-0215736b240e
# ╠═5ab7852b-ea4e-42de-829f-a101e58a0b93
# ╠═175391bb-3cdc-409f-8963-8aa3d613b601
# ╠═5e475f09-3177-45b5-a9a1-f25c1cfeef37
# ╠═2fac7472-ddf7-477e-9142-0640bda99890
# ╟─d2df242d-259c-4935-b71b-e1d444ac6f07
# ╠═b7b657d6-0723-4209-9ecb-efcc72be402f
# ╠═8357bef1-d31c-4fa5-992f-f2a6eb6d9573
# ╠═c33160d1-b12e-42ab-a3c3-b8d6bc58bd9d
# ╠═544fbeae-ca73-487d-92b5-8d619d9b8db7
# ╟─bcee6550-63e7-4ff2-a324-37bef45a28b9
# ╠═50eee498-9548-4ee7-8375-90c17b376c0a
# ╠═90eb146e-9dbf-45ef-a093-437fb938b10d
# ╟─f5cb08ba-814e-4a7b-b676-5967a352d2a7
# ╠═223e5415-0ab4-4016-b506-658adc75c57b
# ╟─dd26958e-55f1-4b2e-b557-1bd24eb43fa0
# ╠═8b87edb5-d8a7-4234-ac8d-f9e0c8f5f8a0
# ╠═f257c2c9-6cdf-4df8-8b53-7a8d3794b7c3
# ╠═ec462215-2451-4e8e-b440-58503378faae
# ╠═4fb9645f-ee4d-410a-aec3-06f13fdec761
# ╠═cc833913-383f-46ac-8ff2-5b6c8c8cbff4
# ╟─c3f1fae0-4422-4023-b43f-de1651a89983
# ╠═5a1bd5be-9329-44e8-848e-efb7bac1e5f4
# ╟─1d8e3b27-8af6-4541-923a-a7744b271464
# ╠═d775d1a6-906c-4013-8ab7-5951cb5fe071
# ╠═2ad38587-5e55-4e24-9c23-a3180dd1bf99
# ╟─e0fa8b35-19fd-4155-959d-51c2817199e6
# ╟─c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks