text:
#!/usr/bin/env bash
    #=
    exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

ENV["GKS_ENCODING"] = "utf-8"

include("random_numbers_core.jl")

function add_plot_info(plt::AbstractPlot; extra::AbstractString=string(), counter::Integer=0)
    x = xlims(Plots.current())[1]
    y = ylims(Plots.current())[2]
    factor1 = y > 1 ? 2 : 1.6
    factor2 = y > 1 ? 4 : 2.5
    # plot!(annotations = (max(x * 2, x+100), (y - 10 * counter) / factor1,
    #        text(extra, :left, 13)))
	annotate!(plt, max(x * 2, x + 100), ((y - 10 * counter) / factor1), extra)
end

function graphing(q::Integer, n::Integer, d::Integer, data::Union{AbstractArray, DataFrame}, stop_at::Integer; m::Integer=10_000)
	save_path = joinpath(dirname(@__DIR__), "other", "random_number_analysis", "random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).pdf")

	if data isa DataFrame
		data = Matrix(data)
	end

	🐖_obtained = ifelse(🐖 ∈ data, true, false)

	## Plot points
	theme(:solarized)

	x_max = 0
	if 🍖 - maximum(data) ∈ [0:upper_bound_adjustment...]
		x_max = 🍖
	end
	if 🐺 - maximum(data) ∈ [0:upper_bound_adjustment...]
		x_max = 🐺
	end

	bin_adjustment = ifelse(🐖_obtained, 0.5, minimum(abs.(extrema(vcat(data, x_max)) .- 🐖)) + 0.5) # ensure the bins encompass the greedy
	bins = minimum(data) - bin_adjustment : max(maximum(data), x_max) + bin_adjustment

    plt = histogram(
					data,
					fontfamily = font("Times"),
					xlabel = "|C|",
					ylabel = "Frequency of Times Found",
					label = "",
					title = "(q, n, d) = ($q, $n, $d): $(format(stop_at, commas = true)) Random Samples, with Selection Size $(format(m, commas = true))",
					# xticks = minimum(data) : max(maximum(data), x_max),
					bins = bins
					)

	## Annotate plot
	annotation_counter = 0
	if isempty(data)
		add_plot_info(plt, extra="🐖 = $(🐖)", counter=annotation_counter); annotation_counter += 1
		add_plot_info(plt, extra="🍖 = $(🍖)", counter=annotation_counter); annotation_counter += 1
		add_plot_info(plt, extra="🐺 = $(🐺)", counter=annotation_counter); annotation_counter += 1
		savefig(plt, save_path)
	    println("Plot saved at $save_path.")
		return nothing
	end

	🐖_idx = searchsortedlast(bins, 🐖) # <- the index of that value
	🐖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐖_idx] # height of bar of choice
	annotate!(plt, 🐖, 🐖_y_top, "🐖")

	annotation_factor = 10
	x_annotation_factor = abs(bins[end] - bins[1]) * 0.05
	y_annotation_factor = abs(bins[end] - bins[1]) * 1.2
	🍖_y, 🐺_y = 0, 0
	🍖_annotated, 🐺_annotated = false, false

	if abs(x_max - 🍖) ≤ upper_bound_adjustment && abs(x_max - 🐺) ≤ upper_bound_adjustment # the singleton AND hamming bounds can fit on the plot
		🍖_idx = searchsortedlast(bins, 🍖) # <- the index of that value
		🍖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🍖_idx] # height of bar of choice
		🐺_idx = searchsortedlast(bins, 🐺) # <- the index of that value
		🐺_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐺_idx] # height of bar of choice

		# adjust emoticon positioning if overlapping values.
		if 🐖 == 🍖
			🍖_y = 🍖_y_top + y_annotation_factor
		end

		if 🐖 == 🐺
			if 🐖 == 🍖 # then all of 🍖, 🐖, and 🐺 are the same
				🍖_y = 🍖_y_top + y_annotation_factor
				🐺_y = 🐺_y_top + 2*y_annotation_factor
			else
				🐺_y = 🐺_y_top + y_annotation_factor
			end
		end

		annotate!(plt, 🍖, 🍖_y, "🍖")
		annotate!(plt, 🐺, 🐺_y, "🐺")

		🍖_annotated = true
		🐺_annotated = true
	end

	if abs(x_max - 🍖) ≤ upper_bound_adjustment && ! 🍖_annotated # only the hamming bound fits on the plot
		🍖_idx = searchsortedlast(bins, 🍖) # <- the index of that value
		🍖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🍖_idx] # height of bar of choice

		# adjust emoticon positioning if overlapping values.
		if 🐖 == 🍖
			🍖_y = 🍖_y_top + y_annotation_factor
		end

		annotate!(plt, 🍖, 🍖_y, "🍖")

		🍖_annotated = true
	end

	if abs(x_max - 🐺) ≤ upper_bound_adjustment && ! 🐺_annotated # only the singleton bound fits on the plot
		🐺_idx = searchsortedlast(bins, 🐺) # <- the index of that value
		🐺_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐺_idx] # height of bar of choice

		# adjust emoticon positioning if overlapping values.
		if 🐖 == 🐺
			🐺_y = 🐺_y_top + y_annotation_factor
		end

		annotate!(plt, 🐺, 🐺_y, "🐺")

		🐺_annotated = true
	end

	if ! 🍖_annotated
		annotate!(plt, maximum(filter(!isnan, plt.series_list[1].plotattributes[:x])) - x_annotation_factor, maximum(filter(!isnan, plt.series_list[1].plotattributes[:y])) + 0, "🍖 = $(🍖)")
		# add_plot_info(plt, extra="🍖 = $(🍖)", counter=annotation_counter); annotation_counter += 1
	end

	if ! 🐺_annotated
		annotate!(plt, maximum(filter(!isnan, plt.series_list[1].plotattributes[:x])) - x_annotation_factor, maximum(filter(!isnan, plt.series_list[1].plotattributes[:y])) - y_annotation_factor, "🐺 = $(🐺)")
		# add_plot_info(plt, extra="🐺 = $(🐺)", counter=annotation_counter); annotation_counter += 1
	end

    savefig(plt, save_path)
    println("Plot saved at $save_path.")
	return nothing
end

graphing(q::Integer, n::Integer, d::Integer, stop_at::Integer; m::Integer=10_000) =
	graphing(q, n, d, obtain_data(q, n, d, stop_at, m = m), stop_at, m = m)


text:
"""
    convert2features(array)

Convert the given MNIST tensor to a feature matrix (or feature
vector in the case of a single image). The purpose of this
function is to drop the spatial dimensions such that traditional
ML algorithms can process the dataset.

```julia
julia> MNIST.convert2features(MNIST.traintensor()) # full training data
784×60000 Array{N0f8,2}:
[...]

julia> MNIST.convert2features(MNIST.traintensor(1)) # first observation
784-element Array{N0f8,1}:
[...]
```
"""
convert2features(array::AbstractMatrix{<:Number}) = vec(array)

function convert2features(array::AbstractArray{<:Number,3})
    nrows, ncols, nimages = size(array)
    reshape(array, (nrows * ncols, nimages))
end

"""
    convert2image(array) -> Array{Gray}

Convert the given MNIST horizontal-major tensor (or feature matrix)
to a vertical-major `Colorant` array. The values are also color
corrected according to the website's description, which means that
the digits are black on a white background.

```julia
julia> MNIST.convert2image(MNIST.traintensor()) # full training dataset
28×28×60000 Array{Gray{N0f8},3}:
[...]

julia> MNIST.convert2image(MNIST.traintensor(1)) # first training image
28×28 Array{Gray{N0f8},2}:
[...]
```
"""
function convert2image(array::AbstractVector{<:Number})
    @assert length(array) % 784 == 0
    if length(array) == 784
        convert2image(reshape(array, 28, 28))
    else
        n = Int(length(array) / 784)
        convert2image(reshape(array, 28, 28, n))
    end
end

function convert2image(array::AbstractMatrix{T}) where {T<:Number}
    if size(array) == (28, 28)
        # simple check to see if values are normalized to [0,1]
        if any(x->x > 50, array)
            _colorview(Gray, T(1) .- transpose(array) ./ T(255))
        else
            _colorview(Gray, T(1) .- transpose(array))
        end
    else # feature matrix
        @assert size(array, 1) == 784
        n = size(array, 2)
        convert2image(reshape(array, 28, 28, n))
    end
end

function convert2image(array::AbstractArray{T,3}) where {T<:Number}
    h, w, n = size(array)
    @assert h == 28 && w == 28
    # simple check to see if values are normalized to [0,1]
    if any(x->x > 50, array)
        _colorview(Gray, permutedims(T(1) .- array ./ T(255), [2,1,3]))
    else
        _colorview(Gray, permutedims(T(1) .- array, [2,1,3]))
    end
end

ext: jl


text:
# This is an extended version of approx_space_time_inference.jl. It combines it with
# Optim + ParameterHandling + Zygote to learn the kernel parameters.
# If you understand how to use Optim + ParameterHandling + Zygote for an AbstractGP,
# e.g. that shown on the README for this package, and how approx_space_time_inference.jl
# works, then you should understand this file.

using AbstractGPs
using KernelFunctions
using TemporalGPs

using TemporalGPs: Separable, approx_posterior_marginals, RegularInTime

# Load standard packages from the Julia ecosystem
using Optim # Standard optimisation algorithms.
using ParameterHandling # Helper functionality for dealing with model parameters.
using Zygote # Algorithmic Differentiation

using ParameterHandling: flatten

# Declare model parameters using `ParameterHandling.jl` types.
flat_initial_params, unflatten = flatten((
    var_kernel = positive(0.6),
    λ_space = positive(0.5),
    λ_time = positive(0.1),
    var_noise = positive(0.1),
));

# Construct a function to unpack flattened parameters and pull out the raw values.
unpack = ParameterHandling.value ∘ unflatten;
params = unpack(flat_initial_params);

function build_gp(params)
    k_space = SEKernel() ∘ ScaleTransform(params.λ_space)
    k_time = Matern52Kernel() ∘ ScaleTransform(params.λ_time)
    k = params.var_kernel * Separable(k_space, k_time)
    return to_sde(GP(k), ArrayStorage(Float64))
end

# Construct inputs. Spatial locations change at each point in time.
# Also works with RectilinearGrids of inputs.
# Times must be increasing, points in space can be anywhere.
N = 50;
T = 1000;
points_in_space = [randn(N) for _ in 1:T];
points_in_time = RegularSpacing(0.0, 0.1, T);
x = RegularInTime(points_in_time, points_in_space);

# Since it's not straightforward to generate samples from this GP at `x`, use a known
# function, under a bit of iid noise.
xs = collect(x);
y = sin.(first.(xs)) .+ cos.(last.(xs)) + sqrt.(params.var_noise) .* randn(length(xs));

# Spatial pseudo-point inputs.
z_r = collect(range(-3.0, 3.0; length=5));

# Specify an objective function for Optim to minimise in terms of x and y.
# We choose the usual negative log marginal likelihood (NLML).
function objective(params)
    f = build_gp(params)
    return -elbo(f(x, params.var_noise), y, z_r)
end

# Optimise using Optim. Takes a little while to compile because Zygote.
training_results = Optim.optimize(
    objective ∘ unpack,
    θ -> only(Zygote.gradient(objective ∘ unpack, θ)),
    flat_initial_params,
    BFGS(
        alphaguess = Optim.LineSearches.InitialStatic(scaled=true),
        linesearch = Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(show_trace = true);
    inplace=false,
);

# Extracting the final values of the parameters.
# Should be close to truth.
final_params = unpack(training_results.minimizer);
# final_params = unpack(flat_initial_params);

# Locations in space at which to make predictions. Assumed to be the same at each point in
# time, but this assumption could easily be relaxed.
N_pr = 150;
x_r_pr = range(-5.0, 5.0; length=N_pr);

# Compute the approximate posterior marginals.
fx_final = build_gp(final_params)(x, final_params.var_noise)
f_post_marginals = approx_posterior_marginals(dtc, fx_final, y, z_r, x_r_pr);
m_post_marginals = mean.(f_post_marginals);
σ_post_marginals = std.(f_post_marginals);

# Visualise the posterior marginals. We don't do this during in CI because it causes
# problems.
if get(ENV, "TESTING", "FALSE") == "FALSE"
    using Plots
    savefig(
        plot(
            heatmap(reshape(m_post_marginals, N_pr, T)),
            heatmap(reshape(σ_post_marginals, N_pr, T));
            layout=(1, 2),
        ),
        "posterior.png",
    );
end


text:
export submodules, minimal_submodules, maximal_submodules, composition_series,
       composition_factors, meataxe

################################################################################
#
#  Tools for MeatAxe
#
################################################################################

#
# Given a matrix $M$ in echelon form and a vector, it returns
# the vector reduced with respect to $M$
#
function cleanvect(M::T, v::T) where {T}
  @assert nrows(v)==1
  w=deepcopy(v)
  if iszero(v)
    return w
  end
  for i=1:nrows(M)
    if iszero_row(M,i)
      continue
    end
    ind=1
    while iszero(M[i,ind])
      ind+=1
    end
    if iszero(w[1,ind])
      continue
    end
    mult=divexact(w[1,ind], M[i,ind])
    w[1,ind] = parent(M[1,1])(0)
    for k=ind+1:ncols(M)
      w[1,k]-= mult*M[i,k]
    end
  end
  return w

end

@doc Markdown.doc"""
    closure(C::T, G::Array{T,1}) where T <: MatElem

Given a matrix C representing a subspace of K^n and a list of matrices G representing endomorphisms of K^n,
the function returns a matrix representing the closure of the subspace under the action, i.e. the smallest
subspace of K^n invariant under the endomorphisms.
"""
function closure(C::T, G::Array{T,1}) where {T}
  rref!(C)
  i=1
  while i <= nrows(C)
    w=view(C, i:i, 1:ncols(C))
    for j=1:length(G)
      res = cleanvect(C, w*G[j])
      if !iszero(res)
        C = vcat(C,res)
        if nrows(C)==ncols(C)
          i=ncols(C)+1
          break
        end
      end
    end
    i+=1
  end
  r = rref!(C)
  if r != nrows(C)
    C = sub(C, 1:r, 1:ncols(C))
  end
  return C
end

#
#  Function to obtain the action of G on the quotient and on the submodule
#

function clean_and_quotient(M::T, N::T, pivotindex::Set{Int}) where {T}

  coeff=zero_matrix(parent(M[1,1]),nrows(N),nrows(M))
  for i=1:nrows(N)
    for j=1:nrows(M)
      if iszero_row(M,j)
        continue
      end
      ind=1
      while iszero(M[j,ind])
        ind+=1
      end
      coeff[i,j]=divexact(N[i,ind], M[j,ind])
      for s=1:ncols(N)
        N[i,s]-=coeff[i,j]*M[j,s]
      end
    end
  end
  vec= zero_matrix(parent(M[1,1]),nrows(N),ncols(M)-length(pivotindex))
  for i=1:nrows(N)
    pos=0
    for s=1:ncols(M)
      if !(s in pivotindex)
        pos+=1
        vec[i,pos]=N[i,s]
      end
    end
  end
  return coeff, vec
end

#  Restriction of the action to the submodule generated by C and the quotient

function __split(C::T, G::Vector{T}) where {T <: MatElem}
# I am assuming that C defines an submodule
  equot=Vector{T}(undef, length(G))
  esub=Vector{T}(undef, length(G))
  pivotindex=Set{Int}()
  for i = 1:nrows(C)
    if iszero_row(C,i)
      continue
    end
    ind = 1
    while iszero(C[i, ind])
      ind += 1
    end
    push!(pivotindex, ind)
  end
  for a = 1:length(G)
    subm, vec=clean_and_quotient(C, C*G[a], pivotindex)
    esub[a] = subm
    s = zero_matrix(base_ring(C),ncols(G[1]) - length(pivotindex), ncols(G[1]) - length(pivotindex))
    pos = 0
    for i= 1:nrows(G[1])
      if !(i in pivotindex)
        m, vec= clean_and_quotient(C, sub(G[a], i:i, 1:nrows(G[1])), pivotindex)
        for j=1:ncols(vec)
          s[i - pos,j] = vec[1, j]
        end
      else
        pos += 1
      end
    end
    equot[a] = s
  end
  return ModAlgAss(esub), ModAlgAss(equot), pivotindex
end

#  Restriction of the action to the submodule generated by C
function _actsub(C::T, G::Vector{T}) where {T <: MatElem}
  esub = Vector{T}(undef, length(G))
  pivotindex = Set{Int}()
  for i=1:nrows(C)
    ind = 1
    while iszero(C[i, ind])
      ind += 1
    end
    push!(pivotindex, ind)
  end
  for a=1:length(G)
    subm, vec = clean_and_quotient(C, C*G[a], pivotindex)
    esub[a] = subm
  end
  return ModAlgAss(esub)
  #return ModAlgAss(esub)
end

#  Restriction of the action to the quotient by the submodule generated by C
function _actquo(C::T,G::Vector{T}) where {T <: MatElem}
  equot = Vector{T}(undef, length(G))
  pivotindex = Set{Int}()
  for i=1:nrows(C)
    ind = 1
    while iszero(C[i,ind])
      ind += 1
    end
    push!(pivotindex,ind)
  end
  for a=1:length(G)
    s = zero_matrix(base_ring(C), ncols(G[1]) - length(pivotindex), ncols(G[1]) - length(pivotindex))
    pos = 0
    for i=1:nrows(G[1])
      if !(i in pivotindex)
        m, vec = clean_and_quotient(C, sub(G[a],i:i,1:nrows(G[1])), pivotindex)
        for j=1:ncols(vec)
          s[i - pos, j]=vec[1, j]
        end
      else
        pos += 1
      end
    end
    equot[a] = s
  end
  return ModAlgAss(equot), pivotindex
end

#
#  Function that determine if two modules are isomorphic, provided that the first is irreducible
#
function isisomorphic(M::ModAlgAss{S, T, V}, N::ModAlgAss{S, T, V}) where {S, T, V}
  @assert M.isirreducible == 1
  @assert base_ring(M) == base_ring(N)
  @assert length(M.action) == length(N.action)
  if dimension(M) != dimension(N)
    return false
  end

  if M.dimension==1
    return M.action==N.action
  end

  K = base_ring(M)
  Kx, x = PolynomialRing(K, "x", cached=false)

  if length(M.action) == 1
    f = charpoly(Kx, M.action[1])
    g = charpoly(Kx, N.action[1])
    if f==g
      return true
    else
      return false
    end
  end
  rel = _relations(M,N)
  return iszero(rel[N.dimension, N.dimension])

end

function _enum_el(K,v,dim)
  if dim == 0
    return [v]
  else
    list=[]
    push!(v,K(0))
    for x in K
      v[length(v)]=x
      push!(list,deepcopy(v))
    end
    list1=[]
    for x in list
      append!(list1,_enum_el(K,x, dim-1))
    end
    return list1
  end
end

function dual_space(M::ModAlgAss{S, T, V}) where {S, T, V}
  G = T[transpose(g) for g in M.action]
  return ModAlgAss(G)
end

function _subst(f::Nemo.PolyElem{T}, a::MatElem{T}) where {T <: Nemo.RingElement}
   n = degree(f)
   if n < 0
      return similar(a)
   elseif n == 0
      return coeff(f, 0) * identity_matrix(base_ring(a), nrows(a))
   elseif n == 1
      return coeff(f, 0) * identity_matrix(base_ring(a), nrows(a)) + coeff(f, 1)*a
   end
   d1 = isqrt(n)
   d = div(n, d1)
   A = powers(a, d)
   s = coeff(f, d1*d)*A[1]
   for j = 1:min(n - d1*d, d - 1)
      c = coeff(f, d1*d + j)
      if !iszero(c)
         s += c*A[j + 1]
      end
   end
   for i = 1:d1
      s *= A[d + 1]
      s += coeff(f, (d1 - i)*d)*A[1]
      for j = 1:min(n - (d1 - i)*d, d - 1)
         c = coeff(f, (d1 - i)*d + j)
         if !iszero(c)
            s += c*A[j + 1]
         end
      end
   end
   return s
end

#################################################################
#
#  MeatAxe, Composition Factors and Composition Series
#
#################################################################

@doc Markdown.doc"""
    meataxe(M::ModAlgAss) -> Bool, MatElem

Given module M, returns true if the module is irreducible (and the identity matrix) and false if the space is reducible, togheter with a basis of a submodule

"""
function meataxe(M::ModAlgAss{S, T, V}) where {S, T, V}

  K=base_ring(M)
  Kx,x=PolynomialRing( K,"x", cached=false)
  n=dimension(M)
  H = M.action
  if n == 1
    M.isirreducible = 1
    return true, identity_matrix(K, n)
  end

  G = deepcopy(H)
  filter!(x -> !iszero(x), G)

  if length(G) == 0
    return false, matrix(K, 1, n, [one(K) for i = 1:n])
  end

  if isone(length(G))
    A = G[1]
    poly = minpoly(Kx, A)
    sq = factor_squarefree(poly)
    lf = factor(first(keys(sq.fac)))
    t = first(keys(lf.fac))
    if degree(t)==n
      M.isirreducible = 1
      return true, identity_matrix(K, n)
    else
      N = _subst(t, A)
      null, kern = kernel(N, side = :left)
      B = closure(sub(kern, 1:1, 1:n), G)
      return false, B
    end
  end

  #
  #  Adding generators to obtain randomness
  #

  Gt = T[transpose(x) for x in M.action]
  cnt = 0
  while true
    cnt += 1
    if cnt > 1000
      error("Too many attempts")
    end
    # At every step, we add a generator to the group.
    new_gen = G[rand(1:length(G))]*G[rand(1:length(G))]
    while iszero(new_gen)
      new_gen = G[rand(1:length(G))]*G[rand(1:length(G))]
    end
    push!(G, new_gen)

    #
    # Choose a random combination of the generators of G
    #
    A = zero_matrix(K, n, n)
    for i=1:length(G)
      add!(A, A, rand(K)*G[i])
    end

    #
    # Compute the characteristic polynomial and, for irreducible factor f, try the Norton test
    #
    poly = minpoly(Kx, A)
    sqfpart = keys(factor_squarefree(poly).fac)
    for el in sqfpart
      sq = el
      i=1
      while !isone(sq)
        f = gcd(powmod(x, order(K)^i, sq)-x,sq)
        sq = divexact(sq, f)
        lf=factor(f)
        for t in keys(lf.fac)
          N = _subst(t, A)
          a, kern = kernel(N, side = :left)
          @assert a > 0
          #  Norton test
          B = closure(sub(kern, 1:1, 1:n), M.action)
          if nrows(B) != n
            M.isirreducible= 2
            return false, B
          end
          aa, kernt = kernel(transpose(N), side = :left)
          @assert aa == a
          Bt = closure(sub(kernt, 1:1, 1:n), Gt)
          if nrows(Bt) != n
            aa, Btnu = kernel(Bt)
            subst = transpose(Btnu)
            #@assert nrows(subst)==nrows(closure(subst,G))
            M.isirreducible = 2
            return false, subst
          end
          if degree(t) == a
            # f is a good factor, irreducibility!
            M.isirreducible = 1
            return true, identity_matrix(K, n)
          end
        end
        i+=1
      end
    end
  end
end

@doc Markdown.doc"""
    composition_series(M::ModAlgAss) -> Array{MatElem,1}

Given a Fq[G]-module M, it returns a composition series for M, i.e. a sequence of submodules such that the quotient of two consecutive element is irreducible.

"""
function composition_series(M::ModAlgAss{S, T, V}) where {S, T, V}

  if M.isirreducible == 1 || M.dimension == 1
    return [identity_matrix(base_ring(M.action[1]), dimension(M))]
  end

  bool, C = meataxe(M)
  #
  #  If the module is irreducible, we return a basis of the space
  #
  if bool
    return [identity_matrix(base_ring(M.action[1]), dimension(M))]
  end
  #
  #  The module is reducible, so we call the algorithm on the quotient and on the subgroup
  #
  G = M.action
  K = M.base_ring

  rref!(C)

  esub,equot,pivotindex=__split(C,G)
  sub_list = composition_series(esub)
  quot_list = composition_series(equot)
  #
  #  Now, we have to write the submodules of the quotient and of the submodule in terms of our basis
  #
  list=Vector{T}(undef, length(sub_list)+length(quot_list))
  for i=1:length(sub_list)
    list[i]=sub_list[i]*C
  end
  for z=1:length(quot_list)
    s=zero_matrix(K,nrows(quot_list[z]), ncols(C))
    for i=1:nrows(quot_list[z])
      pos=0
      for j=1:ncols(C)
        if j in pivotindex
          pos+=1
        else
          s[i,j]=quot_list[z][i,j-pos]
        end
      end
    end
    list[length(sub_list)+z]=vcat(C,s)
  end
  return list
end

@doc Markdown.doc"""
    composition_factors(M::ModAlgAss)

Given a Fq[G]-module M, it returns, up to isomorphism, the composition factors of M with their multiplicity,
i.e. the isomorphism classes of modules appearing in a composition series of M

"""
function composition_factors(M::ModAlgAss{S, T, V}; dimension::Int=-1) where {S, T, V}

  if M.isirreducible == 1 || M.dimension == 1
    if dimension != -1
      if M.dimension == dimension
        return Tuple{ModAlgAss{S, T, V}, Int}[(M,1)]
      else
        return Tuple{ModAlgAss{S, T, V}, Int}[]
      end
    else
      return Tuple{ModAlgAss{S, T, V}, Int}[(M,1)]
    end
  end

  K=M.base_ring
  bool, C = meataxe(M)
  #  If the module is irreducible, we just return a basis of the space
  if bool
    if dimension!= -1
      if M.dimension==dimension
        return Tuple{ModAlgAss{S, T}, Int}[(M,1)]
      else
        return Tuple{ModAlgAss{S, T}, Int}[]
      end
    else
      return Tuple{ModAlgAss{S, T}, Int}[(M,1)]
    end
  end
  G=M.action

  #  The module is reducible, so we call the algorithm on the quotient and on the subgroup
  rref!(C)
  sub, quot, pivotindex = __split(C, G)
  sub_list = composition_factors(sub)
  quot_list = composition_factors(quot)
  #  Now, we check if the factors are isomorphic
  for i=1:length(sub_list)
    for j=1:length(quot_list)
      if isisomorphic(sub_list[i][1], quot_list[j][1])
        sub_list[i]=(sub_list[i][1], sub_list[i][2]+quot_list[j][2])
        deleteat!(quot_list,j)
        break
      end
    end
  end
  return append!(sub_list, quot_list)

end

function _relations(M::ModAlgAss{S, T, V}, N::ModAlgAss{S, T, V}) where {S, T, V}
  @assert M.isirreducible == 1
  G =M.action
  H=N.action
  K=base_ring(M)
  n=dimension(M)

  sys=zero_matrix(K,2*dimension(N),dimension(N))
  matrices=T[]
  first=true
  B=zero_matrix(K,1,dimension(M))
  B[1,1]=K(1)
  X=B
  push!(matrices, identity_matrix(base_ring(B), dimension(N)))
  i=1
  while i<=nrows(B)
    w=sub(B, i:i, 1:n)
    for j=1:length(G)
      v=w*G[j]
      res=cleanvect(X,v)
      if !iszero(res)
        X=rref(vcat(X,v))[2]
        B=vcat(B,v)
        push!(matrices, matrices[i]*H[j])
      else
        fl, x = can_solve(B, v, side = :left)
        @assert fl
        A=sum([x[1,q]*matrices[q] for q=1:ncols(x)])
        A=A-(matrices[i]*H[j])
        if first
          for s=1:N.dimension
            for t=1:N.dimension
              sys[s,t]=A[t,s]
            end
          end
          first=false
        else
          for s=1:N.dimension
            for t=1:N.dimension
              sys[N.dimension+s,t]=A[t,s]
            end
          end
        end
        rref!(sys)
      end
    end
    if sys[N.dimension,N.dimension]!=0
      break
    end
    i=i+1
  end
  return view(sys, 1:N.dimension, 1:N.dimension)
end

function _irrsubs(M::ModAlgAss{S, T}, N::ModAlgAss{S, T}) where {S, T}

  @assert M.isirreducible == 1
  K=M.base_ring
  rel=_relations(M,N)
  if rel[N.dimension, N.dimension]!=0
    return T[]
  end
  a,kern = nullspace(rel)
  # TODO: Remove this once fixed.
  if !(kern isa T)
    a, kern = kern, a
  end
  kern = transpose(kern)
  if a == 1
    return T[closure(kern, N.action)]
  end
  vects = T[sub(kern, i:i, 1:N.dimension) for i=1:a]
  #First sieve: The vectors may be dependent when considering the structure of G-module.
  #So I reduce the number of vectors by looking at the subspace they generate.
  #If there is an inclusion, then I can remove them.
  to_reduce = closure(vects[1], N.action)
  final_vect_list = [vects[1]]
  for i = 2:length(vects)
    w = cleanvect(to_reduce, vects[i])
    if !iszero(w)
      push!(final_vect_list, vects[i])
      to_reduce = vcat(vects[i], to_reduce)
      to_reduce = closure(to_reduce, N.action)
    end
  end
  if isone(length(final_vect_list))
    return [closure(final_vect_list[1], N.action)]
  end
  # Now, I have a list of generators as G-modules.
  # Therefore I need to consider combinations as G-modules
  return _submodules_direct_sum(final_vect_list, N)

end

function _all_combinations(M::MatElem{T}) where T
  K = base_ring(M)
  els = collect(x for x in K)
  @assert fits(Int, fmpz(length(els))^nrows(M))
  res = Vector{typeof(M)}(undef, length(els)^nrows(M))
  ind = 1
  m = zero_matrix(K, 1, nrows(M))
  for i in CartesianIndices(Tuple(1:length(els) for i in 1:nrows(M)))
    for j = 1:nrows(M)
      m[1, j] = els[i[j]]
    end
    res[ind] = m * M
    ind += 1
  end
  return res
end

function _submodules_direct_sum(gens::Vector{T}, N::ModAlgAss{S, T, V}) where {S, T, V}
  K = base_ring(N)
  closures = [closure(x, N.action) for x in gens]
  all_combinations = [_all_combinations(x) for x in closures]
  res = Vector{T}()
  for i = 1:length(gens)-1
    #I have to list all the elements that have 1 in the first component and all the possible elements in the other.
    it = CartesianIndices(Tuple(1:length(all_combinations[j]) for j = i+1:length(gens)))
    m = zero_matrix(K, 0, dimension(N))
    non_zero_el = all_combinations[i][1]
    ind = 2
    while iszero(non_zero_el)
      non_zero_el = all_combinations[i][ind]
      ind += 1
    end
    for I in it
      m = non_zero_el
      for s = 1:length(I)
        m += all_combinations[i+s][I[s]]
      end
      push!(res, closure(m, N.action))
    end
  end
  non_zero_el = all_combinations[length(gens)][1]
  ind = 2
  while iszero(non_zero_el)
    non_zero_el = all_combinations[length(gens)][ind]
    ind += 1
  end
  push!(res, closure(non_zero_el, N.action))
  return res
end

@doc Markdown.doc"""
    minimal_submodules(M::ModAlgAss)

Given a Fq[G]-module M, it returns all the minimal submodules of M
"""
function minimal_submodules(M::ModAlgAss{S, T, V}, dim::Int=M.dimension+1, lf = Tuple{ModAlgAss{S, T, V}, Int}[]) where {S, T, V}

  K = M.base_ring
  n = M.dimension
  if isone(M.isirreducible)
    if dim >= n
      return T[identity_matrix(K, n)]
    else
      return T[]
    end
  end

  list = T[]
  if isempty(lf)
    lf = composition_factors(M)
  end
  if isone(length(lf)) && isone(lf[1][2])
    return T[identity_matrix(K, n)]
  end
  if dim!=n+1
    lf=[x for x in lf if x[1].dimension==dim]
  end
  if isempty(lf)
    return list
  end
  for x in lf
    append!(list, Hecke._irrsubs(x[1], M))
  end
  return list
end

@doc Markdown.doc"""
    maximal_submodules(M::ModAlgAss)

Given a $G$-module $M$, it returns all the maximal submodules of M

"""
function maximal_submodules(M::ModAlgAss{S, T, V}, index::Int=M.dimension, lf = Tuple{ModAlgAss{S, T, V}, Int}[]) where {S, T, V}

  M_dual = dual_space(M)
  minlist = minimal_submodules(M_dual, index+1, lf)
  maxlist = Array{T, 1}(undef, length(minlist))
  for j=1:length(minlist)
    maxlist[j]=transpose(nullspace(minlist[j])[2])
  end
  return maxlist

end

@doc Markdown.doc"""
    submodules(M::ModAlgAss)

Given a $G$-module $M$, it returns all the submodules of M

"""
function submodules(M::ModAlgAss{S, T, V}) where {S, T, V}

  K=M.base_ring
  list = T[]
  if M.dimension == 1
    return [zero_matrix(K, 1, 1), identity_matrix(K, 1)]
  end
  lf=composition_factors(M)
  minlist = minimal_submodules(M, M.dimension+1, lf)
  for x in minlist
    rref!(x)
    N, pivotindex =_actquo(x,M.action)
    ls=submodules(N)
    for a in ls
      s=zero_matrix(K,nrows(a), M.dimension)
      for t=1:nrows(a)
        pos=0
        for j=1:M.dimension
          if j in pivotindex
            pos+=1
          else
            s[t,j]=a[t,j-pos]
          end
        end
      end
      push!(list,vcat(x,s))
    end
  end
  for i = 1:length(list)
    x = list[i]
    rk = rref!(x)
    if rk < nrows(x)
      nx = sub(x, 1:rk, 1:ncols(x))
      list[i] = nx
    end
  end
  push!(list, zero_matrix(K, 0, M.dimension))
  push!(list, identity_matrix(K, M.dimension))
  append!(list, minlist)
  res = T[list[1]]
  for i = 2:length(list)
    x = list[i]
    found = false
    for j = 1:length(res)
      if x == res[j]
        found = true
        break
      end
    end
    if !found
      push!(res, x)
    end
  end
  return res

end

@doc Markdown.doc"""
    submodules(M::ModAlgAss, index::Int)

Given a $G$-module $M$, it returns all the submodules of M of index q^index, where q is the order of the field

"""
function submodules(M::ModAlgAss{S, T}, index::Int; comp_factors=Tuple{ModAlgAss{S, T}, Int}[]) where {S, T}

  K=M.base_ring
  if index==M.dimension
    return T[zero_matrix(K,1,M.dimension)]
  end
  list=T[]
  if index>= M.dimension/2
    if index== M.dimension -1
      if isempty(comp_factors)
        lf=composition_factors(M, dimension=1)
      else
        lf=comp_factors
      end
      list=minimal_submodules(M,1,lf)
      return list
    end
    if isempty(comp_factors)
      lf=composition_factors(M)
    else
      lf=comp_factors
    end
    for i=1: M.dimension-index-1
      minlist=minimal_submodules(M,i,lf)
      for x in minlist
        N, pivotindex= _actquo(x, M.action)
        #
        #  Recover the composition factors of the quotient
        #
        Sub=_actsub(x, M.action)
        lf1=[(x[1], x[2]) for x in lf]
        for j=1:length(lf1)
          if isisomorphic(lf1[j][1], Sub)
            if lf1[j][2]==1
              deleteat!(lf1,j)
            else
              lf1[j]=(lf1[j][1], lf1[j][2]-1)
            end
            break
          end
        end
        #
        #  Recursively ask for submodules and write their bases in terms of the given set of generators
        #
        ls=submodules(N,index, comp_factors=lf1)
        for a in ls
          s=zero_matrix(K,nrows(a)+nrows(x), M.dimension)
          for t=1:nrows(a)
            pos=0
            for j=1:M.dimension
              if j in pivotindex
               pos+=1
             else
               s[t,j]=a[t,j-pos]
              end
            end
          end
          for t=nrows(a)+1:nrows(s)
            for j=1:ncols(s)
              s[t,j]=x[t-nrows(a),j]
            end
          end
          push!(list,s)
        end
      end
    end

  #
  #  Eliminating repeatitions
  #

    for x in list
      rref!(x)
    end
    i=1
    while i<=length(list)
      k=i+1
      while k<=length(list)
        if list[i]==list[k]
          deleteat!(list, k)
        else
          k+=1
        end
      end
      i+=1
    end
    append!(list,minimal_submodules(M,M.dimension-index, lf))
  else
  #
  #  Duality
  #
    M_dual=dual_space(M)
    dlist=submodules(M_dual, M.dimension-index)
    list=T[transpose(nullspace(x)[2]) for x in dlist]
  end
  return list

end

## Make Nmod iteratible

Base.iterate(R::NmodRing) = (zero(R), zero(UInt))

function Base.iterate(R::NmodRing, st::UInt)
  if st == R.n - 1
    return nothing
  end

  return R(st + 1), st + 1
end

Base.eltype(::Type{NmodRing}) = nmod

Base.IteratorSize(::Type{NmodRing}) = Base.HasLength()

Base.length(R::NmodRing) = R.n

Base.iterate(R::GaloisField) = (zero(R), zero(UInt))

function Base.iterate(R::GaloisField, st::UInt)
  if st == R.n - 1
    return nothing
  end

  return R(st + 1), st + 1
end

Base.eltype(::Type{GaloisField}) = gfp_elem

Base.IteratorSize(::Type{GaloisField}) = Base.HasLength()

Base.length(R::GaloisField) = R.n

function powmod(f::Zmodn_poly, e::fmpz, g::Zmodn_poly)
  if nbits(e) <= 63
    return powmod(f, Int(e), g)
  else
    _e = BigInt()
    z = parent(f)()
    ccall((:fmpz_get_mpz, :libflint), Nothing, (Ref{BigInt}, Ref{fmpz}), _e, e)
    ccall((:nmod_poly_powmod_mpz_binexp, :libflint), Nothing,
          (Ref{Zmodn_poly}, Ref{Zmodn_poly}, Ref{BigInt}, Ref{Zmodn_poly}),
           z, f, e, g)
    return z
  end
end


text:
function get_sphere(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    radius=0.5,
    contact=true,
    contact_type=:nonlinear,
    T=Float64)

    origin = Origin{T}(name=:origin)
    mass = 1.0
    bodies = [Sphere(radius, mass, name=:sphere)]
    joints = [JointConstraint(Floating(origin, bodies[1]), name=:floating_joint)]
    mechanism = Mechanism(origin, bodies, joints,
        timestep=timestep,
        gravity=gravity)

    if contact
        contact = [0.0, 0.0, 0.0]
        normal = [0.0, 0.0, 1.0]
        contacts = [contact_constraint(get_body(mechanism, :sphere), normal,
            friction_coefficient=friction_coefficient,
            contact_origin=contact,
            contact_radius=radius,
            contact_type=contact_type)]
        set_minimal_coordinates!(mechanism, get_joint(mechanism, :floating_joint), [0.0; 0.0; radius; zeros(3)])
        mechanism = Mechanism(origin, bodies, joints, contacts,
            gravity=gravity,
            timestep=timestep)
    end
    return mechanism
end

function initialize_sphere!(mechanism::Mechanism{T};
    x=zeros(3),
    q=one(Quaternion),
    v=zeros(3),
    ω=zeros(3)) where T

    r = mechanism.bodies[1].shape.r
    joint = get_joint(mechanism, :floating_joint)
    zero_velocity!(mechanism)
    set_minimal_coordinates!(mechanism, joint, [x + [0.0, 0.0, r] rotation_vector(q)])
    set_minimal_velocities!(mechanism, joint, [v; ω])
end


text:
# script to run inference for the OU SDEMEM model
using Pkg
using LinearAlgebra
using DataFrames
import Statistics.mean
import Statistics.std
using Printf

# load functions
include(pwd()*"/src/SDEMEM OU neuron data/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU neuron data/mcmc.jl")

M_subjects = 100

seed = 100 #parse(Int,ARGS[1])

y,prior_parameters_η,prior_parameters_σ_ϵ = set_up(M_subjects,seed)

# plot data
using PyPlot

# Plot data (this plot should be similar to Fig 3 in the paper)
PyPlot.figure()
for i in 1:M_subjects
    PyPlot.plot(y[i].mSec, y[i].mV, "k", alpha = 0.2)
end
PyPlot.xlabel("Time (msec)")
PyPlot.ylabel("depolarization mV")

nbr_obs = 0
for i in 1:M_subjects
    global nbr_obs = nbr_obs + length(y[i].mSec)
end

nbr_obs

startval_ϕ = ones(M_subjects,3)

for j = 1:3

    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_ϕ[:,j,1] = μ_0_j*ones(M_subjects)
    #startval_ϕ[:,j,1] = μ_0_j*ones(M_subjects) + randn(M_subjects)*0.05

end

startval_log_σ_ϵ = log(0.2)

posterior_mean_kalman = CSV.read("data/SDEMEM OU neuron data/posterior_mean_sigma_epsilon_and_phi.csv")

log_σ_ϵ  = posterior_mean_kalman[1,1]
ϕ  = Matrix(posterior_mean_kalman[2:end,:])

# set time step
dt = diff(y[1].mSec)[1]

# run kalman filter + plot paths

ll1, path1 = @time kalman(y, exp(log_σ_ϵ), ϕ, dt, true)

loglik_kalman = sum(ll1)

# bootstrap filter

function run_bootstrap(particels)

    nbr_particles = particels*ones(Int64, M_subjects)
    ρ = 0.0

    # set number of obs for each data set
    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y[i].mV); end

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    run_time = @elapsed ll1,path1, nbr_resamlping =  cpf(y, exp(log_σ_ϵ), ϕ, dt, u_prop_old, u_resample_old, nbr_particles, true, true)

    return ll1, run_time, path1,nbr_resamlping

end

ll1, run_time, path1, nbr_resamlping = run_bootstrap(500)

sum(ll1)

# bridge pf

function run_bridge_pf(particels)

    nbr_particles = particels*ones(Int64, M_subjects)
    ρ = 0.0

    # set number of obs for each data set
    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y[i].mV); end

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    run_time = @elapsed ll1,path1,nbr_resamlping =  bridge_pf(y, exp(log_σ_ϵ), ϕ, dt, u_prop_old, u_resample_old, nbr_particles, true, true)

    return ll1,run_time,path1,nbr_resamlping
end

ll1b,run_timeb,path1b,nbr_resamlping = run_bridge_pf(500)

sum(ll1b)

subject = rand(1:M_subjects)
PyPlot.figure(figsize=(8,5))
PyPlot.subplot(121)
PyPlot.plot(path1[subject][:,:]', "r")
PyPlot.plot(y[subject].mV, "g")
PyPlot.xlabel("Time (msec)")
PyPlot.ylabel("depolarization mV")
PyPlot.title("Bootstrap")
PyPlot.subplot(122)
PyPlot.plot(path1b[subject][:,:]', "r")
PyPlot.plot(y[subject].mV, "g")
PyPlot.title("Bridge")
PyPlot.xlabel("Time (msec)")
PyPlot.ylabel("depolarization mV")
PyPlot.savefig("figures/boostrap_bridge_tracking.pdf", format="pdf", dpi=1000)
run(`pdftops -eps figures/boostrap_bridge_tracking.pdf figures/boostrap_bridge_tracking.eps`)

# randomly selected subject is subject 60
subject = 60

y[subject].mSec

PyPlot.figure(figsize=(12,5))
PyPlot.subplot(131)
PyPlot.plot(y[subject].mSec, y[subject].mV, "k")
PyPlot.title("Data")
PyPlot.xlabel("Time (msec)")
PyPlot.ylabel("depolarization mV")
PyPlot.subplot(132)
PyPlot.plot(y[subject].mSec,path1[subject][:,:]', "k")
PyPlot.xlabel("Time (msec)")
PyPlot.ylabel("depolarization mV")
PyPlot.title("Bootstrap")
PyPlot.subplot(133)
PyPlot.plot(y[subject].mSec,path1b[subject][:,:]', "k")
PyPlot.title("Bridge")
PyPlot.xlabel("Time (msec)")
PyPlot.ylabel("depolarization mV")
PyPlot.savefig("figures/boostrap_bridge_tracking.pdf", format="pdf", dpi=1000)
run(`pdftops -eps figures/boostrap_bridge_tracking.pdf figures/boostrap_bridge_tracking.eps`)

# test variance of loglik est

N_filters = 100

data = zeros(4, N_filters)
run_time_kalman = zeros(N_filters)
ll_all_boostrap = zeros(M_subjects, N_filters)
ll_all_bridge = zeros(M_subjects, N_filters)

for i in 1:N_filters

    println(i)

    ll_boostrap, run_time_boostrap, path1 = run_bootstrap(500)
    ll1_bridge, run_time_bridge, path1 = run_bridge_pf(500)

    ll_all_boostrap[:,i] = ll_boostrap
    ll_all_bridge[:,i] = ll1_bridge

    data[:,i] = [sum(ll_boostrap);run_time_boostrap; sum(ll1_bridge);run_time_bridge]

end

for i in 1:N_filters
    run_time_kalman[i] = @elapsed ll1, path1 = kalman(y, exp(log_σ_ϵ), ϕ, dt, true)
end

@printf "----------------\n"
@printf "Loglik kalman: %.4f\n" round(loglik_kalman, digits = 0)
@printf "Kalman - runtime (mean): %.4f\n" round(mean(run_time_kalman), digits = 4)
@printf "Bootstrap (mean): %.4f\n" round(mean(data[1,:]), digits = 0)
@printf "Bootstrap (var): %.4f\n" round(var(data[1,:]), digits = 0)
@printf "Bootstrap (std): %.4f\n" round(std(data[1,:]), digits = 0)
@printf "Bootstrap (var) ind ll est:\n"
println(var(ll_all_boostrap, dims = 2))
@printf "Bootstrap - runtime (mean): %.4f\n" round(mean(data[2,:]), digits = 4)
@printf "Bridge (mean): %.4f\n" round(mean(data[3,:]), digits = 0)
@printf "Bridge (var): %.4f\n" round(var(data[3,:]), digits = 0)
@printf "Bridge (std): %.4f\n" round(std(data[3,:]), digits = 0)
@printf "Bridge (var) ind ll est:\n"
println(var(ll_all_bridge, dims = 2))
@printf "Bridge - runtime (mean): %.4f\n" round(mean(data[4,:]), digits = 4)

max_stars_repo_path: src/SDEMEM OU neuron data/analyses_particle_filters.jl


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int)::Int = (n - 5) ÷ 7 * 3 + 2

compute(8)

compute(1000000)

@benchmark compute(1000000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Pkg
Pkg.add( "Interact" )
Pkg.add( "ReadStat" )
Pkg.add( "DataVoyager" )
Pkg.add( "TableTraits" )
Pkg.add( "Juno" )
Pkg.add( "HTTP" )
Pkg.add( "IteratorInterfaceExtensions" )
Pkg.add( "VegaDatasets" )
Pkg.add( "DataStructures" )
Pkg.add( "StatsFuns" )
Pkg.add( "DataFrames" )
Pkg.add( "JuliaDB" )
Pkg.add( "StatsModels" )
Pkg.add( "Distributions" )
Pkg.add( "ZipFile" )
Pkg.add( "VegaLite" )
Pkg.add( "CodeTracking" )
Pkg.add( "Missings" )
Pkg.add( "Plotly" )
Pkg.add( "Formatting" )
Pkg.add( "Queryverse" )
Pkg.add( "JDBC" )
Pkg.add( "HDF5" )
Pkg.add( "Gadfly" )
Pkg.add( "FileIO" )
Pkg.add( "Printf" )
Pkg.add( "MySQL" )
Pkg.add( "IterableTables" )
Pkg.add( "CSV" )
Pkg.add( "Genie" )
Pkg.add( "Revise" )
Pkg.add( "Test" )
Pkg.add( "CSVFiles" )
Pkg.add( "Atom" )
Pkg.add( "Optim" )
Pkg.add( "StatPlots" )
Pkg.add( "Dagger" )
Pkg.add( "Memento" )
Pkg.add( "DataFramesMeta" )
Pkg.add( "JSON" )
Pkg.add( "StatsBase" )
Pkg.add( "IJulia" )
Pkg.add( "Plots" )
Pkg.add( "PyPlot" )
Pkg.add( "Profile" )
Pkg.add( "Pandas" )
Pkg.add( "GLM" )
Pkg.add( "Documenter" )
Pkg.add( "Query" )
Pkg.add( "LineSearches" )
Pkg.add( "DelimitedFiles" )
Pkg.add( "IniFile" )
Pkg.add( "StatFiles" )
Pkg.add( "NLsolve" )
Pkg.add( "DataValues" )
Pkg.add( "LightXML" )


text:
fs2()

@testset "figalt, fig" begin
    mkpath(joinpath(F.PATHS[:site], "assets"))
    write(joinpath(F.PATHS[:site], "assets", "testimg.png"), "png code")
    h = raw"""
        A figure:
        \figalt{fig 1}{/assets/testimg.png}
        \fig{/assets/testimg.png}
        \fig{/assets/testimg}
        Done.
        """ |> seval
    @test isapproxstr(h, """

A figure:

            Done.

            """)
    p = mkpath(joinpath(F.PATHS[:site], "assets", "output"))
    write(joinpath(p, "testimg_2.png"), "png code")
    h = raw"""
        Another figure:
        \figalt{fig blah "hello!"}{/assets/testimg_2.png}
        """ |> seval
    @test isapproxstr(h, """

Another figure:

            """)
    h = raw"""
        No fig:
        \fig{/assets/testimg_3.png}
        """ |> seval
    @test isapproxstr(h, """

No fig: $(F.html_err("Image matching '/assets/testimg_3.png' not found."))

            """)
end

@testset "table" begin
    #
    # has header in source
    #
    testcsv = "h1,h2,h3\nstring1, 1.567, 0\n,,\n l i n e ,.158,99999999"
    write(joinpath(F.PATHS[:site], "assets", "testcsv.csv"), testcsv)
    # no header specified
    h = raw"""
        A table:
        \tableinput{}{/assets/testcsv.csv}
        Done.
        """ |> seval

    # NOTE: in VERSION > 1.4, Markdown has alignment for tables:
    # -- https://github.com/JuliaLang/julia/pull/33849
    if VERSION >= v"1.4.0-"
        shouldbe = """

A table:

h1	h2	h3

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    else
        shouldbe = """

A table:

h1	h2	h3

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    end
    @test isapproxstr(h, shouldbe)
    # header specified
    h = raw"""
        A table:
        \tableinput{A,B,C}{/assets/testcsv.csv}
        Done.
        """ |> seval
    if VERSION >= v"1.4.0-"
        shouldbe = """

A table:

A	B	C

h1	h2	h3

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    else
        shouldbe = """

A table:

A	B	C

h1	h2	h3

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    end
    @test isapproxstr(h, shouldbe)
    # wrong header
    h = raw"""
        A table:
        \tableinput{,}{/assets/testcsv.csv}
        Done.
        """ |> seval
    shouldbe = """
A table:
// In `\\tableinput`: header size (2) and number of columns (3) do not match. //

            Done.
"""
    @test isapproxstr(h, shouldbe)

    #
    # does not have header in source
    #

    testcsv = "string1, 1.567, 0\n,,\n l i n e ,.158,99999999"
    write(joinpath(F.PATHS[:site], "assets", "testcsv.csv"), testcsv)
    # no header specified
    h = raw"""
        A table:
        \tableinput{}{/assets/testcsv.csv}
        Done.
        """ |> seval
    if VERSION >= v"1.4.0-"
        shouldbe = """

A table:

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    else
        shouldbe = """

A table:

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    end
    @test isapproxstr(h, shouldbe)
    # header specified
    h = raw"""
        A table:
        \tableinput{A,B,C}{/assets/testcsv.csv}
        Done.
        """ |> seval

    if VERSION >= v"1.4.0-"
        shouldbe = """

A table:

A	B	C

string1	1.567	0

l i n e	.158	99999999

            Done.
"""
    else
        shouldbe = """
A table:

A	B	C

string1	1.567	0

l i n e	.158	99999999

                Done.
"""
    end
    @test isapproxstr(h, shouldbe)
    # wrong header
    h = raw"""
        A table:
        \tableinput{A,B}{/assets/testcsv.csv}
        Done.
        """ |> seval

    shouldbe = """
A table:
// In `\\tableinput`: header size (2) and number of columns (3) do not match. //

            Done.
"""
    @test isapproxstr(h, shouldbe)
end


text:
using Thebes, Luxor, Random

function main()
    Drawing(800, 800, "test3.png")
    background("azure")
    origin()

    eyepoint(Point3D(500, 500, 300))
    perspective(1500)
    sethue("pink")
    carpet(400)
    setline(0.2)

    for k in 1:4:20
        for i in 1:20
            object = make(Pyramid)
            objects = [object]
            s = rescale(k, 5, 20, 10, 20)
            scaleby!.(objects, s, s, s)
            moveby!.(objects, rescale(k, 1, 20, 50, 400), 1, 0)
            a = rescale(i, 1, 20, 0, 2pi)
            rotateby!.(objects, 0, 0, a)
            sortfaces!(objects, eyepoint = eyepoint)
            # default is stroke at the moment
            pin.(objects)
        end
    end
    finish()
end

main()


text:
module Harlequin

export use_mpi

using Requires

# Defining use_mpi as a function instead of a constant allows the
# Julia compiler to completely optimize out "if" statements like
#
#     if use_mpi() && ...
#         ...
#     end
#
# A "const" would do the same, but we would not be able to redefine it
# in the "@require" statement below.

@doc raw"""
    use_mpi()

Return `true` if Harlequin is taking advantage of MPI, `false` if not.

"""
use_mpi() = false

function __init__()
    @require MPI="da04e1cc-30fd-572f-bb4f-1f8673147195" include("mpi.jl")
end

include("quaternions.jl")
include("genpointings.jl")
include("dipole.jl")
include("beams.jl")
include("mapmaking.jl")

end # module


text:
for p in ("Knet","AutoGrad","ArgParse","Compat")
    Pkg.installed(p) == nothing && Pkg.add(p)
end

"""
charlm.jl: Knet8 version (c) Emre Yolcu, Deniz Yuret, 2016

This example implements an LSTM network for training and testing
character-level language models inspired by ["The Unreasonable
Effectiveness of Recurrent Neural
Networks"](http://karpathy.github.io/2015/05/21/rnn-effectiveness) from
the Andrej Karpathy blog.  The model can be trained with different
genres of text, and can be used to generate original text in the same
style.

Example usage:

* `julia charlm.jl`: trains a model using its own code.

* `julia charlm.jl --data foo.txt`: uses foo.txt to train instead.

* `julia charlm.jl --data foo.txt bar.txt`: uses foo.txt for training
  and bar.txt for validation.  Any number of files can be specified,
  the first two will be used for training and validation, the rest for
  testing.

* `julia charlm.jl --best foo.jld --save bar.jld`: saves the best
  model (according to validation set) to foo.jld, last model to
  bar.jld.

* `julia charlm.jl --load foo.jld --generate 1000`: generates 1000
  characters from the model in foo.jld.

* `julia charlm.jl --help`: describes all available options.

"""
module CharLM

using Knet,AutoGrad,ArgParse,Compat

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="charlm.jl (c) Emre Yolcu, Deniz Yuret, 2016. Character level language model based on http://karpathy.github.io/2015/05/21/rnn-effectiveness."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--datafiles"; nargs='+'; help="If provided, use first file for training, second for dev, others for test.")
        ("--loadfile"; help="Initialize model from file")
        ("--savefile"; help="Save final model to file")
        ("--bestfile"; help="Save best model to file")
        ("--generate"; arg_type=Int; default=0; help="If non-zero generate given number of characters.")
        ("--hidden"; nargs='+'; arg_type=Int; default=[256]; help="Sizes of one or more LSTM layers.")
        ("--embed"; arg_type=Int; default=256; help="Size of the embedding vector.")
        ("--epochs"; arg_type=Int; default=3; help="Number of epochs for training.")
        ("--batchsize"; arg_type=Int; default=128; help="Number of sequences to train on in parallel.")
        ("--seqlength"; arg_type=Int; default=100; help="Number of steps to unroll the network for.")
        ("--decay"; arg_type=Float64; default=0.9; help="Learning rate decay.")
        ("--lr"; arg_type=Float64; default=4.0; help="Initial learning rate.")
        ("--gclip"; arg_type=Float64; default=3.0; help="Value to clip the gradient norm at.")
        ("--winit"; arg_type=Float64; default=0.3; help="Initial weights set to winit*randn().")
        ("--gcheck"; arg_type=Int; default=0; help="Check N random gradients.")
        ("--seed"; arg_type=Int; default=-1; help="Random number seed.")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        #TODO ("--dropout"; arg_type=Float64; default=0.0; help="Dropout probability.")
    end
    println(s.description)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && srand(o[:seed])
    o[:atype] = eval(parse(o[:atype]))
    if any(f->(o[f]!=nothing), (:loadfile, :savefile, :bestfile))
        Pkg.installed("JLD")==nothing && Pkg.add("JLD") # error("Please Pkg.add(\"JLD\") to load or save files.")
        eval(Expr(:using,:JLD))
    end

    # we initialize a model from loadfile, train using datafiles (both optional).
    # if the user specifies neither, train a model using the charlm.jl source code.
    isempty(o[:datafiles]) && o[:loadfile]==nothing && push!(o[:datafiles],@__FILE__) # shakespeare()

    # read text and report lengths
    text = map((@compat readstring), o[:datafiles])
    !isempty(text) && info("Chars read: $(map((f,c)->(basename(f),length(c)),o[:datafiles],text))")

    # vocab (char_to_index) comes from the initial model if there is one, otherwise from the datafiles.
    # if there is an initial model make sure the data has no new vocab
    if o[:loadfile]==nothing
        vocab = Dict{Char,Int}()
        for t in text, c in t; get!(vocab, c, 1+length(vocab)); end
        model = initweights(o[:atype], o[:hidden], length(vocab), o[:embed], o[:winit])
    else
        info("Loading model from $(o[:loadfile])")
        vocab = load(o[:loadfile], "vocab")
        for t in text, c in t; haskey(vocab, c) || error("Unknown char $c"); end
        model = map(p->convert(o[:atype],p), load(o[:loadfile], "model"))
    end
    info("$(length(vocab)) unique chars.")
    if !isempty(text)
        train!(model, text, vocab, o)
    end
    if o[:savefile] != nothing
        info("Saving last model to $(o[:savefile])")
        save(o[:savefile], "model", model, "vocab", vocab)
    end
    if o[:generate] > 0
        state = initstate(o[:atype],o[:hidden],1)
        generate(model, state, vocab, o[:generate])
    end
end

function train!(model, text, vocab, o)
    s0 = initstate(o[:atype], o[:hidden], o[:batchsize])
    data = map(t->minibatch(t, vocab, o[:batchsize]), text)
    lr = o[:lr]
    if o[:fast]
        @time (for epoch=1:o[:epochs]
               train1(model, copy(s0), data[1]; slen=o[:seqlength], lr=lr, gclip=o[:gclip])
               end; gpu()>=0 && Knet.cudaDeviceSynchronize())
        return
    end
    losses = map(d->loss(model,copy(s0),d), data)
    println((:epoch,0,:loss,losses...))
    devset = ifelse(length(data) > 1, 2, 1)
    devlast = devbest = losses[devset]
    for epoch=1:o[:epochs]
        @time train1(model, copy(s0), data[1]; slen=o[:seqlength], lr=lr, gclip=o[:gclip])
        @time losses = map(d->loss(model,copy(s0),d), data)
        println((:epoch,epoch,:loss,losses...))
        if o[:gcheck] > 0
            gradcheck(loss, model, copy(s0), data[1], 1:o[:seqlength]; gcheck=o[:gcheck])
        end
        devloss = losses[devset]
        if devloss < devbest
            devbest = devloss
            if o[:bestfile] != nothing
                info("Saving best model to $(o[:bestfile])")
                save(o[:bestfile], "model", model, "vocab", vocab)
            end
        end
        if devloss > devlast
            lr *= o[:decay]
            info("New learning rate: $lr")
        end
        devlast = devloss
    end
end

# sequence[t]: input token at time t
# state is modified in place
function train1(param, state, sequence; slen=100, lr=1.0, gclip=0.0)
    for t = 1:slen:length(sequence)-slen
        range = t:t+slen-1
        gloss = lossgradient(param, state, sequence, range)
        gscale = lr
        if gclip > 0
            gnorm = sqrt(mapreduce(sumabs2, +, 0, gloss))
            if gnorm > gclip
                gscale *= gclip / gnorm
            end
        end
        for k in 1:length(param)
            # param[k] -= gscale * gloss[k]
            axpy!(-gscale, gloss[k], param[k])
        end
        isa(state,Vector{Any}) || error("State should not be Boxed.")
        # The following is needed in case AutoGrad boxes state values during gradient calculation
        for i = 1:length(state)
            state[i] = AutoGrad.getval(state[i])
        end
    end
end

# param[2k-1,2k]: weight and bias for the k'th lstm layer
# param[end-2]: embedding matrix
# param[end-1,end]: weight and bias for final prediction
function initweights(atype, hidden, vocab, embed, winit)
    param = Array(Any, 2*length(hidden)+3)
    input = embed
    for k = 1:length(hidden)
        param[2k-1] = winit*randn(input+hidden[k], 4*hidden[k])
        param[2k]   = zeros(1, 4*hidden[k])
        param[2k][1:hidden[k]] = 1 # forget gate bias
        input = hidden[k]
    end
    param[end-2] = winit*randn(vocab,embed)
    param[end-1] = winit*randn(hidden[end],vocab)
    param[end] = zeros(1,vocab)
    return map(p->convert(atype,p), param)
end

# state[2k-1,2k]: hidden and cell for the k'th lstm layer
function initstate(atype, hidden, batchsize)
    state = Array(Any, 2*length(hidden))
    for k = 1:length(hidden)
        state[2k-1] = zeros(batchsize,hidden[k])
        state[2k] = zeros(batchsize,hidden[k])
    end
    return map(s->convert(atype,s), state)
end

function lstm(weight,bias,hidden,cell,input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

# s[2k-1,2k]: hidden and cell for the k'th lstm layer
# w[2k-1,2k]: weight and bias for k'th lstm layer
# w[end-2]: embedding matrix
# w[end-1,end]: weight and bias for final prediction
# state is modified in place
function predict(w, s, x)
    x = x * w[end-2]
    for i = 1:2:length(s)
        (s[i],s[i+1]) = lstm(w[i],w[i+1],s[i],s[i+1],x)
        x = s[i]
    end
    return x * w[end-1] .+ w[end]
end

# sequence[t]: input token at time t
# state is modified in place
function loss(param,state,sequence,range=1:length(sequence)-1)
    total = 0.0; count = 0
    atype = typeof(AutoGrad.getval(param[1]))
    input = convert(atype,sequence[first(range)])
    for t in range
        ypred = predict(param,state,input)
        ynorm = logp(ypred,2) # ypred .- log(sum(exp(ypred),2))
        ygold = convert(atype,sequence[t+1])
        total += sum(ygold .* ynorm)
        count += size(ygold,1)
        input = ygold
    end
    return -total / count
end

lossgradient = grad(loss)

function generate(param, state, vocab, nchar)
    index_to_char = Array(Char, length(vocab))
    for (k,v) in vocab; index_to_char[v] = k; end
    input = oftype(param[1], zeros(1,length(vocab)))
    index = 1
    for t in 1:nchar
        ypred = predict(param,state,input)
        input[index] = 0
        index = sample(exp(logp(ypred)))
        print(index_to_char[index])
        input[index] = 1
    end
    println()
end

function sample(p)
    p = convert(Array,p)
    r = rand()
    for c = 1:length(p)
        r -= p[c]
        r < 0 && return c
    end
end

function shakespeare()
    file = Knet.dir("data","100.txt")
    if !isfile(file)
        info("Downloading 'The Complete Works of William Shakespeare'")
        url = "http://www.gutenberg.org/files/100/100.txt"
        download(url,file)
    end
    return file
end

function minibatch(chars, char_to_index, batch_size)
    nbatch = div(length(chars), batch_size)
    vocab_size = length(char_to_index)
    data = [ falses(batch_size, vocab_size) for i=1:nbatch ] # using BitArrays
    cidx = 0
    for c in chars            # safest way to iterate over utf-8 text
        idata = 1 + cidx % nbatch
        row = 1 + div(cidx, nbatch)
        row > batch_size && break
        col = char_to_index[c]
        data[idata][row,col] = 1
        cidx += 1
    end
    return data
end

# To be able to load/save KnetArrays:
if Pkg.installed("JLD") != nothing
    import JLD: writeas, readas
    type KnetJLD; a::Array; end
    writeas(c::KnetArray) = KnetJLD(Array(c))
    readas(d::KnetJLD) = KnetArray(d.a)
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia charlm.jl --epochs 10
# julia> CharLM.main("--epochs 10")
!isinteractive() && !isdefined(Core.Main,:load_only) && main(ARGS)

end  # module

# Note: 10.txt used in the sample runs below was generated using
#   head -10000 100.txt > 10.txt
# where 100.txt is the file downloaded by shakespeare().

### SAMPLE RUN 74a2e6c+ Mon Sep 19 14:03:10 EEST 2016
### Implemented multi-layer.  Removed the keepstate option fixing it to true.
### Note that winit default changed so I specify it below for comparison.
### The slight difference is due to keepstate.

# julia> CharLM.main("--data 10.txt --winit 0.3 --fast")
# charlm.jl (c) Emre Yolcu, Deniz Yuret, 2016. Character level language model based on http://karpathy.github.io/2015/05/21/rnn-effectiveness.
# opts=(:lr,4.0)(:atype,"KnetArray{Float32}")(:winit,0.3)(:savefile,nothing)(:loadfile,nothing)(:generate,0)(:bestfile,nothing)(:gclip,3.0)(:hidden,[256])(:epochs,3)(:decay,0.9)(:gcheck,0)(:seqlength,100)(:seed,42)(:embed,256)(:batchsize,128)(:datafiles,Any["10.txt"])(:fast,true)
# INFO: Chars read: [("10.txt",425808)]
# INFO: 87 unique chars.
#   1.406687 seconds (1.74 M allocations: 196.799 MB, 2.49% gc time)
# (:epoch,0,:loss,6.1075509900258)
#   4.002638 seconds (6.12 M allocations: 374.188 MB, 2.41% gc time)
#   3.990772 seconds (6.11 M allocations: 374.129 MB, 2.44% gc time)
#   4.006249 seconds (6.12 M allocations: 374.240 MB, 2.53% gc time)
#   1.405878 seconds (1.75 M allocations: 197.059 MB, 2.56% gc time)
# (:epoch,3,:loss,1.8713183968407767)

### SAMPLE RUN 4ce58d1+ Fri Sep 16 12:24:00 EEST 2016
### Transposed everything so getindex does not need to copy

# charlm.jl (c) Emre Yolcu, Deniz Yuret, 2016. Character level language model based on http://karpathy.github.io/2015/05/21/rnn-effectiveness.
# opts=(:keepstate,false)(:lr,4.0)(:atype,"KnetArray{Float32}")(:winit,0.3)(:savefile,nothing)(:loadfile,nothing)(:generate,0)(:bestfile,nothing)(:gclip,3.0)(:embedding,256)(:hidden,256)(:epochs,3)(:decay,0.9)(:gcheck,0)(:seqlength,100)(:seed,42)(:batchsize,128)(:datafiles,Any["data/10.txt"])(:fast,true)
# INFO: Chars read: [("10.txt",425808)]
# INFO: 87 unique chars.
#   1.388652 seconds (1.73 M allocations: 196.686 MB, 2.04% gc time)
# (:epoch,0,:loss,6.1075509900258)
#   3.940298 seconds (6.06 M allocations: 373.166 MB, 2.06% gc time)
#   3.935995 seconds (6.06 M allocations: 373.244 MB, 2.07% gc time)
#   3.934983 seconds (6.06 M allocations: 373.245 MB, 2.08% gc time)
#   1.390374 seconds (1.73 M allocations: 196.820 MB, 2.11% gc time)
# (:epoch,3,:loss,1.860654126576015)

### SAMPLE RUN 31136d5+ Wed Sep 14 17:51:44 EEST 2016: using vcat(x,h) and vcat(w...)
### optimized learning parameters: winit=0.3, lr=4.0, gclip=3.0

# charlm.jl (c) Emre Yolcu, Deniz Yuret, 2016. Character level language model based on http://karpathy.github.io/2015/05/21/rnn-effectiveness.
# opts=(:keepstate,false)(:lr,4.0)(:atype,"KnetArray{Float32}")(:winit,0.3)(:savefile,nothing)(:loadfile,nothing)(:generate,0)(:bestfile,nothing)(:gclip,3.0)(:embedding,256)(:hidden,256)(:epochs,3)(:decay,0.9)(:gcheck,0)(:seqlength,100)(:seed,42)(:batchsize,128)(:datafiles,Any["10.txt"])(:fast,true)
# INFO: Chars read: [("10.txt",425808)]
# INFO: 87 unique chars.
#   1.596959 seconds (1.79 M allocations: 197.432 MB, 2.56% gc time)
# (:epoch,0,:loss,5.541976199042528)
#   4.421566 seconds (6.23 M allocations: 375.843 MB, 2.54% gc time)
#   4.418540 seconds (6.25 M allocations: 376.058 MB, 2.51% gc time)
#   4.402317 seconds (6.26 M allocations: 376.297 MB, 2.66% gc time)
#   1.594677 seconds (1.81 M allocations: 197.737 MB, 2.64% gc time)
# (:epoch,3,:loss,1.8484550957572192)

### SAMPLE RUN 80503e7+ Wed Sep 14 17:35:36 EEST 2016: using vcat(x,h)
#
# charlm.jl (c) Emre Yolcu, Deniz Yuret, 2016. Character level language model based on http://karpathy.github.io/2015/05/21/rnn-effectiveness.
# opts=(:keepstate,false)(:lr,1.0)(:atype,"KnetArray{Float32}")(:savefile,nothing)(:loadfile,nothing)(:generate,0)(:bestfile,nothing)(:embedding,256)(:gclip,5.0)(:hidden,256)(:epochs,3)(:decay,0.9)(:gcheck,0)(:seqlength,100)(:seed,42)(:batchsize,128)(:datafiles,Any["10.txt"])(:fast,true)
# INFO: Chars read: [("10.txt",425808)]
# INFO: 87 unique chars.
#   1.930180 seconds (1.95 M allocations: 213.741 MB, 1.82% gc time)
# (:epoch,0,:loss,4.462641664662756)
#   4.968101 seconds (7.47 M allocations: 454.259 MB, 2.24% gc time)
#   4.963733 seconds (7.47 M allocations: 454.363 MB, 2.26% gc time)
#   4.967413 seconds (7.45 M allocations: 454.024 MB, 2.14% gc time)
#   1.945658 seconds (1.98 M allocations: 214.183 MB, 2.02% gc time)
# (:epoch,3,:loss,3.2389672966290237)

### SAMPLE RUN 65f57ff+ Wed Sep 14 10:02:30 EEST 2016: separate x, h, w, b
#
# charlm.jl (c) Emre Yolcu, Deniz Yuret, 2016. Character level language model based on http://karpathy.github.io/2015/05/21/rnn-effectiveness.
# opts=(:keepstate,false)(:lr,1.0)(:atype,"KnetArray{Float32}")(:savefile,nothing)(:loadfile,nothing)(:generate,0)(:bestfile,nothing)(:embedding,256)(:gclip,5.0)(:hidden,256)(:epochs,3)(:decay,0.9)(:gcheck,0)(:seqlength,100)(:seed,42)(:batchsize,128)(:datafiles,Any["10.txt"])(:fast,true)
# INFO: Chars read: [("10.txt",425808)]
# INFO: 87 unique chars.
#   2.156358 seconds (2.31 M allocations: 237.913 MB, 2.30% gc time)
# (:epoch,0,:loss,4.465127425659868)
#   6.287736 seconds (9.54 M allocations: 574.703 MB, 2.84% gc time)
#   6.272144 seconds (9.54 M allocations: 574.633 MB, 2.80% gc time)
#   6.277462 seconds (9.54 M allocations: 574.637 MB, 2.86% gc time)
#   2.165516 seconds (2.34 M allocations: 238.323 MB, 2.56% gc time)
# (:epoch,3,:loss,3.226540256084356)

### SAMPLE OUTPUT (with head -10000 100.txt): first version
# julia> CharLM.main("--gpu --data 10.txt")
# opts=(:lr,1.0)(:savefile,nothing)(:loadfile,nothing)(:dropout,0.0)(:generate,0)(:bestfile,nothing)(:embedding,256)(:gclip,5.0)(:hidden,256)(:epochs,10)(:nlayer,1)(:decay,0.9)(:gpu,true)(:seqlength,100)(:seed,42)(:batchsize,128)(:datafiles,Any["10.txt"])
# INFO: Chars read: [("10.txt",425808)]
# INFO: 87 unique chars.
# (0,4.465127425659868)
#   2.182693 seconds (2.36 M allocations: 240.394 MB, 1.74% gc time)
#   7.861079 seconds (10.12 M allocations: 601.311 MB, 1.84% gc time)
# (1,3.3244698245543285)
#   2.159062 seconds (2.35 M allocations: 239.217 MB, 1.80% gc time)
#   6.200085 seconds (9.55 M allocations: 575.573 MB, 2.28% gc time)
# (2,3.24593908969621)
#   2.352389 seconds (2.34 M allocations: 239.381 MB, 1.51% gc time)
#   6.211946 seconds (9.55 M allocations: 575.568 MB, 2.21% gc time)
# (3,3.226540256084356)


text:
#
# Date created: 2022-03-16
# Author: aradclif
#
#
############################################################################################
function findmax_loops(N::Int, D)
    b4 = preexpr2(:typemin)
    inner = innerloop(N, D)
    push!(inner.args, maxblock3(N, D))
    push!(b4.args, inner)
    push!(b4.args, postexpr3(N, D).args...)
    outer = outerloop(N, D)
    push!(outer.args, b4)
    outer
end
function findmax_loops_innerturbo(N::Int, D, multithreaded::Bool)
    b4 = preexpr2(:typemin)
    inner = innerloop(N, D)
    push!(inner.args, maxblock3(N, D))
    push!(b4.args, Expr(:macrocall, multithreaded ? Symbol("@tturbo") : Symbol("@turbo"), (), inner))
    push!(b4.args, postexpr3(N, D).args...)
    outer = outerloop(N, D)
    push!(outer.args, b4)
    outer
end

function bfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = Expr(:block, ntuple(d -> offsetk(d + 2), N - 2)...)
    b3 = Expr(:block, totaloffsetraw(N), Expr(:(=), :Dstar, Expr(:call, :-, :Dstar)))
    outer = findmax_loops(N, D)
    return quote
        $b1
        $b2
        $b3
        $outer
    end
end
function ivfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = Expr(:block, ntuple(d -> offsetk(d + 2), N - 2)...)
    b3 = Expr(:block, totaloffsetraw(N), Expr(:(=), :Dstar, Expr(:call, :-, :Dstar)))
    outer = findmax_loops_innerturbo(N, D, false)
    return quote
        $b1
        $b2
        $b3
        $outer
    end
end
function ivtfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = Expr(:block, ntuple(d -> offsetk(d + 2), N - 2)...)
    b3 = Expr(:block, totaloffsetraw(N), Expr(:(=), :Dstar, Expr(:call, :-, :Dstar)))
    outer = findmax_loops_innerturbo(N, D, true)
    return quote
        $b1
        $b2
        $b3
        $outer
    end
end

function vfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = Expr(:block, ntuple(d -> offsetk(d + 2), N - 2)...)
    b3 = Expr(:block, totaloffsetraw(N), Expr(:(=), :Dstar, Expr(:call, :-, :Dstar)))
    outer = findmax_loops(N, D)
    return quote
        $b1
        $b2
        $b3
        @turbo $outer
    end
end

function vtfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = Expr(:block, ntuple(d -> offsetk(d + 2), N - 2)...)
    b3 = Expr(:block, totaloffsetraw(N), Expr(:(=), :Dstar, Expr(:call, :-, :Dstar)))
    outer = findmax_loops(N, D)
    return quote
        $b1
        $b2
        $b3
        @tturbo $outer
    end
end

@generated _bfindmax!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                      A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D} = bfindmax_quote(N, D)
function bfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
    Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
    B = similar(A, Dᴮ)
    C = similar(A, Int, Dᴮ)
    _bfindmax!(B, C, A, Dᴮ′)
    B, CartesianIndices(A)[C]
end

@generated _vfindmax!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                      A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D} = vfindmax_quote(N, D)
function vfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        return vfindmax1(A)
    elseif 1 ∈ dims
        return ivfindmax(A, dims)
    else
        Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
        Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
        B = similar(A, Dᴮ)
        C = similar(A, Int, Dᴮ)
        _vfindmax2!(B, C, A, Dᴮ′)
        return B, CartesianIndices(A)[C]
    end
end

@generated _vtfindmax!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                       A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D} = vtfindmax_quote(N, D)
function vtfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        return vtfindmax1(A)
    elseif 1 ∈ dims
        return ivtfindmax(A, dims)
    else
        Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
        Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
        B = similar(A, Dᴮ)
        C = similar(A, Int, Dᴮ)
        _vtfindmax!(B, C, A, Dᴮ′)
        return B, CartesianIndices(A)[C]
    end
end

# For handling cases where the first dimension is being reduced. Alas, still fails
# for some cases, e.g. A ∈ ℝᴵˣᴶˣᴷˣᴸ, dims=(1,2,3)
@generated _ivfindmax!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                       A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D} = ivfindmax_quote(N, D)
function ivfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
    Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
    B = similar(A, Dᴮ)
    C = similar(A, Int, Dᴮ)
    _ivfindmax!(B, C, A, Dᴮ′)
    B, CartesianIndices(A)[C]
end
@generated _ivtfindmax!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                        A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D} = ivtfindmax_quote(N, D)
function ivtfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
    Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
    B = similar(A, Dᴮ)
    C = similar(A, Int, Dᴮ)
    _ivtfindmax!(B, C, A, Dᴮ′)
    B, CartesianIndices(A)[C]
end

function vfindmax1(A::AbstractArray{T, N}) where {T, N}
    m = typemin(T)
    j = 0
    @turbo for i ∈ eachindex(A)
        y = A[i] > m
        m = y ? A[i] : m
        j = y ? i : j
    end
    m, CartesianIndices(A)[j]
end
function vtfindmax1(A::AbstractArray{T, N}) where {T, N}
    m = typemin(T)
    j = 0
    @tturbo for i ∈ eachindex(A)
        y = A[i] > m
        m = y ? A[i] : m
        j = y ? i : j
    end
    m, CartesianIndices(A)[j]
end

################
# Handle scalar dims by wrapping in Tuple
vfindmax(A::AbstractArray{T, N}, dims::Int) where {T, N} = vfindmax(A, (dims,))
# Convenience dispatches to match JuliaBase
vfindmax(A::AbstractArray{T, N}; dims=:) where {T, N} = vfindmax(A, dims)
vfindmax(A::AbstractArray{T, N}) where {T, N} = vfindmax1(A)
vfindmax(A::AbstractArray{T, N}, ::Colon) where {T, N} = vfindmax1(A)
vtfindmax(A::AbstractArray{T, N}, dims::Int) where {T, N} = vtfindmax(A, (dims,))
# Convenience dispatches to match JuliaBase
vtfindmax(A::AbstractArray{T, N}; dims=:) where {T, N} = vtfindmax(A, dims)
vtfindmax(A::AbstractArray{T, N}) where {T, N} = vtfindmax1(A)
vtfindmax(A::AbstractArray{T, N}, ::Colon) where {T, N} = vtfindmax1(A)

# @timev vfindmax1(A)
# @timev vtfindmax1(A)
# @timev vfindmax2(A, (2,));
# @timev vfindmax3(A, (2,));
# @timev vtfindmax2(A, (2,));
# @timev vtfindmax3(A, (2,));


text:
"""
Return average average value for an array of cost vectors.
If the length of the input is less than 1, mean is not well defined,
the key values returned are 'nothing'.
"""
function dict_average(dict_list)
    avcosts = zero_costvector()

    if length(dict_list) == 0
        for cost_type in keys(avcosts)
            avcosts[cost_type] = NaN
        end
        return avcosts
    end

    for cost_type in keys(avcosts)
        costs = collect(map(x->x[cost_type], dict_list))
        avcosts[cost_type] = mean(costs)
    end
    return avcosts
end

"""
Return average standard error for an array of cost vectors.
If the length of the input is less than 2, error is not well defined,
the key values returned are 'nothing'.
"""
function dict_err(dict_list)
    averr = zero_costvector()
    len = length(dict_list)

    if len < 2
        for cost_type in keys(averr)
            averr[cost_type] = NaN
        end
        return averr
    end

    for cost_type in keys(averr)
        costs = collect(map(x->x[cost_type], dict_list))
        averr[cost_type] = std(costs)/(sqrt(length(costs)))
    end
    return averr
end

function percentage_error(main_dict, error_dict)
    perc_dict = Dict()
    for key in keys(main_dict)
        perc_dict[key] = error_dict[key] / main_dict[key]
    end
    return perc_dict
end

# """Generate a list of user_pairs for a QNetwork"""
# function make_user_pairs(net::QNetwork, num_pairs::Int)::Vector{Tuple{Int64, Int64}}
#     num_nodes = length(net.nodes)
#     @assert num_nodes >= num_pairs*2 "Graph space too small for number of pairs"
#     rand_space = Array(collect(1:num_nodes))
#     pairs = Vector{Tuple}()
#     i = 0
#     while i < num_pairs
#         idx = rand(1:length(rand_space))
#         u = rand_space[idx]
#         deleteat!(rand_space, idx)
#         idx = rand(1:length(rand_space))
#         v = rand_space[idx]
#         deleteat!(rand_space, idx)
#         chosen_pair = (u, v)
#         push!(pairs, chosen_pair)
#         i += 1
#     end
#     return pairs
# end

# TODO: This is probably a better function than above
function make_user_pairs(net::QNetwork, num_pairs::Int; node_list=nothing)::Vector{Tuple{Int64, Int64}}
    num_nodes = length(net.nodes)
    @assert num_nodes >= num_pairs*2 "Graph space too small for number of pairs"
    if node_list != nothing
        rand_space = node_list
    else
        rand_space = Array(collect(1:num_nodes))
    end
    pairs = Vector{Tuple}()
    i = 0
    while i < num_pairs
        idx = rand(1:length(rand_space))
        u = rand_space[idx]
        deleteat!(rand_space, idx)
        idx = rand(1:length(rand_space))
        v = rand_space[idx]
        deleteat!(rand_space, idx)
        chosen_pair = (u, v)
        push!(pairs, chosen_pair)
        i += 1
    end
    return pairs
end

"""
Generate a list of user_pairs for a square grid network, ensuring that there is an
edge buffer of a given thickness.
"""
function make_user_pairs_wbuffer(squarenet::QNetwork, num_pairs::Int; buffer::Int=1)::Vector{Tuple{Int64, Int64}}
    num_nodes = length(squarenet.nodes)
    d = sqrt(num_nodes)
    @assert isinteger(d) "Grid network is not square"
    @assert buffer > 0
    @assert d - 2 * buffer > 1 "Gridsize too small for given buffer"
    rand_space = Array(collect(0:num_nodes-1))
    # Filter out nodes in the margin
    # NOTE Cool trick: n % d == x coord, n ÷ d == y coord
    filter!(i -> (buffer <= (i%d)) || (buffer <= (i÷d)) || (i%d < d-buffer) || (i÷d < d-buffer), rand_space)
    # Set first index to 1
    rand_space .+ 1
    pairs = make_user_pairs(squarenet, num_pairs, node_list=rand_space)
    return pairs
end

"""
Generate a list of end-users for a TemporalGraph.
src_layer and dst_layer specify the temporal locations of the source and dst nodes
respectively. The default value for these is -1, which indicates the end-users should
be asynchronus.
"""
function make_user_pairs(net::QuNet.TemporalGraph, num_pairs::Int;
    src_layer::Int64=-1, dst_layer::Int64=-1)::Vector{Tuple{Int64, Int64}}

    num_nodes = net.nv
    @assert num_nodes >= num_pairs*2 "Graph space too small for number of pairs"
    @assert dst_layer <= net.steps "dst_layer must be between 1 and $(net.steps) -- or -1 for async nodes"

    rand_space = Array(collect(1:num_nodes))
    pairs = Vector{Tuple}()
    i = 0
    while i < num_pairs
        # Random source
        idx = rand(1:length(rand_space))
        u = rand_space[idx]
        deleteat!(rand_space, idx)

        # Random dest
        idx = rand(1:length(rand_space))
        v = rand_space[idx]
        deleteat!(rand_space, idx)

        # Update u and v to point to the specified source and dest layers
        # If src_layer == -1, index to async_nodes
        if src_layer == -1
            u += num_nodes * net.steps
        elseif src_layer > 0
            u += (src_layer - 1) * num_nodes
        else
            error("Invalid src_layer. Choose from {1, ..., T.steps} or -1 for asynchronus nodes")
        end

        if dst_layer == -1
            v += num_nodes * net.steps
        elseif dst_layer > 0
            v += (dst_layer - 1) * num_nodes
        else
            error("Invalid dst_layer. Choose from {1, ..., T.steps} or -1 for asynchronus nodes")
        end

        chosen_pair = (u, v)
        push!(pairs, chosen_pair)
        i += 1
    end
    return pairs
end

"""
Given a tally for the number of paths used by each end-user in a greedy_protocol:
(i.e. [3,4,2,1] meaning 3 end-users used no paths, 4, end-users used 1 path, etc.)
This function finds the average number of paths used in the protocol.
"""
function ave_paths_used(pathuse_count::Vector{Float64})
    ave_pathuse = 0.0
    len = length(pathuse_count)
    for i in 1:len
        ave_pathuse += (i-1) * pathuse_count[i]
    end
    ave_pathuse = ave_pathuse / sum(pathuse_count)
    return ave_pathuse
end

"""
Takes a network as input and return greedy_multi_path! performance statistics for some number of
random user pairs. Ensure graph is refreshed before starting.
"""
function net_performance(network::Union{QNetwork, QuNet.TemporalGraph},
    num_trials::Int64, num_pairs::Int64; max_paths=3, src_layer::Int64=-1,
    dst_layer::Int64=-1, edge_perc_rate=0.0, get_apd=false, get_postman=false,
    async_cost=1)

    # Sample of average routing costs between end-users
    ave_cost_data = []
    # Sample of path usage statistics for the algorithm
    pathcount_data = []
    # number of paths found by the last pair in the network
    last_count = []
    # Average distances between userpairs in the network.
    ave_path_distance = []
    ave_path_distance_err = []
    # Average manhattan distances between userpairs post-selected on paths existing
    ave_postman = []
    ave_postman_err = []
    # List of pathlengths > 0 for every trial, and every userpair.
    pathlengths = []
    shortest_avail_path = []

    for i in 1:num_trials
        net = deepcopy(network)

        # Generate random communication pairs
        if typeof(network) == TemporalGraph
            user_pairs = make_user_pairs(network, num_pairs, src_layer=src_layer, dst_layer=dst_layer)
            # Add asynchronus nodes to the network copy
            add_async_nodes!(net, user_pairs, ϵ=async_cost)
        else
            user_pairs = make_user_pairs(network, num_pairs)
        end

        # Percolate edges
        # WARNING: Edge percolation will not be consistant between temporal layers
        # if typeof(net) = QuNet.TemporalGraph
        if edge_perc_rate != 0.0
            @assert 0 <= edge_perc_rate <= 1.0 "edge_perc_rate out of bounds"
            net = QuNet.percolate_edges(net, edge_perc_rate)
            refresh_graph!(net)
        end

        # Get data from greedy_multi_path
        pathset, routing_costs, pathuse_count, lastpair_numpaths = QuNet.greedy_multi_path!(net, purify, user_pairs, max_paths)

        # Add values to ave_path_difference
        if get_apd == true
            add_man_dists!(network, pathset, user_pairs, ave_path_distance)
        end

        # Add values to add_postman_dists!
        if get_postman == true
            add_postman_dists!(network, pathset, user_pairs, ave_postman)
        end

        push!(last_count, lastpair_numpaths)

        # # Collect path lengths from pathset, add them to pathlengths
        # for user in pathset
        #     for path in user
        #         l = length(path)
        #         if l != 0
        #             push!(pathlengths, l)
        #         end
        #     end
        # end

        # Filter out entries where no paths were found and costs are not well defined
        filter!(x->x!=nothing, routing_costs)

        # If the mean is well defined, average the routing costs and push to ave_cost_data
        if length(routing_costs) > 0
            # Average the data
            ave = dict_average(routing_costs)
            push!(ave_cost_data, ave)
        end
        push!(pathcount_data, pathuse_count)
    end

    # Find the mean and standard error of ave_cost_data. Call this the performance
    performance = dict_average(ave_cost_data)
    performance_err = dict_err(ave_cost_data)
    # Get percentage error from standard error:
    perc_err = percentage_error(performance, performance_err)

    # Find the mean and standard error of the path usage statistics:

    # Usage:
        # Each entry in pathcount_data is a vector of Ints of length (max_paths + 1)
        # An example entry is [3, 4, 3, 1]
        # Where 3 end-users found 0 paths, 4 end-users found 1 path, etc.
        # Given a collection of path statistics, ie. [[0, 1, 2, 1], [0, 0, 3, 0], [0, 1, 2, 1]]
        # we want to return vectors of average path usage with associated error:
        # ie. [0, 0.666, 2.333, 0.666] for means
    ave_pathcounts = [0.0 for i in 0:max_paths]
    ave_pathcounts_err = [0.0 for i in 0:max_paths]

    for i in 1:max_paths+1
        data = [pathcount_data[j][i] for j in 1:num_trials]
        ave_pathcounts[i] = mean(data)
        ave_pathcounts_err[i] = std(data)/(sqrt(length(data)))
    end

    # Convert costs from decibelic to metric form
    metric_performance = convert_costs(performance)

    # Get metric error from percentage error:
    metric_err = Dict()
    for key in keys(perc_err)
        metric_err[key] = perc_err[key] * metric_performance[key]
    end

    # Get stats for ave_lastcount
    ave_lastcount = mean(last_count)

    # Get stats for num_zeropaths
    num_zeropaths = count(i->(i==0), last_count)

    # ave_path_distance = mean(pathlengths)
    # ave_path_distance_err = sqrt(var(pathlengths)/num_trials)

    # Get stats for apd
    if get_apd == true
        ave_path_distance_err = sqrt(var(ave_path_distance)/num_trials)
        ave_path_distance = mean(ave_path_distance)
        # sqrt(var(ave_path_distance)/float(num_trials))
    end

    # Get stats for postman
    if get_postman == true
        ave_postman_err = sqrt(var(ave_postman)/num_trials)
        ave_postman = mean(ave_postman)
    end

    # return performance, performance_err, ave_pathcounts, ave_pathcounts_err
    # # TODO:
    return metric_performance, metric_err, ave_pathcounts, ave_pathcounts_err, ave_lastcount, num_zeropaths,
    ave_path_distance, ave_path_distance_err, ave_postman, ave_postman_err
end

# Current version
# Add manhattan distances for given network and userpairs onto man_paths array
# Find L1 distance:
#   For each element in pathset, check if path exists
#   If path exists, find L1 distance
#   Add L1 distance to relevant array
function add_man_dists!(network::QNetwork, pathset, userpairs, man_paths)
    g = network.graph["Z"]
    for (pairidx, paths) in enumerate(pathset)
        if length(paths) != 0
            pair = userpairs[pairidx]
            p = shortest_path(g, pair[1], pair[2])
            l = length(p)
            push!(man_paths, l)
        end
    end
end

# Shortest availible path
# For each element in pathset,
    # Check if path exists
    # If path exists, get length of first path
    # Add path to array
function add_postman_dists!(network, pathset, userpairs, postman_paths)
    for (idx, pair) in enumerate(pathset)
        if length(pair) != 0
            l = length(first(pair))
            push!(postman_paths, l)
        end
    end
end

# # Add manhattan distances for given network and userpairs onto man_paths array
# function add_man_dists!(network::QNetwork, userpairs, man_paths)
#     g = network.graph["Z"]
#     for pair in userpairs
#         path = shortest_path(g, pair[1], pair[2])
#         l = length(path)
#         if l != 0
#             push!(man_paths, l)
#         end
#     end
# end

# Add manhattan distances for given network and userpairs post-selected
# path between userpairs existing in percolated network.
# TODO: Include competition effects.
# function add_postman_dists!(network::QNetwork, perc_network::QNetwork, userpairs, postman_paths)
#     gperc = perc_network.graph["Z"]
#     g = network.graph["Z"]
#     for pair in userpairs
#         path = shortest_path(gperc, pair[1], pair[2])
#         l = length(path)
#         if l != 0
#             # Path exists, therefore find ave manhattan distance
#             path = shortest_path(g, pair[1], pair[2])
#             l = length(path)
#             push!(postman_paths, l)
#         end
#     end
# end

# function add_postman_dists!(network, pathset, userpairs, postman_paths)
#     kill_list = []
#     filtered_pairs = deepcopy(userpairs)
#     for (idx, pair) in enumerate(pathset)
#         if length(pair) != 0
#             push!(kill_list, idx)
#         end
#     end
#     deleteat!(filtered_pairs, kill_list)
#     add_man_dists!(network, filtered_pairs, postman_paths)
# end

"""
Takes a network as input and returns greedy_multi_path! heatmap data for
some number of user pairs. (i.e. a list of efficiency fidelity coordinates
for all end users over all num_trials)
"""
function heat_data(network::Union{QNetwork, QuNet.TemporalGraph},
    num_trials::Int64, num_pairs::Int64; max_paths=3, src_layer::Int64=-1,
    dst_layer::Int64=-1, edge_perc_rate=0.0)

    # e,f coords for all end users over all num_trials
    coord_list = []

    for i in 1:num_trials
        println("Collecting data for trial $i")
        net = deepcopy(network)

        # Generate random communication pairs
        if typeof(network) == TemporalGraph
            user_pairs = make_user_pairs(network, num_pairs, src_layer=src_layer, dst_layer=dst_layer)
            # Add asynchronus nodes to the network copy
            add_async_nodes!(net, user_pairs)
        else
            user_pairs = make_user_pairs(network, num_pairs)
        end

        # Percolate edges
        # WARNING: Edge percolation will not be consistant between temporal layers
        # if typeof(net) = QuNet.TemporalGraph
        if edge_perc_rate != 0.0
            @assert 0 <= edge_perc_rate <= 1.0 "edge_perc_rate out of bounds"
            net = QuNet.percolate_edges(net, edge_perc_rate)
            refresh_graph!(net)
        end

        # Get data from greedy_multi_path
        dummy, routing_costs, pathuse_count = QuNet.greedy_multi_path!(net, purify, user_pairs, max_paths)

        # Filter out entries where no paths were found and costs are not well defined
        filter!(x->x!=nothing, routing_costs)

        # Put end-user costs into tuples (e, f)
        for usercost in routing_costs
            coord = Vector{Float64}()
            push!(coord, usercost["loss"])
            push!(coord, usercost["Z"])
            # Add this new coordinate to the list
            push!(coord_list, coord)
        end
    end
    return coord_list
end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
@testset "Reporting" begin

    p = DSProblem(3, objective = x -> sum(x.^2), initial_point= [10, 10, 10])
    Optimize!(p)

    @testset "report_config" begin
        rc = DS.report_config(p)

        @test rc.title == "Config"

        names = [e.first for e in rc.entries if e != nothing]

        @test "Search" in names
        @test "Poll" in names
        @test "Mesh" in names
        @test "Mesh Scale" in names
        @test "Opportunistic" in names
        @test "Number of processes" in names
        @test "Max simultanious evaluations" in names

        #just check that it doesn't crash
        @test typeof(DS.format(rc)) == String
    end

    @testset "report_config" begin
        rc = DS.report_status(p)

        @test rc.title == "Status"

        names = [e.first for e in rc.entries if e != nothing]

        @test "Function Evaluations" in names
        @test "Iterations" in names
        @test "Optimization Status" in names
        @test "Runtime" in names
        @test "Search Time" in names
        @test "Poll Time" in names
        @test "Blackbox Evaluation Time" in names

        @test typeof(DS.format(rc)) == String
    end

    @testset "report_problem" begin
        rc = DS.report_problem(p)

        @test rc.title == "Optimization Problem"

        names = [e.first for e in rc.entries if e != nothing]

        @test "Variables" in names
        @test "Initial Point" in names
        @test "Sense" in names
        @test "Feasible Solution" in names
        @test "Feasible Cost" in names
        @test "Infeasible Solution" in names
        @test "Infeasible Cost" in names

        @test typeof(DS.format(rc)) == String
    end
end


text:
info(
"======================================================================
                              Running Behavior Tests
      ======================================================================")

@testset "OpenCL Hello World Test" begin

    hello_world_kernel = "
        #pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable

        __constant char hw[] = \"hello world\";

        __kernel void hello(__global char *out) {
            int tid = get_global_id(0);
            out[tid] = hw[tid];
        }"

    hello_world_str = "hello world"

    for device in cl.devices()
        if device[:platform][:name] == "Portable Computing Language"
            warn("Skipping OpenCL.Kernel mem/workgroup size for Portable Computing Language Platform")
            continue
        end

        ctx   = cl.Context(device)
        queue = cl.CmdQueue(ctx)

        str_len  = length(hello_world_str) + 1
        out_buf  = cl.Buffer(Cchar, ctx, :w, sizeof(Cchar) * str_len)

        prg   = cl.Program(ctx, source=hello_world_kernel) |> cl.build!
        kern  = cl.Kernel(prg, "hello")

        queue(kern, str_len, nothing, out_buf)
        h = cl.read(queue, out_buf)

        @test cl.CLString(h) == hello_world_str
    end
end

@testset "OpenCL Low Level Api Test" begin

  test_source = "
    __kernel void sum(__global const float *a,
                      __global const float *b,
                      __global float *c,
                      const unsigned int count)
    {
      unsigned int gid = get_global_id(0);
      if (gid < count) {
          c[gid] = a[gid] + b[gid];
      }
    }
    "

    for device in cl.devices()

        length = 1024
        h_a = Array(cl.CL_float, length)
        h_b = Array(cl.CL_float, length)
        h_c = Array(cl.CL_float, length)
        h_d = Array(cl.CL_float, length)
        h_e = Array(cl.CL_float, length)
        h_f = Array(cl.CL_float, length)
        h_g = Array(cl.CL_float, length)

        for i in 1:length
            h_a[i] = cl.cl_float(rand())
            h_b[i] = cl.cl_float(rand())
            h_e[i] = cl.cl_float(rand())
            h_g[i] = cl.cl_float(rand())
        end

        err_code = Ref{cl.CL_int}()

        # create compute context (TODO: fails if function ptr's not passed...)
        ctx_id = cl.api.clCreateContext(C_NULL, 1, [device.id],
                                        C_NULL,
                                        C_NULL,
                                        err_code)
        if err_code[] != cl.CL_SUCCESS
            throw(cl.CLError(err_code[]))
        end

        q_id = cl.api.clCreateCommandQueue(ctx_id, device.id, 0, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Failed to create command queue")
        end

        # create program
        bytesource = String(test_source)
        prg_id = cl.api.clCreateProgramWithSource(ctx_id, 1, [bytesource], C_NULL, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Failed to create program")
        end

        # build program
        err = cl.api.clBuildProgram(prg_id, 0, C_NULL, C_NULL, C_NULL, C_NULL)
        if err != cl.CL_SUCCESS
            error("Failed to build program")
        end

        # create compute kernel
        k_id = cl.api.clCreateKernel(prg_id, "sum", err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Failed to create compute kernel")
        end

        # create input array in device memory
        Aid = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_READ_ONLY | cl.CL_MEM_COPY_HOST_PTR,
                                    sizeof(cl.CL_float) * length, h_a, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer A")
        end
        Bid = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_READ_ONLY | cl.CL_MEM_COPY_HOST_PTR,
                                    sizeof(cl.CL_float) * length, h_b, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer B")
        end
        Eid = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_WRITE_ONLY | cl.CL_MEM_COPY_HOST_PTR,
                                    sizeof(cl.CL_float) * length, h_e, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer E")
        end
        Gid = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_WRITE_ONLY | cl.CL_MEM_COPY_HOST_PTR,
                                    sizeof(cl.CL_float) * length, h_g, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer G")
        end

        # create output arrays in device memory

        Cid = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_READ_WRITE,
                                    sizeof(cl.CL_float) * length, C_NULL, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer C")
        end
        Did = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_READ_WRITE,
                                    sizeof(cl.CL_float) * length, C_NULL, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer D")
        end
        Fid = cl.api.clCreateBuffer(ctx_id, cl.CL_MEM_WRITE_ONLY,
                                    sizeof(cl.CL_float) * length, C_NULL, err_code)
        if err_code[] != cl.CL_SUCCESS
            error("Error creating buffer F")
        end

        err  = cl.api.clSetKernelArg(k_id, 0, sizeof(cl.CL_mem), [Aid])
        err |= cl.api.clSetKernelArg(k_id, 1, sizeof(cl.CL_mem), [Bid])
        err |= cl.api.clSetKernelArg(k_id, 2, sizeof(cl.CL_mem), [Cid])
        err |= cl.api.clSetKernelArg(k_id, 3, sizeof(cl.CL_uint), cl.CL_uint[length])
        if err != cl.CL_SUCCESS
            error("Error setting kernel 1 args")
        end

        nglobal = Ref{Csize_t}(length)
        err = cl.api.clEnqueueNDRangeKernel(q_id, k_id,  1, C_NULL,
                                            nglobal, C_NULL, 0, C_NULL, C_NULL)
        if err != cl.CL_SUCCESS
            error("Failed to execute kernel 1")
        end

        err  = cl.api.clSetKernelArg(k_id, 0, sizeof(cl.CL_mem), [Eid])
        err |= cl.api.clSetKernelArg(k_id, 1, sizeof(cl.CL_mem), [Cid])
        err |= cl.api.clSetKernelArg(k_id, 2, sizeof(cl.CL_mem), [Did])
        if err != cl.CL_SUCCESS
            error("Error setting kernel 2 args")
        end
        err = cl.api.clEnqueueNDRangeKernel(q_id, k_id,  1, C_NULL,
                                            nglobal, C_NULL, 0, C_NULL, C_NULL)
        if err != cl.CL_SUCCESS
            error("Failed to execute kernel 2")
        end

        err  = cl.api.clSetKernelArg(k_id, 0, sizeof(cl.CL_mem), [Gid])
        err |= cl.api.clSetKernelArg(k_id, 1, sizeof(cl.CL_mem), [Did])
        err |= cl.api.clSetKernelArg(k_id, 2, sizeof(cl.CL_mem), [Fid])
        if err != cl.CL_SUCCESS
            error("Error setting kernel 3 args")
        end
        err = cl.api.clEnqueueNDRangeKernel(q_id, k_id,  1, C_NULL,
                                            nglobal, C_NULL, 0, C_NULL, C_NULL)
        if err != cl.CL_SUCCESS
            error("Failed to execute kernel 3")
        end

        # read back the result from compute device...
        err = cl.api.clEnqueueReadBuffer(q_id, Fid, cl.CL_TRUE, 0,
                                         sizeof(cl.CL_float) * length, h_f, 0, C_NULL, C_NULL)
        if err != cl.CL_SUCCESS
            error("Failed to read output array")
        end

        # test results
        for i in 1:length
            tmp = h_a[i] + h_b[i] + h_e[i] + h_g[i]
            @test tmp ≈ h_f[i]
        end
    end
end

immutable Params
    A::Float32
    B::Float32
    #TODO: fixed size arrays?
    X1::Float32
    X2::Float32
    C::Int32
    Params(a, b, x, c) = begin
        new(Float32(a),
            Float32(b),
            Float32(x[1]),
            Float32(x[2]),
            Int32(c))
    end
end

let test_struct = "
    typedef struct Params
    {
        float A;
        float B;
        float x[2];  //padding
        int C;
    } Params;

    __kernel void part3(__global const float *a,
                        __global const float *b,
                        __global float *c,
                        __constant struct Params* test)
    {
        int gid = get_global_id(0);
        c[gid] = test->A * a[gid] + test->B * b[gid] + test->C;
    }
"

@testset "OpenCL Struct Buffer Test" begin
    for device in cl.devices()

        if device[:platform][:name] == "Portable Computing Language"
            warn("Skipping OpenCL Struct Buffer Test for Portable Computing Language Platform")
            continue
        end

        ctx = cl.Context(device)
        q   = cl.CmdQueue(ctx)
        p   = cl.Program(ctx, source=test_struct) |> cl.build!

        part3 = cl.Kernel(p, "part3")

        X = fill(1f0, 10)
        Y = fill(1f0, 10)

        P = [Params(0.5, 10.0, [0.0, 0.0], 3)]

        #TODO: constructor for single immutable types.., check if passed parameter isbits
        P_buf = cl.Buffer(Params, ctx, :r, length(P))
        cl.write!(q, P_buf, P)

        X_buf = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=X)
        Y_buf = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=Y)
        R_buf = cl.Buffer(Float32, ctx, :w, length(X))

        global_size = size(X)
        q(part3, global_size, nothing, X_buf, Y_buf, R_buf, P_buf)

        r = cl.read(q, R_buf)
        @test all(x -> x == 13.5, r)
    end
end
end