text:
export readMesh, save, load, read_data, write_data, convert_mat, read_strain_stress,
        meshread, load_mesh, psread

"""
    meshread(gmshfile::String)

Reads a gmsh file `gmshfile` and return (nodes, elements) tuple.
"""
function meshread(gmshfile::String)
    cnt = read(gmshfile, String)
    cnt = replace(cnt, "\r"=>"")
    r = r"\$MeshFormat\n(\S*).*\n\$EndMeshFormat"s
    version = match(r, cnt)[1]
    println("Gmsh file version ... $version")

    r = r"\$Nodes\n(.*)\n\$EndNodes"s
    nodes = match(r, cnt)[1]
    nodes = split(nodes,'\n')
    nodes = [parse.(Float64, split(x)) for x in nodes]
    nodes = filter(x->length(x)==3,nodes)
    nodes = hcat(nodes...)'[:,1:2]
    println("Nodes ... $(size(nodes,1))")

    r = r"\$Elements\n(.*)\n\$EndElements"s
    elems = match(r, cnt)[1]
    elems = split(elems,'\n')
    elems = [parse.(Int64, split(x)) for x in elems]
    elems = filter(x->length(x)==5,elems)
    elems = hcat(elems...)'[:,2:5]
    println("Elements ... $(size(elems,1))")

    println("Remove redundant nodes...")
    s = Set(unique(elems[:]))
    k = 1
    node_new_id = Dict()
    new_id = Int64[]
    for i = 1:size(nodes, 1)
        if i in s
            node_new_id[i] = k
            push!(new_id, i)
            k += 1
        end
    end
    for i = 1:size(elems,1)
        for j = 1:4
            elems[i,j] = node_new_id[elems[i,j]]
        end
    end
    nodes = nodes[new_id,:]
    println("Preprocessed nodes ... $(size(nodes,1))")
    println("Preprocessed elements ... $(size(elems,1))")

    return nodes, elems
end

@doc raw"""
    psread(gmshfile::String)

Reads the physical group information from `gmshfile`. Returns a dictionary

```
PhysicalGroupName => Index
```

`Index` is a $n\times 2$ coordinate matrix.
We do not output indices because there may be postprocessing procedures after creating the mesh.
"""
function psread(gmshfile::String)
    cnt = read(gmshfile, String)
    cnt = replace(cnt, "\r"=>"")
    r = r"\$PhysicalNames\n(.*?)\n\$EndPhysicalNames"s
    physicalnames = match(r, cnt)[1]
    physicalnames = split(physicalnames,'\n')[2:end]

    r = r"\$Entities\n(.*?)\n\$EndEntities"s
    r = match(r, cnt)[1]
    info = split(r, '\n')
    n = parse(Int64, split(info[1])[1])
    D = Dict{String, Array{Float64, 2}}()

    r = r"\$Nodes\n(.*?)\n\$EndNodes"s
    r = match(r, cnt)[1]
    nodeinfo = split(r, '\n')[2:end]

    ps = physicalnames
    # convert physical tag to entity tag
    for p in ps
        dim = split(p)[1]
        tag = split(p)[2]
        name = strip(split(p)[3], '"')

        entity_tag = []
        if dim=="0"
            for i = 1:n
                if length(split(info[i+1]))>=6 && split(info[i+1])[6]==tag
                    push!(entity_tag, "0 $(split(info[i+1])[1]) ")
                end
            end
        elseif dim=="1"
            for i = 1:n
                if split(info[n+i+1])[9]==tag
                    push!(entity_tag, "1 $(split(info[n+i+1])[1]) ")
                end
            end
        end

        all_nodes = []
        for et in entity_tag
            # extract entity nodes
            entity_nodes = []
            k = 0
            while k<length(nodeinfo)
                k += 1
                if startswith(nodeinfo[k], et) && length(split(nodeinfo[k]))==4
                    N = parse(Int64, split(nodeinfo[k])[end])
                    for j = 1:N
                        k+=1
                    end
                    for j = 1:N
                        k += 1
                        out = map(x->parse(Float64, x), split(nodeinfo[k])[1:2])
                        push!(entity_nodes, [out[1] out[2]])
                    end
                end
            end
            entity_nodes = vcat(entity_nodes...)
            push!(all_nodes, entity_nodes)
        end
        if haskey(D, name)
            D[name] = vcat(D[name], all_nodes...)
        else
            D[name] = vcat(all_nodes...)
        end
    end
    for k in D
        D[k] = unique(D[k], dims=1)
    end
    return D
end

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

"""
    readMesh(gmshFile::String)

Reads a `gmsh` file and extracts element, coordinates and boundaries.
"""
function readMesh(gmshFile::String)
    fp = open(gmshFile);
    boundaries = Dict{String, Array}()
    physicalnames = Dict{Int64, String}()
    nodes = nothing
    elements = nothing

    line = readline(fp)
    if (line == "\$MeshFormat")
        format = readline(fp)
        @assert(readline(fp) == "\$EndMeshFormat")
    end

    line = readline(fp)
    if (line == "\$PhysicalNames")
        nphys = parse(Int64,readline(fp))
        for i = 1:nphys
            line = readline(fp)
            physicalid = parse(Int64, split(line)[2])
            physicalnames[physicalid] = split(line)[3]
            boundaries[split(line)[3]] = []
        end
        @assert(readline(fp) == "\$EndPhysicalNames")
    end

    line = readline(fp)
    if (line == "\$Nodes")
        nnodes = parse(Int64, readline(fp))
        nodes = zeros(nnodes,2)
        for i = 1:nnodes
            l = readline(fp)|>split
            nodes[i,1] = parse(Float64, l[2])
            nodes[i,2] = parse(Float64, l[3])
        end
        @assert(readline(fp) == "\$EndNodes")
    end
    # println(physicalnames)
    # error()
    line = readline(fp)
    if (line == "\$Elements")
        nelems = readline(fp)
        nelems = parse(Int64,nelems)
        elements = []
        for i = 1:nelems
            l = readline(fp)|>split
            l4 = parse(Int64, l[4])
            physicalname = physicalnames[l4]
            if startswith(physicalname,"\"Dirichlet")
                k1,k2 = parse(Int64, l[6]),parse(Int64, l[7])
                push!(boundaries[physicalname],[k1;k2])
            elseif startswith(physicalname,"\"Neumann")
                k1,k2 = parse(Int64, l[6]),parse(Int64, l[7])
                push!(boundaries[physicalname],[k1;k2])
            else
                # println(physicalnames[l[4]], startswith(physicalnames[l[4]],"Dirichlet"))
                # println(l)
                k = [parse(Int64, l[i]) for i = 6:9]
                push!(elements, k)
            end
        end
        # elem = zeros(length(elements),4)
        # for i = 1:length(elements)
        #     elem[i,:] = elements[i]
        # end
        @assert(readline(fp) == "\$EndElements")
    end
    @warn("Jobs are not finished!!! Users need to compute EBC, g, NBC, f.")
    return elements, nodes, boundaries
end

function save(file::String, domain::Domain, globaldata::GlobalData)
    @save file domain globaldata
end

function load(file::String)
    @load file domain globaldata
    return domain, globaldata
end

# domain.state, domain.fint, fext
function write_data(file::String, domain::Domain)
    file = matopen(file, "w")
    write(file, "state", domain.history["state"])
    write(file, "fext", domain.history["fext"])
    write(file, "fint", domain.history["fint"])
    write(file, "strain", domain.history["strain"])
    write(file, "stress", domain.history["stress"])
    close(file)
end

# state, fext
function read_data(file::String)
    vars = matread(file)
    #use fint, strain, stress, for debugging purpose
    vars["state"], vars["fext"]
end

# state, fext
function read_strain_stress(file::String)
    vars = matread(file)
    #use fint, strain, stress, for debugging purpose
    vars["strain"], vars["stress"]
end

# type == nn2array: read mat convert to 1D array
# type == array2nn: read 1D array convert to mat
function convert_mat(type::String, nnlayers::Array{Int64}, file1::String)
    if type == "nn2array"
       nn = matread(file1)
       nnlayerNum = length(nnlayers) - 1
       theta_indptr = ones(Int64, 2*nnlayerNum+1)
       for i = 1:nnlayerNum
       	   theta_indptr[2*i]   = theta_indptr[2*i-1] + nnlayers[i]*nnlayers[i+1]
	       theta_indptr[2*i+1] = theta_indptr[2*i] + nnlayers[i+1]
       end

       theta = zeros(Float64, theta_indptr[end]-1)
       for i = 1:nnlayerNum
           for name in nn
               if occursin("connected_" * string(i) * "backslashweights", name.first)
                   theta[theta_indptr[2*(i+1)-1]:theta_indptr[2*(i+1)]-1] = transpose(name.second)[:]
               end
  	       if occursin("connected_" * string(i) * "backslashbiases", name.first)
                   theta[theta_indptr[2*(i+1)]:theta_indptr[2*(i+1)+1]-1] = transpose(name.second)[:]
               end
  	       if occursin("connectedbackslashweights", name.first)
                   theta[theta_indptr[1]:theta_indptr[2]-1] = transpose(name.second)[:]
               end
  	       if occursin("connectedbackslashbiases", name.first)
                   theta[theta_indptr[2]:theta_indptr[3]-1] = transpose(name.second)[:]
               end
           end
       end

       return theta
   end
end

"""
    load_mesh(s::String = "holes")

Load built-in meshes from `$(abspath(joinpath(@__DIR__, "../../deps/Data/")))`.
"""
function load_mesh(s::String = "holes")
    meshdata = joinpath(@__DIR__, "../../deps/Data/$s.mat")
    data = matread(meshdata)
    data["node"], data["elem"]
end


text:
module StructuredOptimization

using LinearAlgebra
using RecursiveArrayTools
using AbstractOperators
using ProximalOperators
using ProximalAlgorithms

import ProximalAlgorithms:ForwardBackward, ZeroFPR, PANOC
export ForwardBackward, ZeroFPR, PANOC

include("syntax/syntax.jl")
include("calculus/precomposeNonlinear.jl") # TODO move to ProximalOperators?
include("arraypartition.jl") # TODO move to ProximalOperators?

# problem parsing
include("solvers/terms_extract.jl")
include("solvers/terms_properties.jl")
include("solvers/terms_splitting.jl")

# solver calls
include("solvers/solvers_options.jl")
include("solvers/build_solve.jl")
include("solvers/minimize.jl")

end

ext: jl


text:
#Demo 7: Bode Plots
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)
j = im

#Input data
#-------------------------------------------------------------------------------
fmax = 10e9
fstep = 5e6
f = collect(0:fstep:fmax)
fmin_disp = 10e6
G = 20.0
f_p1 = 50e6
f_p2 = .9e9
f_z = Inf #.8e9
f1 = .8e9

#==Helper functions
===============================================================================#
phase(x) = rad2deg(angle(x))
function roots(a,b,c)
	r = sqrt(b*b-4*a*c)
	result = (-b .+ [r, -r]) / (2a)
	return (result[1], result[2])
end
xf_o2(f, G, f_z, f_p1, f_p2) = #2nd order transfer function
	G*(1+j*(f./f_z))./( (1+j*(f./f_p1)).*(1+j*(f./f_p2)) )
function f_unitygain_f(f, G, f_z, f_p1, f_p2) #Compute unity gain frequency
	invsq(x) = (result = 1/x; return result*result)
	_z² = invsq(f_z); G² = G*G
	_p1² = invsq(f_p1); _p2² = invsq(f_p2)
	r = roots(_p1²*_p2², _p1²+_p2²-G²*_z², 1-G²)
	return sqrt(max(r...)) #Select sqrt(positive root)
end
f_3dB(f, G, f_z, f_p1, f_p2) = #Compute 3dB frequency
	f_unitygain_f(f, sqrt(2), f_z, f_p1, f_p2)

#==Calculations
===============================================================================#
X = xf_o2(f, G, f_z, f_p1, f_p2)
f0 = f_unitygain_f(f, G, f_z, f_p1, f_p2)
fBW = f_3dB(f, G, f_z, f_p1, f_p2)
phase0 = phase(xf_o2(f0, G, f_z, f_p1, f_p2))
pmargin = 180+phase0
@show pmargin

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Bode Plot")
mplot.hplot*=.6 #Bode plot looks better with wider aspect ratio
mplot.ncolumns = 1

splot = plot_mag = add(mplot, InspectDR.Plot2D)
	splot.axes = InspectDR.axes(:log10, :dB20)
	splot.ext_full = InspectDR.PExtents2D(xmin=fmin_disp, ymin=-10)
	wfrm = add(splot, f, X)
		wfrm.line = line(color=blue, width=3)
	a = splot.annotation
#		a.title = ""
		a.xlabel = "Frequency (Hz)"; a.ylabel = "Magnitude (dB)"
	splot.layout.htitle = 8

splot = plot_phase = add(mplot, InspectDR.Plot2D)
	splot.axes = InspectDR.axes(:log10, :lin)
	splot.ext_full = InspectDR.PExtents2D(xmin=fmin_disp, ymin=-180)
	wfrm = add(splot, f, phase(X))
		wfrm.line = line(color=blue, width=3)
	a = splot.annotation
#		a.title = ""
		a.xlabel = "Frequency (Hz)"; a.ylabel = "Phase (°)"
	splot.layout.htitle = 8

#Add specialized annotation
#-------------------------------------------------------------------------------
afont = InspectDR.Font(12) #Font for annotation
markerline = line(style=:dash, width=2.5)
markerline_light = line(style=:dash, width=2.5, color=RGB24(.4,.4,.4))

#Add vertical markers to both plots:
for splot in [plot_mag, plot_phase]
	splot.layout.grid = grid(vmajor=true, vminor=true, hmajor=false)
	add(splot, vmarker(f0, markerline))
	add(splot, vmarker(fBW, markerline))
end

#Add annotation to Magnitude plot:
add(plot_mag, atext("0dB", y=1, xoffset=0.5, yoffset=2/100, font=afont, align=:bc))
add(plot_mag, hmarker(1, markerline_light))

fstr = "f3dB=" * @sprintf("%.1f MHz", fBW/1e6)
add(plot_mag, atext(fstr, x=fBW, xoffset=-1/100, yoffset=.5, font=afont, angle=-90, align=:bc))

fstr = "f0=" * @sprintf("%.1f MHz", f0/1e6)
add(plot_mag, atext(fstr, x=f0, xoffset=-1/100, yoffset=.5, font=afont, angle=-90, align=:bc))

#Add annotation to Phase plot:
fstr = "PM=" * @sprintf("%.1f°", pmargin)
add(plot_phase, atext(fstr, y=(phase0-180)/2, xoffset=0.5, font=afont, align=:cc))
add(plot_phase, hmarker(phase0, markerline))

#==Render plot
===============================================================================#
gplot = display(InspectDR.GtkDisplay(), mplot)

InspectDR.write_png("export_bode.png", mplot)
#InspectDR.write_svg("export_bode.svg", mplot)
#InspectDR.write_eps("export_bode.eps", mplot)
#InspectDR.write_pdf("export_bode.pdf", mplot)

:DONE


text:
include("dmixture.jl")
include("update_Z.jl")
include("update_Z_v2.jl")
# include("repFAM/update_Z_repFAM.jl")  # this is done in Model.jl
include("update_delta.jl")
include("update_alpha.jl")
include("update_v.jl")
include("update_W.jl")
include("update_sig2.jl")
include("update_eta.jl")
include("update_lam.jl")
include("update_eps.jl")
include("update_gam.jl")
include("update_y_imputed.jl")
include("compute_loglike.jl")

function update_state!(s::State, c::Constants, d::Data, tuners::Tuners,
                       ll::Vector{Float64}, fix::Vector{Symbol},
                       use_repulsive::Bool, joint_update_Z::Bool, sb_ibp::Bool)

  # NOTE: `@doIf` is defined in "util.jl"

  # Return true if parameter (sym) is not fixed
  isRandom(sym::Symbol)::Bool = !(sym in fix)

  # Gibbs.
  if isRandom(:Z)
    if joint_update_Z
      update_Z_v2!(s, c, d, tuners, sb_ibp, use_repulsive=use_repulsive)
    else
      # Do regular updates
      update_Z!(s, c, d, sb_ibp, use_repulsive=use_repulsive)
    end
  end

  @doIf isRandom(:v)          update_v!(s, c, d, tuners, sb_ibp)
  @doIf isRandom(:alpha)      update_alpha!(s, c, d, sb_ibp)
  @doIf isRandom(:lam)        update_lam!(s, c, d)
  @doIf isRandom(:W)          update_W!(s, c, d)
  @doIf isRandom(:eps)        update_eps!(s, c, d)

  # gam update must be done after updating Z and before updating delta
  @doIf isRandom(:gam)        update_gam!(s, c, d)
  @doIf isRandom(:eta)        update_eta!(s, c, d)

  @doIf isRandom(:delta)      update_delta!(s, c, d)
  @doIf isRandom(:sig2)       update_sig2!(s, c, d)

  # Metropolis.
  @doIf isRandom(:y_imputed)  update_y_imputed!(s, c, d, tuners)

  # Compute loglikelihood.
  append!(ll, compute_loglike(s, c, d))
end


text:
using BenchmarkTools

using QSimulator

const SUITE = BenchmarkGroup()

########################################## Basic Operators #########################################
SUITE["operators"] = BenchmarkGroup()

for n = 2:5
    q = DuffingTransmon("q0", n, DuffingSpec(5, -0.2))

    SUITE["operators"]["number(q) ($n levels)"] = @benchmarkable number($q)
    SUITE["operators"]["raising(q) ($n levels)"] = @benchmarkable raising($q)
    SUITE["operators"]["raising(q, 0.123) ($n levels)"] = @benchmarkable raising($q, 0.25)
    SUITE["operators"]["lowering(q) ($n levels)"] = @benchmarkable lowering($q)
    SUITE["operators"]["lowering(q, 0.123) ($n levels)"] = @benchmarkable lowering($q, 0.25)
    SUITE["operators"]["X(q) ($n levels)"] = @benchmarkable X($q)
    SUITE["operators"]["X(q, 0.123) ($n levels)"] = @benchmarkable X($q)

end

########################################## PerturbativeTransmon ####################################
SUITE["transmon"] = BenchmarkGroup()

for n in (3, 10, 20)
    SUITE["transmon"]["perturbative_transmon_freq ($n terms)"] =
        @benchmarkable perturbative_transmon_freq(0.165, 16.0, 2.75, 0.123; num_terms = $n)
    SUITE["transmon"]["perturbative_transmon_anharm ($n terms)"] =
        @benchmarkable perturbative_transmon_anharm(0.165, 16.0, 2.75, 0.123; num_terms = $n)
    SUITE["transmon"]["perturbative_transmon_λ ($n terms)"] =
        @benchmarkable perturbative_transmon_λ(0.165, 16.0, 2.75, 0.123; num_terms = $n)
    SUITE["transmon"]["perturbative_transmon_Λ ($n terms)"] =
        @benchmarkable perturbative_transmon_Λ(0.165, 16.0, 2.75, 0.123; num_terms = $n)
end

SUITE["transmon"]["perturbative_transmon_Ueigen"] =
    @benchmarkable perturbative_transmon_Ueigen(0.15)

########################################## CompositeSystems ########################################
SUITE["embedding"] = BenchmarkGroup()

for n = 2:2:8
    r = rand(3,3)
    dims = fill(3, n)
    acting_on = [n]
    SUITE["embedding"]["QSimulator.embed(r, $acting_on, $dims)"] = @benchmarkable QSimulator.embed($r, $acting_on, $dims)
end

########################################## Basis Utilities #########################################
SUITE["basis_utils"] = BenchmarkGroup()

for n = 2:4
    basis_tuple = ((2 for i in 1:n)...,)
    SUITE["basis_utils"]["tensor_product_basis ($n qubits)"] = @benchmarkable TensorProductBasis($basis_tuple);
    b = TensorProductBasis(basis_tuple)
    SUITE["basis_utils"]["basis_states ($n qubits)"] = @benchmarkable basis_states($b);
    bs = basis_states(b)
    SUITE["basis_utils"]["vec ($n qubits)"] = @benchmarkable QSimulator.vec($bs[1]);
    SUITE["basis_utils"]["index ($n qubits)"] = @benchmarkable QSimulator.index($bs[1]);
    SUITE["basis_utils"]["getindex ($n qubits)"] = @benchmarkable getindex($b, 1);
end

######################################### Time Evolution ###########################################

SUITE["unitary"] = BenchmarkGroup()
SUITE["unitary"]["propagator"] = BenchmarkGroup()
SUITE["unitary"]["pure state"] = BenchmarkGroup()

# Free evolution of an n level transmon
prop_suite = SUITE["unitary"]["propagator"]["free evolution"] = BenchmarkGroup()
state_suite = SUITE["unitary"]["pure state"]["free evolution"] = BenchmarkGroup()
for n = 2:4
    q0 = DuffingTransmon("q0", n, DuffingSpec(5, -0.2))
    cqs = CompositeQSystem([q0])
    add_hamiltonian!(cqs, q0)
    times = range(0, stop=1.0, length=201)
    # superposition of all levels
    ψ₀ = (1/sqrt(dimension(cqs))) * ones(ComplexF64, dimension(cqs))
    prop_suite["single transmon ($n levels)"] = @benchmarkable unitary_propagator($cqs, $times);
    state_suite["single transmon ($n levels)"] = @benchmarkable unitary_state($cqs, $times, $ψ₀);
end

# Dipole interaction between a chain of fixed transmons
for n = 2:4
    qs = [DuffingTransmon("q$ct", 3, DuffingSpec(4.0 + 0.1*ct, -(0.2 + 0.01*ct))) for ct = 0:(n-1)]
    cqs = CompositeQSystem(qs)
    for q = qs
        add_hamiltonian!(cqs, hamiltonian(q), q)
    end
    for (qa,qb) = zip(qs[1:end-1], qs[2:end])
        add_hamiltonian!(cqs, 0.005*X([qa, qb]), [qa, qb])
    end
    times = range(0, stop=1.0, length=201)
    # superposition of all levels
    ψ₀ = (1/sqrt(dimension(cqs))) * ones(ComplexF64, dimension(cqs))
    prop_suite["dipole chain ($n transmons)"] = @benchmarkable unitary_propagator($cqs, $times);
    state_suite["dipole chain ($n transmons)"] = @benchmarkable unitary_state($cqs, $times, $ψ₀);
end

# Lab frame Rabi oscillations of an n level transmon
prop_suite = SUITE["unitary"]["propagator"]["rabi flops"] = BenchmarkGroup()
state_suite = SUITE["unitary"]["pure state"]["rabi flops"] = BenchmarkGroup()

qubit_freq = 5.0
nutation_freq = 0.02
for n = 2:4
    q0 = DuffingTransmon("q0", n, DuffingSpec(qubit_freq, -0.2))
    cqs = CompositeQSystem([q0]);
    add_hamiltonian!(cqs, q0)
    add_hamiltonian!(cqs, dipole_drive(q0, t -> nutation_freq*cos(2π*qubit_freq * t)), q0);
    # start in ground state
    ψ₀ = zeros(ComplexF64, n)
    ψ₀[1] = 1
    times = range(0, stop=100, length=101)
    prop_suite["$n level transmon"] = @benchmarkable unitary_propagator($cqs, $times);
    state_suite["$n level transmon"] = @benchmarkable unitary_state($cqs, $times, $ψ₀);
end

# Lab frame parametric interaction between two transmons with spectators
prop_suite = SUITE["unitary"]["propagator"]["lab frame parametric 2Q gate"] = BenchmarkGroup()
state_suite = SUITE["unitary"]["pure state"]["lab frame parametric 2Q gate"] = BenchmarkGroup()

# helper function to add flux drive
for n = 2:3
    q0 = DuffingTransmon("q0", 3, DuffingSpec(3.94015, -0.1807))
    q1 = PerturbativeTransmon("q1",  3, TransmonSpec(0.172, 16.4, 0.55))

    # add fixed frequency spectators
    spectator_qs = [DuffingTransmon("q$ct", 3, DuffingSpec(4.0 + 0.1*ct, -(0.2 + 0.01*ct))) for ct = 2:(n-1)]

    # should get an iSWAP interaction at ≈ 122 MHz with drive amplitude 0.323 Φ₀
    mod_freq = 122.1/1e3
    amp = 0.323
    all_qs = [q0, q1, spectator_qs...]
    cqs = CompositeQSystem(all_qs)
    add_hamiltonian!(cqs, q0)
    add_hamiltonian!(cqs, parametric_drive(q1, t -> sin(2π*mod_freq*t)), q1)
    add_hamiltonian!(cqs, 0.006*X([q0, q1]), [q0, q1])

    # add hamiltoians for spectators coupled to tunable transmon
    for q = all_qs[3:end]
        add_hamiltonian!(cqs, q)
        add_hamiltonian!(cqs, 0.006*X([q, q1]), [q, q1])
    end

    times = 0:200
    ψ₀ = ComplexF64[0.0; 1.0; 0.0] ⊗ ComplexF64[1.0; 0.0; 0.0] # start in 10 state
    for ct = 1:n-2
        ψ₀ = ψ₀ ⊗ ComplexF64[1.0; 0.0; 0.0]
    end

    prop_suite["$n transmons"] = @benchmarkable unitary_propagator($cqs, $times)
    state_suite["$n transmons"] = @benchmarkable unitary_state($cqs, $times, $ψ₀)
end

# Rotating frame parametric interaction between two transmons with spectators
prop_suite = SUITE["unitary"]["propagator"]["rotating frame parametric 2Q gate"] = BenchmarkGroup()
state_suite = SUITE["unitary"]["pure state"]["rotating frame parametric 2Q gate"] = BenchmarkGroup()

for n = 2:4

    q0 = DuffingTransmon("q0", 3, DuffingSpec(3.94015, -0.1807))
    q1 = PerturbativeTransmon("q1",  3, TransmonSpec(0.172, 16.4, 0.55))

    # add fixed frequency spectators
    spectator_qs = [DuffingTransmon("q$ct", 3, DuffingSpec(4.0 + 0.1*ct, -(0.2 + 0.01*ct))) for ct = 2:(n-1)]

    # should get an CZ02 interaction at ≈ 115 MHz
    freq = 115.5/1e3
    amp = 0.245

    all_qs = [q0, q1, spectator_qs...]
    cqs = CompositeQSystem(all_qs)
    add_hamiltonian!(cqs, q0)
    # add rotating frame Hamiltonian shifts
    add_hamiltonian!(cqs, -spec(q0).frequency*number(q0), q0)
    q1_freq = hamiltonian(q1, 0.0)[2,2]
    add_hamiltonian!(cqs, -q1_freq*number(q1), q1)

    diff_freq = spec(q0).frequency - q1_freq
    # time dependent flip flop interaction
    add_hamiltonian!(cqs, t -> 0.006*.5*X_Y([q0, q1], [diff_freq*t, 0.0]), [q0,q1])
    add_hamiltonian!(cqs, parametric_drive(q1, t -> amp*sin(2π*freq*t)), q1)

    # add hamiltoians for spectators coupled to tunable transmon
    for q = all_qs[3:end]
        add_hamiltonian!(cqs, q)
        # add rotating frame Hamiltonian shifts
        add_hamiltonian!(cqs, -spec(q).frequency*number(q), q)
        diff_freq = q1_freq - spec(q).frequency
        add_hamiltonian!(cqs, t -> 0.006*.5*X_Y([q1, q], [diff_freq*t, 0.0]), [q1, q])
    end

    times = 0:200
    ψ₀ = ComplexF64[0.0; 1.0; 0.0] ⊗ ComplexF64[0.0; 1.0; 0.0] # start in 11 state
    for ct = 1:n-2
        ψ₀ = ψ₀ ⊗ ComplexF64[1.0; 0.0; 0.0]
    end

    prop_suite["$n transmons"] = @benchmarkable unitary_propagator($cqs, $times)
    state_suite["$n transmons"] = @benchmarkable unitary_state($cqs, $times, $ψ₀);
end


text:
using Test
import ArchGDAL; const AG = ArchGDAL
import GeoFormatTypes; const GFT = GeoFormatTypes

# Tests high level convert methods

@testset "convert point format" begin
    point = AG.createpoint(100, 70)
    json = convert(GFT.GeoJSON, point)
    kml = convert(GFT.KML, point)
    gml = convert(GFT.GML, point)
    wkb = convert(GFT.WellKnownBinary, point)
    wkt = convert(GFT.WellKnownText, point)
    @test json.val == AG.toJSON(point)
    @test kml.val == AG.toKML(point)
    @test gml.val == AG.toGML(point)
    @test wkb.val == AG.toWKB(point)
    @test wkt.val == AG.toWKT(point)
    @test convert(GFT.GeoJSON, json) == convert(GFT.GeoJSON, wkb) ==
          convert(GFT.GeoJSON, wkt) == convert(GFT.GeoJSON, gml) == json
    @test convert(GFT.KML, gml) == convert(GFT.KML, wkt)
end

@testset "convert crs format" begin
    proj4326 = GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs")
    @test convert(GFT.ProjString, GFT.CRS(),
                  convert(GFT.WellKnownText, GFT.CRS(),
                          convert(GFT.ESRIWellKnownText, GFT.CRS(),
                                  GFT.EPSG(4326)))) == proj4326
end


text:
module FixedPointSinCosApproximations

    using LoopVectorization, DocStringExtensions
    using LoopVectorization: Vec

    export fpsin,
        fpcos,
        fpsincos

    const VInt16 = Union{Int16,Vec{<:Any,Int16}}
    const VInt32 = Union{Int32,Vec{<:Any,Int32}}

    include("second_order.jl")
    include("third_order.jl")
    include("fourth_order.jl")
    include("fifth_order.jl")
    include("sixth_order.jl")
    include("approximations.jl")

end # module


text:
module NeuralNetworks

# TODO: Add exports

using LinearAlgebra, Zygote, ForwardDiff

##
## Data types
##

"""
A single layer of a neural network.
"""
struct Layer
    # Each row corresponds to a neuron in the previous layer, including bias
    # term.
    # Each column corresponds to a neuron in this layer, not including the bias
    # term.
    # Biases for a given layer are implicitly added by the next layer, and so do
    # not show up in dimensionality of the same layer. Thus, the dimensionality
    # of this matrix is (num_prev_outputs + 1) * (num_curr_outputs).
    weights::AbstractMatrix{<:Number}
    # Vector-valued function (f: R^d -> R^d) mapping the neurons in this layer
    # to their activations.
    activation_fn::Function
end

"""
A complete neural network.
"""
struct NeuralNetwork
    # Layers of the network, applied in order.
    layers::AbstractVector{<:Layer}
    # Single loss function f: (R^n * R^n) -> R applied to the network's final
    # layer output.
    loss_fn::Function
end

"""
A data structure containing the values computed by a single layer during the
forward-pass.
"""
struct LayerForwardPass
    # Inputs from the previous layer. Includes bias.
    inputs::AbstractVector{<:Number}
    # Sums of each neuron before activation; same dimensionality as layer
    # output.
    sums::AbstractVector{<:Number}
    # Activation of each neuron; same dimensionality as layer output.
    activations::AbstractVector{<:Number}
end

"""
A data structure containing the partial derivatives of the loss with respect to
a given layer during the backward-pass.
"""
struct LayerBackwardPass
    # Partial derivatives w.r.t. pre-activation sums.
    dL_dS::AbstractVector{<:Number}
    # Partial derivatives w.r.t. each weight.
    dL_dW::AbstractMatrix{<:Number}
    # Partial derivatives w.r.t. inputs. This gets used by the previous layer in
    # backprop (these inputs are exactly the outputs of the previous layer).
    dL_dI::AbstractVector{<:Number}
end

function Base.show(io::IO, layer::Layer)
    print("  $(size(layer.weights, 2)) neurons\n")
    print("  Weights " *
          "[$(size(layer.weights, 1))x$(size(layer.weights, 2))]: " *
          "$(layer.weights)\n")
end

function Base.show(io::IO, nn::NeuralNetwork)
    num_layers = length(nn.layers)
    print("Neural network with $(num_layers) layers.\n")
    print("Input dimensionality: $(in_ndims_sans_bias(nn.layers[1]))\n")
    print("Output dimensionality: $(out_ndims(nn.layers[num_layers]))\n")
    for l in 1:num_layers
        print("Layer $(l):\n")
        print(nn.layers[l])
    end
end

function Base.show(io::IO, layer_forward_pass::LayerForwardPass)
    print("Input vector: $(layer_forward_pass.inputs), " *
          "$(size(layer_forward_pass.inputs))\n")
    print("Intermediate sums: $(layer_forward_pass.sums), " *
          "$(size(layer_forward_pass.sums))\n")
    print("Output vector: $(layer_forward_pass.activations), " *
          "$(size(layer_forward_pass.activations))\n")
end

function Base.show(io::IO, layer_backward_pass::LayerBackwardPass)
    print("dL_dS: $(layer_backward_pass.dL_dS), " *
          "$(size(layer_backward_pass.dL_dS))\n")
    print("dL_dW: $(layer_backward_pass.dL_dW), " *
          "$(size(layer_backward_pass.dL_dW))\n")
    print("dL_dI: $(layer_backward_pass.dL_dI), " *
          "$(size(layer_backward_pass.dL_dI))\n")
end

"""
Data structure containing loss stats for each epoch of training.
"""
struct TrainingStats
    losses::AbstractVector
end

##
## Constructors
##

"""
    create_layer(in_ndims_sans_bias, out_ndims, activation_fn)

Creates a single layer of a neural network with the given dimensions and
activation function, and randomly assigned weights.

In this package, a "layer" is composed of (1) weights, which are applied to the
previous layer, and (2) activation functions, which are applied after weights.
"""
function create_layer(in_ndims_sans_bias::Integer,
                      out_ndims::Integer,
                      activation_fn::Function)::Layer
    weights = randn((in_ndims_sans_bias+1, out_ndims)) ./ sqrt(in_ndims_sans_bias)
    return Layer(weights, activation_fn)
end

"""
    validate_layer(layer, num_inputs_sans_bias, num_outputs)

Asserts that the layer parameters have the desired dimensionality, not including
the bias.
"""
function validate_layer(layer::Layer,
                        num_inputs_sans_bias::Integer,
                        num_outputs::Integer)
    @assert in_ndims_sans_bias(layer) == num_inputs_sans_bias
    @assert out_ndims(layer) == num_outputs
end

"""
    compose_layers(input_ndims, output_ndims, loss_fn, layers...)

Creates a neural network from the supplied layers.
"""
function compose_layers(input_ndims::Integer,
                        output_ndims::Integer,
                        loss_fn::Function,
                        layers...)::NeuralNetwork
    num_layers = length(layers)
    nn_layers = Vector{Layer}(undef, num_layers)

    # Verify layers are dimensionally compatible.
    for l in 1:num_layers
        expected_in_ndims_sans_bias = (l == 1 ?
                                       input_ndims :
                                       out_ndims(layers[l-1]))
        expected_out_ndims = (l == num_layers ?
                              output_ndims :
                              in_ndims_sans_bias(layers[l+1]))
        validate_layer(layers[l],
                       expected_in_ndims_sans_bias,
                       expected_out_ndims)
    end

    return NeuralNetwork([layers...], loss_fn)
end

"""
    create_nn(layer_dims, activation_fn, loss_fn,
              network_input_dims, network_output_dims)

Creates a neural network by building layers with the specified dimensionality,
`activation_fn`, and `loss_fn`.

- `layer_dims::AbstractVector`: a vector of integers, where the `n`-th value
corresponds to the number of neurons in the `n`-th layer of the network.
- `activation_fn`: the activation function used for all layers in the network.
- `loss_fn`: a single loss function applied at the end of a network, producing a
  scalar loss value.
- `network_input_dims`: the dimensionality of inputs to the network.
- `network_output_dims`: the dimensionality of outputs of the network.

TODO: This function is convenient, but it's probably better to construct and
compose layers individually.
"""
function create_nn(layer_dims::AbstractVector{<:Integer},
                   activation_fn::Function,
                   loss_fn::Function,
                   network_input_ndims::Integer,
                   network_output_ndims::Integer)::NeuralNetwork
    num_layers = length(layer_dims)
    layers = Vector{Layer}(undef, num_layers)
    for l in 1:num_layers
        in_ndims_sans_bias = (l == 1 ?
                              network_input_ndims :
                              layer_dims[l-1])
        out_ndims = layer_dims[l]
        layers[l] = create_layer(in_ndims_sans_bias, out_ndims,
                                 activation_fn)
    end

    return compose_layers(network_input_ndims, network_output_ndims,
                          loss_fn, layers...)
end

##
## Training and inference
##

"""
    run_forward_pass(layer, inputs)

Computes the activations of `layer` on `inputs`, returning the cached results.
"""
function run_forward_pass(layer::Layer,
                          inputs::AbstractVector{<:Number})::LayerForwardPass
    @assert ndims(inputs) == 1
    @assert length(inputs) == in_ndims_sans_bias(layer)

    # Append the bias term.
    inputs = copy(inputs)
    push!(inputs, 1)

    # Perform matrix multiplication `weights`^T * `input`.
    # For the `weights` matrix, the rows index the inputs and the columns index
    # the outputs. `weights` is transposed before multiplication so that the
    # rows index outputs, and the result is a column vector with an entry for
    # each neuron output in this layer.
    sums = transpose(layer.weights) * inputs

    # Apply each neuron's activation function. The activation function is
    # vector-valued, so the result is a vector of the same dimensionality.
    activations = layer.activation_fn(sums)

    return LayerForwardPass(inputs, sums, activations)
end

"""
    compute_loss_gradient(net_output, label, loss_fn)

Computes the gradient of `loss_fn` parameterized by `label`, with respect to
`net_output`, evaluated at `net_output`.

- `loss_fn`: a function f: (R^n * R^n) -> R

Returns a vector with the size of the output layer.
"""
function compute_loss_gradient(net_output::AbstractVector{<:Number},
                               label::AbstractVector{<:Number},
                               loss_fn::Function)::AbstractVector{<:Number}
    @assert length(net_output) == length(label)

    partial_loss(x) = loss_fn(x, label)
    dL_dO = gradient(partial_loss, net_output)[1]

    return dL_dO
end

"""
    run_backward_pass(layer, forward_pass, dL_dO)

Does backprop to compute the loss gradients with respect to the components of a
single layer, using the values from the `forward_pass`.

Returns a LayerBackwardPass containing the gradients, without mutating the
original layer.
"""
function run_backward_pass(layer::Layer,
                           forward_pass::LayerForwardPass,
                           dL_dO::AbstractVector{<:Number})::LayerBackwardPass
    @assert length(dL_dO) == out_ndims(forward_pass)

    # dL_dS|S = dL_dO|O * dO_dS|S
    # dO_dS is just the derivative (Jacobian) of the activation function.
    # (d_out * d_out)
    dO_dS = ForwardDiff.jacobian(layer.activation_fn, forward_pass.sums)

    # dL_dS = dL_dO * dO_dS
    # For activation functions which are elementwise-only, this is actually
    # equivalent to a vector-vector Hadamard product, because the Jacobian is
    # diagonal.
    # (d_out)
    dL_dS = transpose(dO_dS) * dL_dO

    # dL_dW|W = dL_dS|S * dS_dW|W
    # dS_dW for each weight is the value of the output which is scaled by that
    # weight.
    # (d_in+1 * d_out)
    dL_dW = forward_pass.inputs * transpose(dL_dS)

    # dL_dI|I = sum(dL_dS|S * dS_dI|I)
    # dS_dI is just the applied weight; the matrix multiplication applies the
    # vectorized sum.
    # (d_in+1)
    dL_dI = layer.weights * dL_dS

    # Note that the last element in dL_dW corresponds to bias, and so should be
    # omitted by the previous layer during backprop.
    return LayerBackwardPass(dL_dS, dL_dW, dL_dI)
end

"""
    update_layer_weights!(layer, backward_pass, num_samples, learning_rate)

Updates weights according to backprop gradient calculation by applying values
from a previously computed `backward_pass`.
"""
function update_layer_weights!(layer::Layer,
                               backward_pass::LayerBackwardPass,
                               num_samples::Number,
                               learning_rate::Number)
    @assert size(layer.weights) == size(backward_pass.dL_dW)

    layer.weights .-= (learning_rate .* backward_pass.dL_dW ./ num_samples)
end

"""
    train!(nn, samples, labels, learning_rate, batch_size)
"""
function train!(nn::NeuralNetwork,
                samples::AbstractVector,
                labels::AbstractVector,
                learning_rate::Number,
                batch_size::Integer)::TrainingStats
    @assert length(samples) == length(labels)

    num_layers = length(nn.layers)

    # Divide the data into batch_size chunks.
    sample_batches = Iterators.partition(samples, batch_size) |> collect
    label_batches = Iterators.partition(labels, batch_size) |> collect

    num_passes = length(sample_batches)

    losses = Vector(undef, length(sample_batches))

    for (batch_index, input_batch) in enumerate(sample_batches)
        label_batch = label_batches[batch_index]
        num_samples = length(input_batch)

        # Record the results of forward and backward passes for each layer.
        # These are only used for backprop and are overwritten on every batch.
        # Rows index samples; columns index layers.
        forward_passes = Matrix{LayerForwardPass}(undef, (num_samples, num_layers))
        backward_passes = Matrix{LayerBackwardPass}(undef, (num_samples, num_layers))

        batch_loss = 0.0

        # Compute updates for the batch.
        for (point_index, input_point) in enumerate(input_batch)
            point_label = label_batch[point_index]

            # Run forward-pass.
            for l in 1:num_layers
                forward_passes[point_index, l] = run_forward_pass(nn.layers[l], input_point)
                input_point = forward_passes[point_index, l].activations
            end

            # Compute loss.
            batch_loss += nn.loss_fn(input_point, point_label)

            # Backprop.
            loss_grad = compute_loss_gradient(input_point, point_label, nn.loss_fn)

            for l in reverse(1:num_layers)
                backward_passes[point_index, l] =
                    run_backward_pass(nn.layers[l],
                                      forward_passes[point_index, l],
                                      l == num_layers ?
                                      loss_grad :
                                      backward_passes[point_index, l+1].dL_dI[1:end-1])
            end
        end

        # Update weights. The weight updates have already been computed, so this
        # can be done in any order.
        for i in 1:num_samples, l in 1:num_layers
            update_layer_weights!(nn.layers[l],
                                  backward_passes[i, l],
                                  num_samples,
                                  learning_rate)
        end

        batch_loss /= num_samples

        losses[batch_index] = batch_loss
        if (batch_index % 1000 == 0)
            print("Iteration $batch_index / $(length(sample_batches)) : " *
                  "loss=$(batch_loss)\n")
        end
    end

    return TrainingStats(losses)
end

"""
    predict(nn, inputs)

Computes the end-to-end application of `nn` on `inputs`.
"""
function predict(nn::NeuralNetwork, inputs::AbstractVector{<:Number})
    for layer in nn.layers
        inputs = run_forward_pass(layer, inputs).activations
    end
    return inputs
end

##
## Vectorized implementations
##

# It's not ideal that the implementation is forked, but it's rather challenging
# to have a single implementation that handles both vectorized and
# non-vectorized updates. In particular, the vectorized code is significantly
# more complex, so is not optimized for clarity.

"""
Vectorized forward-pass data structure.
"""
struct LayerForwardPassVectorized
    # (n * d_in+1)
    inputs::AbstractMatrix{<:Number}
    # (n * d_out)
    sums::AbstractMatrix{<:Number}
    # (n * d_out)
    activations::AbstractMatrix{<:Number}
end

"""
Vectorized backward-pass data structure.
"""
struct LayerBackwardPassVectorized
    # Allocate space for the Jacobian of the activation function. The partial
    # derivatives can't be hard-coded since activation functions may vary, and
    # each sample requires a (d * d) matrix, since the activation function can
    # map a vector to a vector.
    # (n * d_out * d_out)
    dO_dS::AbstractArray{<:Number, 3}
    # (n * d_out)
    dL_dS::AbstractMatrix{<:Number}
    # Since only the total gradient update is needed, the n-dimension is collapsed.
    # Otherwise, this would be a 3D (n * d_in+1 + d_out) array.
    # (d_in+1 * d_out)
    dL_dW::AbstractMatrix{<:Number}
    # (n * d_in+1)
    dL_dI::AbstractMatrix{<:Number}
end

"""
    run_forward_pass_vectorized!(layer, inputs, num_samples, batch_size,
                                 forward_pass_out)

Vectorized forward-pass.

Content of extra rows when `num_samples` < `batch_size` is undefined.

- `inputs`: (n * d_in) matrix of samples.
- `num_samples`: n_s, number of input samples.
- `batch_size`: n, batch size, corresponding to the number of rows preallocated
  for various matrices; may be larger than `num_samples` for incomplete batches.
"""
function run_forward_pass_vectorized!(layer::Layer,
                                      inputs::AbstractMatrix{<:Number},
                                      num_samples::Integer,
                                      batch_size::Integer,
                                      forward_pass_out::LayerForwardPassVectorized)
    @assert ndims(inputs) == 2
    @assert size(inputs, 2) == in_ndims_sans_bias(layer)

    # Append the bias column (extra dimension for each sample).
    forward_pass_out.inputs .= [inputs ones(batch_size,1)]

    # sums = biased_inputs * layer.weights
    # (n * d_in+1) * (d_in+1 * d_out)
    mul!(forward_pass_out.sums, forward_pass_out.inputs, layer.weights)

    # 1/ One-liner alternative for multiple activation functions:
    #
    # activations = mapslices(row -> map.(layer.activation_fns, row), sums, dims=2)
    #
    # `mapslices(f, m, dims=2)` broadcasts f to each row of m.
    # `map.(f, v)` performs elementwise function application for equal sized
    # vectors of functions and elements.
    # The net result is that the dth function is applied to every element in the
    # dth column.
    #
    # 2/ Straightforward, but optimized, version:
    #
    # Threads.@threads for j in 1:out_ndims(layer)
    #     map!(layer.activation_fns[j],
    #          view(forward_pass_out.activations, 1:num_samples, j),
    #          view(forward_pass_out.sums, 1:num_samples, j))
    # end

    # Skip activation computation for the last (batch_size - num_samples) rows.
    # Poor cache access pattern.
    Threads.@threads for i in 1:num_samples
        view(forward_pass_out.activations, i, :) .=
            layer.activation_fn(view(forward_pass_out.sums, i, :))
    end
end

"""
    compute_loss_gradient_vectorized!(net_outputs, labels, loss_fn,
                                      num_samples, batch_size,
                                      loss_grad_out)

Vectorized gradient computation.

Content of extra rows when n_s < n is undefined.

- `net_outputs`: (n * d) matrix of output predictions (rows are samples).
- `labels`: (n_s * d) matrix of labels (rows are samples); note that there may
  be a size mismatch between n and n_s if the batch is incomplete.
- `loss_fn`: a function that computes a scalar loss from two matrices.
- `num_samples`: n_s, number of input samples.
- `batch_size`: n, batch size, corresponding to the number of rows preallocated
  for various matrices; may be larger than `num_samples` for incomplete batches.

Returns (n * d) matrix of gradients.

This function may return a different (scaled) result compared to
`compute_loss_gradient`. When the values disagree, the non-vectorized version
should be treated as the reference implementation.
"""
function compute_loss_gradient_vectorized!(net_outputs::AbstractMatrix{<:Number},
                                           labels::AbstractMatrix{<:Number},
                                           loss_fn::Function,
                                           num_samples::Integer,
                                           batch_size::Integer,
                                           loss_grad_out::AbstractMatrix{<:Number})
    @assert size(net_outputs, 2) == size(labels, 2)

    partial_loss(xs) = loss_fn(xs, labels)

    # Compute dL_dO.
    if num_samples == batch_size
        loss_grad_out .= gradient(partial_loss, net_outputs)[1]
    else
        view(loss_grad_out, 1:num_samples, :) .=
            gradient(partial_loss, view(net_outputs, 1:num_samples, :))[1]
    end
end

"""
    run_backward_pass_vectorized!(layer, forward_pass, dL_dO, num_samples, batch_size,
                                  backward_pass_out)

Vectorized backward-pass.

Content of extra rows when `num_samples` < `batch_size` is undefined.

- `dL_dO`: (n * d_out) matrix of gradients.
- `num_samples`: n_s, number of input samples.
- `batch_size`: n, batch size, corresponding to the number of rows preallocated
  for various matrices; may be larger than `num_samples` for incomplete batches.
"""
function run_backward_pass_vectorized!(layer::Layer,
                                       forward_pass::LayerForwardPassVectorized,
                                       dL_dO::AbstractMatrix{<:Number},
                                       num_samples::Integer,
                                       batch_size::Integer,
                                       backward_pass_out::LayerBackwardPassVectorized)
    @assert size(dL_dO, 2) == out_ndims(forward_pass)

    # 1/ One-liner alternative for multiple activation functions:
    #
    # Slice by row, and for each, do element-wise activation function
    # application.
    #
    # dO_dS = mapslices(row -> map(gradient, layer.activation_fns, row),
    #                   forward_pass.sums,
    #                   dims=2) |>
    #     vec -> map(t -> t[1], vec)
    #
    # 2/ Optimized version:
    #
    # Threads.@threads for j in 1:length(layer.activation_fns)
    #     for i in 1:num_samples
    #         backward_pass_out.dL_dS[i,j] = gradient(layer.activation_fns[j],
    #                                                 forward_pass.sums[i,j])[1]
    #     end
    # end

    # Skip the last (batch_size - num_samples) rows.
    # Poor cache access pattern.
    # (n * d_out * d_out)
    Threads.@threads for i in 1:num_samples
        ForwardDiff.jacobian!(view(backward_pass_out.dO_dS, i, :, :),
                              layer.activation_fn,
                              view(forward_pass.sums, i, :))
    end

    # dL_dS = dL_dO * dO_dS
    # (n * d_out)
    # Poor cache access pattern.
    Threads.@threads for i in 1:num_samples
        mul!(view(backward_pass_out.dL_dS, i, :),
             transpose(view(backward_pass_out.dO_dS, i, :, :)),
             view(dL_dO, i, :))
    end

    # 1/ One-liner for 3D solution (keep each sample, (n * d_in+1 * d_out))
    # Slice by row (each row is a (d_in+1 * d_out) matrix!), and for each
    # compute the outer product which computes the weight gradients for the
    # given sample.
    #
    # dL_dW = permutedims(
    #     cat([forward_pass.inputs[i,:] * transpose(dL_dS[i,:]) for i in 1:num_samples]...,
    #         dims=3),
    #     (3,1,2))

    # 2/ Vectorized 3D, (n * d_in * d_out)
    #
    # dL_dW = Array{Float64, 3}(undef, (num_samples, # n
    #                                   size(forward_pass.inputs, 2), # d_in+1
    #                                   size(dL_dS, 2))) # d_out
    # for i in 1:num_samples
    #     # Poor cache locality
    #     dL_dW[i,:,:] .= forward_pass.inputs[i,:] * transpose(dL_dS[i,:])
    # end

    # inputs are (n * d_in+1), dL_dS is (n * d_out).
    # (d_in * d_out)
    fill!(backward_pass_out.dL_dW, 0.0)
    # TODO: Parallelize this with thread-local memory, and combine at end
    # Skip the last (batch_size - num_samples) rows.
    for i in 1:num_samples
        # dL_dW .+= forward_pass.inputs[i,:] * transpose(dL_dS[i,:])
        BLAS.axpy!(1,
                   view(forward_pass.inputs,i,:) *
                   transpose(view(backward_pass_out.dL_dS,i,:)),
                   backward_pass_out.dL_dW)
    end

    # Weights are (d_in+1, d_out)
    # (n * d_in+1)
    # dL_dI = dL_dS * transpose(layer.weights)
    # Uses BLAS to avoid allocation for the transpose.
    BLAS.gemm!('N',
               'T',
               1.0,
               backward_pass_out.dL_dS,
               layer.weights,
               0.0,
               backward_pass_out.dL_dI)
end

"""
    update_layer_weights_vectorized!(layer, backward_pass, learning_rate)

Vectorized layer weight update
"""
function update_layer_weights_vectorized!(layer::Layer,
                                          backward_pass::LayerBackwardPassVectorized,
                                          learning_rate::Number)
    @assert size(layer.weights) == size(backward_pass.dL_dW)

    # One-liner for n-dimension collapse:
    # layer.weights .-= (learning_rate .*
    #                    dropdims(sum(backward_pass.dL_dW, dims=1), dims=1))

    BLAS.axpy!(-learning_rate, backward_pass.dL_dW, layer.weights)
end

"""
    train_vectorized!(nn, samples, labels, learning_rate, batch_size)

Vectorized gradient descent.

TODO: The conditional logic involved in properly handling partial batches causes
a noticeable performance regression.
"""
function train_vectorized!(nn::NeuralNetwork,
                           samples::AbstractVector,
                           labels::AbstractVector,
                           learning_rate::Number,
                           batch_size::Integer)::TrainingStats
    @assert length(samples) == length(labels)

    num_layers = length(nn.layers)

    # Divide the data into batch_size chunks.
    sample_batches = Iterators.partition(samples, batch_size) |>
        batches -> map(batch -> transpose(hcat(batch...)), batches) |>
        collect
    label_batches = Iterators.partition(labels, batch_size) |>
        batches -> map(batch -> transpose(hcat(batch...)), batches) |>
        collect

    num_passes = length(sample_batches)

    # Preallocate memory for each batch iteration.
    forward_passes = Vector{LayerForwardPassVectorized}(undef, num_layers)
    for l in 1:num_layers
        inputs = zeros(batch_size, in_ndims_with_bias(nn.layers[l]))
        sums = zeros(batch_size, out_ndims(nn.layers[l]))
        activations = zeros(batch_size, out_ndims(nn.layers[l]))
        forward_passes[l] = LayerForwardPassVectorized(inputs, sums, activations)
    end

    backward_passes = Vector{LayerBackwardPassVectorized}(undef, num_layers)
    for l in 1:num_layers
        dO_dS = zeros(batch_size, out_ndims(nn.layers[l]), out_ndims(nn.layers[l]))
        dL_dS = zeros(batch_size, out_ndims(nn.layers[l]))
        dL_dW = zeros(in_ndims_with_bias(nn.layers[l]),
                      out_ndims(nn.layers[l]))
        dL_dI = zeros(batch_size, in_ndims_with_bias(nn.layers[l]))
        backward_passes[l] = LayerBackwardPassVectorized(dO_dS, dL_dS, dL_dW, dL_dI)
    end

    loss_grad = zeros(batch_size, out_ndims(nn.layers[end]))

    losses = Vector(undef, num_passes)

    for (index, input_batch) in enumerate(sample_batches)
        num_samples = size(input_batch, 1)

        # Individual subroutines expect matrices with size batch_size, but
        # otherwise correctly handle mismatches between num_samples and
        # batch_size.
        if (num_samples != batch_size)
            print("Incomplete batch; $num_samples / $batch_size\n")
            input_batch = vcat(input_batch,
                               zeros(batch_size - num_samples,
                                     size(input_batch, 2)))
        end

        # Run forward-pass.
        for l in 1:num_layers
            run_forward_pass_vectorized!(nn.layers[l],
                                         l == 1 ?
                                             input_batch :
                                             forward_passes[l-1].activations,
                                         num_samples,
                                         batch_size,
                                         forward_passes[l])
        end

        # Compute loss.
        # Avoid taking a view unless there isn't a full batch.
        loss = (num_samples == batch_size) ?
            nn.loss_fn(forward_passes[num_layers].activations,
                       label_batches[index]) :
            nn.loss_fn(
                view(forward_passes[num_layers].activations,1:num_samples),
                label_batches[index])

        # Backprop.
        compute_loss_gradient_vectorized!(forward_passes[num_layers].activations,
                                          label_batches[index],
                                          nn.loss_fn,
                                          num_samples,
                                          batch_size,
                                          loss_grad)

        for l in reverse(1:num_layers)
            run_backward_pass_vectorized!(nn.layers[l],
                                          forward_passes[l],
                                          l == num_layers ?
                                          loss_grad :
                                          # Chop off the bias column.
                                          view(backward_passes[l+1].dL_dI,
                                               :, 1:out_ndims(nn.layers[l])),
                                          num_samples,
                                          batch_size,
                                          backward_passes[l])
        end

        # Update weights. The weight updates have already been computed, so this
        # can be done in any order.
        # This could be multithreaded, but the network would need to be very
        # deep for the parallelism to outweigh the thread overhead, so it's not
        # recommended.
        for l in 1:num_layers
            update_layer_weights_vectorized!(nn.layers[l],
                                             backward_passes[l],
                                             learning_rate)
        end

        losses[index] = loss
        if (index % 1000 == 0)
            print("Iteration $index / $(num_passes) : loss=$(loss)\n")
        end
    end

    return TrainingStats(losses)
end

##
## Miscellaneous
##

"""
    in_ndims_sans_bias(layer)

Returns the number of inputs to `layer`, not including the bias term (in other
words, only the inputs that come from the previous layer).
"""
function in_ndims_sans_bias(layer::Layer)
    return size(layer.weights, 1) - 1
end

"""
    in_ndims_sans_bias(layer)

Returns the number of inputs to `layer`, not including the bias term (in other
words, only the inputs that come from the previous layer).
"""
function in_ndims_sans_bias(layer::LayerBackwardPass)
    @assert length(dL_dI) == size(layer.dL_dW, 1) - 1
    return length(dL_dI)
end

"""
    in_ndims_with_bias(layer)

Returns the number of inputs to `layer`, including the bias term.
"""
function in_ndims_with_bias(layer::Layer)
    return size(layer.weights, 1)
end

"""
    out_ndims(layer)

Returns the number of outputs from `layer`.
"""
function out_ndims(layer::Layer)
    return size(layer.weights, 2)
end

"""
    out_ndims(layer)

Returns the number of outputs from `layer`.
"""
function out_ndims(layer::LayerForwardPass)
    @assert length(layer.sums) == length(layer.activations)
    return length(layer.activations)
end

"""
    out_ndims(layer)

Returns the number of outputs from `layer`.
"""
function out_ndims(layer::LayerForwardPassVectorized)
    @assert size(layer.sums) == size(layer.activations)

    return size(layer.activations, 2)
end

##
## Buyer beware
##

"""
DEPRECATED

    compute_loss(net_output, label, loss_fn)

Computes the loss between `label` and `net_output` using `loss_fn`. Despite the
loss being a scalar, this transforms it into a singleton array for consistency
with the multi-objective function.

- `loss_fn`: a function f: (R^n * R^n) -> R
"""
function compute_loss(net_output::AbstractVector{<:Number},
                      label::AbstractVector{<:Number},
                      loss_fn::Function)::AbstractVector{<:Number}
    return [loss_fn(label, net_output)]
end

"""
DEPRECATED

    compute_multi_objective_loss_gradient(net_output, label, loss_fn)

Computes the gradient of `loss_fn` parameterized by `label`, with respect to
`net_output`, evaluated at `net_output`.

You probably don't want to do this (this package does not provide a good
solution to multi-objective optimization problems).

- `loss_fn`: a function f: (R^n * R^n) -> R^n
"""
function compute_multi_objective_loss_gradient(net_output::AbstractVector{<:Number},
                                               label::AbstractVector{<:Number},
                                               loss_fn::Function)::AbstractVector{<:Number}
    @assert length(net_output) == length(label)

    # Vector of partially-applied loss functions, where each entry lossmap[i]
    # contains f(x) = loss_fn(label[i], x)
    lossmap = map(y -> (x -> loss_fn(y, x)), label)

    # Take the gradient of each loss function with respect to each output,
    # evaluated at the outputs.
    # dL_dO is just the derivative of the loss function.
    dL_dO = map(gradient, lossmap, net_output) |>
        vec -> map(t -> t[1], vec)

    return dL_dO
end

"""
DEPRECATED

    compute_multi_objective_loss(net_output, label, loss_fn)

Computes the loss between `label` and `net_output` using `loss_fn`.

You probably don't want to do this (this package does not provide a good
solution to multi-objective optimization problems).

- `loss_fn`: a function f: (R^n * R^n) -> R^n
"""
function compute_multi_objective_loss(net_output::AbstractVector{<:Number},
                                      label::AbstractVector{<:Number},
                                      loss_fn::Function)::AbstractVector{<:Number}
    # Vector of partially-applied loss functions, where each entry lossmap[i]
    # contains f(x) = loss_fn(label[i], x)
    lossmap = map(y -> (x -> loss_fn(y, x)), label)

    # Returns a vector loss.
    return map.(lossmap, net_output)
end

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
# ---
# title: 460. LFU Cache
# id: problem460
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Design
# link: <https://leetcode.com/problems/lfu-cache/description/>
# hidden: true
# ---
#
# Design and implement a data structure for [Least Frequently Used
# (LFU)](https://en.wikipedia.org/wiki/Least_frequently_used) cache.
#
# Implement the `LFUCache` class:
#
#   * `LFUCache(int capacity)` Initializes the object with the `capacity` of the data structure.
#   * `int get(int key)` Gets the value of the `key` if the `key` exists in the cache. Otherwise, returns `-1`.
#   * `void put(int key, int value)` Sets or inserts the value if the `key` is not already present. When the cache reaches its `capacity`, it should invalidate the least frequently used item before inserting a new item. For this problem, when there is a tie (i.e., two or more keys with the same frequency), **the least recently** used `key` would be evicted.
#
# **Notice  that** the number of times an item is used is the number of calls to
# the `get` and `put` functions for that item since it was inserted. This number
# is set to zero when the item is removed.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["LFUCache", "put", "put", "get", "put", "get", "get", "put", "get", "get", "get"]
#     [[2], [1, 1], [2, 2], [1], [3, 3], [2], [3], [4, 4], [1], [3], [4]]
#     **Output**
#     [null, null, null, 1, null, -1, 3, null, -1, 3, 4]
#
#     **Explanation**
#     LFUCache lfu = new LFUCache(2);
#     lfu.put(1, 1);
#     lfu.put(2, 2);
#     lfu.get(1);      // return 1
#     lfu.put(3, 3);   // evicts key 2
#     lfu.get(2);      // return -1 (not found)
#     lfu.get(3);      // return 3
#     lfu.put(4, 4);   // evicts key 1.
#     lfu.get(1);      // return -1 (not found)
#     lfu.get(3);      // return 3
#     lfu.get(4);      // return 4
#
#
#
#
# **Constraints:**
#
#   * `0 <= capacity, key, value <= 104`
#   * At most `105` calls will be made to `get` and `put`.
#
#
#
# **Follow up:** Could you do both operations in `O(1)` time complexity?
#
#
## @lc code=start
using LeetCode

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


text:
"""
The entrypoint to run a set of experiments.

    function run_experiments!(
        results_table::AbstractResultTable,
        instances::Dict{String},
        approaches::Dict{String},
        result_type::Type{<:AbstractResult};
        config::Dict{String,Any}=Dict();
        parallelize::Bool = true,
    )

This will consider every pair of `AbstractInstance` from `keys(instances)`
with each `AbstractApproach` from `keys(approaches)`, build the corresponding
 `AbstractModel`, solve it, and record the resulting result of type
`result_type` in `record_table`. The keys of the dicts `instances` and
`approaches` are String names used to identify each in `results_table`. If
`parallize=true`, the experiments will be parallelized across the available
threads, as determined by the `threads` flag used to launch Julia.
"""
function run_experiments!(
    results_table::AbstractResultTable,
    instances::Dict{String},
    approaches::Dict{String},
    result_type::Type{<:AbstractResult};
    config::Dict{String} = Dict{String,Any}(),
    parallelize::Bool = true,
)
    approach_keys = collect(keys(approaches))
    for (instance_name, instance) in instances
        @info "Prepping instance $instance_name"
        prep_instance!(instance)
        if parallelize
            Threads.@threads for approach_name in approach_keys
                @info "Current experiment: $instance_name with $approach_name"
                approach = approaches[approach_name]
                model = build_model(approach, instance, config)
                optimize!(model)
                result = tear_down(model, instance, result_type)
                record_result!(results_table, result, instance_name, approach_name)
            end
        else
            for (approach_name, approach) in approaches
                @info "Current experiment: $instance_name with $approach_name"
                model = build_model(approach, instance, config)
                optimize!(model)
                result = tear_down(model, instance, result_type)
                record_result!(results_table, result, instance_name, approach_name)
            end
        end
    end
    return nothing
end


text:
const table="accession2taxid"

function create!(db::SQLite.DB, source::String; overwrite::Bool=false, header::Bool=true, delim::Union{Nothing, Char, String}=nothing , accession_col::Int=2, taxid_col::Int=3)
    overwrite ? DBInterface.execute(db,"DROP TABLE IF EXISTS $table") : nothing
    col_width = max(accession_col, taxid_col)
    temp_colname= [ "Column_$i" for i in 1:col_width ]
    temp_colname[accession_col] = "accession"
    temp_colname[taxid_col] = "taxid"
    CSV.File(source; header=temp_colname, select=[accession_col,taxid_col], delim=delim) |> SQLite.load!(db, table)
    SQLite.createindex!(db, table, "accessionindex", "accession")
    return db
end

function Base.get(db::SQLite.DB, accession::String, default::Any)
    result = DBInterface.execute(db, "SELECT * FROM $table WHERE accession = \"$accession\"")
    id = [row.taxid for row in result]
    @assert length(id) < 2
    length(id) == 0 && return default
    return id[1]
end


text:
@timedtestset "gaussian_process_probabilistic_programme" begin

    @timedtestset "split" begin
        x = BlockData(randn(5), randn(4))
        @testset "Vector" begin
            x1, x2 = split(x, randn(9))
            @test length(x1) == 5
            @test length(x2) == 4
        end
        @testset "Matrix" begin
            x1, x2 = split(x, randn(9, 3))
            @test size(x1) == (5, 3)
            @test size(x2) == (4, 3)
        end
    end

    # Build a toy collection of processes.
    gpc = GPC()
    f1 = wrap(GP(sin, SEKernel()), gpc)
    f2 = wrap(GP(cos, Matern52Kernel()), gpc)
    f3 = f1 + 3 * f2

    # Use them to build a programme.
    f = Stheno.GPPP((f1 = f1, f2 = f2, f3 = f3), gpc)

    # The same answers should be obtained manually or via the GPPP.
    @timedtestset "External Consistency" begin

        x0 = GPPPInput(:f1, randn(4))
        x1 = GPPPInput(:f3, randn(3))

        @test mean(f1, x0.x) == mean(f, x0)
        @test mean(f3, x1.x) == mean(f, x1)

        @test cov(f1, x0.x) == cov(f, x0)
        @test cov(f3, x1.x) == cov(f, x1)

        @test cov(f1, f3, x0.x, x1.x) == cov(f, x0, x1)
        @test var(f3, f1, x1.x, x0.x) == var(f, x1, x0)

        y = rand(f(x1))
        @test cov(posterior(f3(x1.x), y)(x1.x)) == cov(posterior(f(x1), y)(x1))
    end

    # The GPPP must be self-consistent like any other AbstractGP.
    # This should hold for all of the various permutations of applicable input types.
    @testset "Internal Conistency ($(typeof(x0)), $(typeof(x1))" for (x0, x1) in [
        (
            GPPPInput(:f1, randn(4)),
            GPPPInput(:f3, randn(3)),
        ),
        (
            GPPPInput(:f1, randn(4)),
            BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))]),
        ),
        (
            BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))]),
            GPPPInput(:f1, randn(4)),
        ),
        (
            BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))]),
            BlockData([GPPPInput(:f1, randn(6))]),
        ),
        (
            collect(GPPPInput(:f1, randn(4))),
            collect(GPPPInput(:f3, randn(3))),
        ),
        (
            GPPPInput(:f1, randn(4)),
            collect(GPPPInput(:f3, randn(3))),
        ),
        (
            collect(BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))])),
            collect(GPPPInput(:f1, randn(4))),
        ),
        (
            collect(BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))])),
            GPPPInput(:f1, randn(4)),
        ),
        (
            BlockData([collect(GPPPInput(:f2, randn(3))), GPPPInput(:f3, randn(2))]),
            GPPPInput(:f1, randn(4)),
        ),
    ]
        test_internal_abstractgps_interface(MersenneTwister(123456), f, x0, x1)
    end

    @timedtestset "gppp macro" begin

        # Declare a GPPP using the helper functionality.
        f = @gppp let
            f1 = GP(SEKernel())
            f2 = GP(Matern52Kernel())
            f3 = f1 + f2
        end
    end

    # No custom rules to worry about, just need to make sure that nothing errors.
    @timedtestset "Zygote" begin
        x = GPPPInput(:f3, randn(5))
        s = 0.1
        y = rand(f(x, s))
        Zygote.gradient((x, y, f, s) -> logpdf(f(x, s), y), x, y, f, s)
    end

    # Check that we can use one GPPP inside another.
    @timedtestset "nested gppp" begin

        gpc_outer = GPC()
        f1_outer = Stheno.wrap(f, gpc_outer)
        f2_outer = 5 * f1_outer
        f_outer = Stheno.GPPP((f1=f1_outer, f2=f2_outer), gpc_outer)

        x0 = GPPPInput(:f1, randn(5))
        x1 = GPPPInput(:f2, randn(4))
        x0_outer = GPPPInput(:f1, x0)
        x1_outer = GPPPInput(:f2, x1)
        rng = MersenneTwister(123456)
        test_internal_abstractgps_interface(rng, f_outer, x0_outer, x1_outer)
    end
end


text:
@assert begin
    intp=InterpolationFunction(
        :(a0+a1*x+a2*y+a11*x*y),
        [
            :x, :y
        ],
        [
            :a0, :a1, :a2, :a11
        ],
        [
            [0, 0]
        ]
    )

    mat=get_interp_linmat(
        intp,
        [
            [0.0, 0.0],
            [1.0, 0.0],
            [1.0, 1.0],
            [0.0, 1.0]
        ]
    )

    isapprox(
        [1.0 0.0 0.0 0.0; 1.0 1.0 0.0 0.0; 1.0 1.0 1.0 1.0; 1.0 0.0 1.0 0.0],
        mat
    )
end


text:
using Dates

abstract type Object <: OrgComponent end

include("../data/entities.jl")
struct Entity <: Object
    name::AbstractString
end

struct LaTeXFragment <: Object
    contents::AbstractString
    delimiters::Union{Tuple{AbstractString, AbstractString}, Nothing}
end

struct ExportSnippet <: Object
    backend::AbstractString
    snippet::AbstractString
end

mutable struct FootnoteReference{L <: Union{<:AbstractString, Nothing},
                                 D <: Union{Vector{Object}, Nothing}} <: Object
    label::L
    definition::D
end

mutable struct CitationReference <: Object
    prefix::Vector{Object}
    key::AbstractString
    suffix::Vector{Object}
end

mutable struct Citation <: Object
    style::Tuple{Union{AbstractString, Nothing},
                 Union{AbstractString, Nothing}}
    globalprefix::Vector{Object}
    citerefs::Vector{CitationReference}
    globalsuffix::Vector{Object}
end

struct InlineBabelCall <: Object
    name::AbstractString
    header::Union{AbstractString, Nothing}
    arguments::Union{AbstractString, Nothing}
end

struct InlineSourceBlock <: Object
    lang::AbstractString
    options::Union{AbstractString, Nothing}
    body::AbstractString
end

struct LineBreak <: Object end

abstract type Link <: Object end

struct LinkPath <: Object
    protocol::Union{Symbol, AbstractString}
    path::AbstractString
end

struct PlainLink <: Link
    path::LinkPath
end

struct AngleLink <: Link
    path::LinkPath
end

struct RegularLink <: Link
    path::LinkPath
    description::Union{Vector{Object}, Nothing}
end

struct Macro <: Object
    name::AbstractString
    arguments::Vector{AbstractString}
end

mutable struct RadioTarget <: Object
    contents::Vector{Object}
end

mutable struct RadioLink <: Link
    radio::RadioTarget
end

mutable struct Target <: Object
    target::AbstractString
end

abstract type StatisticsCookie <: Object end
struct StatisticsCookiePercent <: StatisticsCookie
    percentage::AbstractString
end
struct StatisticsCookieFraction <: StatisticsCookie
    complete::Union{Integer, Nothing}
    total::Union{Integer, Nothing}
end

abstract type Script <: Object end
struct Subscript <: Script
    char::Char
    script::AbstractString
end
struct Superscript <: Script
    char::Char
    script::AbstractString
end

mutable struct TableCell <: Object
    contents::Vector{Object}
end

abstract type Timestamp <: Object end
mutable struct TimestampRepeaterOrDelay
    type::Symbol
    value::Real
    unit::Char
end
mutable struct TimestampDiary <: Timestamp
    sexp::AbstractString
end
abstract type TimestampInstant <: Timestamp end
mutable struct TimestampActive <: TimestampInstant
    date::Date
    time::Union{Time, Nothing}
    repeater::Union{TimestampRepeaterOrDelay, Nothing}
    warning::Union{TimestampRepeaterOrDelay, Nothing}
end
mutable struct TimestampInactive <: TimestampInstant
    date::Date
    time::Union{Time, Nothing}
    repeater::Union{TimestampRepeaterOrDelay, Nothing}
    warning::Union{TimestampRepeaterOrDelay, Nothing}
end
abstract type TimestampRange <: Timestamp end
mutable struct TimestampActiveRange <: TimestampRange
    start::TimestampActive
    stop::TimestampActive
end
mutable struct TimestampInactiveRange <: TimestampRange
    start::TimestampInactive
    stop::TimestampInactive
end

mutable struct TextMarkup{C <: Union{Vector{Object}, <:AbstractString}} <: Object
    formatting::Symbol
    contents::C
end

struct TextPlain{S <: AbstractString} <: Object
    text::S
end


text:
# ---
# title: 832. Flipping an Image
# id: problem832
# author: Qling
# date: 2021-02-24
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/flipping-an-image/description/>
# hidden: true
# ---
#
# Given a binary matrix `A`, we want to flip the image horizontally, then invert
# it, and return the resulting image.
#
# To flip an image horizontally means that each row of the image is reversed.
# For example, flipping `[1, 1, 0]` horizontally results in `[0, 1, 1]`.
#
# To invert an image means that each `0` is replaced by `1`, and each `1` is
# replaced by `0`. For example, inverting `[0, 1, 1]` results in `[1, 0, 0]`.
#
# **Example 1:**
#
#
#
#     Input: [[1,1,0],[1,0,1],[0,0,0]]
#     Output: [[1,0,0],[0,1,0],[1,1,1]]
#     Explanation: First reverse each row: [[0,1,1],[1,0,1],[0,0,0]].
#     Then, invert the image: [[1,0,0],[0,1,0],[1,1,1]]
#
#
# **Example 2:**
#
#
#
#     Input: [[1,1,0,0],[1,0,0,1],[0,1,1,1],[1,0,1,0]]
#     Output: [[1,1,0,0],[0,1,1,0],[0,0,0,1],[1,0,1,0]]
#     Explanation: First reverse each row: [[0,0,1,1],[1,0,0,1],[1,1,1,0],[0,1,0,1]].
#     Then invert the image: [[1,1,0,0],[0,1,1,0],[0,0,0,1],[1,0,1,0]]
#
#
# **Notes:**
#
#   * `1 <= A.length = A[0].length <= 20`
#   * `0 <= A[i][j] <= 1`
#
#
## @lc code=start
using LeetCode

function flip_and_invert_image(A::Vector{Vector{Int}})::Vector{Vector{Int}}
    m = length(A)
    for i = 1:m
        reverse!(A[i]) .⊻= 1
    end

    return A
end

## @lc code=end


text:
using StatsBase
using Statistics

function fuel2pos(positions,finalposition)
    sum(abs.(positions .- finalposition))
end
function answer1()
    input = "input"
    crap_positions0 = parse.(Int,split(readline(input),","))
    test_final_position = median(crap_positions0)
    println("range: ", minimum(crap_positions0)," - ", maximum(crap_positions0))
    min = Inf
    minpos = -1
    for pos in minimum(crap_positions0):1:maximum(crap_positions0)
        distance = fuel2pos(crap_positions0,pos)
        if distance < min
            min = distance
            minpos = pos
        end
    end
    min
    # crap_positions0
end

@show answer1()

triangulardistance(n) = Int(n*(n+1)/2)
function fuel2pos2(positions,finalposition)
    sum(triangulardistance.(abs.(positions .- finalposition)))
end
function answer2()
    input = "input"
    crap_positions0 = parse.(Int,split(readline(input),","))
    println("range: ", minimum(crap_positions0)," - ", maximum(crap_positions0))
    min = Inf
    minpos = -1
    for pos in minimum(crap_positions0):1:maximum(crap_positions0)
        distance = fuel2pos2(crap_positions0,pos)
        if distance < min
            min = distance
            minpos = pos
        end
    end
    min
    # crap_positions0
end

@show answer2()


text:
module metrics

using ..operators, ..operators_sparse

export tracedistance

function tracedistance(rho::DenseOperator, sigma::DenseOperator)
    delta = (rho - sigma).data
    @assert size(delta, 1) == size(delta, 2)
    for i=1:size(delta,1)
        delta[i,i] = real(delta[i,i])
    end
    s = eigvals(Hermitian(delta))
    return 0.5*sum(abs(s))
end

function tracedistance{T<:Operator}(rho::T, sigma::T)
    throw(ArgumentError("tracedistance not implemented for $(T). Use dense operators instead."))
end

end # module


text:
using Receive: Break, receive
using Select

include("select.jl")
include("csp.jl")

cases = [
("Select.jl", receive_select),
("CSP", receive_csp),
("Receive.jl", receive)
]

for (str, receive) in cases

    state = let
        quit = Channel()
        results = Channel()

        t3 = @async begin
            for i in 1:1000000
                put!(results, i)
            end
        end

        t4 = @async begin
            sleep(1)
            put!(quit, nothing)
        end

        state = 0
        curvalue = 0

        # Select.jl has starvation issue if instead
        # `receive_select(results, quit)` is used!
        receive(quit, results) do ch, value

            if ch==quit
                close(quit)
                return Break()
            end

            state += 1
            curvalue = value

        end

        @assert state == curvalue
        state
    end

    println("$str: $state")
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

module Tesseract
using Tesseract_jll
using Leptonica_jll

const TESSERACT = Tesseract_jll.libtesseract
const LEPTONICA = Leptonica_jll.liblept
const TESS_VERSION = "5.1.0"
const LEPT_VERSION = "leptonica-1.82.0"
const TESS_DATA = "tessdata"

include("leptonica/_package.jl")

include("sample.jl")

include("tesseract/_package.jl")

#include("tess_result_renderer.jl")

export tess_version
export TessInst, TessParam, Tsv
export tess_delete!, tess_image, tess_recognize, tess_text, tess_hocr, tess_alto,
       tess_tsv, tess_text_box, tess_word_box, tess_lstm_box, tess_unlv, tess_unlv_latin1,
       tess_confidences, tess_parsed_tsv, tess_resolution
export tess_init, tess_initialized_languages, tess_loaded_languages,
       tess_available_languages, tess_params, tess_read_config,
       tess_read_debug_config, tess_params_parsed,
       tess_set_input_name, tess_get_input_name
export tess_get_param, tess_set_param, tess_set_debug_param
export sample_tiff, sample_pix

export Pix, PixBox
export IFF, IFF_UNKNOWN, IFF_BMP, IFF_JFIF_JPEG, IFF_PNG, IFF_TIFF, IFF_TIFF_PACKBITS,
       IFF_TIFF_RLE, IFF_TIFF_G3, IFF_TIFF_G4, IFF_TIFF_LZW, IFF_TIFF_ZIP, IFF_PNM, IFF_PS,
       IFF_GIF, IFF_JP2, IFF_WEBP, IFF_LPDF, IFF_TIFF_JPEG, IFF_DEFAULT, IFF_SPIX
export pix_read_bmp, pix_write_bmp, pix_read_gif, pix_write_gif, pix_read_jp2k,
       pix_write_jp2k, pix_read_jpeg, pix_write_jpeg, pix_write_pdf, pix_write_pam,
       pix_read_png, pix_write_png, pix_read_pnm, pix_write_pnm, pix_write_ps,
       pix_read_spix, pix_write_spix, pix_read_tiff, pix_write_tiff, pix_read_webp,
       pix_write_webp, pix_read, pix_write, pix_write_implied_format
export pix_get_depth, pix_get_dimensions
export pix_delete!, lept_version

export GitHubProject
export update_languages, download_languages, update_pdf_font, download_pdf_font

#include("tess_result_renderer.jl")
export TessPipeline, TessOutput

export tess_run_pipeline, is_available, tess_pipeline_alto, tess_pipeline_hocr,
       tess_pipeline_lstm_box, tess_pipeline_pdf, tess_pipeline_text, tess_pipeline_tsv,
       tess_pipeline_unlv, tess_pipeline_unlv_latin1, tess_pipeline_word_box
end


text:
struct FlattenedArray{V,L,M,N,P<:NestedArray{V,M,N},InAx<:Anys{M}} <: AbstractArray{V,L}
    parent::P
    inner_axes::InAx
    @inline function FlattenedArray{V,L,M,N,P,InAx}(parent, inner_axes) where {V,L,M,N,P<:NestedArray{V,M,N},InAx<:Anys{M}}
        # TODO check type parameters
        new{V,L,M,N,P,InAx}(parent, inner_axes)
    end
end

@inline function FlattenedArray(parent::NestedArray{V,M,N}, inner_axes::Anys{M}) where {V,M,N}
    FlattenedArray{V,M+N,M,N,typeof(parent),typeof(inner_axes)}(parent, inner_axes)
end

FlattenedArray(parent::NestedArray) = FlattenedArray(parent, inner_axes(parent))

####
#### Core Array Interface
####

@inline Base.axes(F::FlattenedArray) = (F.inner_axes..., axes(F.parent)...)

@inline Base.size(F::FlattenedArray) = map(Base.unsafe_length, axes(F))

@propagate_inbounds function Base.getindex(F::FlattenedArray{<:Any,L,M}, I::Vararg{Int,L}) where {L,M}
    Iin, Iout = tuple_split(I, Val(M))
    return F.parent[Iout...][Iin...]
end

@propagate_inbounds function Base.setindex!(
    F::FlattenedArray{<:Any,L,M},
    v,
    I::Vararg{Int,L},
) where {L,M}
    Iin, Iout = tuple_split(I, Val(M))
    F.parent[Iout...][Iin...] = v
    return F
end

####
#### Misc
####

Base.copy(F::FlattenedArray) = FlattenedArray(deepcopy(F.parent), F.inner_axes)

Base.parent(F::FlattenedArray) = F.parent

function Base.showarg(io::IO, A::FlattenedArray, toplevel)
    print(io, "flatten(")
    Base.showarg(io, parent(A), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(A))
    return nothing
end

# TODO verify that dataids should be the union of the child arrays
# Because F::FlattenedArray is a flattened view of an AbstractArray{<:AbstractArray},
# Base.dataids should return the equivalent of union(Base.dataids(x), Base.dataids.(x)...).
# This works, but the default recursive implementation is slow when length(parent(F)) is large.
function Base.mightalias(A::AbstractArray, F::FlattenedArray)
    !isbits(A) && !isbits(F) && !isdisjoint(Base.dataids(A), _dataids(F))
end
function Base.mightalias(F::FlattenedArray, A::AbstractArray)
    !isbits(F) && !isbits(A) && !isdisjoint(_dataids(F), Base.dataids(A))
end
function Base.mightalias(F1::FlattenedArray, F2::FlattenedArray)
    !isbits(F1) && !isbits(F2) && !isdisjoint(_dataids(F1), _dataids(F2))
end

Base.dataids(F::FlattenedArray) = Tuple(_dataids(F))

function _dataids(F::FlattenedArray)
    ids = collect(Base.dataids(F.parent))
    for A in F.parent
        append!(ids, Base.dataids(A))
    end
    return ids
end

####
#### Extra
####

"""
    flatten(A::AbstractArray{<:AbstractArray{V,M},N})

Return a `M+N`-dimensional flattened view of `A`. If `A` is not a subtype of
`AbstractArray{<:AbstractArray{V,M},N}`, the return value is `A` itself.
Throws an error if the elements of `A` do not have equal size.
Equivalent to `align(A, 1:M...)`.

See also: [`align`](@ref).

# Examples

```jldoctest
julia> A = [[1 2; 3 4], [5 6; 7 8]]
2-element Array{Array{Int64,2},1}:
 [1 2; 3 4]
 [5 6; 7 8]

julia> flatten(A)
2×2×2 flatten(::Array{Array{Int64,2},1}) with eltype Int64:
[:, :, 1] =
 1  2
 3  4

[:, :, 2] =
 5  6
 7  8

julia> flatten(A) == reshape(reduce(hcat, A), 2, 2, 2)
true
```
"""
flatten(A::NestedArray) = FlattenedArray(A)
flatten(A::AbstractArray) = A

"""
    $(SIGNATURES)

A recursive version of [`flatten`](@ref).

# Examples

```jldoctest
julia> x = [ [[1,2], [3,4]], [[5,6], [7,8]]]
2-element Array{Array{Array{Int64,1},1},1}:
 [[1, 2], [3, 4]]
 [[5, 6], [7, 8]]

julia> deep_flatten(x)
2×2×2 flatten(flatten(::Array{Array{Array{Int64,1},1},1})) with eltype Int64:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> deep_flatten(x) == flatten(flatten(x))
true
```
"""
deep_flatten(A::NestedArray) = deep_flatten(flatten(A))
deep_flatten(A::AbstractArray) = A


text:
@info "Defining spyCliqMat(..) for visualizing association matrix of a clique in the Bayes (Junction) tree"

exportimg(pl) = Gadfly.PNG(pl)

export spyCliqMat

"""
    $SIGNATURES

Draw the clique association matrix, with keyword arguments for more or less console print outs.

Notes
* Columns are variables, rows are factors.
* Drawn from up message passing perspective.
* Blue color implies no factor association.
* Frontal, separator, and upmessages are all drawn at different intensity of red.
* Downward messages not shown, as they would just be singletons of the full separator set.
"""
function spyCliqMat(cliq::TreeClique; showmsg=true, suppressprint::Bool=false)
  mat = deepcopy(getCliqMat(cliq, showmsg=showmsg))
  # TODO -- add improved visualization here, iter vs skip
  mat = map(Float64, mat)*2.0.-1.0
  numlcl = size(getCliqAssocMat(cliq),1)
  mat[(numlcl+1):end,:] *= 0.9
  mat[(numlcl+1):end,:] .-= 0.1
  numfrtl1 = floor(Int,length(getCliqueData(cliq).frontalIDs) + 1)
  mat[:,numfrtl1:end] *= 0.9
  mat[:,numfrtl1:end] .-= 0.1
  if !suppressprint
    @show getCliqueData(cliq).itervarIDs
    @show getCliqueData(cliq).directvarIDs
    @show getCliqueData(cliq).msgskipIDs
    @show getCliqueData(cliq).directFrtlMsgIDs
    @show getCliqueData(cliq).directPriorMsgIDs
  end
  if size(mat,1) == 1
    mat = [mat; -ones(size(mat,2))']
  end
  sp = Gadfly.spy(mat)
  push!(sp.guides, Gadfly.Guide.title("$(getLabel(cliq)) || $(getCliqueData(cliq).frontalIDs) :$(getCliqueData(cliq).separatorIDs)"))
  push!(sp.guides, Gadfly.Guide.xlabel("fmcmcs $(getCliqueData(cliq).itervarIDs)"))
  push!(sp.guides, Gadfly.Guide.ylabel("lcl=$(numlcl) || msg=$(size(getCliqMsgMat(cliq),1))" ))
  return sp
end
function spyCliqMat(bt::AbstractBayesTree, lbl::Symbol; showmsg=true, suppressprint::Bool=false)
  spyCliqMat(getClique(bt,lbl), showmsg=showmsg, suppressprint=suppressprint)
end


text:
struct ParticleList{T, ST <: AbstractMatrix{T}, VT <: NamedTuple, PT <: NamedTuple, PART <: AbstractVector, VART <: NamedTuple}
    list::ST
    views::VT
    params::PT
    particles::PART
    variables::VART

    function ParticleList(list, views, params, particles, variables)
        @assert isdisjoint(fieldnames(ParticleList), keys(views))
        @assert isdisjoint(fieldnames(ParticleList), keys(params))
        @assert isdisjoint(keys(views), keys(params))
        new{eltype(list), typeof(list), typeof(views), typeof(params), typeof(particles), typeof(variables)}(list, views, params, particles, variables)
    end
end

function ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple())
    views = map(idx_range -> view(list, idx_range, :), variables)
    vars = map(idx_range -> [view(list, idx_range, i) for i in axes(list, 2)], variables)
    particles = [Particle(p; variables = variables, parameters = parameters) for p in eachcol(list)]
    ParticleList(list, views, parameters, particles, vars)
end

function ParticleList(DT, np, nd; kwargs...)
    ParticleList(zeros(DT, nd, np); kwargs...)
end

@inline function Base.hasproperty(::ParticleList{T,ST,VT,PT}, s::Symbol) where {T,ST,VT,PT}
    hasfield(VT, s) || hasfield(PT, s) || hasfield(Particle, s)
end

@inline function Base.getproperty(p::ParticleList{T,ST,VT,PT}, s::Symbol) where {T, ST, VT, PT}
    if hasfield(VT, s)
        return getfield(p, :views)[s]
    elseif hasfield(PT, s)
        return getfield(p, :params)[s]
    else
        return getfield(p, s)
    end
end

Base.eltype(::ParticleList{T}) where {T} = T

Base.length(pl::ParticleList) = length(pl.particles)

Base.size(pl::ParticleList) = size(pl.list)

Base.getindex(pl::ParticleList, i::Int) = getindex(pl.particles, i)

Base.getindex(pl::ParticleList, I::Vararg{Int}) = getindex(pl.list, I...)

Base.getindex(pl::ParticleList, I...) = getindex(pl.list, I...)

Base.setindex!(pl::ParticleList, v, i::Int) = setindex!(pl.particles, v, i)

Base.setindex!(pl::ParticleList, v, I::Vararg{Int}) = setindex!(pl.list, v, I...)

Base.setindex!(pl::ParticleList, X, I...) = setindex!(pl.list, X, I...)

Base.eachindex(pl::ParticleList) = eachindex(pl.particles)

Base.iterate(pl::ParticleList) = (pl[1], 1)

Base.iterate(pl::ParticleList, i::Int) = i < length(pl) ? (pl[i+1], i+1) : nothing

# function ParticleList(h5::H5DataStore, path::AbstractString = "/")
#     group = h5[path]

#     x = read(group["x"])
#     v = read(group["v"])
#     w = read(group["w"])

#     ParticleList(x, v, w)
# end

# function ParticleList(fpath::AbstractString, path::AbstractString = "/")
#     h5open(fpath, "r") do file
#         ParticleList(file, path)
#     end
# end

# function h5save(h5::H5DataStore, p::ParticleList; path::AbstractString = "/")
#     group = _create_group(h5, path)

#     group["x"] = p.x
#     group["v"] = p.v
#     group["w"] = p.w
# end

# function h5load(::Type{ParticleList}, h5::H5DataStore; path::AbstractString = "/")
#     ParticleList(h5, path)
# end

# Base.length(pl::ParticleList) = length(pl.w)