text:
"""
    mutable struct MDP{T,R}
        observationspace::DiscreteSpace
        actionspace::DiscreteSpace
        state::Int64
        trans_probs::Array{T, 2}
        reward::R
        initialstates::Array{Int64, 1}
        isterminal::Array{Int64, 1}

A Markov Decision Process with `ns` states, `na` actions, current `state`,
`na`x`ns` - array of transition probabilites `trans_props` which consists for
every (action, state) pair of a (potentially sparse) array that sums to 1 (see
[`getprobvecrandom`](@ref), [`getprobvecuniform`](@ref),
[`getprobvecdeterministic`](@ref) for helpers to constract the transition
probabilities) `na`x`ns` - array of `reward`, array of initial states
`initialstates`, and `ns` - array of 0/1 indicating if a state is terminal.
"""
mutable struct MDP{T,R}
    observationspace::DiscreteSpace
    actionspace::DiscreteSpace
    state::Int64
    trans_probs::Array{T, 2}
    reward::R
    initialstates::Array{Int64, 1}
    isterminal::Array{Int64, 1}
end
function MDP(ospace, aspace, state, trans_probs::AbstractArray{T, 2},
             reward::R, initialstates, isterminal) where {T, R}
    if R <: AbstractMatrix
        reward = DeterministicStateActionReward(reward)
    end
    MDP{T,typeof(reward)}(ospace, aspace, state, trans_probs, reward, initialstates, isterminal)
end

function interact!(env::MDP, action)
    oldstate = env.state
    run!(env, action)
    r = reward(env.reward, oldstate, action, env.state)
    (observation = env.state, reward = r, isdone = env.isterminal[env.state] == 1)
end
function getstate(env::MDP)
    (observation = env.state, isdone = env.isterminal[env.state] == 1)
end
function reset!(env::MDP)
    env.state = rand(ENV_RNG, env.initialstates)
    (observation = env.state, )
end

"""
    struct DeterministicNextStateReward
        value::Vector{Float64}
"""
struct DeterministicNextStateReward
    value::Vector{Float64}
end
reward(r::DeterministicNextStateReward, s, a, s′) = r.value[s′]
"""
    struct DeterministicStateActionReward
        value::Array{Float64, 2}

`value` should be a `na × ns`-matrix.
"""
struct DeterministicStateActionReward
    value::Array{Float64, 2}
end
reward(r::DeterministicStateActionReward, s, a, s′) = r.value[a, s]
"""
    struct NormalNextStateReward
        mean::Vector{Float64}
        std::Vector{Float64}
"""
struct NormalNextStateReward
    mean::Vector{Float64}
    std::Vector{Float64}
end
reward(r::NormalNextStateReward, s, a, s′) = r.mean[s′] + randn(ENV_RNG) * r.std[s′]
"""
    struct NormalStateActionReward
        mean::Array{Float64, 2}
        std::Array{Float64, 2}

`mean` and `std` should be `na × ns`-matrices.
"""
struct NormalStateActionReward
    mean::Array{Float64, 2}
    std::Array{Float64, 2}
end
reward(r::NormalStateActionReward, s, a, s′) = r.mean[a, s] + randn(ENV_RNG) * r.std[a, s]

"""
    getprobvecrandom(n)

Returns an array of length `n` that sums to 1. More precisely, the array is a
sample of a [Dirichlet
distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution) with `n`
categories and ``α_1 = ⋯  = α_n = 1``.
"""
getprobvecrandom(n) = normalize(rand(ENV_RNG, n), 1)
"""
    getprobvecrandom(n, min, max)

Returns an array of length `n` that sums to 1 where all elements outside of
`min`:`max` are zero.
"""
getprobvecrandom(n, min, max) = sparsevec(collect(min:max),
                                          getprobvecrandom(max - min + 1), n)
"""
    getprobvecuniform(n)  = fill(1/n, n)
"""
getprobvecuniform(n) = fill(1/n, n)
"""
    getprobvecdeterministic(n, min = 1, max = n)

Returns a `SparseVector` of length `n` where one element in `min`:`max` has
value 1.
"""
getprobvecdeterministic(n, min = 1, max = n) = sparsevec([rand(ENV_RNG, min:max)], [1.], n)
# constructors
"""
    MDP(ns, na; init = "random")
    MDP(; ns = 10, na = 4, init = "random")

Return MDP with `init in ("random", "uniform", "deterministic")`, where the
keyword init determines how to construct the transition probabilites (see also
[`getprobvecrandom`](@ref), [`getprobvecuniform`](@ref),
[`getprobvecdeterministic`](@ref)).
"""
function MDP(ns, na; init = "random")
    r = randn(ENV_RNG, na, ns)
    func = eval(Meta.parse("getprobvec" * init))
    T = [func(ns) for a in 1:na, s in 1:ns]
    MDP(DiscreteSpace(ns, 1), DiscreteSpace(na, 1), rand(ENV_RNG, 1:ns), T, r,
        1:ns, zeros(ns))
end
MDP(; ns = 10, na = 4, init = "random") = MDP(ns, na, init = init)

actionspace(env::MDP) = env.actionspace

"""
    treeMDP(na, depth; init = "random", branchingfactor = 3)

Returns a tree structured MDP with na actions and `depth` of the tree.
If `init` is random, the `branchingfactor` determines how many possible states a
(action, state) pair has. If `init = "deterministic"` the `branchingfactor =
na`.
"""
function treeMDP(na, depth;
                 init = "random",
                 branchingfactor = 3)
    isdet = (init == "deterministic")
    if isdet
        branchingfactor = na
        ns = na.^(0:depth - 1)
    else
        ns = branchingfactor.^(0:depth - 1)
    end
    cns = cumsum(ns)
    func = eval(Meta.parse("getprobvec" * init))
    T = Array{SparseVector, 2}(undef, na, cns[end])
    for i in 1:depth - 1
        for s in 1:ns[i]
            for a in 1:na
                lb = cns[i] + (s - 1) * branchingfactor + (a - 1) * isdet + 1
                ub = isdet ? lb : lb + branchingfactor - 1
                T[a, (i == 1 ? 0 : cns[i-1]) + s] = func(cns[end] + 1, lb, ub)
            end
        end
    end
    r = zeros(na, cns[end] + 1)
    isterminal = [zeros(cns[end]); 1]
    for s in cns[end-1]+1:cns[end]
        for a in 1:na
            r[a, s] = -rand(ENV_RNG)
            T[a, s] = getprobvecdeterministic(cns[end] + 1, cns[end] + 1,
                                              cns[end] + 1)
        end
    end
    MDP(DiscreteSpace(cns[end] + 1, 1), DiscreteSpace(na, 1), 1, T, r, 1:1, isterminal)
end

function emptytransprob!(v::SparseVector)
    empty!(v.nzind); empty!(v.nzval)
end
emptytransprob!(v::Array{Float64, 1}) = v[:] .*= 0.

"""
    setterminalstates!(mdp, range)

Sets `mdp.isterminal[range] .= 1`, empties the table of transition probabilities
for terminal states and sets the reward for all actions in the terminal state to
the same value.
"""
function setterminalstates!(mdp, range)
    mdp.isterminal[range] .= 1
    for s in findall(mdp.isterminal)
        mdp.reward[:, s] .= mean(mdp.reward[:, s])
        for a in 1:mdp.na
            emptytransprob!(mdp.trans_probs[a, s])
        end
    end
end

# run MDP

"""
    run!(mdp::MDP, action::Int64)

Transition to a new state given `action`. Returns the new state.
"""
function run!(mdp::MDP, action::Int64)
    if mdp.isterminal[mdp.state] == 1
        reset!(mdp)
    else
        mdp.state = wsample(ENV_RNG, mdp.trans_probs[action, mdp.state])
        (observation = mdp.state,)
    end
end

"""
    run!(mdp::MDP, policy::Array{Int64, 1}) = run!(mdp, policy[mdp.state])

"""
run!(mdp::MDP, policy::Array{Int64, 1}) = run!(mdp, policy[mdp.state])


text:
function julia_main()::Cint
    try
        processargs()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function processargs()

    if length(ARGS) >= 3

        ignorelane = false
        indexfile = ""
        metafile = ""
        outfolder = "."
        index = 1
        readfile = ""
        while index <= length(ARGS)
            if ARGS[index] == "-ignorelane"
                ignorelane = true
            elseif ARGS[index] == "-indexfile"
                index += 1
                indexfile = ARGS[index]
            elseif ARGS[index] == "-meta"
                index += 1
                metafile = ARGS[index]
            elseif ARGS[index] == "-outfolder"
                index += 1
                outfolder = ARGS[index]
            else
                readfile = ARGS[index]
            end
            index += 1
        end
        rundemultiplex(readfile, indexfile, metafile, outfolder, ignorelane)
        # if index_file == ""
        #     split_fastq(file, meta, folder, ignore_lane)
        # else
        #     split_fastq_annot(file, index_file, meta, folder, ignore_lane)
        # end

    else

        println("Usage: DemultiplexFASTQ [-ignorelane] -meta <meta> -outfolder <outfolder> -indexfile <index_file> <fastqfile>")
        printboldln("\t-ignorelane\tSet if all lanes contain same pool of barcodes, note sample meta file must contain column *Lane* if *ignorelane* not set")
        printboldln("\t-meta\t\tSample meta file, a tab-delimited file containing columns: *SampleName*, *Index*, *[Lane]*")
        printboldln("\t-outfolder\tOutput folder")
        printboldln("\t-indexfile\tFASTQ file containing index read")
    end

end

function printboldln(s)
    f = split(s, "*")
    for i = 1:2:length(f)
        print(f[i])
        (i < length(f)) && printstyled(f[i+1], bold=true)
    end
    println("")
end

ext: jl


text:
using Clang.Generators
using snappy_jll

cd(@__DIR__)

include_dir = normpath(snappy_jll.artifact_dir, "include")

opts = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

headers = filter(file -> endswith(file, ".h"), readdir(include_dir, join=true))

ctx = create_context(headers, args, opts)

build!(ctx)


text:
"""
SamplerState for the stochastic gradient Nosé-Poincaré Hamiltonian Dynamcics algorithm following

Roychowdhury, A., Kulis, B., & Parthasarathy, S. (2016). Robust Monte Carlo Sampling using Riemannian Nos’{e}-Poincar’{e} Hamiltonian Dynamics. In ICML (pp. 2673–2681).

Fields:

    - x: the current state
    - p: current momentum
    - xi: thermostat variable
    - t: iterations counter
    - stepsize: s.stepsize(s.t) gives the stepsize to be used in iteration s.t
    - mu: Constant in the thermostat equation
    - D: friction parameters
"""
type SGNHTState <: SamplerState
    x::Array{Float64,1}
    p::Array{Float64,1}
    xi::Float64

    t::Int
    stepsize::Function
    mu::Float64
    D::Float64
    function SGNHTState(x,stepsize::Function;mu=1.,D=1.0)
        p = randn(length(x))
        xi = D
        new(x,p,xi,0,stepsize,mu,D)
    end
    function SGNHTState(x,stepsize::Float64;mu=1.,D=1.0)
        p = randn(length(x))
        xi = D
        eps = stepsize
        stepsize = niters::Int -> eps
        new(x,p,xi,0,stepsize,mu,D)
    end

end

"""
sample! performs one SGLD update and returns the updated SGLDState.
"""
function sample!(s::SGNHTState,grad::Function)
    nparams = length(s.x)
    eps = s.stepsize(s.t)
    Base.LinAlg.axpy!(.5*eps,grad(s.x),s.p)
    s.t += 1
    Base.LinAlg.axpy!(0.5*eps,s.p,s.x)
    s.xi += .5*eps/s.mu * (sum(s.p.*s.p) - nparams)
    if abs(eps*s.xi) > 1e-6
    s.p = exp(-eps*s.xi)*s.p + sqrt(.5*(1-exp(-2.0*eps*s.xi))/s.xi * s.D) .* randn(nparams)
    else
    Base.LinAlg.axpy!(sqrt(eps*s.D), randn(nparams),s.p)
    end
    s.xi += .5*eps/s.mu * (sum(s.p.*s.p) - nparams)
    Base.LinAlg.axpy!(.5*eps,s.p,s.x)
    Base.LinAlg.axpy!(.5*eps,grad(s.x),s.p)
    s
end


text:
using PyCall
using StatsBase

pushfirst!(PyVector(pyimport("sys")["path"]), "")
#ushfirst!(PyVector(pyimport("sys")["path"]), @__DIR__)
@info "torchmodels compiled. path added -> $(pwd())"

abstract type TorchModel end

"""
	fAnoGAN_GP(;idim=(1,2,2), zdim=100, kernelsizes=(9,7,5,3), channels=(128,64,32,16), scalings=(2,2,2,1),
		activation="relu", batchnorm=false, usegpu=true))

The fAnoGAN model with gradient penalization.
"""
mutable struct fAnoGAN_GP <: TorchModel
	model
end

function fAnoGAN_GP(;kwargs...)
pushfirst!(PyVector(pyimport("sys")["path"]), "")

py"""
import numpy
import torch
import fanogan_gp # <- this module must be located in directory where you want to run it from.

def c_fanogan_gp(kwargs):
	return fanogan_gp.fAnoGAN(**kwargs)
"""
	return fAnoGAN_GP(py"c_fanogan_gp"(kwargs))
end

function StatsBase.fit!(model::fAnoGAN_GP, X::Array{T, 4}; max_iters=10000, lr_gan=1e-4, lr_enc=1e-4,
		batch_size=64, n_critic=5) where T<:Real
	# transposition since Python models are row major
	X = Array(permutedims(X, [4,3,2,1]))
	_, history = model.model.fit(X, max_iters, lr_gan, lr_enc, batch_size, n_critic)
	return model, history
end

function StatsBase.predict(model::fAnoGAN_GP, X::Array{T, 4}) where T<:Real
	X = Array(permutedims(X, [4,3,2,1]))
	return Array(model.model.predict(X))
end


text:
using Revise
using ADSeismic
using ADCME
using MAT
using PyPlot
using PyCall
using DelimitedFiles
# matplotlib.use("Agg")
close("all")
if has_gpu()
  gpu = true
else
  gpu = false
end

data_dir = "data/acoustic"
if !ispath(data_dir)
  mkpath(data_dir)
end
figure_dir = "figure/PhysGNN_V2/"
if !ispath(figure_dir)
  mkpath(figure_dir)
end
result_dir = "result/PhysGNN_V2/"
if !ispath(result_dir)
  mkpath(result_dir)
end
# if isfile(joinpath(result_dir, "loss.txt"))
#   rm(joinpath(result_dir, "loss.txt"))
# end

################### Inversion using Automatic Differentiation #####################
reset_default_graph()
batch_size = 8 ## number of models per source
sample_size = 8 ## number of sampled y per source
model_name = "models/marmousi2-model-smooth.mat"

## load model
params = load_params(model_name)
src = load_acoustic_source(model_name)
rcv = load_acoustic_receiver(model_name)
nsrc = length(src)
rcv = Array{Any}(undef, nsrc, batch_size)
for i = 1:nsrc
  for j = 1:batch_size
    rcv[i,j] = load_acoustic_receiver(model_name)[i]
  end
end

vp0 = matread(model_name)["vp"]
mean_vp0 = mean(vp0)
std_vp0 = std(vp0)
vp0 = constant(vp0)

## load data
Rs = Array{Array{Float64,2}}(undef, length(src))
for i = 1:length(src)
    Rs[i] = readdlm(joinpath(data_dir, "marmousi-r$i.txt"))
end

## original fwi
# vp_ = Variable(vp0 / mean_vp0)
# vp = vp_ * mean_vp0
# vars = vp_

## add NN
x, isTrain, y, z = sampleUQ(batch_size, sample_size, (size(src[1].srcv,1)+1,length(rcv[1].rcvi)),
                            z_size=8, base=4, ratio=size(vp0)[1]/size(vp0)[2])
x = x * std_vp0
vp = vp0 + tf.slice(x, (size(x).-(size(x)[1], size(vp0)...)).÷2, (size(x)[1], size(vp0)...))

## assemble acoustic propagator model
# model = x->AcousticPropagatorSolver(params, x, vp^2)
models = Array{Any}(undef, nsrc, batch_size)
for i = 1:nsrc
  for j = 1:batch_size
    vp = vp0 + tf.slice(x[j], (size(x[j]).-size(vp0)).÷2, size(vp0))
    models[i,j] = x->AcousticPropagatorSolver(params, x, vp)
  end
end

@pywith tf.device("/cpu:0") begin
  global variable_rcv_ = placeholder(Float64, shape=(nsrc, size(Rs[1])...))
  global Rs_ = Array{PyObject}(undef, nsrc, sample_size)
  for i = 1:nsrc
    for j = 1:sample_size
        Rs_[i,j] = variable_rcv_[i] + y[i, j]
    end
  end

  global variable_src = Array{Any}(undef, nsrc)
  global si_ = placeholder(Int32, shape=[nsrc,length(src[1].srci)])
  global sj_ = placeholder(Int32, shape=[nsrc,length(src[1].srcj)])
  global sv_ = placeholder(Float64, shape=[nsrc,size(src[1].srcv)...])
  for i = 1:nsrc
    variable_src[i] = AcousticSource(si_[i], sj_[i], sv_[i])
  end

  global loss = sampling_compute_loss_and_grads_GPU_v2(models, variable_src, rcv, Rs_, reg=200.0, method="lp")
end

lr = 0.001
opt = AdamOptimizer(lr).minimize(loss, colocate_gradients_with_ops=true)

si = zeros(size(si_)...)
sj = zeros(size(sj_)...)
sv = zeros(size(sv_)...)
for i = 1:nsrc
  si[i,:] = src[i].srci
  sj[i,:] = src[i].srcj
  sv[i,:,:] = src[i].srcv
end
variable_rcv = zeros(size(variable_rcv_)...)
for i = 1:nsrc
  variable_rcv[i,:,:] = Rs[i]
end

dic = Dict(
  isTrain=>true,
  z=>rand(Float32, size(z)...),
  y=>randn(size(y)...),
  si_=>si,
  sj_=>sj,
  sv_=>sv,
  variable_rcv_=>variable_rcv
)

sess = Session(); init(sess)
@info "Initial loss: ", run(sess, loss, feed_dict=dic)

## run inversion
losses = []
σ = 0.01
fixed_z = randn(Float32, size(z)...)
time = 0
for iter = 1:1000
    if iter%10==1
      dic = Dict(
        isTrain=>true,
        z=>fixed_z
      )
      plot_result(sess, x, dic, iter, dirs=figure_dir, var_name="x")
      plot_result(sess, vp, dic, iter, dirs=figure_dir, var_name="vp")
    end

    dic = Dict(
      isTrain=>true,
      z=>rand(Float32, size(z)...),
      y=>randn(size(y)...) .* σ .* [std(Rs[i]) for i = 1:nsrc],
      si_=>si,
      sj_=>sj,
      sv_=>sv,
      variable_rcv_=>variable_rcv
    )

    global time += @elapsed ls, _ = run(sess, [loss, opt], feed_dict=dic)
    if iter == 1
      global mean_ls = ls
    else
      global mean_ls
      mean_ls += 1/(iter+1) * (ls - mean_ls)
    end
    # println("[#$iter] loss = $loss_, mean loss = $mean_ls")
    println("   $iter\t$ls\t$mean_ls")
    println(" * time: $time")
    push!(losses, ls)
end


text:
using Clipper
using Test

function test(run::Function, name::AbstractString; verbose=true)
    test_name = "Test $(name)"
    if verbose
        println(stderr, test_name)
    end
    try
        run()
    catch e
        println("Error running: ", test_name)
        rethrow(e)
    end
end

include("clipper_test.jl")
include("clipper_offset_test.jl")


text:
function worker_size()
    ws = VectorizationBase.nextpow2(num_threads())
    IfElse.ifelse(Static.lt(ws,StaticInt{8}()), StaticInt{8}(), ws)
end
worker_type() = VectorizationBase.mask_type(worker_size())
worker_pointer_type() = Ptr{worker_type()}

const WORKERS = Ref(zero(UInt128)) # 0 = unavailable, 1 = available
const STATES = UInt8[]

worker_pointer() = Base.unsafe_convert(worker_pointer_type(), pointer_from_objref(WORKERS))

function reserved(id)
    p = Base.unsafe_convert(worker_pointer_type(), STATES)
    __vload(p, VectorizationBase.lazymul(cache_linesize(), id), False(), register_size())
end
function reserve_threads!(id, reserve)
    p = Base.unsafe_convert(worker_pointer_type(), STATES)
    __vstore!(p, reserve % worker_type(), VectorizationBase.lazymul(cache_linesize(), id), False(), False(), False(), register_size())
    nothing
end
function free_threads!(freed_threads)
    ThreadingUtilities._atomic_or!(worker_pointer(), freed_threads)
    nothing
end
function free_local_threads!()
    tid = Base.Threads.threadid()
    tmask = one(worker_type()) << (tid - one(tid))
    r = reserved(tid) | tmask
    reserve_threads!(tid, zero(worker_type()))
    free_threads!(r)
end

@inline function _request_threads(id::UInt32, num_requested::UInt32)
    reserved_threads = reserved(id)
    reserved_count = count_ones(reserved_threads)%UInt32
    no_threads = zero(worker_type())
    # reserved_count ≥ num_requested && return reserved_threads, no_threads
    reserved_count%Int32 ≥ num_requested%Int32 && return UnsignedIteratorEarlyStop(reserved_threads, num_requested), no_threads
    # to get more, we xchng, setting all to `0`
    # then see which we need, and free those we aren't using.
    wp = worker_pointer()
    _all_threads = all_threads = ThreadingUtilities._atomic_xchg!(wp, no_threads)
    additional_threads = count_ones(all_threads)
    # num_requested === StaticInt{-1}() && return reserved_threads, all_threads
    num_requested === StaticInt{-1}() && return UnsignedIteratorEarlyStop(reserved_threads | all_threads), all_threads
    excess = additional_threads + reserved_count - num_requested
    # signed(excess) ≤ 0 && return reserved_threads, all_threads
    signed(excess) ≤ 0 && return UnsignedIteratorEarlyStop(reserved_threads | all_threads), all_threads
    # we need to return the `excess` to the pool.
    lz = leading_zeros(all_threads) % UInt32
    # i = 8
    while true
        # start by trying to trim off excess from lz
        lz += excess%UInt32
        m = (one(worker_type()) << (UInt32(worker_size()) - lz)) - one(worker_type())
        masked = (all_threads & m) ⊻ all_threads
        excess -= count_ones(masked)
        all_threads &= (~masked)
        # @show bitstring(masked), count_ones(masked), bitstring(unused_threads), excess, lz, bitstring(all_threads)
        excess == 0 && break
        # i -= 1
        # @assert i > 0
    end
    ThreadingUtilities._atomic_store!(wp, _all_threads & (~all_threads))
    return UnsignedIteratorEarlyStop(reserved_threads | all_threads, num_requested), all_threads
end
@inline function request_threads(id, num_requested)
    _request_threads(id % UInt32, num_requested % UInt32)
end
reserved_threads(id) = UnsignedIteratorEarlyStop(reserved(id))
reserved_threads(id, count) = UnsignedIteratorEarlyStop(reserved(id), count % UInt32)


text:
using SymbolicUtils: istree

function nterms(t)
    if istree(t)
        return reduce(+, map(nterms, arguments(t)), init=0)
    else
        return 1
    end
end
# Soft pivoted
# Note: we call this function with a matrix of Union{SymbolicUtils.Symbolic, Any}
# It should work as-is with Operation type too.
function sym_lu(A)
    m, n = size(A)
    L = fill!(Array{Any}(undef, size(A)),0) # TODO: make sparse?
    for i=1:min(m, n)
        L[i,i] = 1
    end
    U = copy!(Array{Any}(undef, size(A)),A)
    p = BlasInt[1:m;]
    for k = 1:m-1
        _, i = findmin(map(ii->iszero(U[ii, k]) ? Inf : nterms(U[ii,k]), k:n))
        i += k - 1
        # swap
        U[k, k:end], U[i, k:end] = U[i, k:end], U[k, k:end]
        L[k, 1:k-1], L[i, 1:k-1] = L[i, 1:k-1], L[k, 1:k-1]
        p[k] = i

        for j = k+1:m
            L[j,k] = U[j, k] / U[k, k]
            U[j,k:m] .= U[j,k:m] .- L[j,k] .* U[k,k:m]
        end
    end
    for j=1:m
        for i=j+1:n
            U[i,j] = 0
        end
    end

    (L, U, LinearAlgebra.ipiv2perm(p, m))
end

# Given a vector of equations and a
# list of independent variables,
# return the coefficient matrix `A` and a
# vector of constants (possibly symbolic) `b` such that
# A \ b will solve the equations for the vars
function A_b(eqs, vars)
    exprs = rhss(eqs) .- lhss(eqs)
    for ex in exprs
        @assert islinear(ex, vars)
    end
    A = jacobian(exprs, vars)
    b = A * vars - exprs
    A, b
end

"""
    solve_for(eqs::Vector, vars::Vector)

Solve the vector of equations `eqs` for a set of variables `vars`.

Assumes `length(eqs) == length(vars)`

Currently only works if all equations are linear.
"""
function solve_for(eqs, vars)
    A, b = A_b(eqs, vars)
    _solve(A, b)
end

function _solve(A, b)
    A = SymbolicUtils.simplify.(to_symbolic.(A), polynorm=true)
    b = SymbolicUtils.simplify.(to_symbolic.(b), polynorm=true)
    SymbolicUtils.simplify.(ldiv(sym_lu(A), b))
end

# ldiv below

_iszero(x::Number) = iszero(x)
_isone(x::Number) = isone(x)
_iszero(::Term) = false
_isone(::Term) = false

function simplifying_dot(x,y)
    isempty(x) && return 0
    muls = map(x,y) do xi,yi
        _isone(xi) ? yi : _isone(yi) ? xi : _iszero(xi) ? 0 : _iszero(yi) ? 0 : xi * yi
    end

    reduce(muls) do acc, x
        _iszero(acc) ? x : _iszero(x) ? acc : acc + x
    end
end

function ldiv((L,U,p), b)
    m, n = size(L)
    x = Vector{Any}(undef, length(b))
    b = b[p]

    for i=n:-1:1
        sub = simplifying_dot(x[i+1:end], U[i,i+1:end])
        den = U[i,i]
        x[i] = _iszero(sub) ? b[i] : b[i] - sub
        x[i] = _isone(den) ? x[i] : _isone(-den) ? -x[i] : x[i] / den
    end

    # unit lower triangular solve first:
    for i=1:n
        sub = simplifying_dot(b[1:i-1], L[i, 1:i-1]) # this should be `b` not x
        x[i] = _iszero(sub) ? x[i] : x[i] - sub
    end
    x
end


text:
# Autogenerated wrapper script for libfdk_aac_jll for aarch64-apple-darwin
export libfdk

JLLWrappers.@generate_wrapper_header("libfdk_aac")
JLLWrappers.@declare_library_product(libfdk, "@rpath/libfdk-aac.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfdk,
        "lib/libfdk-aac.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "ExtremelySimpleNL" begin
# input
s1 = lorenzx + 0.1noise # input
s2 = roesslerz + 0.1noise # input
k = 5 # backward dim
ℓ = 5 # forward dim
ε = 0.3 # range
w = 1 # theiler window

for (name, s) in zip(("lorenz", "roessler"), (s1, s2))

    τ = estimate_delay(s, "mi_min") #  5 # delaytime
    method = ExtremelySimpleNL(k, ℓ, τ, w, ε)

    x, r = decompose(s, method)

    err = nrmse(x, s)
    @test err < 0.2

    # figure()
    # title(name)
    # plot(s; alpha = 0.75, label = "original")
    # plot(x; alpha = 0.75, ls = "dashed", label = "noiseless, err=$err")
    # legend()
end
end

@testset "ManifoldProjection" begin
    @testset "Henon" begin
        he = Systems.henon()
        tr = trajectory(he, 10000; Ttr = 100)
        Random.seed!(151521)
        z = tr[:, 1]
        s = z .+ randn(10001)*0.1*std(z)
        # s = lorenzx + 0.1noise
        m = 5
        w = 1 # theiler window
        metric = Euclidean()
        searchtype = 30
        Q = 2
        Q = [2, 2, 2, 3, 3, 3, 3]

        x, r = decompose(s, ManifoldProjection(m, Q, searchtype))

        err = nrmse(x, tr[:, 1])
        @test err < 0.1

        # figure()
        # plot(s, label = "input")
        # plot(z, color = "k", ls = "dashed", label = "real")
        # plot(x, alpha = 0.5, label = "output")

        # figure()
        # scatter(s[1:end-1], s[2:end], s = 0.1, label = "input")
        # scatter(z[1:end-1], z[2:end], s = 0.1, label = "real", color = "k")
        # scatter(x[1:end-1], x[2:end], s = 0.1, label = "output", alpha = 0.5)
        # legend()
    end

    @testset "lorenz" begin
        lo = Systems.lorenz()
        tr = trajectory(lo, 1000; Ttr = 100, dt = 0.1)
        Random.seed!(151521)
        z = tr[:, 1]
        s = z .+ randn(10001)*0.2*std(z)
        m = 13
        w = 1 # theiler window
        metric = Euclidean()
        k = 30
        Q = 2
        Q = [2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]

        x, r = decompose(s, ManifoldProjection(m, Q, 30))
        err = nrmse(x, z)
        @test err < 0.2

        # figure()
        # plot(s, label = "input")
        # plot(z, color = "k", ls = "dashed", label = "real")
        # plot(x, alpha = 0.5, label = "output")
        # estimate_delay(z, "mi_min")
        # xlim(0, 200)
        #
        # τ = 1
        # d = 2
        # figure()
        # plot(columns(embed(s, d, τ))..., marker = "o", ms = 1, lw = 0.1, label = "input")
        # plot(columns(embed(z, d, τ))..., marker = "o", ms = 1, lw = 0.1, label = "real", color = "k")
        # plot(columns(embed(x, d, τ))..., marker = "o", ms = 1, lw = 0.1, label = "real", alpha = 0.5)
    end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

using ExecutableSpecifications.Gherkin: Given, When, Then
using AttributeType
using TypeDBClientException
using ThingSteps: get, put

###### When Steps ######################
#                                      #
########################################

@when("\$x = attribute(is-alive) as(boolean) put: true") do context
    @fail "Implement me"
end

#maybe like this? :D
@when("{var:Var} = attribute({type_label}) as(boolean) put: {value:Bool}") do context
    context[Symbol(var)] = [:tx(), :concepts(), get_attribute_type(:type_label::String), :as_remote(context[:tx()]), :as_boolean(), :put(value::bool)]
end

#=
@step("{var:Var} = attribute({type_label}) as(boolean) put: {value:Bool}")
def step_impl(context: Context, var: str, type_label: str, value: bool):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_boolean().put(value))
=#

@when("\$x = attribute(age) as(long) put: 21") do context
    @fail "Implement me"
end

@when("\$x = attribute(score) as(double) put: 123.456") do context
    @fail "Implement me"
end

@when("\$x = attribute(name) as(string) put: alice") do context
    @fail "Implement me"
end

@when("\$x = attribute(email) as(string) put: alice@email.com") do context
    @fail "Implement me"
end

@when("attribute(email) as(string) put: alice-email-com; throws exception") do context
    @fail "Implement me"
end

@when("\$x = attribute(birth-date) as(datetime) put: 1990-01-01 11:22:33") do context
    @fail "Implement me"
end

@when("\$x = attribute(score) as(double) put: 123") do context
    @fail "Implement me"
end

@when("delete attribute: \$x") do context
    @fail "Implement me"
end

###### Then Steps ######################
#                                      #
########################################

@then("attribute \$x is null: false") do context
    @fail "Implement me"
end

@then("attribute \$x has type: is-alive") do context
    @fail "Implement me"
end

@then("attribute \$x has value type: boolean") do context
    @fail "Implement me"
end

@then("attribute \$x has boolean value: true") do context
    @fail "Implement me"
end

@then("attribute \$x has type: age") do context
    @fail "Implement me"
end

@then("attribute \$x has value type: long") do context
    @fail "Implement me"
end

@then("attribute \$x has long value: 21") do context
    @fail "Implement me"
end

@then("attribute \$x has type: score") do context
    @fail "Implement me"
end

@then("attribute \$x has value type: double") do context
    @fail "Implement me"
end

@then("attribute \$x has double value: 123.456") do context
    @fail "Implement me"
end

@then("attribute \$x has type: name") do context
    @fail "Implement me"
end

@then("attribute \$x has value type: string") do context
    @fail "Implement me"
end

@then("attribute \$x has string value: alice") do context
    @fail "Implement me"
end

@then("attribute \$x has type: email") do context
    @fail "Implement me"
end

@then("attribute \$x has string value: alice@email.com") do context
    @fail "Implement me"
end

@then("attribute \$x has type: birth-date") do context
    @fail "Implement me"
end

@then("attribute \$x has value type: datetime") do context
    @fail "Implement me"
end

@then("attribute \$x has datetime value: 1990-01-01 11:22:33") do context
    @fail "Implement me"
end

@then("attribute(is-alive) get instances contain: \$x") do context
    @fail "Implement me"
end

@then("attribute(age) get instances contain: \$x") do context
    @fail "Implement me"
end

@then("attribute(score) get instances contain: \$x") do context
    @fail "Implement me"
end

@then("attribute(name) get instances contain: \$x") do context
    @fail "Implement me"
end

@then("attribute(birth-date) get instances contain: \$x") do context
    @fail "Implement me"
end

@then("attribute \$x is deleted: true") do context
    @fail "Implement me"
end

@then("attribute \$x is null: true") do context
    @fail "Implement me"
end

@then("attribute \$x has double value: 123") do context
    @fail "Implement me"
end

#=
from datetime import datetime

from behave import *
from hamcrest import *

from typedb.api.concept.type.attribute_type import AttributeType
from typedb.common.exception import TypeDBClientException
from tests.behaviour.context import Context

@step("attribute({type_label}) get instances contain: {var:Var}")
def step_impl(context: Context, type_label: str, var: str):
    assert_that(context.get(var), is_in(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).get_instances()))

@step("attribute({type_label}) get instances is empty")
def step_impl(context: Context, type_label: str):
    assert_that(calling(next).with_args(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).get_instances()), raises(StopIteration))

@step("attribute {var1:Var} get owners contain: {var2:Var}")
def step_impl(context: Context, var1: str, var2: str):
    assert_that(context.get(var2), is_in(context.get(var1).as_remote(context.tx()).get_owners()))

@step("attribute {var1:Var} get owners do not contain: {var2:Var}")
def step_impl(context: Context, var1: str, var2: str):
    assert_that(context.get(var2), not_(is_in(context.get(var1).as_remote(context.tx()).get_owners())))

@step("attribute {var:Var} has value type: {value_type:ValueType}")
def step_impl(context: Context, var: str, value_type: AttributeType.ValueType):
    assert_that(context.get(var).get_type().get_value_type(), is_(value_type))

@step("attribute({type_label}) as(boolean) put: {value:Bool}; throws exception")
def step_impl(context: Context, type_label: str, value: bool):
    assert_that(calling(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_boolean().put).with_args(value), raises(TypeDBClientException))

@step("{var:Var} = attribute({type_label}) as(boolean) put: {value:Bool}")
def step_impl(context: Context, var: str, type_label: str, value: bool):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_boolean().put(value))

@step("attribute({type_label}) as(long) put: {value:Int}; throws exception")
def step_impl(context: Context, type_label: str, value: int):
    assert_that(calling(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_long().put).with_args(value), raises(TypeDBClientException))

@step("{var:Var} = attribute({type_label}) as(long) put: {value:Int}")
def step_impl(context: Context, var: str, type_label: str, value: int):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_long().put(value))

@step("attribute({type_label}) as(double) put: {value:Float}; throws exception")
def step_impl(context: Context, type_label: str, value: float):
    assert_that(calling(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_double().put).with_args(value), raises(TypeDBClientException))

@step("{var:Var} = attribute({type_label}) as(double) put: {value:Float}")
def step_impl(context: Context, var: str, type_label: str, value: float):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_double().put(value))

@step("attribute({type_label}) as(string) put: {value}; throws exception")
def step_impl(context: Context, type_label: str, value: str):
    assert_that(calling(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_string().put).with_args(value), raises(TypeDBClientException))

@step("{var:Var} = attribute({type_label}) as(string) put: {value}")
def step_impl(context: Context, var: str, type_label: str, value: str):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_string().put(value))

@step("attribute({type_label}) as(datetime) put: {value:DateTime}; throws exception")
def step_impl(context: Context, type_label: str, value: datetime):
    assert_that(calling(context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_datetime().put).with_args(value), raises(TypeDBClientException))

@step("{var:Var} = attribute({type_label}) as(datetime) put: {value:DateTime}")
def step_impl(context: Context, var: str, type_label: str, value: datetime):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_datetime().put(value))

@step("{var:Var} = attribute({type_label}) as(boolean) get: {value:Bool}")
def step_impl(context: Context, var: str, type_label: str, value: bool):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_boolean().get(value))

@step("{var:Var} = attribute({type_label}) as(long) get: {value:Int}")
def step_impl(context: Context, var: str, type_label: str, value: int):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_long().get(value))

@step("{var:Var} = attribute({type_label}) as(double) get: {value:Float}")
def step_impl(context: Context, var: str, type_label: str, value: float):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_double().get(value))

@step("{var:Var} = attribute({type_label}) as(string) get: {value}")
def step_impl(context: Context, var: str, type_label: str, value: str):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_string().get(value))

@step("{var:Var} = attribute({type_label}) as(datetime) get: {value:DateTime}")
def step_impl(context: Context, var: str, type_label: str, value: datetime):
    context.put(var, context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).as_datetime().get(value))

@step("attribute {var:Var} has boolean value: {value:Bool}")
def step_impl(context: Context, var: str, value: bool):
    assert_that(context.get(var).get_value(), is_(value))

@step("attribute {var:Var} has long value: {value:Int}")
def step_impl(context: Context, var: str, value: int):
    assert_that(context.get(var).get_value(), is_(value))

@step("attribute {var:Var} has double value: {value:Float}")
def step_impl(context: Context, var: str, value: float):
    assert_that(context.get(var).get_value(), is_(value))

@step("attribute {var:Var} has string value: {value}")
def step_impl(context: Context, var: str, value: str):
    assert_that(context.get(var).get_value(), is_(value))

@step("attribute {var:Var} has datetime value: {value:DateTime}")
def step_impl(context: Context, var: str, value: datetime):
    assert_that(context.get(var).get_value(), is_(value))
=#


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

using Markdown
using InteractiveUtils

# ╔═╡ d9ed9070-64c2-11eb-316d-9f82842705be
using PlutoUI

# ╔═╡ b40446e0-64a6-11eb-3652-85fdb1612e5e
md"""# The Way of the Program

The goal of this book is to teach you to think like a computer scientist. This way of thinking combines some of the best features of mathematics, engineering, and natural science. Like mathematicians, computer scientists use formal languages to denote ideas (specifically computations). Like engineers, they design things, assembling components into systems and evaluating trade-offs among alternatives. Like scientists, they observe the behavior of complex systems, form hypotheses, and test predictions.

The single most important skill for a computer scientist is problem solving. Problem solving means the ability to formulate problems, think creatively about solutions, and express a solution clearly and accurately. As it turns out, the process of learning to program is an excellent opportunity to practice problem-solving skills. That’s why this chapter is called “The Way of the Program.”

On one level, you will be learning to program, a useful skill by itself. On another level, you will use programming as a means to an end. As we go along, that end will become clearer."""

# ╔═╡ 6b1e8390-64a7-11eb-2249-adcf7f79fc61
md"""## What Is a Program?

A _program_ is a sequence of instructions that specifies how to perform a computation. The computation might be something mathematical, such as solving a system of equations or finding the roots of a polynomial, but it can also be a symbolic computation, such as searching for and replacing text in a document, or something graphical, like processing an image or playing a video.

The details look different in different languages, but a few basic instructions appear in just about every language:

* *Input*: Get data from the keyboard, a file, the network, or some other device.
* *Output*: Display data on the screen, save it in a file, send it over the network, etc.
* *Math*: Perform basic mathematical operations like addition and multiplication.
* *Conditional execution*: Check for certain conditions and run the appropriate code.
* *Repetition*: Perform some action repeatedly, usually with some variation.

Believe it or not, that’s pretty much all there is to it. Every program you’ve ever used, no matter how complicated, is made up of instructions that look pretty much like these. So you can think of programming as the process of breaking a large, complex task into smaller and smaller subtasks until the subtasks are simple enough to be performed with one of these basic instructions."""

# ╔═╡ b2560330-64bd-11eb-0a30-57ada26f6ad1
md"""## Running Julia

One of the challenges of getting started with Julia is that you might have to install it and related software on your computer. Fortunately, the installation of Julia is straightforward. Julia is available for various operating systems at the [Julia](https://julialang.org/downloads/) website.

We will start out running Julia in the terminal and in a browser. Later, when you are comfortable with Julia, I’ll make suggestions for installing a dedicated programming environment on your computer.

The Julia *REPL* (Read–Eval–Print Loop) is a program that reads and executes Julia code. You can start the REPL by opening a terminal on your computer and typing **`julia`** on the command line. When it starts, you should see output like this:

```
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.5.3 (2020-11-09)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia>
```

The first lines contain information about the REPL, so it might be different for you. But you should check that the version number is at least 1.2.0.

The last line is a *prompt* that indicates that the REPL is ready for you to enter code. If you type a line of code and hit Enter, the REPL displays the result:

```julia
julia> 1 + 1
2
```

Code snippets can be copied and pasted verbatim, including the `julia>` prompt and any output.

Now you’re ready to get started. From here on, I assume that you know how to start the Julia REPL and run code.
"""

# ╔═╡ cdb192f0-64bf-11eb-0a21-2987e4b81296
md"""## The First Program
Traditionally, the first program you write in a new language is called “Hello, World!” because all it does is display the words “Hello, World!” In Julia, it looks like this:

```julia
julia> println("Hello, World!")
Hello, World!
```

This is an example of a *print statement*, although it doesn’t actually print anything on paper. It displays a result on the screen.

The quotation marks in the program mark the beginning and end of the text to be displayed; they don’t appear in the result.

The parentheses indicate that `println` is a function. We’ll get to functions in Chapter 3.
"""

# ╔═╡ 224a71b0-64e3-11eb-0f83-4d671426a415
md"""
!!! languages
    Other programming languages have a similar first program.

    In the C programming language, the following lines have to be put in a file and compiled before the program can run:
    ```c
    #include <stdio.h>

    int main(void) {
        printf("hello, world\n");
    }
    ```

    Python provides a REPL:
    ```python
    >>> print('Hello, World!')
    ```
"""

# ╔═╡ 5a40a0b0-64c2-11eb-00e4-6fb54d1a31bb
md"""## Arithmetic Operators

After “Hello, World!” the next step is arithmetic. Julia provides *operators*, which are symbols that represent computations like addition and multiplication.

The operators `+`, `-`, and `*` perform addition, subtraction, and multiplication, as in the following examples:

```julia
julia> 40 + 2
42
julia> 43 - 1
42
julia> 6 * 7
42
```

The operator `/` performs division:

```julia
julia> 84 / 2
42.0
```

You might wonder why the result is `42.0` instead of `42`. I’ll explain in the next section.

Finally, the operator `^` performs exponentiation; that is, it raises a number to a power:

```julia
julia> 6^2 + 6
42
```
"""

# ╔═╡ 97df5b90-64c3-11eb-1b21-31020fae745d
md"""## Values and Types

A `value` is one of the basic things a program works with, like a letter or a number. Some values we have seen so far are `2`, `42.0`, and `"Hello, World!"`.

These values belong to different types: `2` is an *integer*, `42.0` is a *floating-point number*, and `"Hello, World!"` is a *string*, so called because the letters it contains are strung together.

If you are not sure what type a value has, the REPL can tell you:

```julia
julia> typeof(2)
Int64
julia> typeof(42.0)
Float64
julia> typeof("Hello, World!")
String
```

Integers belong to the type `Int64`, strings belong to `String`, and floating-point numbers belong to `Float64`.

What about values like `"2"` and `"42.0"`? They look like numbers, but they are in quotation marks like strings. These are strings too:

```juliz
julia> typeof("2")
String
julia> typeof("42.0")
String
```

When you type a large integer, you might be tempted to use commas between groups of digits, as in `1,000,000`. This is not a legal integer in Julia, but it is legal:

```julia
julia> 1,000,000
(1, 0, 0)
```

That’s not what we expected at all! Julia parses `1,000,000` as a comma-separated sequence of integers. We’ll learn more about this kind of sequence later.

You can get the expected result using `1_000_000`, however.
"""

# ╔═╡ 437336b0-64e3-11eb-39de-3bec9939b989
md"""
!!! languages
    Python has a similar feature:
    ```python
    >>> type(2)
    <class 'int'>
    >>> type(42.0)
    <class 'float'>
    >>> type('Hello, World!')
    <class 'str'>
    ```
"""

# ╔═╡ 92130620-64c4-11eb-1790-df61511ecf87
md"""## Formal and Natural Languages

*Natural languages* are the languages people speak, such as English, Spanish, and French. They were not designed by people (although people try to impose some order on them); they evolved naturally.

*Formal languages* are languages that are designed by people for specific applications. For example, the notation that mathematicians use is a formal language that is particularly good at denoting relationships among numbers and symbols. Chemists use a formal language to represent the chemical structure of molecules. And most importantly, programming languages are formal languages that have been designed to express computations.

Formal languages tend to have strict *syntax* rules that govern the structure of statements. For example, in mathematics the statement ``3 + 3 = 6`` has correct syntax, but ``3 + = 3\$6`` does not. In chemistry, ``H_2O`` is a syntactically correct formula, but ``_2Zz`` is not.

Syntax rules come in two flavors: tokens and structure. Tokens are the basic elements of the language, such as words, numbers, and chemical elements. One of the problems with ``3 + = 3\$6`` is that ``\$`` is not a legal token in mathematics (at least as far as I know). Similarly, ``_2Zz`` is not legal because there is no element with the abbreviation ``Zz``.

The second type of syntax rule pertains to the way tokens are combined. The equation ``3 + = 3`` is illegal because even though ``+`` and ``=`` are legal tokens, you can’t have one right after the other. Similarly, in a chemical formula the subscript comes after the element name, not before.

This is ``@`` well-structured ``Engli\$h`` sentence with invalid ``t*kens`` in it. This sentence all valid tokens has, but invalid structure with.

When you read a sentence in English or a statement in a formal language, you have to figure out the structure (although in a natural language you do this subconsciously). This process is called *parsing*.

Although formal and natural languages have many features in common—tokens, structure, and syntax—there are some differences:

* *Ambiguity*:
  Natural languages are full of ambiguity, which people deal with by using contex! tual clues and other information. Formal languages are designed to be nearly or completely unambiguous, which means that any statement has exactly one meaning, regardless of context.
* *Redundancy*:
  In order to make up for ambiguity and reduce misunderstandings, natural languages employ lots of redundancy. As a result, they are often verbose. Formal languages are less redundant and more concise.
* *Literalness*:
  Natural languages are full of idiom and metaphor. If I say, “The penny dropped,” there is probably no penny and nothing dropping (this idiom means that someone understood something after a period of confusion). Formal languages mean exactly what they say.

Because we all grow up speaking natural languages, it is sometimes hard to adjust to formal languages. The difference between formal and natural language is like the difference between poetry and prose, but more so:

* *Poetry*:
  Words are used for their sounds as well as for their meaning, and the whole poem together creates an effect or emotional response. Ambiguity is not only common but often deliberate.
* *Prose*:
  The literal meaning of words is more important, and the structure contributes more meaning. Prose is more amenable to analysis than poetry but still often ambiguous.
* *Programs*:
  The meaning of a computer program is unambiguous and literal, and can be understood entirely by analysis of the tokens and structure.

Formal languages are more dense than natural languages, so it takes longer to read them. Also, the structure is important, so it is not always best to read from top to bottom, left to right. Instead, you’ll learn to parse the program in your head, identifying the tokens and interpreting the structure. Finally, the details matter. Small errors in spelling and punctuation, which you can get away with in natural languages, can make a big difference in a formal language.
"""

# ╔═╡ 75792d30-64c6-11eb-2b86-b13f35f42768
md"""## Debugging

Programmers make mistakes. For whimsical reasons, programming errors are called *bugs* and the process of tracking them down is called *debugging*.

Programming, and especially debugging, sometimes brings out strong emotions. If you are struggling with a difficult bug, you might feel angry, despondent, or embarrassed.

There is evidence that people naturally respond to computers as if they were people. When they work well, we think of them as teammates, and when they are obstinate or rude, we respond to them the same way we respond to rude, obstinate people.[^1]

Preparing for these reactions might help you deal with them. One approach is to think of the computer as an employee with certain strengths, like speed and precision, and particular weaknesses, like lack of empathy and inability to grasp the big picture.

Your job is to be a good manager: find ways to take advantage of the strengths and mitigate the weaknesses. And find ways to use your emotions to engage with the problem, without letting your reactions interfere with your ability to work effectively.

Learning to debug can be frustrating, but it is a valuable skill that is useful for many activities beyond programming. At the end of each chapter there is a section, like this one, with my suggestions for debugging. I hope they help!

[^1]: Reeves, Byron, and Clifford Ivar Nass. 1996. “The Media Equation: How People Treat Computers, Television, and New Media Like Real People and Places.” Chicago, IL: Center for the Study of Language and Information; New York: Cambridge University Press.
"""

# ╔═╡ a7270080-64c8-11eb-175a-192fbae1322b
md"""## Glossary

*problem solving*:
The process of formulating a problem, finding a solution, and expressing it.

*program*:
A sequence of instructions that specifies a computation.

*REPL*:
A program that repeatedly reads input, executes it, and outputs results.

*prompt*:
Characters displayed by the REPL to indicate that it is ready to take input from the user.

*print statement*:
An instruction that causes the Julia REPL to display a value on the screen.

*operator*:
A symbol that represents a simple computation like addition, multiplication, or string concatenation.

*value*:
A basic unit of data, like a number or string, that a program manipulates.

*type*:
A category of values. The types we have seen so far are integers (Int64), floatingpoint numbers (Float64), and strings (String).

*integer*:
A type that represents whole numbers.

*floating-point*:
A type that represents numbers with a decimal point.

*string*:
A type that represents sequences of characters.

*natural language*:
Any one of the languages that people speak that evolved naturally.

*formal language*:
Any one of the languages that people have designed for specific purposes, such as representing mathematical ideas or computer programs. All programming lan guages are formal languages.

*syntax*:
The rules that govern the structure of a program.

*token*:
One of the basic elements of the syntactic structure of a program, analogous to a word in a natural language.

*structure*:
The way tokens are combined.

*parse*:
To examine a program and analyze the syntactic structure.

*bug*:
An error in a program.

*debugging*:
The process of finding and correcting bugs.
"""

# ╔═╡ ce5393c0-64c9-11eb-25ec-c55227798aa8
md"""## Exercises

!!! tip
    It is a good idea to read this book in front of a computer so you can try out the examples as you go.

#### Exercise 1-1

Whenever you are experimenting with a new feature, you should try to make mistakes. For example, in the “Hello, World!” program, what happens if you leave out one of the quotation marks? What if you leave out both? What if you spell println wrong?

This kind of experiment helps you remember what you read; it also helps when you are programming, because you get to know what the error messages mean. It is better to make mistakes now and on purpose rather than later and accidentally.

1. In a print statement, what happens if you leave out one of the parentheses, or both?
2. If you are trying to print a string, what happens if you leave out one of the quotation marks, or both?
3. You can use a minus sign to make a negative number like ``-2``. What happens if you put a plus sign before a number? What about ``2++2``?
4. In math notation, leading zeros are okay, as in 02. What happens if you try this in Julia?
5. What happens if you have two values with no operator between them?
"""

# ╔═╡ 6472f940-64ca-11eb-0492-0953cf692c54
md"""#### Exercise 1-2

Start the Julia REPL and use it as a calculator.

1. How many seconds are there in 42 minutes 42 seconds?
2. How many miles are there in 10 kilometers? Note that there are 1.61 kilometers in a mile.
3. If you run a 10-kilometer race in 37 minutes 48 seconds, what is your average pace (time per mile in minutes and seconds)? What is your average speed in miles per hour?
"""

# ╔═╡ Cell order:
# ╟─d9ed9070-64c2-11eb-316d-9f82842705be
# ╟─b40446e0-64a6-11eb-3652-85fdb1612e5e
# ╟─6b1e8390-64a7-11eb-2249-adcf7f79fc61
# ╟─b2560330-64bd-11eb-0a30-57ada26f6ad1
# ╟─cdb192f0-64bf-11eb-0a21-2987e4b81296
# ╟─224a71b0-64e3-11eb-0f83-4d671426a415
# ╟─5a40a0b0-64c2-11eb-00e4-6fb54d1a31bb
# ╟─97df5b90-64c3-11eb-1b21-31020fae745d
# ╟─437336b0-64e3-11eb-39de-3bec9939b989
# ╟─92130620-64c4-11eb-1790-df61511ecf87
# ╟─75792d30-64c6-11eb-2b86-b13f35f42768
# ╟─a7270080-64c8-11eb-175a-192fbae1322b
# ╟─ce5393c0-64c9-11eb-25ec-c55227798aa8
# ╟─6472f940-64ca-11eb-0492-0953cf692c54


text:
# wrappers for functionality provided by the CUDA toolkit

# extensions to the C language
include(joinpath("cuda", "memory_shared.jl"))
include(joinpath("cuda", "indexing.jl"))
include(joinpath("cuda", "synchronization.jl"))
include(joinpath("cuda", "warp_vote.jl"))
include(joinpath("cuda", "warp_shuffle.jl"))
include(joinpath("cuda", "output.jl"))
include(joinpath("cuda", "assertion.jl"))
include(joinpath("cuda", "memory_dynamic.jl"))
include(joinpath("cuda", "atomics.jl"))
include(joinpath("cuda", "misc.jl"))

# functionality from libdevice
#
# > The libdevice library is a collection of NVVM bitcode functions that implement common
# > functions for NVIDIA GPU devices, including math primitives and bit-manipulation
# > functions. These functions are optimized for particular GPU architectures, and are
# > intended to be linked with an NVVM IR module during compilation to PTX.
include(joinpath("cuda", "math.jl"))
# TODO: native mathematical functions, CUDA C programming guide" > "C language extensions"
#       https://docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__INTRINSIC__DOUBLE.html
#       see /path/to/cuda/include/sm_20_intrinsics.h

# functionality from libcudadevrt
#
# The libcudadevrt library is a collection of PTX bitcode functions that implement parts of
# the CUDA API for execution on the device, such as device synchronization primitives,
# dynamic kernel APIs, etc.
const cudaError_t = Cint
include(joinpath("cuda", "cooperative_groups.jl"))
include(joinpath("cuda", "dynamic_parallelism.jl"))


text:
import Base: getindex, setindex!
import Base.Broadcast: materialize!
using Base.Broadcast: Broadcasted
using CUDA: CuArray, CuPtr, cuMemcpyDtoD_v2
using Knet.LibKnet8: @knet8
# include("kptr.jl") ## KnetPtr, Cptr
# include("karray.jl") ## KnetArray, KnetMatrix, KnetVector

### INDEXING
## Indexing with Real
## Indexing with Tuple{Real}
## Indexing with CartesianIndex: calls Tuple{Real}
## Indexing with AbstractUnitRange
## Indexing with Colon
## Indexing with KnetArray{Int32}: low level, only Int32 supported, no bound checking
## Indexing with (Colon,KnetArray{Int32})
## Indexing with (KnetArray{Int32},Colon)
## Indexing with AbstractArray{Real} calls KnetArray{Int32} after boundchecking
## Indexing with AbstractArray{CartesianIndex} calls AbstractArray{Real}
## Indexing with Empty Array or other unrecognized AbstractArray calls AbstractArray{Real}
## Indexing with (Colon,AbstractVector{Real}) calls (Colon,KnetArray{Int32}) after bound checking
## Indexing with (AbstractVector{Real},Colon) calls (KnetArray{Int32},Colon) after bound checking
## Indexing with StepRange calls AbstractArray{Real}
## Indexing with (StepRange,Colon) calls (AbstractArray{Real},Colon)
## Indexing with (Colon,StepRange) calls (Colon,AbstractArray{Real})
## Indexing with AbstractArray{Bool} calls KnetArray{Int32}; no need for bound checking
## Indexing with (Colon,AbstractVector{Bool}) calls (Colon,KnetArray{Int32}); no need for bound checking
## Indexing with (AbstractVector{Bool},Colon) calls (KnetArray{Int32},Colon); no need for bound checking
## Indexing with KnetArray{T} for logicals calls KnetArray{Int32}
## Indexing with Pair{Union{Real,AbstractUnitRange,Colon}}

# CuArray fallbacks: these rely on KnetArray<->CuArray conversions having shared pointers
# and cover the cases which are not explicitly defined using kernels below.

# Fallback for A[I...]  => getindex(A, I...)
function getindex(A::KnetArray, I...)
    _A = CuArray(A)
    _B = getindex(_A, I...)
    B = (_B isa CuArray ? KnetArray(_B) : _B)
    return B
end

# Fallback for A[I...] = B  => setindex!(A, B, I...)
function setindex!(A::KnetArray, B, I...)
    _A = CuArray(A)
    _B = (B isa KnetArray || B isa AbstractArray ? CuArray(B) : B)
    setindex!(_A, _B, I...)  ## This only works for x[I...] = y but not for x[I...] .= y
    return A
end

# Fallback for A[I...] .= B  => materialize!(dotview(A, I...), broadcasted(identity, B))
function materialize!(A::SubArray{T,N,<:KnetArray}, B) where {T,N}
    _A = view(CuArray(A.parent), A.indices...)
    _B = (B isa KnetArray || B isa AbstractArray ? CuArray(B) : B)
    materialize!(_A, _B)
    return A
end

# For contiguous I, dotview(A, I...) gives a shared-memory KnetArray rather than a view
function materialize!(A::KnetArray, B) where {T,N}
    _A = CuArray(A)
    _B = (B isa KnetArray || B isa AbstractArray ? CuArray(B) : B)
    materialize!(_A, _B)
    return A
end

# Ambiguity fixes:
function materialize!(A::SubArray{T,N,<:KnetArray}, B::Broadcasted{S}) where {T,N,S}
    _A = view(CuArray(A.parent), A.indices...)
    materialize!(_A, B)
    return A
end

function materialize!(A::KnetArray{T,N}, B::Broadcasted{S}) where {T,N,S}
    _A = CuArray(A)
    materialize!(_A, B)
    return A
end

# The following fallback version tried to do all allocations using KnetArrays but was recently broken (Issue 618).
# This only makes a difference if we are using the Knet allocator (i.e. cuallocator[] = false) which is not the default any more.
# Based on julia-1.4.2/base: getindex@abstractarray.jl:980, _getindex@multidimensional.jl:726, _unsafe_getindex!@multidimensional.jl:738
# function getindex_broken(A::KnetArray, I...)
#     _A = CuArray(A)
#     I = Base.to_indices(_A, I)
#     checkbounds(_A, I...)
#     shape = Base.index_shape(I...)
#     B = similar(A, length.(shape))
#     _B = CuArray(B)
#     Base._unsafe_getindex!(_B, _A, I...)
#     return B
# end

## Indexing with Real

function getindex(A::KnetArray{T}, I::Real) where {T}
    J = Int(I)
    if !(1 <= J <= length(A)); throw(BoundsError(A,J)); end
    copyto!(T[0], 1, A, J, 1)[1]
end

function setindex!(A::KnetArray{T}, v, I::Real) where {T}
    J = Int(I)
    if !(1 <= J <= length(A)); throw(BoundsError(A,J)); end
    copyto!(A, J, T[v], 1, 1)
    return A
end

## Indexing with Tuple{Real}
# Julia #14770
# If I is shorter than ndims(A) but longer than 1 the remaining indices assumed =1
# Also extra 1's at the end of I are ignored

function getindex(A::KnetArray{T}, I::Real...) where {T}
    J = ntuple(i->Int(I[i]), length(I)) # deprecated: Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        if !(1 <= J[j] <= size(A,j)); throw(BoundsError(A,J)); end
    end
    i = (LinearIndices(size(A)))[J...]
    copyto!(T[0], 1, A, i, 1)[1]
end

function setindex!(A::KnetArray{T}, v, I::Real...) where {T}
    J = ntuple(i->Int(I[i]), length(I)) # deprecated: Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        if !(1 <= J[j] <= size(A,j)); throw(BoundsError(A,J)); end
    end
    i = (LinearIndices(size(A)))[J...]
    copyto!(A, i, T[v], 1, 1)
    return A
end

## Indexing with CartesianIndex: calls Tuple{Real}

function getindex(A::KnetArray{T}, c::CartesianIndex) where {T}
    getindex(A, c.I...)
end

function setindex!(A::KnetArray{T}, v, c::CartesianIndex) where {T}
    setindex!(A, v, c.I...)
end

## Indexing with AbstractUnitRange
# We will implement indexing ranges as views not copies, if possible (when contiguous).
# For contiguous memory without stride all but the last >1 dimension must be full
# The original getindex(a,i:j...) for AbstractArray copies:
# function _getindex(l::LinearIndexing, A::AbstractArray, I::Union{Real, AbstractArray, Colon}...)
# in abstractarray.jl:487,multidimensional.jl:184.

function getindex(A::KnetArray{T}, I::AbstractUnitRange) where {T}
    if !(1 <= first(I) <= last(I) <= length(A)); throw(BoundsError(A,I)); end
    off = 1+(first(I)-1)*sizeof(T)
    len = length(I)*sizeof(T)
    ptr = KnetPtr(A.ptr, off, len)
    KnetArray{T,1}(ptr, (length(I),))
end

# Efficient fill:
for S in (32,64); T = Symbol("Float$S"); F = "fill_$S"
    @eval function unsafe_setindex!(a::KnetArray{$T},v::$T,I::AbstractUnitRange)
        @knet8($F,(Cint,$T,Ptr{$T}),length(I),v,pointer(a,first(I)))
        return a
    end
end

function setindex!(A::KnetArray{T}, v::Real, I::AbstractUnitRange) where {T}
    if !(1 <= first(I) <= last(I) <= length(A)); throw(BoundsError(A,I)); end
    if length(I)==0; return A; end
    unsafe_setindex!(A,T(v),I)
    return A
end

function setindex!(A::KnetArray{T}, v::Real, I::AbstractUnitRange{Bool}) where {T} # julia4 ambig fix
    if !(1 <= first(I) <= last(I) <= length(A)); throw(BoundsError(A,I)); end
    if length(I)==0; return A; end
    unsafe_setindex!(A,T(v),I)
    return A
end

function setindex!(A::KnetArray{T}, v, I::AbstractUnitRange) where {T}
    if !(1 <= first(I) <= last(I) <= length(A)); throw(BoundsError(A,I)); end
    if length(v)!=length(I); throw(DimensionMismatch()); end
    if length(I)==0; return A; end
    if eltype(v)!=T; v = convert(Array{T},v); end
    copyto!(A,first(I),v,1,length(I))
    return A
end

## Indexing with Colon
# Note that getindex(a,:) returns a view not a copy

function getindex(A::KnetArray, I::Colon)
    reshape(A,length(A))
end

function setindex!(A::KnetArray{T}, v::Real, I::Colon) where {T}
    if length(A)==0; return A; end
    unsafe_setindex!(A, T(v), 1:length(A))
    return A
end

function setindex!(A::KnetArray{T}, v, I::Colon) where {T}
    if length(v)!=length(A); throw(DimensionMismatch()); end
    if length(v)==0; return A; end
    if eltype(v)!=T; v = convert(Array{T},v); end
    copyto!(A,1,v,1,length(A))
    return A
end

for F in (32,64); T=Symbol("Float$F"); @eval begin

## Indexing with KnetArray{Int32}: low level, only Int32 supported, no bounds checking

    function unsafe_getindex!(x::KnetArray{$T}, y::KnetArray{$T}, i::KnetArray{Int32})
        @knet8($("getents_$F"),(Cint,Ptr{Int},Ptr{$T},Ptr{$T}), length(i), i, x, y)
        return y
    end

    function unsafe_setindex!(x::KnetArray{$T}, y::$T, i::KnetArray{Int32})
        @knet8($("setent1_$F"),(Cint,Ptr{Int},Ptr{$T},$T), length(i), i, x, y)
        return x
    end

    function unsafe_setindex!(x::KnetArray{$T}, y::KnetArray{$T}, i::KnetArray{Int32})
        @knet8($("setents_$F"),(Cint,Ptr{Int},Ptr{$T},Ptr{$T}), length(i), i, x, y)
        return x
    end

## Indexing with (Colon,KnetArray{Int32})
# TODO: Just special case rows and columns in matrices until we have a more general solution

    function unsafe_getindex!(x::KnetMatrix{$T}, y::KnetMatrix{$T}, ::Colon, i::KnetVector{Int32})
        @knet8($("getcols_$F"),(Cint,Cint,Cint,Ptr{Int},Ptr{$T},Ptr{$T}),
               size(x,1), size(x,2), length(i), i, x, y)
        return y
    end

    function unsafe_setindex!(x::KnetMatrix{$T}, y::$T, ::Colon, i::KnetVector{Int32})
        @knet8($("setcol1_$F"),(Cint,Cint,Cint,Ptr{Int},Ptr{$T},$T),
               size(x,1), size(x,2), length(i), i, x, y)
        return x
    end

    function unsafe_setindex!(x::KnetMatrix{$T}, y::KnetMatrix{$T}, ::Colon, i::KnetVector{Int32})
        @knet8($("setcols_$F"),(Cint,Cint,Cint,Ptr{Int},Ptr{$T},Ptr{$T}),
               size(x,1), size(x,2), length(i), i, x, y)
        return x
    end

## Indexing with (KnetArray{Int32},Colon)

    function unsafe_getindex!(x::KnetMatrix{$T}, y::KnetMatrix{$T}, i::KnetVector{Int32}, ::Colon)
        @knet8($("getrows_$F"),(Cint,Cint,Cint,Ptr{Int},Ptr{$T},Ptr{$T}),
               size(x,1), size(x,2), length(i), i, x, y)
        return y
    end

    function unsafe_setindex!(x::KnetMatrix{$T}, y::KnetMatrix{$T}, i::KnetVector{Int32}, ::Colon)
        @knet8($("setrows_$F"),(Cint,Cint,Cint,Ptr{Int},Ptr{$T},Ptr{$T}),
               size(x,1), size(x,2), length(i), i, x, y)
        return x
    end

    function unsafe_setindex!(x::KnetMatrix{$T}, y::$T, i::KnetVector{Int32}, ::Colon)
        @knet8($("setrow1_$F"),(Cint,Cint,Cint,Ptr{Int},Ptr{$T},$T),
               size(x,1), size(x,2), length(i), i, x, y)
        return x
    end

end; end

# bound checking

function checkbetween(i::AbstractArray{I},lo::L,hi::H) where {I<:Integer,L<:Integer,H<:Integer}
    checkbetween(Array{Int32}(i),Int32(lo),Int32(hi))
end

function checkbetween(i::Array{Int32},lo::Int32,hi::Int32)
    @inbounds for ii in i
        if !(lo <= ii <= hi)
            throw(BoundsError(lo:hi, ii))
        end
    end
end

## Indexing with AbstractArray{Real} calls KnetArray{Int32} after boundchecking

function getindex(x::KnetArray{T}, i::AbstractArray{I}) where {T,I<:Real}
    y = similar(x, size(i))
    if isempty(y); return y; end
    i = Array{Int32}(i)
    checkbetween(i, 1, length(x))
    unsafe_getindex!(x,y,KnetArray{Int32}(i))
    return y
end

function setindex!(x::KnetArray{T}, y::Real, i::AbstractArray{I}) where {T,I<:Real}
    if isempty(i); return x; end
    i = Array{Int32}(i)
    checkbetween(i, 1, length(x))
    unsafe_setindex!(x,T(y),KnetArray{Int32}(i))
    return x
end

function setindex!(x::KnetArray{T}, y, i::AbstractArray{I}) where {T,I<:Real}
    if length(y) != length(i); throw(DimensionMismatch()); end
    if isempty(i); return x; end
    i = Array{Int32}(i)
    checkbetween(i, 1, length(x))
    unsafe_setindex!(x,KnetArray{T}(y),KnetArray{Int32}(i))
    return x
end

## Indexing with (Colon,AbstractVector{Real}) calls (Colon,KnetArray{Int32}) after bound checking

function getindex(x::KnetMatrix{T}, c::Colon, i::AbstractVector{I}) where {T,I<:Real}
    xrows,xcols = size(x); ycols = length(i)
    y = similar(x, xrows, ycols)
    if isempty(y); return y; end
    i = Array{Int32}(i)
    checkbetween(i,1,xcols)
    unsafe_getindex!(x,y,c,KnetArray{Int32}(i))
    return y
end

function setindex!(x::KnetMatrix{T}, y::Real, c::Colon, i::AbstractVector{I}) where {T,I<:Real}
    if isempty(i); return x; end
    xrows,xcols = size(x); ycols=length(i)
    i = Array{Int32}(i)
    checkbetween(i,1,xcols)
    unsafe_setindex!(x,T(y),c,KnetArray{Int32}(i))
    return x
end

function setindex!(x::KnetMatrix{T}, y, c::Colon, i::AbstractVector{I}) where {T,I<:Real}
    if ndims(y) != 2; throw(DimensionMismatch()); end
    xrows,xcols = size(x); yrows,ycols=size(y)
    if yrows != xrows; throw(DimensionMismatch()); end
    if ycols != length(i); throw(DimensionMismatch()); end
    if isempty(y); return x; end
    i = Array{Int32}(i)
    checkbetween(i,1,xcols)
    unsafe_setindex!(x,KnetArray{T}(y),c,KnetArray{Int32}(i))
    return x
end

# function setindex!(x::KnetMatrix{T}, y, c::AbstractUnitRange, i::AbstractVector{I}) where {T,I<:Real}
#     if c == 1:size(x,1)
#         setindex!(x, y, :, i)
#     else
#         throw(MethodError(setindex!,x,y,c,i))
#     end
# end

## Indexing with (AbstractVector{Real},Colon) calls (KnetArray{Int32},Colon) after bound checking

function getindex(x::KnetMatrix{T}, i::AbstractVector{I}, c::Colon) where {T,I<:Real}
    xrows,xcols = size(x); yrows = length(i)
    y = similar(x, yrows, xcols)
    if isempty(y); return y; end
    i = Array{Int32}(i)
    checkbetween(i,1,xrows)
    unsafe_getindex!(x,y,KnetArray{Int32}(i),c)
    return y
end

function setindex!(x::KnetMatrix{T}, y::Real, i::AbstractVector{I}, c::Colon) where {T,I<:Real}
    if isempty(i); return x; end
    xrows,xcols = size(x); yrows=length(i)
    i = Array{Int32}(i)
    checkbetween(i,1,xrows)
    unsafe_setindex!(x,T(y),KnetArray{Int32}(i),c)
    return x
end

function setindex!(x::KnetMatrix{T}, y, i::AbstractVector{I}, c::Colon) where {T,I<:Real}
    if ndims(y) != 2; throw(DimensionMismatch()); end
    xrows,xcols = size(x); yrows,ycols=size(y)
    if ycols != xcols; throw(DimensionMismatch()); end
    if yrows != length(i); throw(DimensionMismatch()); end
    if isempty(y); return x; end
    i = Array{Int32}(i)
    checkbetween(i,1,xrows)
    unsafe_setindex!(x,KnetArray{T}(y),KnetArray{Int32}(i),c)
    return x
end

# function setindex!(x::KnetMatrix{T}, y, i::AbstractVector{I}, c::AbstractUnitRange) where {T,I<:Real}
#     if c == 1:size(x,2)
#         setindex!(x, y, i, :)
#     else
#         throw(MethodError(setindex!,x,y,i,c))
#     end
# end

## Indexing with AbstractArray{CartesianIndex} calls AbstractArray{Real}

c2i(d::Dims,i::AbstractArray{I}) where {I<:CartesianIndex} = Int32[(LinearIndices(d))[c.I...] for c in i]
getindex(x::KnetArray{T}, i::AbstractArray{I}) where {T,I<:CartesianIndex} = getindex(x, c2i(size(x),i))
setindex!(x::KnetArray{T}, y, i::AbstractArray{I}) where {T,I<:CartesianIndex} = setindex!(x, y, c2i(size(x),i))

## Indexing with Empty Array or other unrecognized AbstractArray calls AbstractArray{Real}

getindex(x::KnetArray{T}, i::AbstractArray) where {T}=getindex(x, Array{Int32}(i))
setindex!(x::KnetArray{T}, y, i::AbstractArray) where {T}=setindex!(x, y, Array{Int32}(i))

## Indexing with StepRange calls AbstractArray{Real}

function getindex(A::KnetArray{T}, I::StepRange) where {T}
    getindex(A, collect(I))
end

function setindex!(A::KnetArray{T}, v::Real, I::StepRange{R}) where {T,R<:Real} # julia4 ambiguity fix
    setindex!(A, v, collect(I))
end

function setindex!(A::KnetArray{T}, v::Real, I::StepRange{Bool}) where {T} # julia4 ambiguity fix
    setindex!(A, v, collect(I))
end

function setindex!(A::KnetArray{T}, v, I::StepRange) where {T}
    setindex!(A, v, collect(I))
end

## Indexing with (StepRange,Colon) calls (AbstractArray{Real},Colon)

function getindex(A::KnetMatrix{T}, I::StepRange, c::Colon) where {T}
    getindex(A, collect(I), c)
end

function setindex!(A::KnetMatrix{T}, v::Real, I::StepRange{R}, c::Colon) where {T,R<:Real} # julia4 ambig fix
    setindex!(A, v, collect(I), c)
end

function setindex!(A::KnetMatrix{T}, v::Real, I::StepRange{Bool}, c::Colon) where {T} # julia4 ambig fix
    setindex!(A, v, collect(I), c)
end

function setindex!(A::KnetMatrix{T}, v, I::StepRange, c::Colon) where {T}
    setindex!(A, v, collect(I), c)
end

# function setindex!(A::KnetMatrix, v, I::StepRange, r::AbstractUnitRange)
#     if r == 1:size(A,2)
#         setindex!(A,v,I,:)
#     else
#         throw(MethodError(setindex!, A, v, I, r))
#     end
# end

## Indexing with (Colon,StepRange) calls (Colon,AbstractArray{Real})

function getindex(A::KnetMatrix{T}, c::Colon, I::StepRange) where {T}
    getindex(A, c, collect(I))
end

function setindex!(A::KnetMatrix{T}, v::Real, c::Colon, I::StepRange{R}) where {T,R<:Real} # julia4 ambig fix
    setindex!(A, v, c, collect(I))
end

function setindex!(A::KnetMatrix{T}, v::Real, c::Colon, I::StepRange{Bool}) where {T} # julia4 ambig fix
    setindex!(A, v, c, collect(I))
end

function setindex!(A::KnetMatrix{T}, v, c::Colon, I::StepRange) where {T}
    setindex!(A, v, c, collect(I))
end

# function setindex!(A::KnetMatrix, v, r::AbstractUnitRange, I::StepRange)
#     if r == 1:size(A,1)
#         setindex!(A,v,:,I)
#     else
#         throw(MethodError(setindex!, A, v, r, I))
#     end
# end

## Indexing with AbstractArray{Bool} calls KnetArray{Int32}; no need for bound checking

function getindex(x::KnetArray{T}, i::AbstractArray{Bool}) where {T}
    if length(i) != length(x); throw(BoundsError(x,i)); end
    j = findall(vec(i))
    y = similar(x, length(j))
    if !isempty(y); unsafe_getindex!(x,y,KnetArray{Int32}(j)); end
    return y
end

function setindex!(x::KnetArray{T}, y::Real, i::AbstractArray{Bool}) where {T}
    if length(i) != length(x); throw(DimensionMismatch()); end
    j = findall(vec(i))
    if !isempty(j); unsafe_setindex!(x,T(y),KnetArray{Int32}(j)); end
    return x
end

function setindex!(x::KnetArray{T}, y, i::AbstractArray{Bool}) where {T}
    if length(i) != length(x); throw(DimensionMismatch()); end
    j = findall(vec(i))
    if length(j) != length(y); throw(BoundsError(y,j)); end
    if !isempty(j); unsafe_setindex!(x,KnetArray{T}(y),KnetArray{Int32}(j)); end
    return x
end

## Indexing with (Colon,AbstractVector{Bool}) calls (Colon,KnetArray{Int32}); no need for bound checking

function getindex(x::KnetMatrix{T}, c::Colon, i::AbstractVector{Bool}) where {T}
    xrows,xcols = size(x)
    if length(i) != xcols; throw(BoundsError(x,(:,i))); end
    j = findall(vec(i)); ycols = length(j)
    y = similar(x, xrows, ycols)
    if !isempty(y); unsafe_getindex!(x,y,c,KnetArray{Int32}(j)); end
    return y
end

function setindex!(x::KnetMatrix{T}, y::Real, c::Colon, i::AbstractVector{Bool}) where {T}
    xrows,xcols = size(x)
    if length(i) != xcols; throw(BoundsError(x,(:,i))); end
    j = findall(vec(i))
    if !isempty(j); unsafe_setindex!(x,T(y),c,KnetArray{Int32}(j)); end
    return x
end

function setindex!(x::KnetMatrix{T}, y, c::Colon, i::AbstractVector{Bool}) where {T}
    if ndims(y) != 2; throw(DimensionMismatch()); end
    xrows,xcols = size(x); yrows,ycols=size(y)
    if yrows != xrows; throw(DimensionMismatch()); end
    if length(i) != xcols; throw(BoundsError(x,(:,i))); end
    j = findall(vec(i))
    if ycols != length(j); throw(DimensionMismatch()); end
    if !isempty(y); unsafe_setindex!(x,KnetArray{T}(y),c,KnetArray{Int32}(j)); end
    return x
end

# function setindex!(x::KnetMatrix, y, c::AbstractUnitRange, i::AbstractVector{Bool})
#     if c == 1:size(x,1)
#         setindex!(x,y,:,i)
#     else
#         throw(MethodError(setindex!,x,y,c,i))
#     end
# end

## Indexing with (AbstractVector{Bool},Colon) calls (KnetArray{Int32},Colon); no need for bound checking

function getindex(x::KnetMatrix{T}, i::AbstractVector{Bool}, c::Colon) where {T}
    xrows,xcols = size(x)
    if length(i) != xrows; throw(BoundsError(x,(i,:))); end
    j = findall(vec(i)); yrows = length(j)
    y = similar(x, yrows, xcols)
    if !isempty(y); unsafe_getindex!(x,y,KnetArray{Int32}(j),c); end
    return y
end

function setindex!(x::KnetMatrix{T}, y::Real, i::AbstractVector{Bool}, c::Colon) where {T}
    xrows,xcols = size(x)
    if length(i) != xrows; throw(BoundsError(x,(i,:))); end
    j = findall(vec(i))
    if !isempty(j); unsafe_setindex!(x,T(y),KnetArray{Int32}(j),c); end
    return x
end

function setindex!(x::KnetMatrix{T}, y, i::AbstractVector{Bool}, c::Colon) where {T}
    if ndims(y) != 2; throw(DimensionMismatch()); end
    xrows,xcols = size(x); yrows,ycols=size(y)
    if ycols != xcols; throw(DimensionMismatch()); end
    if length(i) != xrows; throw(BoundsError(x,(i,:))); end
    j = findall(vec(i))
    if yrows != length(j); throw(DimensionMismatch()); end
    if !isempty(y); unsafe_setindex!(x,KnetArray{T}(y),KnetArray{Int32}(j),c); end
    return x
end

# function setindex!(x::KnetMatrix, y, i::AbstractVector{Bool}, c::AbstractUnitRange)
#     if c == 1:size(x,2)
#         setindex!(x,y,i,:)
#     else
#         throw(MethodError(setindex!,x,y,i,c))
#     end
# end

## Indexing with KnetArray{T} for logicals calls KnetArray{Int32}
# Need this because (k.<0) returns KnetArray{T} instead of BitArray

function getindex(x::KnetArray{T}, i::KnetArray{T}) where {T}
    if length(i) != length(x); throw(BoundsError(x,i)); end
    j = findall(Array(i))
    y = similar(x, length(j))
    if !isempty(y); unsafe_getindex!(x,y,KnetArray{Int32}(j)); end
    return y
end

function setindex!(x::KnetArray{T}, y::Real, i::KnetArray{T}) where {T}
    if length(i) != length(x); throw(DimensionMismatch()); end
    j = findall(Array(i))
    if !isempty(j); unsafe_setindex!(x,T(y),KnetArray{Int32}(j)); end
    return x
end

function setindex!(x::KnetArray{T}, y, i::KnetArray{T}) where {T}
    if length(i) != length(x); throw(DimensionMismatch()); end
    j = findall(Array(i))
    if length(j) != length(y); throw(BoundsError(y,j)); end
    if !isempty(j); unsafe_setindex!(x,KnetArray{T}(y),KnetArray{Int32}(j)); end
    return x
end

# To avoid ambiguity with previous definitions we have these
# TODO: clean this up...
getindex(A::KnetMatrix, I1::Colon, I2::Colon)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B, I1::Colon, I2::Colon)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::Colon, I2::Colon)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::AbstractUnitRange, I2::AbstractUnitRange)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B, I1::AbstractUnitRange, I2::AbstractUnitRange)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::AbstractUnitRange, I2::AbstractUnitRange)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::Colon, I2::AbstractUnitRange)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B::Real, I1::Colon, I2::AbstractUnitRange{Bool})=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Real, I1::Colon, I2::AbstractUnitRange{T}) where {T<:Real}=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B, I1::Colon, I2::AbstractUnitRange)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::Colon, I2::AbstractUnitRange)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::AbstractUnitRange, I2::Colon)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B::Real, I1::AbstractUnitRange{Bool}, I2::Colon)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Real, I1::AbstractUnitRange{R}, I2::Colon) where {R<:Real}=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B, I1::AbstractUnitRange, I2::Colon)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::AbstractUnitRange, I2::Colon)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::Real, I2::AbstractUnitRange)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B, I1::Real, I2::AbstractUnitRange)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::Real, I2::AbstractUnitRange)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::AbstractUnitRange, I2::Real)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B, I1::AbstractUnitRange, I2::Real)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::AbstractUnitRange, I2::Real)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::Colon, I2::Real)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B, I1::Colon, I2::Real)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::Colon, I2::Real)=setindex2!(A,B,I1,I2)
getindex(A::KnetMatrix, I1::Real, I2::Colon)=getindex2(A,I1,I2)
setindex!(A::KnetMatrix, B, I1::Real, I2::Colon)=setindex2!(A,B,I1,I2)
setindex!(A::KnetMatrix, B::Number, I1::Real, I2::Colon)=setindex2!(A,B,I1,I2)

## Indexing with Pair{Union{Real,AbstractUnitRange,Colon}}
# TODO: the following getindex, setindex! work for 1 and 2 dimensions only, write general versions.

const Index3 = Union{Real,AbstractUnitRange,Colon}

function getindex2(A::KnetMatrix{T}, I1::Index3, I2::Index3) where {T}
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I1,I2)
    B1 = isa(I1,Colon) ? size(A,1) : length(I1)
    B2 = isa(I2,Colon) ? size(A,2) : length(I2)
    Bsize = isa(I1,Real) ? (B2,) : isa(I2,Real) ? (B1,) : (B1,B2)
    Bdims = length(Bsize)
    if ncols == 1
        off = 1+(firstindex-1)*sizeof(T)
        len = nrows*sizeof(T)
        ptr = KnetPtr(A.ptr, off, len)
        KnetArray{T,Bdims}(ptr, Bsize)
    else
        B = similar(A, Bsize)
        if isempty(B); return B; end
        nrows *= sizeof(T); astep *= sizeof(T)
        @knet8(xcopy,(Cint,Cint,Cptr,Cint,Cptr,Cint),
               nrows, ncols, pointer(A,firstindex), astep, B, nrows)
        return B
    end
end

function setindex2!(A::KnetMatrix{T}, B, I1::Index3, I2::Index3) where {T}
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I1,I2)
    aptr0 = pointer(A, firstindex)
    if isa(B,Number)
        B = T(B)
        if ncols == 1
            if nelts > 0
                if T <: Float32
                    @knet8(fill_32,(Cint,Cfloat, Ptr{Cfloat}), nelts,B,aptr0)
                elseif T<: Float64
                    @knet8(fill_64,(Cint,Cdouble,Ptr{Cdouble}),nelts,B,aptr0)
                else
                    error("$T not supported")
                end
            end
        elseif nrows > 0 && ncols > 0
            if T <: Float32
                @knet8(xfill_32,(Cint,Cint,Cfloat, Ptr{Cfloat}, Cint),nrows,ncols,B,aptr0,astep)
            elseif T<: Float64
                @knet8(xfill_64,(Cint,Cint,Cdouble,Ptr{Cdouble},Cint),nrows,ncols,B,aptr0,astep)
            else
                error("$T not supported")
            end
        end
    else
        length(B) == nelts || throw(DimensionMismatch())
        B = convert(KnetArray{T},B)
        if ncols == 1
            if nelts > 0
                # @cudart(cudaMemcpy,(Cptr,Cptr,Csize_t,UInt32),
                #         aptr0, B, nelts*sizeof(T), cudadir(A,B))
                cuMemcpyDtoD_v2(CuPtr{Nothing}(UInt(aptr0)), CuPtr{Nothing}(UInt(pointer(B))), nelts*sizeof(T))
            end
        elseif nrows > 0 && ncols > 0
            nrows *= sizeof(T); astep *= sizeof(T)
            @knet8(xcopy,(Cint,Cint,Cptr,Cint,Cptr,Cint), nrows, ncols, B, nrows, aptr0, astep)
        end
    end
    return A
end

function indexparams(A::KnetArray{T,N}, I::Index3...) where {T,N}
    N > 2 && error("setindex for ndims > 2 not implemented yet")
    skipped = false
    nrows = nelts = 1
    subs1 = ones(Int,2)
    astep = length(A)
    for i=1:length(I)
        Ii = I[i]
        Ai = size(A,i)
        if isa(Ii, Colon)
            Li = Ai
            subs1[i] = 1
        elseif isa(Ii, Real)
            1 <= Ii <= Ai || throw(DimensionMismatch())
            Li = 1
            subs1[i] = Int(Ii)
        else
            1 <= first(Ii) <= last(Ii) <= Ai || throw(DimensionMismatch())
            Li = length(Ii)
            subs1[i] = first(Ii)
        end
        nelts *= Li
        if !skipped
            nrows *= Li
            if Li < Ai
                skipped = true
                astep = stride(A,i+1)
            end
        end
    end
    ncols = div(nelts, nrows)
    firstindex = (LinearIndices(size(A)))[subs1...]
    return (nelts,nrows,ncols,firstindex,astep)
end

# These two are not sufficient in spite of what the documentation says:
# display goes into an infinite loop!
# getindex{T}(A::KnetArray{T}, i::Int)=copyto!(T[0], 1, A, i, 1)[1]
# setindex!{T}(A::KnetArray{T}, v, i::Int)=copyto!(A, i, T[v], 1, 1)

# We need x[:,:,t] and hx[:,:,l] for RNNs
function getindex(A::KnetArray, ::Colon, ::Colon, I::Real)
    B = reshape(A, stride(A,3), size(A,3))
    reshape(B[:,I], size(A,1), size(A,2))
end

function getindex(A::KnetArray, ::Colon, ::Colon, ::Colon)
    A
end

function getindex(A::KnetArray, ::Colon, ::Colon, I::Index3)
    B = reshape(A, stride(A,3), size(A,3))
    reshape(B[:,I], size(A,1), size(A,2), length(I))
end

function setindex!(x::KnetArray, y, ::Colon, ::Colon, I::Index3)
    reshape(x, stride(x,3), size(x,3))[:,I] = y
    return x
end

function setindex!(x::KnetArray, y, ::Colon, ::Colon, ::Colon)
    copyto!(x,y)
    return x
end

# function setindex!(x::KnetArray, y, i::AbstractUnitRange, j::AbstractUnitRange, k::Index3)
#     if first(i) == 1 && last(i) == size(x,1) && first(j) == 1 && last(j) == size(x,2)
#         setindex!(x, y, :, :, k)
#     else
#         throw(MethodError(setindex!, (x,y,i,j,k)))
#     end
# end

function getindex(x::KnetArray{T,2}, ::Colon, m::AbstractArray{I,2}) where {T,I<:Integer}
    reshape(x[:,vec(m)], size(x,1), size(m,1), size(m,2))
end


text:
mutable struct SNF
    x::Vector{Int}
    minor::Bool
    left::Int
    right::Int
end

SNF() = SNF(fill(-1, 31), true, 0, 0)

function SNF(s::AbstractString)
    x = fill(-1, 31)
    i = 1
    for c in s
        if c == '['
            i <<= 1
        elseif c == ']'
            i >>= 1
        elseif c != ','
            x[i] = c - '0'
            i += 1
        end
    end
    return SNF(x, false, 0, 0)
end

function add!(z::SNF, x::SNF, y::SNF)
    explode!(x)
    explode!(y)
    z.minor = false
    z.x[31] = y.x[15]
    z.x[30] = y.x[14]
    z.x[29] = y.x[13]
    z.x[28] = y.x[12]
    z.x[27] = y.x[11]
    z.x[26] = y.x[10]
    z.x[25] = y.x[9]
    z.x[24] = y.x[8]
    z.x[23] = x.x[15]
    z.x[22] = x.x[14]
    z.x[21] = x.x[13]
    z.x[20] = x.x[12]
    z.x[19] = x.x[11]
    z.x[18] = x.x[10]
    z.x[17] = x.x[9]
    z.x[16] = x.x[8]
    z.x[15] = y.x[7]
    z.x[14] = y.x[6]
    z.x[13] = y.x[5]
    z.x[12] = y.x[4]
    z.x[11] = x.x[7]
    z.x[10] = x.x[6]
    z.x[9] = x.x[5]
    z.x[8] = x.x[4]
    z.x[7] = y.x[3]
    z.x[6] = y.x[2]
    z.x[5] = x.x[3]
    z.x[4] = x.x[2]
    z.x[3] = y.x[1]
    z.x[2] = x.x[1]
    z.x[1] = -1
    if x.right >= 0
        if y.left >= 0
            z.x[23] = x.right + y.left
        else
            distribute!(z.x, 24, x.right)
        end
    else
        if y.left >= 0
            distribute!(z.x, 23, y.left)
        end
    end
    i = 1
    while true
        while z.x[i] < 0
            i <<= 1
        end
        if z.x[i] > 9
            if i < 16
                a = z.x[i] >> 1
                b = z.x[i] - a
                z.x[i] = -1
                z.x[i << 1] = a
                z.x[(i << 1) | 1] = b
                if z.x[i << 1] > 9
                    i <<= 1
                    continue
                elseif z.x[(i << 1) | 1] > 9
                    i <<= 1
                    i |= 1
                    continue
                end
            else
                a = z.x[i] >> 1
                b = z.x[i] - a
                z.x[i] = 0
                if (i & (i + 1)) != 0
                    distribute!(z.x, i + 1, b)
                end
                if (i & (i - 1)) != 0
                    j = distribute!(z.x, i - 1, a)
                    if z.x[j] > 9
                        i = j
                        continue
                    end
                end
            end
        end
        while (i & 1) != 0
            i >>= 1
        end
        i == 0 && break
        i += 1
    end
end

function explode!(x::SNF)
    x.minor && return
    x.left = -1
    x.right = -1
    for i = 16:2:31
        if x.x[i] >= 0
            a = x.x[i]
            b = x.x[i + 1]
            x.x[i] = -1
            x.x[i + 1] = -1
            x.x[i >> 1] = 0
            if i == 16
                x.left = a
            else
                distribute!(x.x, i - 1, a)
            end
            if i == 30
                x.right = b
            else
                distribute!(x.x, i + 2, b)
            end
        end
    end
    x.minor = true
    return
end

function distribute!(x, i, v)
    while i < 16
        i <<= 1
    end
    while i > 0
        if x[i] >= 0
            x[i] += v
            return i
        end
        i >>= 1
    end
    @assert false
end

function magnitude(x::SNF, i = 1)
    if x.x[i] >= 0
        return x.x[i]
    else
        return 3 * magnitude(x, i << 1) + 2 * magnitude(x, (i << 1) | 1)
    end
end

function part1(input)
    numbers = SNF.(readlines(input))
    s = numbers[1]
    for i = 2:length(numbers)
        add!(s, s, numbers[i])
    end
    return magnitude(s)
end

function part2(input)
    numbers = SNF.(readlines(input))
    s = SNF()
    max_magnitude = 0
    for i = 1:length(numbers)
        for j = 1:length(numbers)
            i == j && continue
            add!(s, numbers[i], numbers[j])
            max_magnitude = max(max_magnitude, magnitude(s))
        end
    end
    return max_magnitude
end