text:
using PyPlot
using CSV
using DataFrames
using SparseArrays
using Statistics
import StatsBase: sample
using LinearAlgebra
using GPUMatFac

showln(x) = (show(x); println())

################################
# LOAD DATA
data = CSV.read("mnist_matrix.tsv", DataFrame; header=false, delim="\t");
#data = convert(Matrix, data);
data = Matrix(data);
full_X = data[:,1:end-1];
full_labels = data[:,end];
classes = sort(unique(full_labels));

################################
# TRAIN/TEST SPLIT
M = size(data,1)

train_idx = sort(sample(1:M, Int(0.9*M), replace=false))
train_idx_set = Set(train_idx)
test_idx = Int[i for i=1:M if !in(i, train_idx_set)]

train_X = full_X[train_idx,:]
train_labels = full_labels[train_idx]
showln("TRAIN SET: ")
showln(size(train_X))

test_X = full_X[test_idx,:]
test_labels = full_labels[test_idx]
showln("TEST SET: ")
showln(size(test_X))

full_X = vcat(train_X, test_X)
full_labels = vcat(train_labels, test_labels)
showln("FULL DATA: ")
showln(size(full_X))

################################
# SOME USEFUL FUNCTIONS
unflatten_mnist(vec) = transpose(reshape(vec, (28,28)))

function coord_to_idx(x, y, H)
    return (x-1)*H + y
end

function idx_to_pixel_coord(idx, H)
    return (div((idx-1),H), ((idx-1) % H))
end

################################
# BUILD THE PIXEL GRAPH
function build_pixel_graph(H, W)
    pixel_graph = [Int64[] for idx=1:(H*W)]

    for x=1:W
        for y=1:H
            u_idx = coord_to_idx(x,y, H)

            # Left
            if x > 1
                v_idx = coord_to_idx(x-1, y, H)
                push!(pixel_graph[u_idx], v_idx)
                push!(pixel_graph[v_idx], u_idx)
            end
            # Up
            if y > 1
                v_idx = coord_to_idx(x, y-1, H)
                push!(pixel_graph[u_idx], v_idx)
                push!(pixel_graph[v_idx], u_idx)
            end

        end
    end

    return pixel_graph
end

function plot_graph(g, idx_to_coord_func)
    for (u, neighbors) in enumerate(g)
        u_xy = idx_to_coord_func(u)
        for v in neighbors
            if v < u
                v_xy = idx_to_coord_func(v)
                plot([u_xy[1], v_xy[1]], [u_xy[2], v_xy[2]])
            end
        end
    end
end

pixel_graph = build_pixel_graph(28,28);

mnist_idx_to_coord_func(idx) = idx_to_pixel_coord(idx,28)

####################################
# BUILD THE TRAIN INSTANCE GRAPH
function build_instance_graph(labels)

    classes = sort(unique(labels))
    class_to_idx = Dict(cls => length(labels)+i for (i, cls) in enumerate(classes))
    n_classes = length(classes)
    n_bk = n_classes #+ 1
                                      # Leaf nodes;  # class nodes;  # root node
    instance_graph = [Int64[] for i=1:(length(labels) + n_bk) ];

    for (leaf_idx, label) in enumerate(labels)
        cls_idx = class_to_idx[label]
        push!(instance_graph[leaf_idx], cls_idx)
        push!(instance_graph[cls_idx], leaf_idx)
    end

    #root_idx = length(instance_graph)
    #for cls in classes
    #    cls_idx = class_to_idx[cls]
    #    push!(instance_graph[cls_idx], root_idx)
    #    push!(instance_graph[root_idx], cls_idx)
    #end

    return instance_graph
end

####################################
# BUILD THE TRAIN INSTANCE MATRIX
train_instance_graph = build_instance_graph(train_labels);
#aug_X = vcat( train_X, fill(NaN, length(classes)+1, 28*28) );
aug_X = vcat( train_X, fill(NaN, length(classes), 28*28) );

showln("AUGMENTED X: ")
showln(size(aug_X))

showln("BUILT PIXEL AND INSTANCE GRAPHS")

function graph_to_spmat(graph; epsilon=0.001, standardize=false)

    I = Int64[]
    J = Int64[]
    V = Float64[]
    N = length(graph)

    diag_entries = fill(epsilon, N)

    for (u, neighbors) in enumerate(graph)
        for v in neighbors
            push!(I, u)
            push!(J, v)
            push!(V, -1.0)
        end
        diag_entries[u] += length(neighbors)
    end
    for (i, v) in enumerate(diag_entries)
        push!(I, i)
        push!(J, i)
        push!(V, v)
    end

    mat = sparse(I, J, V, N, N)

    if standardize
        standardizer = sparse(1:N, 1:N, 1.0./sqrt.(diag_entries), N, N)
        mat = standardizer * mat * standardizer
    end

    return mat
end

pixel_spmat = graph_to_spmat(pixel_graph; standardize=true);
train_inst_spmat = graph_to_spmat(train_instance_graph; standardize=true);

#M_train = length(train_labels)
#train_inst_spmat = sparse(I, M_train, M_train)

####################################
# BUILD THE TEST SET INSTANCE MATRIX
function build_test_graph(train_graph, train_labels, test_labels)

    M_train = length(train_labels)
    M_train_aug = length(train_graph)
    M_test = length(test_labels)
    M_full = M_train_aug + M_test

    classes = sort(unique(train_labels))
    n_classes = length(classes)
    n_bk = n_classes + 1
    @assert n_bk == M_train_aug - M_train

    label_to_idx = Dict(l => (M_train + i) for (i,l) in enumerate(classes))

    test_graph = copy(train_graph)

    for i=1:M_test
        test_idx = M_train_aug + i
        parent_idx = label_to_idx[test_labels[i]]

        push!(test_graph, [parent_idx])
        push!(test_graph[parent_idx], test_idx)
    end

    return test_graph
end

#test_instance_graph = build_test_graph(train_instance_graph, train_labels, test_labels)
#test_spmat = graph_to_spmat(test_instance_graph)
#println("TEST SPARSE MATRIX:")
#showln(test_spmat)

showln("BUILT PIXEL AND INSTANCE MATRICES")

############################################
# Make a matrix of instance covariates -- just
# a one-hot encoding of the label
train_covariates = zeros(length(train_labels) + 11, 10)
for i=1:length(train_labels)
    train_covariates[i,Integer(train_labels[i])+1] = 1.0
end
for i=1:10
    train_covariates[length(train_labels)+i,i] = 1.0
end
#train_covariates[end,:] .= 0.1

############################################
# CONSTRUCT THE MODEL
k = 10

feature_spmat_vec = [copy(pixel_spmat) for i=1:k-1]
train_instance_spmat_vec = [copy(train_inst_spmat) for i=1:k]
#test_instance_spmat_vec = [copy(test_spmat) for i=1:k]

losses = [LogisticLoss(1.0) for i=1:28*28]
model = MatFacModel(train_instance_spmat_vec, feature_spmat_vec, losses;
                    #instance_covariate_coeff_reg=pixel_spmat,
                    feature_offset_reg=pixel_spmat, K=k);

showln("INITIALIZED MODEL")

############################################
# TRAIN MODEL
showln("ABOUT TO FIT")
#fit!(model, aug_X; instance_covariates=train_covariates,
#     method="nesterov", inst_reg_weight=0.0, feat_reg_weight=1.0,
#     max_iter=2000, loss_iter=10, lr=0.5, momentum=0.1,
#     rel_tol=1e-9, a_0_tau=1.0, b_0_tau=1.0)
fit!(model, aug_X;
     #instance_covariates=train_covariates,
     inst_reg_weight=160.0, feat_reg_weight=10.0,
     max_iter=1000, loss_iter=10, lr=0.1,
     rel_tol=1e-9, a_0_tau=1.0, b_0_tau=1.0)
#fit!(model, aug_X; inst_reg_weight=0.1, feat_reg_weight=0.1, max_iter=1000, lr=0.5, rel_tol=1e-9)
#fit!(model, aug_X; method="line_search", inst_reg_weight=0.1, feat_reg_weight=0.1, alpha=0.1, line_search_max_iter=5, grow=1.5, shrink=0.5, c1=1e-4, c2=0.5, max_iter=1000, rel_tol=1e-6)

println("FACTOR IMPORTANCES:")
println(model.factor_importances)

#####################
# SAVE MODEL
save_hdf("mnist_model.hdf", model)

#####################
# RELOAD MODEL
model = load_hdf("mnist_model.hdf")

train_labels = convert(Vector{Int64}, train_labels)
test_labels = convert(Vector{Int64}, test_labels)
colors = collect(keys(PyPlot.colorsm.TABLEAU_COLORS))

function embedding_scatter(X, labels; F=nothing)
    projected = nothing
    if F == nothing
        F = LinearAlgebra.svd(X)
        projected = F.Vt
    else
        projected = diagm(1.0./F.S) * (transpose(F.U)* X)
    end
    M = length(labels)
    for lab in sort(unique(labels))
        lab_idx = findall(labels .== lab)
        scatter3D(projected[1, lab_idx], projected[2,lab_idx], projected[3,lab_idx], color=colors[lab+1], label=string(lab))
    end
    legend()

    return F
end

embedding_F = embedding_scatter(model.X[1:(end-1),:], train_labels)
savefig("embedding_scatter.png", dpi=200)

for i=1:k
    matshow(unflatten_mnist(model.Y[i,:]))
    colorbar()
    savefig(string("embedding_basis_",i,".png"), dpi=200)
end

matshow(unflatten_mnist(model.feature_offset))
colorbar()
savefig(string("embedding_offset.png"), dpi=200)

for i=1:size(model.instance_covariate_coeff,1)
    matshow(unflatten_mnist(model.instance_covariate_coeff[i,:]))
    colorbar()
    savefig(string("covariate_coeff_",i,".png"), dpi=200)
end

######################################
# Use the model to impute data
Z = impute_values(model)
println(size(Z))

#####################################
# TRANSFORM THE HELDOUT DATA

new_factor = 0.001*randn(k, length(test_labels))
new_factor[k,:] .= 1.0

showln("NEW FACTOR:")
showln(size(new_factor))

# TODO UPDATE THE TRANSFORM METHOD
#transformed = GPUMatFac.transform(model, test_X; inst_reg_weight=0.1, max_iter=1000, loss_iter=1, lr=0.1, momentum=0.8, rel_tol=1e-9, X_new=new_factor) #, new_inst_reg_mats=test_instance_spmat_vec)

#println("TRANSFORMED:")
#println(size(transformed))

#embedding_scatter(transformed[1:(end-1),:], test_labels; F=embedding_F)
#savefig("transformed_scatter.png", dpi=200)

ext: jl


text:
using TransmuteDims
using Test, LinearAlgebra, Random

@testset "utils" begin

    using TransmuteDims: sanitise_zero, genperm_zero, invperm_zero, is_off_diag

    @test sanitise_zero((1,2,3), Val(2)) === (1,2,0)
    @test sanitise_zero([1,2,nothing,4], Val(3)) === (1,2,0,0)

    @test genperm_zero((:a,:b,:c,:d), (3,2,4,1)) == Base.PermutedDimsArrays.genperm((:a,:b,:c,:d), (3,2,4,1))
    @test genperm_zero((:a,:b,:c,:d), (3,0,4,1)) == (:c, 1, :d, :a)
    @test genperm_zero((:a,:b,:c,:d), (3,0,4,1), :z) == (:c, :z, :d, :a)

    @test invperm_zero((5,1,3,4,2),5,:A)[1] == invperm((5,1,3,4,2))
    @test invperm_zero((1,2,3),4,:A)[1] == (1,2,3,0)
    @test invperm_zero((1,3),4,:A)[2][2] isa Expr  # two assertions, size(A,2)==1 && size(A,4)==1

    using TransmuteDims: increasing_or_zero, unique_or_zero, may_reshape

    @test increasing_or_zero((1,2,0,3)) == true
    @test increasing_or_zero((1,2,0,2)) == false

    @test unique_or_zero((1,3,0,2)) == true
    @test unique_or_zero((0,1,3)) == true
    @test unique_or_zero((0,1,0,3,0)) == true
    @test unique_or_zero((0,1,1)) == false
    @test unique_or_zero((0,1,2,0,2)) == false

    @test may_reshape(typeof(rand(3))) == true
    @test may_reshape(typeof(rand(3)')) == false
    @test may_reshape(typeof(reshape(1:4,2,2))) == true

end
@testset "eager" begin

    m = rand(1:99, 3,4)
    @test transmutedims(m, (1,0,2)) == reshape(m, 3,1,4)
    @test transmutedims(m, (2,99,1)) == reshape(transpose(m), 4,1,3)
    @test transmutedims(m, (2,nothing,1)) == reshape(permutedims(m), 4,1,3)

    @test transmutedims(m, (3,2,1)) == reshape(m', 1,4,3)
    @test transmutedims(transmutedims(m, (3,2,1)), (3,2)) == m

    o314 = rand(Int, 3,1,4);
    @test transmutedims!(o314, m, (1,0,2)) == reshape(m, 3,1,4)
    o413 = rand(Int, 4,1,3);
    @test transmutedims!(o413, m, (2,99,1)) == reshape(permutedims(m), 4,1,3)
    @test o413 == reshape(permutedims(m), 4,1,3)

    # default perm
    @test transmutedims(m) == m'
    @test transmutedims(m[:,1]) == m[:,1]'

    # errors
    @test_throws ArgumentError transmutedims(m, (1,))    # too few
    @test_throws ArgumentError transmutedims(m, (1,0,0)) # 2 doesn't appear

    @test_throws ArgumentError transmutedims!(o314, m, (1,))
    @test_throws Exception transmutedims!(o314, m, (1,0,0))
    @test_throws Exception transmutedims!(o314, m, (1,2,2))

    @test_throws Exception transmutedims!(o314, m, (2,0,1)) # wrong size

end
@testset "lazy" begin

    m = rand(1:99, 3,4)
    @test transmute(m, (1,0,2)) == reshape(m, 3,1,4)
    @test transmute(m, (2,99,1)) == reshape(transpose(m), 4,1,3)
    @test transmute(m, (2,nothing,1)) == reshape(transpose(m), 4,1,3)
    @test transmute(m, Val((1,0,2))) == reshape(m, 3,1,4)
    @test transmute(m, Val((2,99,1))) == reshape(transpose(m), 4,1,3)
    @test transmute(m, Val((2,nothing,1))) == reshape(transpose(m), 4,1,3)
    @test TransmutedDimsArray(m, (1,0,2)) == reshape(m, 3,1,4)
    @test TransmutedDimsArray(m, (2,99,1)) == reshape(transpose(m), 4,1,3)
    @test TransmutedDimsArray(m, (2,nothing,1)) == reshape(transpose(m), 4,1,3)

    v = m[:,1]
    @test transmute(v, (1,1)) isa Diagonal
    @test transmute(v, Val((1,1))) isa Diagonal

    g = TransmutedDimsArray(m, (1,0,2)) # could be an Array
    t = transmute(m, (2,0,1))

    # setindex!
    g[3,1,4] = 222
    @test m[3,4] == 222
    t[2,1,3] = 333
    @test m[3,2] == 333

    # linear indexing
    @test Base.IndexStyle(t) == IndexCartesian()
    g[4] = 444
    @test m[4] == 444
    t[5] = 555
    @test transpose(m)[5] == 555

    # linear indexing, for more constructors
    y = zeros(1,2,3);
    @test IndexStyle(transmute(y, (1,0,2,0,3))) == IndexLinear() # it's an Array
    @test IndexStyle(transmute(y, (3,0,2,0,1))) == IndexCartesian()
    @test IndexStyle(TransmutedDimsArray(y, (3,0,2,0,1))) == IndexCartesian()

    @test transmute(ones(3) .+ im, (1,))[1] == 1 + im # was an ambiguity error

    # reshape
    @test vec(TransmutedDimsArray(m, (1,0,2))) isa Vector
    @test vec(TransmutedDimsArray(m, (2,0,1))) isa Base.ReshapedArray
    @test transmute(reshape(1:6, 2,3), (1,0,2)) isa Base.ReshapedArray

    # errors
    @test_throws ArgumentError transmute(m, (2,))
    @test_throws ArgumentError transmute(m, (2,0,3,0,2))
    @test_throws ArgumentError transmute(m, Val((2,)))
    @test_throws ArgumentError transmute(m, Val((2,0,3,0,2)))
    @test_throws Exception TransmutedDimsArray(m, (2,))
    @test_throws Exception TransmutedDimsArray(m, (2,0,3,0,2))

end
@testset "unwrapping" begin

    m = rand(1:99, 3,4)
    v = m[:,1]

    # unwrapping: LinearAlgebra
    @test transmute(m', (2,0,1)) == transmute(m, (1,0,2)) # both are reshapes
    @test transmute(m', (2,0,1)) isa Array
    @test transmute(m', Val((2,0,1))) isa Array
    @test transmute(m, (1,0,2)) isa Array
    @test transmute(m, Val((1,0,2))) isa Array
    @test transmute(v |> transpose, (2,3,1)) == transmute(v, (1,0,0))
    @test transmute(v |> transpose, (2,3,1)) isa Array
    @test transmute(v |> transpose, Val((2,3,1))) isa Array

    @test transmute((1:3)',(2,2)) === Diagonal(1:3)
    @test_throws ArgumentError transmute((1:3)',(1,1)) # 1 is trivial dim. Confusing message though!
    @test transmute((1:3)',(1,2)) == (1:3)' # but not ===, unwrapped before noticing.

    @test transmute(Diagonal(1:10), (3,1)) === TransmutedDimsArray(1:10, (0,1))
    @test transmute(Diagonal(rand(10)), (3,1)) isa Matrix
    @test transmute(Diagonal(rand(10)), Val((3,1))) isa Matrix
    @test transmute(TransmutedDimsArray(1:10, (1,1)), (3,1)) === TransmutedDimsArray(1:10, (0,1))
    @test transmute(TransmutedDimsArray(1:10, (1,1)), Val((3,1))) === TransmutedDimsArray(1:10, (0,1))

    # unwrapping: permutations
    @test transmute(PermutedDimsArray(m,(2,1)), (2,nothing,1)) == transmute(m, (1,0,2))
    @test transmute(PermutedDimsArray(m,(2,1)), (2,nothing,1)) isa Array
    @test transmute(PermutedDimsArray(m,(2,1)), Val((2,nothing,1))) isa Array
    @test transmute(TransmutedDimsArray(m,(2,1)), (2,nothing,1)) isa Array
    @test transmute(TransmutedDimsArray(m,(2,1)), Val((2,nothing,1))) isa Array

    @test transmute(PermutedDimsArray(m,(2,1)), (1,3,2)) === transmute(m, (2,3,1))

    x = rand(1:99, 5,4,3,2);
    x1 = permutedims(x, (4,1,2,3))
    x2 = PermutedDimsArray(x, (4,1,2,3))
    @test x1 == x2
    @test TransmutedDimsArray(x1, (2,4,1,3)) == TransmutedDimsArray(x2, (2,4,1,3))
    @test TransmutedDimsArray(x1, (2,0,1,4,3)) == TransmutedDimsArray(x2, (2,0,1,4,3))

    x3 = TransmutedDimsArray(x, (4,1,2,3));
    @test TransmutedDimsArray(x1, (2,4,1,3)) == TransmutedDimsArray(x3, (2,4,1,3))
    @test TransmutedDimsArray(x1, (2,0,1,4,3)) == TransmutedDimsArray(x3, (2,0,1,4,3))

    # unwrapping caused by transpose etc
    @test transpose(TransmutedDimsArray(m, (2,1))) isa Matrix
    @test adjoint(TransmutedDimsArray(v, (0,1))) isa Matrix

end
@testset "view" begin

    m = rand(1:99, 3,4)
    g = TransmutedDimsArray(m, (1,0,2)) # could be an Array
    t = transmute(m, (2,0,1))

    @test view(g, :,1,2) isa SubArray{Int64, 1, <:Matrix}
    @test view(g, :,1,2) == g[:,1,2]

    @test view(t, :,1,2) isa SubArray{Int64, 1, <:Matrix}
    @test view(t, :,1,2) == t[:,1,2]
    @test view(t, 2,1,:) == t[2,1,:]

    @test first(eachcol(TransmutedDimsArray(m,(2,1)))) isa SubArray{Int64, 1, <:Matrix}

    # not unwrapped
    @test view(t, :,1,:) == m'
    @test view(t, 2,:,3) == t[2,:,3]

end
@testset "reductions" begin

    m = rand(1:99, 3,4)

    # reductions
    @test sum(transmute(m,(2,1,3))) == sum(m)
    @test sum(transmute(m,(1,1,2))) == sum(m)
    @test prod(transmute(m,(1,1,2))) == 0
    @test sum(x->x+1, transmute(m,(1,1,2))) == sum(m) + 3*3*4

    @test sum(TransmutedDimsArray(m,(1,2)); dims=1) == sum(m, dims=1)
    @test sum(TransmutedDimsArray(m,(2,1)); dims=1) == sum(m', dims=1)
    @test sum(TransmutedDimsArray(m,(2,1)), dims=2) == sum(m', dims=2)
    @test sum(x->x+1, TransmutedDimsArray(m,(2,1)), dims=2) == sum(x->x+1, m', dims=2)

    @test sum(transmute(m,(1,1,2)), dims=1) == sum(collect(transmute(m,(1,1,2))), dims=1)
    @test sum(transmute(m,(1,1,1,2)), dims=1) == sum(collect(transmute(m,(1,1,1,2))), dims=1)
    @test sum(transmute(m,(1,1,2)), dims=(1,2)) == sum(collect(transmute(m,(1,1,2))), dims=(1,2))
    @test prod(transmute(m,(1,1,2)), dims=(1,2)) == prod(collect(transmute(m,(1,1,2))), dims=(1,2))

    g = TransmutedDimsArray(m, (1,0,2)) # could be an Array
    t = transmute(m, (2,0,1))

    @test sum(g, dims=2) == g
    @test sum(t, dims=2) == t
    @test sum(x->x+1, t, dims=2) == sum(x->x+1, collect(t), dims=2)

    x = rand(1:99, 5,4,3,2);
    x1 = permutedims(x, (4,1,2,3))
    x3 = TransmutedDimsArray(x, (4,1,2,3))
    @test sum(x1, dims=3) == sum(x3, dims=3)
    @test sum(x1, dims=(2,4)) == sum(x3, dims=(2,4))
    @test sum(x->x+1, x1, dims=(2,4)) == sum(x->x+1, x3, dims=(2,4))

end
@testset "dropdims" begin

    m = rand(1:99, 3,4)
    g = TransmutedDimsArray(m, (1,0,2)) # could be an Array
    r = collect(g)
    t = transmute(m, (2,0,1))

    # dropdims
    @test dropdims(g, dims=2) == m
    @test dropdims(g, dims=2) isa Matrix # now unwraps!
    @test dropdims(t, dims=2) == m'
    @test dropdims(t, dims=2) isa Transpose

    h = reshape(1:9, 3,1,3)
    f = transmute(h, (0,3,2,1))

    @test size(f) == (1,3,1,3)
    @test dropdims(f, dims=1) == permutedims(h, (3,2,1))
    @test dropdims(f, dims=3).parent == reshape(h,3,3)
    @test dropdims(f, dims=(1,3)) == reshape(h,3,3)'
    @test dropdims(f, dims=(3,1)) == reshape(h,3,3)'

    # without dropdims
    @test transmute(g, (1,3)) === m # unwraps
    @test transmute(r, (1,3)) == m # reshapes
    @test transmute(r, (1,3)) isa Matrix
    @test transmute(t, (3,1)) === m # unwraps
    @test transmute(t, (1,3)) === transpose(m)

    @test transmute(g, Val((1,3))) === m
    @test transmute(r, Val((1,3))) == m
    @test transmute(r, Val((1,3))) isa Matrix
    @test transmute(t, Val((3,1))) === m
    @test transmute(t, Val((1,3))) === transpose(m)

    @test transmutedims(g, (1,3)) === m
    @test transmutedims(r, (1,3)) == m
    @test transmutedims(r, (1,3)) isa Matrix
    @test transmutedims(t, (3,1)) === m
    @test transmutedims(t, (1,3)) == transpose(m)

    # errors
    @test_throws ArgumentError transmute(g, (1,2))
    @test_throws ArgumentError transmute(g, Val((1,2)))
    @test_throws ArgumentError transmutedims(g, (1,2))

end
@testset "diagonal" begin

    for d in [
        transmute(ones(Int,3), (1,1)),
        transmute(ones(Int,3), Val((1,1))),
        TransmutedDimsArray(ones(Int,3), (1,1)),
        ]
        @test d == [1 0 0; 0 1 0; 0 0 1]
        @test d[2] == 0
        @test (d[3,3] = 33) == 33
        @test d[3,3] == 33
        @test (d[2,1] = 0) == 0
        @test_throws ArgumentError d[1,2] = 99
        @test IndexStyle(d) == IndexCartesian()

        @test sum(d .+ 10) == 90 + 2 + 33
    end

    r = rand(3)
    q = TransmutedDimsArray(r, (1,0,1,1))
    @test q[2,1,2,2] == r[2]
    @test q[2,1,2,3] == 0
    @test_throws ArgumentError q[1,1,1,3] = 99

    # eager
    @test q == transmutedims(r, (1,0,1,1))

    # dropdims
    @test transmute(rand(3,1), (1,1)) isa Diagonal
    @test transmute(rand(1,3), (2,2)) isa Diagonal
    @test_throws Exception transmute(rand(3,2), (1,1))
    @test transmute(rand(3,1), Val((1,1))) isa TransmutedDimsArray{Float64, 2, (1,1), (2,0)}

end
@testset "zero dims" begin

    @test transmute(fill(3), (1,2)) isa Matrix
    @test transmute(fill(3), ()) isa Array{Int,0}
    @test transmute([3], ()) isa Array{Int,0}

    @test transmute(fill(3), Val((1,2))) isa Matrix
    @test transmute(fill(3), Val(())) isa Array{Int,0}
    @test transmute([3], Val(())) isa Array{Int,0}

    @test TransmutedDimsArray(fill(3), (1,2)) isa AbstractMatrix
    @test TransmutedDimsArray(fill(3), ()) isa AbstractArray{Int,0}
    @test TransmutedDimsArray([3], ()) isa AbstractArray{Int,0}

    @test transmutedims(fill(3), (1,2)) isa Matrix
    @test transmutedims(fill(3), ()) isa Array{Int,0}
    @test transmutedims([3], ()) isa Array{Int,0}

end
@testset "tuples" begin

    @test transmute((1,2,3), (2,1)) == [1 2 3]
    @test transmute((1,2,3), (1,1)) == Diagonal(1:3)

    @test transmute((1,2,3), (1,)) isa AbstractVector  # always wraps
    @test transmute((1,2,3), (1,nothing)) isa AbstractMatrix
    @test transmute(transmute((1,2,3), (2,1)), (2,)) isa TupleVector  # unwraps once

    # eager
    @test transmutedims((1,2,3), (2,1)) == [1 2 3]
    @test transmutedims((1,2,3)) == [1 2 3]
    @test transmutedims((1,2,3), (1,1)) isa Matrix

    # promotion
    @test eltype(transmute((1,2,3.0),(2,1)) .+ [1 2 3]) == eltype((1,2,3.0) .+ [1,2,3]) # Real
    @test eltype(transmute((1,2,3+im),(2,1)) .+ [1 2 3]) == eltype((1,2,3+im) .+ [1,2,3]) # Number
    @test transmutedims((1,2,3.0)) isa Matrix{Real}

end
@testset "from Base" begin

    # keeps the num of dim
    p = randperm(5)
    q = randperm(5)
    a = rand(p...)
    b = transmutedims(a,q)
    @test isequal(size(b), tuple(p[q]...))

    # hand made case
    y = zeros(1,2,3)
    for i = 1:6
        y[i]=i
    end

    z = zeros(3,1,2)
    for i = 1:3
        z[i] = i*2-1
        z[i+3] = i*2
    end

    # permutes correctly
    @test isequal(z,transmutedims(y,[3,1,2]))
    @test isequal(z,transmutedims(y,(3,1,2)))

    # of a subarray
    a = rand(5,5)
    s = view(a,2:3,2:3)
    p = transmutedims(s, [2,1])
    @test p[1,1]==a[2,2] && p[1,2]==a[3,2]
    @test p[2,1]==a[2,3] && p[2,2]==a[3,3]

    # of a non-strided subarray
    a = reshape(1:60, 3, 4, 5)
    s = view(a,:,[1,2,4],[1,5])
    c = convert(Array, s)
    for p in ([1,2,3], [1,3,2], [2,1,3], [2,3,1], [3,1,2], [3,2,1])
        @test transmutedims(s, p) == transmutedims(c, p)
        @test TransmutedDimsArray(s, p) == transmutedims(c, p)
    end
    @test_throws ArgumentError transmutedims(a, (1,1,1))
    @test_throws ArgumentError transmutedims(s, (1,1,1))
    @test_throws ArgumentError TransmutedDimsArray(a, (1,1,1))
    @test_throws ArgumentError TransmutedDimsArray(s, (1,1,1))
    cp = TransmutedDimsArray(c, (3,2,1))
    @test pointer(cp) == pointer(c)
    @test_throws ArgumentError pointer(cp, 2)
    @test strides(cp) == (9,3,1)
    ap = TransmutedDimsArray(Array(a), (2,1,3))
    @test strides(ap) == (3,1,12)

    for A in [rand(1,2,3,4),rand(2,2,2,2),rand(5,6,5,6),rand(1,1,1,1)]
        perm = randperm(4)
        @test isequal(A,transmutedims(transmutedims(A,perm),invperm(perm)))
        @test isequal(A,transmutedims(transmutedims(A,invperm(perm)),perm))
    end

    # m = [1 2; 3 4]
    # @test transmutedims(m) == [1 3; 2 4]

    # v = [1,2,3]
    # @test transmutedims(v) == [1 2 3]
end

using TransmuteDims: _float_transmutedims  # exists to test permutedims! path
@testset "sans Strided" begin

    m = rand(1:99, 3,4)
    @test _float_transmutedims(m, (1,0,2)) == reshape(m, 3,1,4)
    @test _float_transmutedims(m, (2,99,1)) == reshape(transpose(m), 4,1,3)
    @test _float_transmutedims(m, (2,nothing,1)) == reshape(permutedims(m), 4,1,3)

    g = transmutedims(m, (1,0,2)) # an Array

    @test _float_transmutedims(g, (3,1)) == m'
    @test _float_transmutedims(g, (3,0,1)) == transmute(m, (2,0,1))

end

using BenchmarkTools
@testset "allocations" begin

    y = ones(1,2,3);

    # wrap
    @test 97 > @ballocated transmute($y, (3,2,0,1)) # zero on 1.7, maybe 1.5?
    @test 17 > @ballocated transmute($y, Val((3,2,0,1)))  # nonzero only on 1.3
    # @code_warntype (y -> transmute(y, (3,2,0,1)))(y)

    # reshape
    @test 129 > @ballocated transmute($y, (1,2,0,3))
    @test 129 > @ballocated transmute($y, Val((1,2,0,3)))
    # @code_warntype (y -> transmute(y, (1,2,0,3)))(y)

end

using OffsetArrays, Random
@testset "offset" begin

    o34 = ones(11:13, 21:24)
    rand!(o34)

    @test transmute(o34, (2,1)) == permutedims(o34)
    @test axes(transmute(o34, (2,0,1))) == (21:24, 1:1, 11:13)

    @test transmute(o34, (2,0,1))[21, 1, 13] == o34[13, 21]

    # dropdims
    @test axes(transmute(sum(o34, dims = 1), (2,))) == (21:24,)
    @test transmute(sum(o34, dims = 1), (2,)) isa OffsetArray # issue 23
    o34 .= 1
    @test all(==(3), transmute(sum(o34, dims = 1), (2,))) # issue 24

end

using Strided
@testset "Strided" begin

    m = rand(1:99, 3,4)
    @test @strided(transmute(m, (0,2,1)) .+ transmute(m, (1,2,0))) isa StridedView

end

using StaticArrays
@testset "StaticArrays" begin

    m = SMatrix{3,4}(rand(1:99, 3,4))
    v = SVector{4}(rand(4))

    permutedims(v) isa Base.ReshapedArray
    @test transmute(v, (2,1)) isa TransmutedDimsArray  # no worse than Base

    permutedims(m) isa MMatrix
    @test transmute(m, (2,1)) isa SMatrix
    @test transmutedims(m, (2,1)) isa MMatrix  # this is permutedims?

    transmute(m, (2,3,1)) isa TransmutedDimsArray
end

if VERSION < v"1.5"
    @warn "skipping tests of GPUArrays"
else
    using GPUArrays, Adapt
    GPUArrays.allowscalar(false)

    jl_file = normpath(joinpath(pathof(GPUArrays), "..", "..", "test", "jlarray.jl"))
    include(jl_file)
    using .JLArrays # a fake GPU array, for testing

    @testset "GPUArrays" begin

        m = rand(4,4)
        jm = JLArray(m)
        @test_throws Exception jm[1]  # scalar getindex is disallowed

        tjm = TransmutedDimsArray(jm, (2,1))  # simple wrapper
        j2 = jm .* log.(tjm) ./ 2
        @test j2 isa JLArray
        @test collect(j2) ≈ m .* log.(m') ./ 2

        jmd = transmute(jm, (1,1,2))  # Diagonal-like
        j3 = jmd .+ 2 .* jmd .+ 1
        @test j3 isa JLArray
        @test collect(j3) ≈ 3 .* transmute(m, (1,1,2)) .+ 1

        # printing, without scalar indexing
        @test sprint(show, jm) isa String
        io = IOBuffer()
        Base.print_array(io, jm)
        @test contains(String(take!(io)), string(m[1]))

        @test jm isa DenseArray  # but reshaping it does not preserve this
        TransmuteDims.may_reshape(::Type{<:JLArray}) = false

        # eager
        @test transmutedims(jm, (0,1,2)) isa JLArray
        @test transmutedims(jm, (2,1)) isa JLArray
        @test_broken transmutedims(jm, (2,0,1)) isa JLArray # permutedims!(dest::Base.ReshapedArray{Float64, 2, JLArray{Float64, 3}, Tuple{}}, src::JLArray{Float64, 2}, perm::Tuple{Int64, Int64})

    end
end

using Tracker
@testset "Tracker: $func" for func in [transmute, transmutedims]
    # subset of the Zygote tests

    @test size(Tracker.gradient(x -> sum(sin, func(x, (2,1))), rand(2,3))[1]) == (2,3)
    @test size(Tracker.gradient(x -> sum(sin, func(x, (2,3,1))), rand(2,3))[1]) == (2,3)
    @test size(Tracker.gradient(x -> sum(sin, func(x, (2,1,1))), rand(2,3))[1]) == (2,3)

    v, m, t = rand(1:99, 3), rand(1:99, 3,3), rand(1:99, 3,3,3)

    fwd, back = Tracker.forward(x -> func(x, (2,1)), v)
    @test fwd == v'
    @test back(ones(1,3))[1] == ones(3)
    @test back(v')[1] == v  # awkward reshape(adjoint)

end

print("... loading Zygote, 🐌 ... ")
using Zygote
println("done")
@testset "Zygote: $func" for func in [transmute, transmutedims]
    NEW = VERSION >= v"1.6-"

    # sizes, and no errors!
    @test size(Zygote.gradient(x -> sum(sin, func(x, (2,1))), rand(2,3))[1]) == (2,3)
    @test size(Zygote.gradient(x -> sum(sin, func(x, (2,3,1))), rand(2,3))[1]) == (2,3)
    @test size(Zygote.gradient(x -> sum(sin, func(x, (2,1,1))), rand(2,3))[1]) == (2,3)
    @test_skip  NEW && func==transmute && @test size(Zygote.gradient(x -> sum(sin, func(x, (2,2,3,1,1))), rand(2,3))[1]) == (2,3)

    @test size(Zygote.gradient(x -> sum(sin, func(x, (1,))), rand(3,1))[1]) == (3,1)
    @test size(Zygote.gradient(x -> sum(sin, func(x, (2,1))), rand(3,1))[1]) == (3,1)
    @test size(Zygote.gradient(x -> sum(sin, func(x, (1,1))), rand(3,1))[1]) == (3,1)
    @test_skip  NEW && func==transmute && @test size(Zygote.gradient(x -> sum(sin, func(x, (1,3,1,3))), rand(3,1))[1]) == (3,1)

    @test size(Zygote.gradient(x -> sum(sin, func(x, (2,3,1))), rand(2,3,4))[1]) == (2,3,4)
    @test size(Zygote.gradient(x -> sum(sin, func(x, (2,4,3,1))), rand(2,3,4))[1]) == (2,3,4)

    # values
    v, m, t = rand(1:99, 3), rand(1:99, 3,3), rand(1:99, 3,3,3)

    fwd, back = Zygote.pullback(x -> func(x, (2,1)), v)
    @test fwd == v'
    @test back(ones(1,3))[1] == ones(3)
    @test back(v')[1] == v  # awkward reshape(adjoint)

    # extracting diagonals
    fwd, back = pullback(x -> func(x, (1,1)), v)
    @test fwd == Diagonal(v)
    @test back(m)[1] == diag(m)
    @test_skip back(m[:,:,:,:])[1] == diag(m)  # trivial extra dimensions, different path

    fwd, back = pullback(x -> func(x, (1,1,2)), m)
    @test fwd[:,:,1] == Diagonal(m[:,1])
    @test back(t)[1] == [t[i,i,k] for i in 1:3, k in 1:3]

end

if VERSION < v"1.6-" || VERSION > v"1.7-"
    @warn "skipping doctests, on Julia $VERSION"  # as these have random numbers
else
    println("... starting Documenter 🗣")
    using Documenter
    @testset "doctests" begin

        doctest(TransmuteDims, manual=false)

    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "BayesOptBuilder"
version = v"0.1.0"

# Collection of sources required to build BayesOptBuilder
sources = [
    "https://github.com/rmcantin/bayesopt.git" =>
    "1326549b68f27801b79984ca227f4183d984107d",

    "https://dl.bintray.com/boostorg/release/1.68.0/source/boost_1_68_0.tar.bz2" =>
    "7f6130bc3cf65f56a618888ce9d5ea704fa10b462be126ad053e80e553d6d8b7",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DBAYESOPT_BUILD_EXAMPLES=OFF -DBAYESOPT_BUILD_SHARED=ON -DNLOPT_BUILD_SHARED=ON bayesopt/ -DBOOST_ROOT=../boost_1_68_0
make -j$nprocs
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libbayesopt", :libbayesopt)
    LibraryProduct(prefix, "libnlopt", :libnlopt)
]

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

]

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


text:
# top type - unifying a type with top yeilds the type
struct Top end

# get the type from a named tuple, given a name
function get_type(t::Type{NamedTuple{N,T}}, k::Symbol) where {N,T<:Tuple}
    if k in N
        return T.types[findfirst(isequal(k), N)]
    else
        return Nothing
    end
end

# unify two types to a single type
unify(a, b) = unify(b, a)
unify(a::Type{T}, b::Type{S}) where {T,S} = Base.promote_typejoin(T, S)
unify(a::Type{T}, b::Type{S}) where {T,S<:T} = T
unify(a::Type{Top}, b::Type{T}) where {T} = T

function unify(
    a::Type{NamedTuple{A,T}},
    b::Type{NamedTuple{B,S}},
) where {A,T<:Tuple,B,S<:Tuple}
    c = Dict()
    for (k, v) in zip(A, T.types)
        c[k] = unify(v, get_type(b, k))
    end

    for (k, v) in zip(B, S.types)
        if k ∉ keys(c)
            c[k] = unify(v, Nothing)
        end
    end

    return NamedTuple{tuple(keys(c)...),Tuple{values(c)...}}
end

function unify(a::Type{Array{T,1}}, b::Type{Array{S,1}}) where {T,S}
    return Array{unify(T, S),1}
end

# parse json into a type
function build_type(o::JSON3.Object)
    d = Dict()
    for (k, v) in o
        d[k] = build_type(v)
    end

    return NamedTuple{tuple(keys(d)...),Tuple{values(d)...}}
end

function build_type(a::JSON3.Array)
    t = Set([])
    nt = Top
    for item in a
        it = build_type(item)
        if it <: NamedTuple
            nt = unify(nt, it)
        else
            push!(t, build_type(item))
        end
    end

    return Array{foldl(unify, t; init = Union{nt}), 1}
end

build_type(x::T) where {T} = T


text:
@doc raw"""
Midpoint extrapolation method with arbitrary order p.

Solves the ordinary differential equation
```math
\begin{aligned}
\dot{x} &= v(t, x) , &
x(t_0) &= x_0 ,
\end{aligned}
```
for $x_1 = x(t_1)$.

Call with
```julia
_midpoint_extrapolation_ode!(v, t₀, t₁, x₀, x₁, s)
```

where

* `v`:  function to compute vector field with signature `v(t,x,ẋ)`
* `t₀`: initial time
* `t₁`: final   time
* `x₀`: initial value $x_0 = x(t_0)$
* `x₁`: final   value $x_1 = x(t_1)$
* `s`:  number of interpolations (order $p=2s+2$)

"""
function _midpoint_extrapolation_ode!(v::Function, t₀::TT, t₁::TT,
                                      x₀::AbstractVector{DT}, x₁::AbstractVector{DT}, s::Int) where {DT,TT}
    @assert size(x₀) == size(x₁)

    local F   = [2i*one(TT) for i in 1:(s+1)]
    local Δt  = t₁ - t₀
    local σ   = Δt ./ F
    local σ²  = σ.^2
    local pts = zeros(DT, length(x₀), s+1)

    local xᵢ₁ = zero(x₀)
    local xᵢ₂ = zero(x₀)
    local xᵢₜ = zero(x₀)
    local vᵢ  = zero(x₀)
    local v₀  = zero(x₀)

    v(t₀, x₀, v₀)

    for i in 1:s+1
        tᵢ   = t₀ + σ[i]
        xᵢ₁ .= x₀
        xᵢ₂ .= x₀ .+ σ[i] .* v₀
        for _ in 1:(F[i]-1)
            v(tᵢ, xᵢ₂, vᵢ)
            xᵢₜ .= xᵢ₁ .+ 2σ[i] .* vᵢ
            xᵢ₁ .= xᵢ₂
            xᵢ₂ .= xᵢₜ
        end
        for k in axes(pts,1)
            pts[k,i] += xᵢ₂[k]
        end
    end

    aitken_neville!(σ², pts, zero(TT), x₁)
    return x₁
end

struct MidpointExtrapolationODE{VT} <: Extrapolation
    s::Int
    v::VT

    function MidpointExtrapolationODE{VT}(v, s) where {VT}
        new(s,v)
    end
end

function MidpointExtrapolationODE(v::VT, s::Int) where {VT}
    MidpointExtrapolationODE{VT}(v, s)
end

function MidpointExtrapolation(equ::ODE, s::Int)
    MidpointExtrapolationODE(_get_v(equ), s)
end

function GeometricBase.evaluate!(extrap::MidpointExtrapolationODE, t₀::TT, t₁::TT,
                          x₀::AbstractVector{DT}, x₁::AbstractVector{DT}) where {DT,TT}
    _midpoint_extrapolation_ode!(extrap.v, t₀, t₁, x₀, x₁, extrap.s)
end

@doc raw"""
Midpoint extrapolation method with arbitrary order p.

Solves the implicit differential equation
```math
\begin{aligned}
\dot{q} &= v(t, q, p) , &
q(t_0) &= q_0 , \\
\dot{p} &= f(t, q, v) , &
p(t_0) &= p_0 ,
\end{aligned}
```
for $q_1 = q(t_1)$ and $p_1 = p(t_1)$.

Call with
```julia
_midpoint_extrapolation_iode!(v, t₀, t₁, q₀, q₁, p₀, p₁, s)
```

where

* `v`:  function to compute vector field with signature `v(t,x,p,ẋ)`
* `f`:  function to compute force  field with signature `f(t,x,v,ṗ)`
* `t₀`: initial time
* `t₁`: final   time
* `q₀`: initial position $q_0 = q(t_0)$
* `q₁`: final   position $q_1 = q(t_1)$
* `p₀`: initial momentum $p_0 = p(t_0)$
* `p₁`: final   momentum $p_1 = p(t_1)$
* `s`:  number of interpolations (order $p=2s+2$)

"""
function _midpoint_extrapolation_iode!(v::Function, f::Function, t₀::TT, t₁::TT,
                                       q₀::AbstractVector{DT}, q₁::AbstractVector{DT},
                                       p₀::AbstractVector{DT}, p₁::AbstractVector{DT}, s::Int) where {DT,TT}
    @assert size(q₀) == size(q₁) == size(p₀) == size(p₁)

    local F   = [2i*one(TT) for i in 1:(s+1)]
    local Δt  = t₁ - t₀
    local σ   = Δt ./ F
    local σ2  = σ.^2

    local qts = zeros(DT, length(q₀), s+1)
    local pts = zeros(DT, length(p₀), s+1)

    local qᵢ₁= zero(q₀)
    local qᵢ₂= zero(q₀)
    local qᵢₜ= zero(q₀)

    local pᵢ₁= zero(p₀)
    local pᵢ₂= zero(p₀)
    local pᵢₜ= zero(p₀)

    local v₀ = zero(q₀)
    local vᵢ = zero(q₀)

    local f₀ = zero(p₀)
    local fᵢ = zero(p₀)

    v(t₀, q₀, v₀)
    f(t₀, q₀, v₀, f₀)

    for i in 1:(s+1)
        tᵢ   = t₀ + σ[i]
        qᵢ₁ .= q₀
        qᵢ₂ .= q₀ .+ σ[i] .* v₀
        pᵢ₁ .= p₀
        pᵢ₂ .= p₀ .+ σ[i] .* f₀
        for _ in 1:(F[i]-1)
            v(tᵢ, qᵢ₂, vᵢ)
            f(tᵢ, qᵢ₂, vᵢ,  fᵢ)
            qᵢₜ .= qᵢ₁ .+ 2σ[i] .* vᵢ
            qᵢ₁ .= qᵢ₂
            qᵢ₂ .= qᵢₜ
            pᵢₜ .= pᵢ₁ .+ 2σ[i] .* fᵢ
            pᵢ₁ .= pᵢ₂
            pᵢ₂ .= pᵢₜ
        end
        for k in axes(qts,1)
            qts[k,i] += qᵢ₂[k]
        end
        for k in axes(pts,1)
            pts[k,i] += pᵢ₂[k]
        end
    end

    aitken_neville!(σ2, qts, zero(TT), q₁)
    aitken_neville!(σ2, pts, zero(TT), p₁)
    return (q₁, p₁)
end

struct MidpointExtrapolationIODE{VT,FT} <: Extrapolation
    s::Int
    v::VT
    f::FT

    function MidpointExtrapolationIODE{VT,FT}(v, f, s) where {VT,FT}
        new(s,v,f)
    end
end

function MidpointExtrapolationIODE(v::VT, f::FT, s::Int) where {VT,FT}
    MidpointExtrapolationIODE{VT,FT}(v, f, s)
end

function MidpointExtrapolation(equ::IODE, s::Int)
    MidpointExtrapolationIODE(_get_v̄(equ), _get_f̄(equ), s)
end

function GeometricBase.evaluate!(extrap::MidpointExtrapolationIODE, t₀::TT, t₁::TT,
                          q₀::AbstractVector{DT}, q₁::AbstractVector{DT},
                          p₀::AbstractVector{DT}, p₁::AbstractVector{DT}) where {DT,TT}
    _midpoint_extrapolation_iode!(extrap.v, extrap.f, t₀, t₁, q₀, q₁, p₀, p₁, extrap.s)
end

@doc raw"""
Midpoint extrapolation method with arbitrary order p.

Solves the partitioned ordinary differential equation
```math
\begin{aligned}
\dot{q} &= v(t, q, p) , &
q(t_0) &= q_0 , \\
\dot{p} &= f(t, q, p) , &
p(t_0) &= p_0 ,
\end{aligned}
```
for $q_1 = q(t_1)$ and $p_1 = p(t_1)$.

Call with
```julia
_midpoint_extrapolation_pode!(v, t₀, t₁, q₀, q₁, p₀, p₁, s)
```

where

* `v`:  function to compute velocity field with signature `v(t,x,p,ẋ)`
* `f`:  function to compute force    field with signature `f(t,x,p,ṗ)`
* `t₀`: initial time
* `t₁`: final   time
* `q₀`: initial position $q_0 = q(t_0)$
* `q₁`: final   position $q_1 = q(t_1)$
* `p₀`: initial momentum $p_0 = p(t_0)$
* `p₁`: final   momentum $p_1 = p(t_1)$
* `s`:  number of interpolations (order $p=2s+2$)

"""
function _midpoint_extrapolation_pode!(v::Function, f::Function, t₀::TT, t₁::TT,
                                       q₀::AbstractVector{DT}, q₁::AbstractVector{DT},
                                       p₀::AbstractVector{DT}, p₁::AbstractVector{DT}, s::Int) where {DT,TT}
    @assert size(q₀) == size(q₁) == size(p₀) == size(p₁)

    local F   = [2i*one(TT) for i in 1:(s+1)]
    local Δt  = t₁ - t₀
    local σ   = Δt ./ F
    local σ2  = σ.^2

    local qts = zeros(DT, length(q₀), s+1)
    local pts = zeros(DT, length(p₀), s+1)

    local qᵢ₁= zero(q₀)
    local qᵢ₂= zero(q₀)
    local qᵢₜ= zero(q₀)

    local pᵢ₁= zero(p₀)
    local pᵢ₂= zero(p₀)
    local pᵢₜ= zero(p₀)

    local v₀ = zero(q₀)
    local vᵢ = zero(q₀)

    local f₀ = zero(p₀)
    local fᵢ = zero(p₀)

    v(t₀, q₀, p₀, v₀)
    f(t₀, q₀, p₀, f₀)

    for i in 1:(s+1)
        tᵢ   = t₀ + σ[i]
        qᵢ₁ .= q₀
        qᵢ₂ .= q₀ .+ σ[i] .* v₀
        pᵢ₁ .= p₀
        pᵢ₂ .= p₀ .+ σ[i] .* f₀
        for _ in 1:(F[i]-1)
            v(tᵢ, qᵢ₂, pᵢ₂, vᵢ)
            f(tᵢ, qᵢ₂, pᵢ₂, fᵢ)
            qᵢₜ .= qᵢ₁ .+ 2σ[i] .* vᵢ
            qᵢ₁ .= qᵢ₂
            qᵢ₂ .= qᵢₜ
            pᵢₜ .= pᵢ₁ .+ 2σ[i] .* fᵢ
            pᵢ₁ .= pᵢ₂
            pᵢ₂ .= pᵢₜ
        end
        for k in axes(qts,1)
            qts[k,i] += qᵢ₂[k]
        end
        for k in axes(pts,1)
            pts[k,i] += pᵢ₂[k]
        end
    end

    aitken_neville!(σ2, qts, zero(TT), q₁)
    aitken_neville!(σ2, pts, zero(TT), p₁)
    return (q₁, p₁)
end

struct MidpointExtrapolationPODE{VT,FT} <: Extrapolation
    s::Int
    v::VT
    f::FT

    function MidpointExtrapolationPODE{VT,FT}(v, f, s) where {VT,FT}
        new(s,v,f)
    end
end

function MidpointExtrapolationPODE(v::VT, f::FT, s::Int) where {VT,FT}
    MidpointExtrapolationPODE{VT,FT}(v, f, s)
end

function MidpointExtrapolation(equ::PODE, s::Int)
    MidpointExtrapolationPODE(_get_v(equ), _get_f(equ), s)
end

function GeometricBase.evaluate!(extrap::MidpointExtrapolationPODE, t₀::TT, t₁::TT,
                          q₀::AbstractVector{DT}, q₁::AbstractVector{DT},
                          p₀::AbstractVector{DT}, p₁::AbstractVector{DT}) where {DT,TT}
    _midpoint_extrapolation_pode!(extrap.v, extrap.f, t₀, t₁, q₀, q₁, p₀, p₁, extrap.s)
end


text:
using Unitful
t = 1.0u"s"

t2 = 1.02u"s"
t+t2

t*t2

sqrt(t)

t + sqrt(t)

using DifferentialEquations
f = (y,p,t) -> 0.5*y
u0 = 1.5u"N"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))
sol = solve(prob,Tsit5())

f = (y,p,t) -> 0.5*y/3.0u"s"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))
sol = solve(prob,Tsit5())

print(sol[:])

using Plots
gr()
plot(ustrip(sol.t),ustrip(sol[:]),lw=3)

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
"""
    AbstractVectorTransportOperator

Abstract type for vector transport operators in the frozen coefficients formulation.
"""
abstract type AbstractVectorTransportOperator end

struct DefaultVectorTransportOperator{TVT<:AbstractVectorTransportMethod} <:
       AbstractVectorTransportOperator
    vtm::TVT
end

"""
    (vto::DefaultVectorTransportOperator)(M::AbstractManifold, p, X, q, params, t_from, t_to)

In the frozen coefficient formulation, transport tangent vector `X` such that
`X = f(p, params, t_from)` to point `q` at time `t_to`. This provides a sort of estimation
of `f(q, params, t_to)`.
"""
function (vto::DefaultVectorTransportOperator)(
    M::AbstractManifold,
    p,
    X,
    q,
    params,
    t_from,
    t_to,
)
    # default implementation, may be customized as needed.
    return vector_transport_to(M, p, X, q)
end

"""
    FrozenManifoldDiffEqOperator{T<:Number,TM<:AbstractManifold,TF,TVT} <: SciMLBase.AbstractDiffEqOperator{T}

DiffEq operator on manifolds in the frozen vector field formulation.
"""
struct FrozenManifoldDiffEqOperator{T<:Number,TF,TVT<:AbstractVectorTransportOperator} <:
       SciMLBase.AbstractDiffEqOperator{T}
    func::TF
    operator_vector_transport::TVT
end

function FrozenManifoldDiffEqOperator{T}(f, ovt) where {T<:Number}
    return FrozenManifoldDiffEqOperator{T,typeof(f),typeof(ovt)}(f, ovt)
end
function FrozenManifoldDiffEqOperator{T}(f) where {T<:Number}
    return FrozenManifoldDiffEqOperator{T}(
        f,
        DefaultVectorTransportOperator(ParallelTransport()),
    )
end

function (L::FrozenManifoldDiffEqOperator)(du, u, p, t)
    return copyto!(du, L.func(u, p, t))
end
function (L::FrozenManifoldDiffEqOperator)(u, p, t)
    return L.func(u, p, t)
end

"""
    LieManifoldDiffEqOperator{T<:Number,TF} <: AbstractDiffEqOperator{T}

DiffEq operator on manifolds in the Lie group action formulation.
"""
struct LieManifoldDiffEqOperator{T<:Number,TF} <: SciMLBase.AbstractDiffEqOperator{T}
    func::TF
end

function LieManifoldDiffEqOperator{T}(f) where {T<:Number}
    return LieManifoldDiffEqOperator{T,typeof(f)}(f)
end

function (L::LieManifoldDiffEqOperator)(du, u, p, t)
    return copyto!(du, L.func(u, p, t))
end
function (L::LieManifoldDiffEqOperator)(u, p, t)
    return L.func(u, p, t)
end


text:
@testset "354.russian-doll-envelopes.jl" begin
    @test max_envelopes([[5,4],[6,4],[6,7],[2,3]]) == 3
end


text:
@testset "d_separated" begin
    # Examples from <https://doi.org/10.1073/PNAS.1510507113>.
    nodes = [:Z1, :W1, :X, :Z3, :W3, :Y, :Z2, :W2]
    mapping = Dict(zip(nodes, 1:length(nodes)))
    m = mapping
    E = [
        :Z1 => :W1,
        :W1 => :X,
        :Z1 => :Z3,
        :Z3 => :X,
        :Z2 => :Z3,
        :Z2 => :W2,
        :Z3 => :Y,
        :W2 => :Y,
        :W3 => :Y,
        :X => :W3
    ]
    edges_num = [m[pair.first] => m[pair.second] for pair in E]
    G = SimpleDiGraph(Edge.(edges_num))
    X = m[:Z1]
    Y = m[:Y]
    Z = Set([m[:X], m[:Z3], m[:W2]])
    @test M.d_separated(G, X, Y, Z)
    Z = Set([m[:X], m[:Z3], m[:W3]])
    @test !(M.d_separated(G, X, Y, Z))
end

@testset "without_in_out" begin
    # From https://youtu.be/pZkCecwE-xE.
    # See also https://youtu.be/9X4pR4jvKmM?t=1025 for the same example with more
    # explanation.
    s = 1  # smoking
    t = 2  # tar
    c = 3  # cancer
    g = 4  # genotype
    E = [
        s => t,
        t => c,
        g => s,
        g => c
    ]
    E = Edge.(E)
    G = SimpleDiGraphFromIterator(E)
    @test Set(M.nodes(G)) == Set([1, 2, 3, 4])

    G_without_in = M.without_incoming(G, Set([s, t]))
    @test G_without_in == SimpleDiGraph(Edge.([t => c, g => c]))

    G_without_out = M.without_outgoing(G, Set([s, t]))
    @test G_without_out == SimpleDiGraph(Edge.([g => s, g => c]))
end


text:
struct Forecast{Fl}
    expected_value::Vector{Vector{Fl}}
    covariance::Vector{Matrix{Fl}}
end

function forecast_expected_value(forec::Forecast)
    return permutedims(cat(forec.expected_value...; dims=2))
end

"""
    forecast(model::StateSpaceModel, steps_ahead::Int; kwargs...)
    forecast(model::StateSpaceModel, exogenous::Matrix{Fl}; kwargs...) where Fl

Forecast the mean and covariance for future observations from a StateSpaceModel.
"""
function forecast end

function forecast(
    model::StateSpaceModel, steps_ahead::Int; filter::KalmanFilter=default_filter(model)
)
    if has_exogenous(model)
        error(
            "The model has exogenous variables, you should use the" *
            "forecast(model::SSM, exogenous::Matrix{Fl}; kwargs...) method",
        )
    end
    # Query the type of model elements
    Fl = typeof_model_elements(model)
    # Observations to forecast
    forecasting_y = if isunivariate(model)
        [model.system.y; fill(NaN, steps_ahead)]
    else
        [model.system.y; fill(NaN, steps_ahead, size(model.system.y, 2))]
    end
    # Copy hyperparameters
    model_hyperparameters = deepcopy(model.hyperparameters)
    # Instantiate a new model
    forecasting_model = reinstantiate(model, forecasting_y)
    # Associate with the model hyperparameters
    forecasting_model.hyperparameters = model_hyperparameters
    # Perform the kalman filter
    fo = kalman_filter(forecasting_model)
    # fill forecast matrices
    expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
    covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
    for i in 1:steps_ahead
        if isunivariate(model)
            expected_value[i] = [
                dot(model.system.Z, fo.a[end - steps_ahead + i - 1]) + model.system.d
            ]
        else
            expected_value[i] =
                model.system.Z * fo.a[end - steps_ahead + i - 1] .+ model.system.d
        end
        covariance[i] = fo.F[end - steps_ahead + i]
    end
    return Forecast{Fl}(expected_value, covariance)
end

function forecast(
    model::StateSpaceModel,
    new_exogenous::Matrix{Fl};
    filter::KalmanFilter=default_filter(model),
) where Fl
    if !has_exogenous(model)
        error(
            "The model does not support exogenous variables, you should use the" *
            "forecast(model::SSM, steps_ahead::Int; kwargs...) where SSM",
        )
    end
    steps_ahead = size(new_exogenous, 1)
    forecasting_y = if isunivariate(model)
        [model.system.y; fill(NaN, steps_ahead)]
    else
        [model.system.y; fill(NaN, steps_ahead, size(model.system.y, 2))]
    end
    forecasting_X = [model.exogenous; new_exogenous]
    # Copy hyperparameters
    model_hyperparameters = deepcopy(model.hyperparameters)
    # Instantiate a new model
    forecasting_model = reinstantiate(model, forecasting_y, forecasting_X)
    # Associate with the model hyperparameters
    forecasting_model.hyperparameters = model_hyperparameters
    # Perform the kalman filter
    fo = kalman_filter(forecasting_model)
    # fill forecast matrices
    expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
    covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
    for i in 1:steps_ahead
        if isunivariate(model)
            expected_value[i] = [
                dot(model.system.Z[end - steps_ahead + i - 1], fo.a[end - steps_ahead + i - 1]) +
                model.system.d[end - steps_ahead + i - 1],
            ]
        else
            expected_value[i] =
                model.system.Z[end - steps_ahead + i - 1] * fo.a[end - steps_ahead + i - 1] +
                model.system.d[end - steps_ahead + i - 1]
        end
        covariance[i] = fo.F[end - steps_ahead + i]
    end
    return Forecast{Fl}(expected_value, covariance)
end

"""
    simulate_scenarios(
        model::StateSpaceModel, steps_ahead::Int, n_scenarios::Int;
        filter::KalmanFilter=default_filter(model)
    ) -> Array{<:AbstractFloat, 3}

Samples `n_scenarios` future scenarios via Monte Carlo simulation for `steps_ahead`
using the desired `filter`.
"""
function simulate_scenarios(
    model::StateSpaceModel,
    steps_ahead::Int,
    n_scenarios::Int;
    filter::KalmanFilter=default_filter(model),
)
    # Query the type of model elements
    Fl = typeof_model_elements(model)
    fo = kalman_filter(model)
    last_state = fo.a[end]
    num_series = size(model.system.y, 2)

    scenarios = Array{Fl,3}(undef, steps_ahead, num_series, n_scenarios)
    for s in 1:n_scenarios
        scenarios[:, :, s] = simulate(model.system, last_state, steps_ahead)
    end
    return scenarios
end


text:
"""
Given a `Modular` object, returns the same object with no empty modules (this
essentially reduces the size of `M.S` if needed).
"""
function no_empty_modules!(M::Modular)
  S = M.S
  S = S[:,[1:size(S)[2]]'[sum(S,1).>0]]
  M.S = S
  return M
end

"""
Check that there are no empty rows in a matrix.
"""
function no_empty_rows(A::Array{Int64, 2})
  return prod(sum(A, 2)) > 0
end

"""
Check that there are no empty columns in a matrix.
"""
function no_empty_columns(A::Array{Int64, 2})
  return prod(sum(A, 1)) > 0
end

"""
Check that two matrices have the same degree distributions.
"""
function same_row_marginals(A::Array{Int64, 2}, B::Array{Int64, 2})
  return sum(A, 2) == sum(B, 2)
end


text:
using Seismic
using Base.Test

# test that FFTOp passes the dot product test
m_rand = randn(20,20,20,20) + 1im*randn(20,20,20,20);
d_rand = randn(20,20,20,20) + 1im*randn(20,20,20,20);
T = rand(20,20,20,20); a = find(T.<=0.5); b = find(T.>0.5); T[a] = 0.; T[b] = 1.;
operators = [WeightingOp,FFTOp,WeightingOp]
parameters = [Dict(:w=>T),Dict(:normalize=>true),Dict(:w=>ones(T))]
a,b = Seismic.DotTest(m_rand,d_rand,operators,parameters)
@test abs((a-b)/a) < 0.01

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
######################################################################################################
# SetCategorical
using LinearAlgebra
using GridInterpolations

struct SetCategorical{S}
    elements::Vector{S} # Set elements (could be repeated)
    distr::Categorical # Categorical distribution over set elements

    function SetCategorical(elements::AbstractVector{S}) where S
        weights = ones(length(elements))
        return new{S}(elements, Categorical(normalize(weights, 1)))
    end

    function SetCategorical(elements::AbstractVector{S}, weights::AbstractVector{Float64}) where S
        ℓ₁ = norm(weights,1)
        if ℓ₁ < 1e-6 || isinf(ℓ₁)
            return SetCategorical(elements)
        end
        distr = Categorical(normalize(weights, 1))
        return new{S}(elements, distr)
    end
end

Distributions.rand(D::SetCategorical) = D.elements[rand(D.distr)]
Distributions.rand(D::SetCategorical, n::Int) = D.elements[rand(D.distr, n)]
function Distributions.pdf(D::SetCategorical, x)
    sum(e == x ? w : 0.0 for (e,w) in zip(D.elements, D.distr.p))
end


text:
# what is it: joGaussian

## helper module
module joGaussian_etc
    using JOLI: jo_convert, LocalVector
    using Random
    using LinearAlgebra
    function fI(v::LocalVector{vdt},m::Integer,n::Integer,rng::AbstractRNG,rns::Array{UInt32,1},rdt::DataType) where vdt<:Number
        rv = zeros(rdt,m)
        lrng=Random.seed!(rng,rns)
        for i=1:n
            rv = rv + randn(lrng,rdt,m) * v[i]
        end
        Random.seed!(rng,rns)
        rv = jo_convert(rdt,rv,false)
        return rv
    end
    function aI(v::LocalVector{vdt},m::Integer,n::Integer,rng::AbstractRNG,rns::Array{UInt32,1},rdt::DataType) where vdt<:Number
        rv = Vector{rdt}(undef,n)
        lrng=Random.seed!(rng,rns)
        for i=1:n
            rv[i] = (randn(lrng,rdt,1,m)*v)[1]
        end
        Random.seed!(rng,rns)
        rv = jo_convert(rdt,rv,false)
        return rv
    end
    function scale(m::Integer,n::Integer,rng::AbstractRNG,rns::Array{UInt32,1},rdt::DataType)
        rv = Vector{rdt}(undef,n)
        lrng=Random.seed!(rng,rns)
        for i=1:n
            #v = rdt(2.) .* randn(lrng,rdt,m)
            v = randn(lrng,rdt,m)
            rv[i] = (one(rdt)/sqrt(v'*v))
        end
        Random.seed!(rng,rns)
        rv = jo_convert(rdt,rv,false)
        return rv
    end
    function fIN(v::LocalVector{vdt},m::Integer,n::Integer,scale::LocalVector{<:Number},rng::AbstractRNG,rns::Array{UInt32,1},rdt::DataType) where vdt<:Number
        rv = zeros(rdt,m)
        lrng=Random.seed!(rng,rns)
        for i=1:n
            rv = rv + randn(lrng,rdt,m) * (scale[i]*v[i])
        end
        Random.seed!(rng,rns)
        rv = jo_convert(rdt,rv,false)
        return rv
    end
    function aIN(v::LocalVector{vdt},m::Integer,n::Integer,scale::LocalVector{<:Number},rng::AbstractRNG,rns::Array{UInt32,1},rdt::DataType) where vdt<:Number
        rv = Vector{rdt}(undef,n)
        lrng=Random.seed!(rng,rns)
        for i=1:n
            rv[i] = scale[i]*((randn(lrng,rdt,1,m)*v)[1])
        end
        Random.seed!(rng,rns)
        rv = jo_convert(rdt,rv,false)
        return rv
    end
end
using .joGaussian_etc

export joGaussian
"""
    julia> op = joGaussian(M,[N];
                [implicit=...,][normalized=...,][orthonormal=...,][RNG=...,]
                [DDT=...,][RDT=...,][name=...])

Gausian matrix

# Signature

    joGaussian(M::Integer,N::Integer=M;
        implicit::Bool=false,normalized::Bool=false,orthonormal::Bool=false,
        RNG::AbstractRNG=Random.seed!(),
        DDT::DataType=joFloat,RDT::DataType=DDT,
        name::String="joGaussian")

# Arguments

- `M[,N]`: sizes
- keywords
    - `implicit`: keyword element-free operator if true
    - `normalized`: keyword normalized (unit global-norm for explicit or unit-column norm for implicit) operator if true
    - `orthonormal`: keyword explict orthonormal operator if true (above implicit/normalized is ignorred); requires M<=N
    - `RNG`: keyword random number generator function with explicit seeding
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Notes

- AbstractRNG has to support seeding (RNG.seed atribute to type); e.g. RandomDevice() will not work

# Examples

not-normalized and explict dense matrix

    A=joGaussian(M,N)

explicit orthonormal dense matrix

    A=joGaussian(M,orthonormal=true)

not-normalized and element-free operator

    A=joGaussian(M,implicit=true)

normalized and explict dense matrix

    A=joGaussian(M,normalized=true)

normalized and element-free operator

    A=joGaussian(M,implicit=true,normalized=true)

examples with DDT/RDT

    A=joGaussian(M,N; DDT=Float32)
    A=joGaussian(M,N; DDT=Float32,RDT=Float64)

"""
function joGaussian(M::Integer,N::Integer=M;
         implicit::Bool=false,normalized::Bool=false,orthonormal::Bool=false,
         RNG::AbstractRNG=Random.seed!(),
         DDT::DataType=joFloat,RDT::DataType=DDT,
         name::String="joGaussian")

    if orthonormal
        M>N && throw(joLinearFunctionException("Orthonormal mode is not supported when M > N"))
        A = qr(randn(RNG,DDT,N,M))
        #a = DDT(0.5)*(A.Q)'
        a = (A.Q)'
        return joMatrix(a;DDT=DDT,RDT=RDT,name=name*"_o")
    else
        if implicit
            rngs=copy(RNG.seed)
            if normalized
                fscale = joGaussian_etc.scale(M,N,RNG,rngs,DDT)
                ascale = joGaussian_etc.scale(M,N,RNG,rngs,RDT)
                return joLinearFunctionFwd_A(M,N,
                    v1->joGaussian_etc.fIN(v1,M,N,fscale,RNG,rngs,RDT),
                    v2->joGaussian_etc.aIN(v2,M,N,ascale,RNG,rngs,DDT),
                    DDT,RDT;
                    name=name*"_in")
            else
                return joLinearFunctionFwd_A(M,N,
                    v1->joGaussian_etc.fI(v1,M,N,RNG,rngs,RDT),
                    v2->joGaussian_etc.aI(v2,M,N,RNG,rngs,DDT),
                    DDT,RDT;
                    name=name*"_i")
            end
        else
            a = randn(RNG,DDT,M,N)
            if normalized
                nf=one(DDT)/sqrt(sum(a.*a))
                a = a * sparse(nf*LinearAlgebra.I,N,N)
                return joMatrix(a;DDT=DDT,RDT=RDT,name=name*"_en")
            else
                return joMatrix(a;DDT=DDT,RDT=RDT,name=name*"_e")
            end
        end
    end
end


text:
# -----------------------------------------------------------
#
# Using boundary terms and integration (summation) by parts
# to construct more accurate derivative matrices for time
# evolution using Lax-Friedricks-type fluxes
#
# All derivatives are computed here in pos basis first
# Then conjugated into hier basis
#
# -----------------------------------------------------------

# -----------------------------------------------------
# Generates the coefficients for a traveling wave
# A cos(\vec{k} \cdot \vec{x} + \phi), \vec{k} = 2 \pi \vec{m}
# with sparse interpolation of type (k,n)
# using periodic boundary conditions in D-dimensionss
# -----------------------------------------------------

function cos_coeffs(k::Int, n::Int, m::Array{Int,1};
					scheme="sparse", phase = 0.0, A = 1.0)
	D = length(m)
	wavenumber = 2*pi*m

	# Begin with writing cos(\sum_i k_i x_i) as a sum of products
	# of sines/cosines of individual k_j x_j
	sines   = [i==1?x->sin(wavenumber[i]*x[1]+phase):x->sin(wavenumber[i]*x[1]) for i in 1:D]
	cosines = [i==1?x->cos(wavenumber[i]*x[1]+phase):x->cos(wavenumber[i]*x[1]) for i in 1:D]

	sine_dicts   = [coeffs_DG(1, k, n, sines[i]) for i in 1:D]
	cosine_dicts = [coeffs_DG(1, k, n, cosines[i]) for i in 1:D]

	ansVect = zeros(get_size(D, k, n; scheme=scheme))

	for SCs in CartesianRange(ntuple(q->2, D))
		num_sines = sum([SCs[i]-1 for i in 1:D])
		if num_sines % 2 == 1
			continue
		end
		sign = num_sines%4==0?1:-1
		# may want to add an if-statement for if any n[i] == 0

		coeff_array = [SCs[i]==1?cosine_dicts[i]:sine_dicts[i] for i in 1:D]
		productDict = tensor_construct(D, k, n, coeff_array; scheme=scheme)
		productVect = D2V(D, k, n, productDict; scheme=scheme)

		ansVect += sign * productVect
	end

	return A * ansVect
end

# -----------------------------------------------------
# The same as above, but using sin
# -----------------------------------------------------

function sin_coeffs(k::Int, n::Int, m::Array{Int,1};
					scheme="sparse", phase=0.0, A=1.0)
	return cos_coeffs(k, n, m; scheme=scheme, phase=phase-pi/2, A=A)
end

# -----------------------------------------------------
# Returns the data for a traveling wave
# using the above methods
# -----------------------------------------------------

function traveling_wave(k::Int, n::Int, m::Array{Int,1};
						scheme="sparse", phase=0.0, A=1.0)
	wavenumber = 2*pi*m
	frequency = sqrt(vecdot(wavenumber,wavenumber))
	u0 = x -> A * cos(vecdot(wavenumber,x) + phase)
	v0 = x -> A * frequency * sin(vecdot(k,x) + phase)
	u0_coeffs = cos_coeffs(k, n, m; scheme=scheme, phase=phase, A=A)
	v0_coeffs = sin_coeffs(k, n, m; scheme=scheme, phase=phase, A=A*frequency)

	return (u0_coeffs, v0_coeffs, u0, v0)
end

# -----------------------------------------------------
# Evolves a traveling wave using the above methods
# between time0 and time1
# using an ODE solver of type 'order' (default 45)
# -----------------------------------------------------
function traveling_wave_solver(k::Int, n::Int, m::Array{Int,1}, time0::Real, time1::Real;
								scheme="sparse", phase=0.0, A=1.0, order="45", kwargs...)
	D = length(m)
	f0coeffs, v0coeffs = traveling_wave(k, n, m; scheme=scheme, phase=phase, A=A)
	laplac = laplacian_matrix(D, k, n; scheme=scheme)

	RHS, y0 = wave_data(laplac, f0coeffs, v0coeffs)

	if order == "78"
		soln=ode78((t,x)->*(RHS,x), y0, [time0,time1]; kwargs...)
	elseif order == "45"
		soln=ode45((t,x)->*(RHS,x), y0, [time0,time1]; kwargs...)
	else
		throw(ArgumentError)
	end
	return soln
end


text:
    safe_email_tlds = ("org", "com", "net")
    free_email_domains = ("gmail.com", "yahoo.com", "hotmail.com")
    tlds = ("com", "com", "com", "com", "com", "com", "biz", "info", "net", "org")

    uri_pages = (
    "index", "home", "search", "main", "post", "homepage", "category", "register", "login", "faq", "about", "terms",
    "privacy", "author")
    uri_paths = (
    "app", "main", "wp-content", "search", "category", "tag", "categories", "tags", "blog", "posts", "list", "explore")
    uri_extensions = (".html", ".html", ".html", ".htm", ".htm", ".php", ".php", ".jsp", ".asp")

    user_name_formats = (
        "$(last_name()).$(first_name())",
        "$(first_name()).$(last_name())",
        "$(first_name())##",
        "?$(last_name())",
    )

    email()= email_formats[rand(1:length(email_formats))]

    user_name()= bothify(user_name_formats[rand(1:length(user_name_formats))])

    safe_email()= "$(user_name())@example.$(safe_email_tlds[rand(1:length(safe_email_tlds))])"

    free_email()= "$(user_name())@$(free_email_domain())"

    domain_suffix()= tlds[rand(1:length(tlds))]

    domain_word()=
       ( company_elements = split(company(), [',', ' ','-']);
        sal = "";
        for i in company_elements
                   if length(i)>0
                       sal*="$(i[1])"
                   end
              end;

        lowercase(sal) )

    domain_name()="$(domain_word()).$(domain_suffix())"

    company_email()= "$(user_name())@$(domain_name())"

    free_email_domain() = free_email_domains[rand(1:length(free_email_domains))]

    url()=url_formats[rand(1:length(url_formats))]

    ipv4()= "$(rand(172:255)).$(rand(0:255)).$(rand(0:255)).$(rand(0:255))"

    ipv6()= (res = "";
        for i = 1:8
            res*= "$(hex(rand(4100:65535)))"
            if i < 8
                res*= ":"
            end
        end;res)

    mac_address()=
        (res = "";
        for i = 1:6
            res*= "$(hex(rand(0x00:0xff)))"
            i < 6 && (res*= ":")
        end;res)

    uri_page()= uri_pages[rand(1:length(uri_pages))]

    uri_path(deep="None")=
        (if deep =="None"
            deep = rand(1:3)
        end; sa="";
            for i =0:deep
               sa*= uri_paths[rand(1:length(uri_paths))]
                i < deep && (sa*= "/")
            end; sa )

    uri_extension()= uri_extensions[rand(1:length(uri_extensions))]

    uri()= uri_formats[rand(1:length(uri_formats))]

    #=slug(value="None")=
        if value == "None"
            value = Lorem.text(20)
            value=#

    image_url(width="None", height="None")=
       image_placeholder_services[rand(1:length(image_placeholder_services))]

    email_formats = (
        "$(user_name())@$(domain_name())",
        "$(user_name())@$(free_email_domain())",
    )
    url_formats = (
        "http://www.$(domain_name())/",
        "http://$(domain_name())/",
    )
    uri_formats = (
        "$(url())",
        "$(url())$(uri_page())/",
        "$(url())$(uri_page())$(uri_extension())",
        "$(url())$(uri_path())/$(uri_page())/",
        "$(url())$(uri_path())/$(uri_page())$(uri_extension())",
    )
    image_placeholder_services = (
        "http://placekitten.com/$(rand(1:1024))/$(rand(1:1024))",
        "http://placehold.it/$(rand(1:1024))x$(rand(1:1024))",
        "http://www.lorempixum.com/$(rand(1:1024))/$(rand(1:1024))",
        "http://dummyimage.com/$(rand(1:1024))x$(rand(1:1024))",
     )


text:
# ---
# title: 781. Rabbits in Forest
# id: problem781
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/rabbits-in-forest/description/>
# hidden: true
# ---
#
# In a forest, each rabbit has some color. Some subset of rabbits (possibly all
# of them) tell you how many other rabbits have the same color as them. Those
# `answers` are placed in an array.
#
# Return the minimum number of rabbits that could be in the forest.
#
#
#
#     **Examples:**
#     Input: answers = [1, 1, 2]
#     Output: 5
#     Explanation:
#     The two rabbits that answered "1" could both be the same color, say red.
#     The rabbit than answered "2" can't be red or the answers would be inconsistent.
#     Say the rabbit that answered "2" was blue.
#     Then there should be 2 other blue rabbits in the forest that didn't answer into the array.
#     The smallest possible number of rabbits in the forest is therefore 5: 3 that answered plus 2 that didn't.
#
#     Input: answers = [10, 10, 10]
#     Output: 11
#
#     Input: answers = []
#     Output: 0
#
#
# **Note:**
#
#   1. `answers` will have length at most `1000`.
#   2. Each `answers[i]` will be an integer in the range `[0, 999]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using FourierTransforms
using OffsetArrays
using Test

include("util/symbolic.jl")

using FourierTransforms: @nontwiddle, @twiddle, CTPlan
FourierTransforms.ωpow(::Type{SymbolicExpr}, n, i) = RootOfUnity(n)^i

@eval @nontwiddle($SymbolicExpr, true, 4)
@eval @nontwiddle($SymbolicExpr, true, 2)
@eval @twiddle($SymbolicExpr, true, 4)
@eval @twiddle($SymbolicExpr, true, 2)

# Input length 8
let inputs = [Var(Symbol(string('x', subscripts[i]))) for i = 0:7]
    plan = CTPlan(SymbolicExpr, true, 8)
    y = plan * inputs
    y = map(simplify, map(expand, y))

    naive = map(simplify, [sum(inputs[i+1]*RootOfUnity(8)^(i*j) for i = 0:7) for j = 0:7])

    # Test that the fft is symbolically equal to the naive expression.
    @test all(map(simplify, map(expand, map(-, y, naive)))) do x
        x == Plus(())
    end
end

# Input length 7 (excercising bluestein fft)
# TODO: This requires more support than the SymbolicExpr type currently has
#=
let inputs = [Var(Symbol(string('x', subscripts[i]))) for i = 0:5]
    plan = CTPlan(SymbolicExpr, true, length(inputs))
    y = plan * inputs
    y = map(simplify, map(expand, y))

    naive = map(simplify, [sum(inputs[i+1]*RootOfUnity(length(inputs))^(i*j) for i = 0:length(inputs)-1) for j = 0:length(inputs)-1])

    # Test that the fft is symbolically equal to the naive expression.
    @test all(map(simplify, map(expand, map(-, y, naive)))) do x
        x == Plus(())
    end
end
=#


text:
using CRBM
using Base.Test

# write your own tests here
@test 1 == 1

# test basic functions - integer to binary vector
@test [0.0, 0.0, 0.0, 0.0] == i2b(0 ,4)
@test [0.0, 0.0, 0.0, 1.0] == i2b(1 ,4)
@test [0.0, 0.0, 1.0, 0.0] == i2b(2 ,4)
@test [0.0, 0.0, 1.0, 1.0] == i2b(3 ,4)
@test [0.0, 1.0, 0.0, 0.0] == i2b(4 ,4)
@test [0.0, 1.0, 0.0, 1.0] == i2b(5 ,4)
@test [0.0, 1.0, 1.0, 0.0] == i2b(6 ,4)
@test [0.0, 1.0, 1.0, 1.0] == i2b(7 ,4)
@test [1.0, 0.0, 0.0, 0.0] == i2b(8 ,4)
@test [1.0, 0.0, 0.0, 1.0] == i2b(9 ,4)
@test [1.0, 0.0, 1.0, 0.0] == i2b(10,4)
@test [1.0, 0.0, 1.0, 1.0] == i2b(11,4)
@test [1.0, 1.0, 0.0, 0.0] == i2b(12,4)
@test [1.0, 1.0, 0.0, 1.0] == i2b(13,4)
@test [1.0, 1.0, 1.0, 0.0] == i2b(14,4)
@test [1.0, 1.0, 1.0, 1.0] == i2b(15,4)

# test basic functions - integer vector to binary vector
@test [0.0, 0.0, 0.0, 0.0] == iv2b([0, 0], 2)
@test [0.0, 1.0, 1.0, 0.0] == iv2b([1, 2], 2)
@test [0.0, 1.0, 0.0, 0.0] == iv2b([4],    4)

# test basic functions - float matrix to binary vector
@test [[0.0 0.0 0.0 0.0], [0.0 0.0 0.0 0.0]] == binarise_matrix([[-1.0 -1.0], [-1.0 -1.0]], 4)

# test inverse binning
@test 0 == b2i([0.0, 0.0])
@test 1 == b2i([0.0, 1.0])
@test 2 == b2i([1.0, 0.0])
@test 3 == b2i([1.0, 1.0])

# test inverse binning - vector
@test [0, 1] == b2iv([0.0, 0.0, 0.0, 1.0], 2)
@test [2, 3] == b2iv([1.0, 0.0, 1.0, 1.0], 2)

@test [-0.75, -0.25] == bv2dv([0.0, 0.0, 0.0, 1.0], 4)
@test [0.25, 0.75]   == bv2dv([1.0, 0.0, 1.0, 1.0], 4)

# test inverse binning - matrix
@test [[-0.75 -0.25], [0.25 0.75]] == unbinarise_matrix([[0.0 0.0 0.0 1.0], [1.0 0.0 1.0 1.0]], 4)

# test binary_draw
@test [1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0] == binary_draw([1.1 1.1 1.1 1.1 -0.1 -0.1 -0.1 -0.1])

# test against old code

function undiscretise(v, nr_of_bins)
  n = int(log2(nr_of_bins))
  l = int(length(v)/log2(nr_of_bins))
  r = zeros(l)
  for i=1:l
    r[i] = binary2integer(v[(i-1)*n+1:i*n])
  end
  r = 2.0 .* (r ./ (nr_of_bins) .+ 1 / (2*nr_of_bins)) .- 1.0
end

function binary2integer(v) # tested
  d = convert(Vector{Int64}, v)
  r = 0
  n = length(d)
  for i=1:n
    r += (d[i]>0)?(1<<(n-i)):0
  end
  return r
end

for i=0:200
  j = float(i-100)/200.0
  k = bv(j, 16, -1.0, 1.0)
  l = i2b(k,4)
  @test undiscretise(l, 16)[1] == bv2dv(l, 16)
end

function int2binary(v::Int64, n::Int64) # checked
  r=zeros(n)
  for i=1:n
    r[i] = (((1 << (n-i)) & v)>0)?1.0:0.0
  end
  return r
end

function discretise(v, nr_of_bins) # checked
  d = convert(Vector{Int64},min(floor(nr_of_bins * (v .+ 1) ./ 2.0), nr_of_bins-1))
  r = Float64[]
  u = int(log2(nr_of_bins))
  for i = 1:length(d)
    r = append!(r,int2binary(d[i], u))
  end
  return r
end

for i=0:200
  j = float(i-100)/200.0
  k = discretise([j], 16)
  l = iv2b(bin_vector([j], -1.0, 1.0, 16), 4)
  @test k == l
end

for i=1:100
  r = rand(12)
  #= @test discretise(r, 16) == binarise_vector(r, 16) =#
  @test discretise(r, 16) == iv2b(bin_vector(r, -1.0, 1.0, 16), 4)
end


text:
module DataTypes

    type Problem
        reynolds    ::Int
        bc_left     ::Int
        bc_bottom   ::Int
        bc_right         ::Int
        bc_top           ::Int
        flow_magnitude   ::Real
        flow_angle       ::Real

        Problem() = new(0,0,0,0,0,0.0,0.0)
    end

    type Grid
        nx          ::Int
        ny          ::Int
        width       ::Real
        height      ::Real

        Grid() = new(0,0,0.0,0.0)
    end

    type Solver
        cg_eps      ::Real
        timestep    ::Real
        nt          ::Int

        Solver() = new(0.0,0.0,0)
    end

    type Parameters
        problem::Problem
        grid::Grid
        solver::Solver

        Parameters() = new(Problem(), Grid(), Solver())
    end

    # Data structure for instantaneous flow data
    type FlowSnapshot
        u::Vector{Float64}
        p::Vector{Float64}
        t::Float64

        FlowSnapshot() = new()
    end

end


text:
""" Here we define types to handle molecular absorption. """

abstract type AbstractAbsorber end

struct NoAbsorption <: AbstractAbsorber end

# The only type that a subtype of ContinuumAbsorber must implement is
# transmittance, which gives the transmitance between two points

""" Transmittance between points a and b.  The vector μ is passed for performance,as it is often already computed in the calling environment.  It must be
μ = (b - a) / |b - a|. """
transmittance(::NoAbsorption, a, b, μ) = 1.0

# Now we define a more useful absorption for a species with a single absorption
# cross-section and a stratified density
struct StratifiedAbsorption{I<:AbstractInterpolation} <: AbstractAbsorber
    "The cross-section of the absorption."
    σ::Float64

    "An interpolator for the cumulative density of the absorbing species."
    cuminterp::I

    # Using the cumulative density allows us to avoid integrals.  We define
    # cuminterp(z) = ∫(from z0 to z) n(z') dz' so the integral from z1 to z2 is
    # cuminterp(z2) - cuminterp(z1).
end

""" Initialize a StratifiedAbsorption with the given cross-section and
    number density evaluated at the values of z.
"""
function StratifiedAbsorption(σ, z::AbstractVector, n::AbstractVector)
    cumn = cumtrapz(z, n)
    cuminterp = LinearInterpolation(z, cumn)
    StratifiedAbsorption(σ, cuminterp)
end

function transmittance(str::StratifiedAbsorption, a, b, μ)
    vertint = str.cuminterp(b[3]) - str.cuminterp(a[3])

    if μ[3] != 0
        return @fastmath exp(-str.σ * abs(vertint / μ[3]))
    else
        # This is the unlikely situation where the ray propagates exaclty
        # horizontally.  This complicates things a bit.
        # Since we only store the cumulative integral we obtain the density
        # by evaluating the derivative.
        n = Interpolations.gradient(a[3])[1]
        return @fastmath exp(-str.σ * n * norm(a - b))
    end
end

struct MultipleAbsorption{T<:Tuple} <: AbstractAbsorber
    absorbers::T
end

# See geometry.jl for a similar construction and why we use generated functions
@generated function transmittance(u::MultipleAbsorption{T}, a, b, μ) where {T}
    L = fieldcount(T)
    out = quote end

    push!(out.args, :(total = 1.0))

    for i in 1:L
        push!(out.args,
              quote
              total *= transmittance(u.absorbers[$i], a, b, μ)
              end
              )
    end
    push!(out.args, :(return total))

    out
end

""" Cumulative trapezoidal integration of a function with values `f` evaluated
    at `x`.  The result is stored in `r`.
"""
function cumtrapz!(r, x, f)
    @assert length(x) == length(f)
    @assert length(r) == length(r)
    cum = zero(eltype(r))

    for i in eachindex(x)
        r[i] = cum
        if i != lastindex(x)
            cum += 0.5 * (f[i + 1] + f[i]) * (x[i + 1] - x[i])
        end
    end
    r
end

""" Like `cumtrapz!` but allocates the output vector. """
function cumtrapz(x, f)
    r = similar(f)
    cumtrapz!(r, x, f)
end


text:
module HuntTests

using GridapMHD: hunt

hunt(
  nc=(4,4),
  L=1.0,
  B=(0.,50.,0.),
  debug=true,
  vtk=true,
  title="hunt",
  solver=:julia,
 )

hunt(
  nc=(4,4),
  L=1.0,
  B=(0.,50.,0.),
  debug=false,
  vtk=true,
  title="hunt",
  solver=:julia,
 )

hunt(
  nc=(4,4),
  np=(1,1),
  backend=:mpi,
  L=1.0,
  B=(0.,50.,0.),
  debug=false,
  vtk=true,
  title="hunt",
  solver=:julia,
 )

hunt(
  nc=(4,4),
  np=(2,2),
  backend=:sequential,
  L=2.0,
  B=(0.,50.,0.),
  debug=false,
  vtk=true,
  title="hunt",
  solver=:julia,
 )

end # module


text: include("vec.jl")


text:
using DReal
reset_ctx!()

a = Var(Float64,"aexpop",0.,1.)
b2 = (a < 0.1) | (a > 0.9)
push_ctx!()
add!(a>0.0)
add!(a<1.0)
add!(b2)
is_satisfiable()
pop_ctx!()
push_ctx!()
add!(a>0.0)
add!(a<1.0)
add!(!b2)
is_satisfiable()


text:
import UUIDs: uuid1

import TableIOInterface: get_example_code, is_extension_supported

"Will hold all 'response handlers': functions that respond to a WebSocket request from the client."
const responses = Dict{Symbol,Function}()

Base.@kwdef struct ClientRequest
    session::ServerSession
    notebook::Union{Nothing,Notebook}
    body::Any=nothing
    initiator::Union{Initiator,Nothing}=nothing
end

require_notebook(r::ClientRequest) = if r.notebook === nothing
    throw(ArgumentError("Notebook request called without a notebook 😗"))
end

###
# RESPONDING TO A NOTEBOOK STATE UPDATE
###

"""
## State management in Pluto

*Aka: how do the server and clients stay in sync?*

A Pluto notebook session has *state*: with this, we mean:

1. The input and ouput of each cell, the cell order, and more metadata about the notebook and cells [^state]

This state needs to be **synchronised between the server and all clients** (we support multiple synchronised clients), and note that:

- Either side wants to update the state. Generally, a client will update cell inputs, the server will update cell outputs.
- Both sides want to *react* to state updates
- The server is in Julia, the clients are in JS
- This is built on top of our websocket+msgpack connection, but that doesn't matter too much

We do this by implementing something similar to how you use Google Firebase: there is **one shared state object, any party can mutate it, and it will synchronise to all others automatically**. The state object is a nested structure of mutable `Dict`s, with immutable ints, strings, bools, arrays, etc at the endpoints.

Some cool things are:

- Our system uses object diffing, so only *changes to the state* are actually tranferred over the network. But you can use it as if the entire state is sent around constantly.
- In the frontend, the *shared state* is part of the *react state*, i.e. shared state updates automatically trigger visual updates.
- Within the client, state changes take effect instantly, without waiting for a round trip to the server. This means that when you add a cell, it shows up instantly.

Diffing is done using `immer.js` (frontend) and `src/webserver/Firebasey.jl` (server). We wrote Firebasey ourselves to match immer's functionality, and the cool thing is: **it is a Pluto notebook**! Since Pluto notebooks are `.jl` files, we can just `include` it in our module.

The shared state object is generated by [`notebook_to_js`](@ref). Take a look! The Julia server orchestrates this firebasey stuff. For this, we keep a **copy** of the latest state of each client on the server (see [`current_state_for_clients`](@ref)). When anything changes to the Julia state (e.g. when a cell finished running), we call [`send_notebook_changes!`](@ref), which will call [`notebook_to_js`](@ref) to compute the new desired state object. For each client, we diff the new state to their last known state, and send them the difference.

### Responding to changes made by a client

When a client updates the shared state object, we want the server to *react* to that change by taking an action. Which action to take depends on which field changes. For example, when `state["path"]` changes, we should rename the notebook file. When `state["cell_inputs"][a_cell_id]["code"]` changes, we should reparse and analyze that cel, etc. This location of the change, e.g.  `"cell_inputs/<a_cell_id>/code"` is called the *path* of the change.

[`effects_of_changed_state`](@ref) define these pattern-matchers. We use a `Wildcard()` to take the place of *any* key, see [`Wildcard`](@ref), and we use the change/update/patch inside the given function.

### Not everything uses the shared state (yet)

Besides `:update_notebook`, you will find more functions in [`responses`](@ref) that respond to classic 'client requests', such as `:reshow_cell` and `:shutdown_notebook`. Some of these requests get a direct response, like the list of autocomplete options to a `:complete` request (in `src/webserver/REPLTools.jl`). On the javascript side, these direct responses can be `awaited`, because every message has a unique ID.

[^state]:
    Two other meanings of _state_ could be:
    2. The reactivity data: the parsed AST (`Expr`) of each cell, which variables are defined or referenced by which cells, in what order will cells run?
    3. The state of the Julia process: i.e. which variables are defined, which packages are imported, etc.

    The first two (1 & 2) are stored in a [`Notebook`](@ref) struct, remembered by the server process (Julia). (In fact, (2) is entirely described by (1), but we store it for performance reasons.) I included (3) for completeness, but it is not stored by us, we hope to control and minimize (3) by keeping track of (1) and (2).

"""
module Firebasey include("./Firebasey.jl") end

# All of the arrays in the notebook_to_js object are 'immutable' (we write code as if they are), so we can enable this optimization:
Firebasey.use_triple_equals_for_arrays[] = true

# the only possible Arrays are:
# - cell_order
# - cell_result > * > output > body
# - bonds > * > value > *

function notebook_to_js(notebook::Notebook)
    Dict{String,Any}(
        "notebook_id" => notebook.notebook_id,
        "path" => notebook.path,
        "in_temp_dir" => startswith(notebook.path, new_notebooks_directory()),
        "shortpath" => basename(notebook.path),
        "process_status" => notebook.process_status,
        "cell_inputs" => Dict{UUID,Dict{String,Any}}(
            id => Dict{String,Any}(
                "cell_id" => cell.cell_id,
                "code" => cell.code,
                "code_folded" => cell.code_folded,
            )
        for (id, cell) in notebook.cells_dict),
        "cell_results" => Dict{UUID,Dict{String,Any}}(
            id => Dict{String,Any}(
                "cell_id" => cell.cell_id,
                "queued" => cell.queued,
                "running" => cell.running,
                "errored" => cell.errored,
                "runtime" => cell.runtime,
                "output" => Dict(
                    "last_run_timestamp" => cell.last_run_timestamp,
                    "persist_js_state" => cell.persist_js_state,
                    "mime" => cell.repr_mime,
                    "body" => cell.output_repr,
                    "rootassignee" => cell.rootassignee,
                ),
            )
        for (id, cell) in notebook.cells_dict),
        "cell_order" => notebook.cell_order,
        "bonds" => Dict{String,Dict{String,Any}}(
            String(key) => Dict("value" => bondvalue.value)
        for (key, bondvalue) in notebook.bonds),
    )
end

"""
For each connected client, we keep a copy of their current state. This way we know exactly which updates to send when the server-side state changes.
"""
const current_state_for_clients = WeakKeyDict{ClientSession,Any}()

"""
Update the local state of all clients connected to this notebook.
"""
function send_notebook_changes!(🙋::ClientRequest; commentary::Any=nothing)
    notebook_dict = notebook_to_js(🙋.notebook)
    for (_, client) in 🙋.session.connected_clients
        if client.connected_notebook !== nothing && client.connected_notebook.notebook_id == 🙋.notebook.notebook_id
            current_dict = get(current_state_for_clients, client, :empty)
            patches = Firebasey.diff(current_dict, notebook_dict)
            patches_as_dicts::Array{Dict} = patches
            current_state_for_clients[client] = deep_enough_copy(notebook_dict)

            # Make sure we do send a confirmation to the client who made the request, even without changes
            is_response = 🙋.initiator !== nothing && client == 🙋.initiator.client

            if !isempty(patches) || is_response
                response = Dict(
                    :patches => patches_as_dicts,
                    :response => is_response ? commentary : nothing
                )
                putclientupdates!(client, UpdateMessage(:notebook_diff, response, 🙋.notebook, nothing, 🙋.initiator))
            end
        end
    end
end

"Like `deepcopy`, but anything onther than `Dict` gets a shallow (reference) copy."
function deep_enough_copy(d::Dict{A,B}) where {A, B}
    Dict{A,B}(
        k => deep_enough_copy(v)
        for (k, v) in d
    )
end
deep_enough_copy(x) = x

"""
A placeholder path. The path elements that it replaced will be given to the function as arguments.
"""
struct Wildcard end

@enum Changed begin
    CodeChanged
    FileChanged
end

# to support push!(x, y...) # with y = []
Base.push!(x::Set{Changed}) = x

const no_changes = Changed[]

const effects_of_changed_state = Dict(
    "path" => function(; request::ClientRequest, patch::Firebasey.ReplacePatch)
        newpath = tamepath(patch.value)
        # SessionActions.move(request.session, request.notebook, newpath)

        if isfile(newpath)
            throw(UserError("File exists already - you need to delete the old file manually."))
        else
            move_notebook!(request.notebook, newpath)
            putplutoupdates!(request.session, clientupdate_notebook_list(request.session.notebooks))
            WorkspaceManager.cd_workspace((request.session, request.notebook), newpath)
        end
        return no_changes
    end,
    "process_status" => function(; request::ClientRequest, patch::Firebasey.ReplacePatch)
        newstatus = patch.value

        @info "Process status set by client" newstatus
    end,
    "in_temp_dir" => function(; _...) no_changes end,
    "cell_inputs" => Dict(
        Wildcard() => function(cell_id, rest...; request::ClientRequest, patch::Firebasey.JSONPatch)
            Firebasey.applypatch!(request.notebook, patch)

            if length(rest) == 0
                [CodeChanged, FileChanged]
            elseif length(rest) == 1 && Symbol(rest[1]) == :code
                request.notebook.cells_dict[UUID(cell_id)].parsedcode = nothing
                [CodeChanged, FileChanged]
            else
                [FileChanged]
            end
        end,
    ),
    "cell_order" => function(; request::ClientRequest, patch::Firebasey.ReplacePatch)
        Firebasey.applypatch!(request.notebook, patch)
        [FileChanged]
    end,
    "bonds" => Dict(
        Wildcard() => function(name; request::ClientRequest, patch::Firebasey.JSONPatch)
            name = Symbol(name)
            Firebasey.applypatch!(request.notebook, patch)
            set_bond_value_reactive(
                session=request.session,
                notebook=request.notebook,
                name=name,
                is_first_value=patch isa Firebasey.AddPatch,
                run_async=true,
            )
            # [BondChanged]
            return no_changes
        end,
    )
)

responses[:update_notebook] = function response_update_notebook(🙋::ClientRequest)
    require_notebook(🙋)
    try
        notebook = 🙋.notebook
        patches = (Base.convert(Firebasey.JSONPatch, update) for update in 🙋.body["updates"])

        if length(patches) == 0
            send_notebook_changes!(🙋)
            return nothing
        end

        if !haskey(current_state_for_clients, 🙋.initiator.client)
            throw(ErrorException("Updating without having a first version of the notebook??"))
        end

        # TODO Immutable ??
        for patch in patches
            Firebasey.applypatch!(current_state_for_clients[🙋.initiator.client], patch)
        end

        changes = Set{Changed}()

        for patch in patches
            (mutator, matches, rest) = trigger_resolver(effects_of_changed_state, patch.path)

            current_changes = if isempty(rest) && applicable(mutator, matches...)
                mutator(matches...; request=🙋, patch=patch)
            else
                mutator(matches..., rest...; request=🙋, patch=patch)
            end

            push!(changes, current_changes...)
        end

        # if CodeChanged ∈ changes
        #     update_caches!(notebook, cells)
        #     old = notebook.topology
        #     new = notebook.topology = updated_topology(old, notebook, cells)
        # end

        # If CodeChanged ∈ changes, then the client will also send a request like run_multiple_cells, which will trigger a file save _before_ running the cells.
        # In the future, we should get rid of that request, and save the file here. For now, we don't save the file here, to prevent unnecessary file IO.
        # (You can put a log in save_notebook to track how often the file is saved)
        if FileChanged ∈ changes && CodeChanged ∉ changes
            save_notebook(notebook)
        end

        send_notebook_changes!(🙋; commentary=Dict(:update_went_well => :👍))
    catch ex
        @error "Update notebook failed"  🙋.body["updates"] exception=(ex, stacktrace(catch_backtrace()))
        response = Dict(
            :update_went_well => :👎,
            :why_not => sprint(showerror, ex),
            :should_i_tell_the_user => ex isa SessionActions.UserError,
        )
        send_notebook_changes!(🙋; commentary=response)
    end
end

function trigger_resolver(anything, path, values=[])
	(value=anything, matches=values, rest=path)
end
function trigger_resolver(resolvers::Dict, path, values=[])
	if isempty(path)
		throw(BoundsError("resolver path ends at Dict with keys $(keys(resolver))"))
	end

	segment = first(path)
	rest = path[firstindex(path)+1:end]
	for (key, resolver) in resolvers
		if key isa Wildcard
			continue
		end
		if key == segment
			return trigger_resolver(resolver, rest, values)
		end
	end

	if haskey(resolvers, Wildcard())
		return trigger_resolver(resolvers[Wildcard()], rest, (values..., segment))
    else
        throw(BoundsError("failed to match path $(path), possible keys $(keys(resolver))"))
	end
end

###
# MISC RESPONSES
###

responses[:connect] = function response_connect(🙋::ClientRequest)
    putclientupdates!(🙋.session, 🙋.initiator, UpdateMessage(:👋, Dict(
        :notebook_exists => (🙋.notebook !== nothing),
        :options => 🙋.session.options,
        :version_info => Dict(
            :pluto => PLUTO_VERSION_STR,
            :julia => JULIA_VERSION_STR,
        ),
    ), nothing, nothing, 🙋.initiator))
end

responses[:ping] = function response_ping(🙋::ClientRequest)
    putclientupdates!(🙋.session, 🙋.initiator, UpdateMessage(:pong, Dict(), nothing, nothing, 🙋.initiator))
end

responses[:reset_shared_state] = function response_reset_shared_state(🙋::ClientRequest)
    delete!(current_state_for_clients, 🙋.initiator.client)
    send_notebook_changes!(🙋; commentary=Dict(:from_reset =>  true))
end

responses[:run_multiple_cells] = function response_run_multiple_cells(🙋::ClientRequest)
    require_notebook(🙋)
    uuids = UUID.(🙋.body["cells"])
    cells = map(uuids) do uuid
        🙋.notebook.cells_dict[uuid]
    end

    if will_run_code(🙋.notebook)
        foreach(c -> c.queued = true, cells)
        send_notebook_changes!(🙋)
    end

    # save=true fixes the issue where "Submit all changes" or `Ctrl+S` has no effect.
    update_save_run!(🙋.session, 🙋.notebook, cells; run_async=true, save=true)
end

responses[:get_all_notebooks] = function response_get_all_notebooks(🙋::ClientRequest)
    putplutoupdates!(🙋.session, clientupdate_notebook_list(🙋.session.notebooks, initiator=🙋.initiator))
end

responses[:interrupt_all] = function response_interrupt_all(🙋::ClientRequest)
    require_notebook(🙋)

    session_notebook = (🙋.session, 🙋.notebook)
    workspace = WorkspaceManager.get_workspace(session_notebook)

    already_interrupting = 🙋.notebook.wants_to_interrupt
    anything_running = !isready(workspace.dowork_token)
    if !already_interrupting && anything_running
        🙋.notebook.wants_to_interrupt = true
        WorkspaceManager.interrupt_workspace(session_notebook)
    end
    # TODO: notify user whether interrupt was successful
end

responses[:shutdown_notebook] = function response_shutdown_notebook(🙋::ClientRequest)
    require_notebook(🙋)
    SessionActions.shutdown(🙋.session, 🙋.notebook; keep_in_session=🙋.body["keep_in_session"])
end

without_initiator(🙋::ClientRequest) = ClientRequest(session=🙋.session, notebook=🙋.notebook)

responses[:restart_process] = function response_restrart_process(🙋::ClientRequest)
    require_notebook(🙋)

    if 🙋.notebook.process_status != ProcessStatus.waiting_to_restart
        🙋.notebook.process_status = ProcessStatus.waiting_to_restart
        send_notebook_changes!(🙋 |> without_initiator)

        SessionActions.shutdown(🙋.session, 🙋.notebook; keep_in_session=true, async=true)

        🙋.notebook.process_status = ProcessStatus.starting
        send_notebook_changes!(🙋 |> without_initiator)

        update_save_run!(🙋.session, 🙋.notebook, 🙋.notebook.cells; run_async=true, save=true)
    end
end

responses[:reshow_cell] = function response_reshow_cell(🙋::ClientRequest)
    require_notebook(🙋)
    cell = let
        cell_id = UUID(🙋.body["cell_id"])
        🙋.notebook.cells_dict[cell_id]
    end
    run = WorkspaceManager.format_fetch_in_workspace((🙋.session, 🙋.notebook), cell.cell_id, ends_with_semicolon(cell.code), (parse(PlutoRunner.ObjectID, 🙋.body["objectid"], base=16), convert(Int64, 🙋.body["dim"])))
    set_output!(cell, run)
    # send to all clients, why not
    send_notebook_changes!(🙋 |> without_initiator)
end

###
# HANDLE NEW BOND VALUES
###

function set_bond_value_reactive(; session::ServerSession, notebook::Notebook, name::Symbol, is_first_value::Bool=false, kwargs...)
    bound_sym = name
    new_value = notebook.bonds[name].value

    variable_exists = is_assigned_anywhere(notebook, notebook.topology, bound_sym)
    if !variable_exists
        # a bond was set while the cell is in limbo state
        # we don't need to do anything
        return
    end

    # TODO: Not checking for any dependents now
    # any_dependents = is_referenced_anywhere(notebook, notebook.topology, bound_sym)

    # fix for https://github.com/fonsp/Pluto.jl/issues/275
    # if `Base.get` was defined to give an initial value (read more about this in the Interactivity sample notebook), then we want to skip the first value sent back from the bond. (if `Base.get` was not defined, then the variable has value `missing`)
    # Check if the variable does not already have that value.
    # because if the initial value is already set, then we don't want to run dependent cells again.
    eq_tester = :(try !ismissing($bound_sym) && ($bound_sym == $new_value) catch; false end) # not just a === comparison because JS might send back the same value but with a different type (Float64 becomes Int64 in JS when it's an integer.)
    if is_first_value && WorkspaceManager.eval_fetch_in_workspace((session, notebook), eq_tester)
        return
    end

    function custom_deletion_hook((session, notebook)::Tuple{ServerSession,Notebook}, to_delete_vars::Set{Symbol}, funcs_to_delete::Set{Tuple{UUID,FunctionName}}, to_reimport::Set{Expr}; to_run::AbstractVector{Cell})
        to_delete_vars = Set([to_delete_vars..., bound_sym]) # also delete the bound symbol
        WorkspaceManager.delete_vars((session, notebook), to_delete_vars, funcs_to_delete, to_reimport)
        WorkspaceManager.eval_in_workspace((session, notebook), :($(bound_sym) = $(new_value)))
    end
    to_reeval = where_referenced(notebook, notebook.topology, Set{Symbol}([bound_sym]))

    update_save_run!(session, notebook, to_reeval; deletion_hook=custom_deletion_hook, save=false, persist_js_state=true, kwargs...)
end

responses[:write_file] = function (🙋::ClientRequest)
    path = 🙋.notebook.path
    reldir = "$(path |> basename).assets"
    dir = joinpath(path |> dirname, reldir)
    file_noext = reduce(*, split(🙋.body["name"], ".")[1:end - 1])
    extension = split(🙋.body["name"], ".")[end]
    save_path = numbered_until_new(joinpath(dir, file_noext); sep=" ", suffix=".$(extension)", create_file=false)

    if !ispath(dir)
        mkpath(dir)
    end
    success = try
        io = open(save_path, "w")
        write(io, 🙋.body["file"])
        close(io)
        true
    catch e
        false
    end

    code = template_code(basename(save_path), reldir, 🙋.body["file"])

    msg = UpdateMessage(:write_file_reply,
        Dict(
            :success => success,
            :code => code
        ), 🙋.notebook, nothing, 🙋.initiator)

    putclientupdates!(🙋.session, 🙋.initiator, msg)
end

# helpers

function template_code(filename, directory, iofilecontents)
    path = """joinpath(split(@__FILE__, '#')[1] * ".assets", "$(filename)")"""
    extension = split(filename, ".")[end]
    varname = replace(basename(path), r"[\"\-,\.#@!\%\s+\;()\$&*\[\]\{\}'^]" => "")

    if extension ∈ ["jpg", "jpeg", "png", "svg", "webp", "tiff", "bmp", "gif", "wav", "aac", "mp3", "mpeg", "mp4", "webm", "ogg"]
        req_code = "import PlutoUI"
        code = """$(extension)_$(varname) = let
    $(req_code)
    PlutoUI.LocalResource($(path))
end"""

    elseif extension ∈ ["txt", "text"]
        code = """txt_$(varname) = let
    $(varname) = open($(path))
    read($(varname), String)
end"""

    elseif extension ∈ ["jl"]
        io = IOBuffer();
        write(io, iofilecontents)
        code = String(take!(io))

    elseif is_extension_supported(extension)
        code = get_example_code(directory, filename)

    else
        code = missing
    end
end