text:
### Stochastic relaxations for LSQ

function train_sr(
  X::Matrix,                  # d-by-n matrix of data points to train on.
  m::Integer,                 # number of codebooks
  h::Integer,                 # number of entries per codebook
  R::Matrix{Float32},         # init rotation
  B::Matrix{Int16},     # init codes
  C::Vector{Matrix{Float32}}, # init codebooks
  niter::Integer,             # number of optimization iterations
  ilsiter::Integer,           # number of ILS iterations to use during encoding
  icmiter::Integer,           # number of iterations in local search
  randord::Bool,              # whether to use random order
  npert::Integer,             # The number of codes to perturb
  method::AbstractString,     # The SR method to use. Either SR_C or SR_D
  p::Float32,                 # SR-D power parameter
  cpp::Bool=true,             # whether to use ICM's cpp implementation
  V::Bool=false)              # whether to print progress

  if V
  println("**********************************************************************************************")
  println("Doing local search with $m codebooks, $npert perturbations, $icmiter icm iterations and random order = $randord");
  println("**********************************************************************************************")
  end

  if !(method in ["SR_C", "SR_D"]); error("SR method unknown"); end
  d, n = size( X )

  RX = R' * X
  if V; @printf("Random error: %e\n", qerror( RX, B, C )); end

  if method == "SR_C"
    # In SR-C we add noise to X
    RX_noisy = SR_C_perturb( RX, 0, niter, p )
    C = update_codebooks_fast_bin( RX_noisy, B, h, V )
  else
    # In SR-D we add noise to C
    C = update_codebooks_fast_bin( RX, B, h, V )
    obj = qerror( RX, B, C )
    if V; @printf("%3d %e \n", -1, obj); end
    C = SR_D_perturb( C, 1, niter, p )
  end

  obj = qerror( RX, B, C );
  if V; @printf("%3d %e \n", -1, obj); end

  # Initialize B
  @time B = encoding_icm( RX, B, C, ilsiter, icmiter, randord, npert, cpp, V )
  obj = qerror( RX, B, C )
  @printf("%3d %e \n", -1, obj)

  obj     = Inf
  objlast = Inf
  objarray = zeros( Float32, niter+1 )

  for iter = 1:niter

    objlast = obj
    obj = qerror( RX, B, C  )
    objarray[iter] = obj
    if V; @printf("%3d %e (%e better) \n", iter, obj, objlast - obj); end

    if method == "SR_C"
      # In SR-C we add noise to X
      RX_noisy = SR_C_perturb( RX, iter, niter, p )
      C = update_codebooks_fast_bin( RX_noisy, B, h, V )
    else
      C = update_codebooks_fast_bin( RX, B, h, V )
      C = SR_D_perturb( C, iter, niter, p )
    end

    # Update the codes with local search
    @time B = encoding_icm( RX, B, C, ilsiter, icmiter, randord, npert, cpp, V )

  end

  objarray[niter+1] = qerror( RX, B, C )

  # Apply the rotation to the codebooks
  for i = 1:m; C[i] = R * C[i]; end

  return C, B, objarray
end

# TODO this should be doable with metaprogramming

function train_sr_cuda(
  X::Matrix,                  # d-by-n matrix of data points to train on.
  m::Integer,                 # number of codebooks
  h::Integer,                 # number of entries per codebook
  R::Matrix{Float32},         # init rotation
  B::Matrix{Int16},     # init codes
  C::Vector{Matrix{Float32}}, # init codebooks
  niter::Integer,             # number of optimization iterations
  ilsiter::Integer,           # number of ILS iterations to use during encoding
  icmiter::Integer,           # number of iterations in local search
  randord::Bool,              # whether to use random order
  npert::Integer,             # The number of codes to perturb
  method::AbstractString,     # The SR method to use. Either SR_C or SR_D
  schedule::Integer,          # 1-3 Schedule to use
  p::AbstractFloat=0.5,       # Temperature decay parameter
  nsplits::Integer=1,         # The number of splits for icm encoding (for limited memory GPUs)
  V::Bool=false)              # whether to print progress

  if V
  println("**********************************************************************************************");
  println("Doing local search with $m codebooks, $npert perturbations, $icmiter icm iterations and random order = $randord");
  println("**********************************************************************************************");
  end

  if !(method in ["SR_C", "SR_D"]); error("SR method unknown"); end
  d, n = size( X )

  RX = R' * X
  if V; @printf("Random error: %e\n", qerror(RX, B, C)); end

  if method == "SR_C"
    # In SR-C we add noise to X
    RX_noisy = SR_C_perturb(RX, 0, niter, schedule, p)
    C = update_codebooks(RX_noisy, B, h, V, "fastbin")
  else
    # In SR-D we add noise to C
    C = update_codebooks(RX, B, h, V, "fastbin")
    obj = qerror(RX, B, C)
    if V; @printf("%3d %e \n", -1, obj); end
    C = SR_D_perturb(C, 1, niter, schedule, p)
  end

  obj = qerror(RX, B, C)
  if V; @printf("%3d %e \n", -1, obj); end

  # Initialize B
  @time B, _ = encode_icm_cuda(RX, B, C, [ilsiter], icmiter, npert, randord, nsplits, V)
  # gc()
  B = B[end]
  obj = qerror(RX, B, C)
  @printf("%3d %e \n", -1, obj)

  obj     = Inf
  objlast = Inf
  objarray = zeros(Float32, niter+1)

  for iter = 1:niter

    objlast = obj
    obj = qerror(RX, B, C)
    objarray[iter] = obj
    if V; @printf("%3d %e (%e better) \n", iter, obj, objlast - obj); end

    if method == "SR_C"
      # In SR-C we add noise to X before updating C
      RX_noisy = SR_C_perturb(RX, iter, niter, schedule, p)
      C = update_codebooks(RX_noisy, B, h, V, "fastbin")
    else
      # In SR-D we add noise to C before updating B
      C = update_codebooks(RX, B, h, V, "fastbin")
      C = SR_D_perturb(C, iter, niter, schedule, p)
    end

    # Update the codes with local search
    @time B, _ = encode_icm_cuda(RX, B, C, [ilsiter], icmiter, npert, randord, nsplits, V)
    # gc()
    B = B[end]

    C = update_codebooks(RX, B, h, V, "fastbin")
  end

  objarray[niter+1] = qerror(RX, B, C)

  # Apply the rotation to the codebooks
  for i = 1:m; C[i] = R * C[i]; end

  return C, B, objarray
end

function experiment_sr(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  R::Matrix{T}, # rotation
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer=8, # number of ILS iterations
  icmiter::Integer=4, # number of ICM iterations
  randord::Bool=true, # whether to explore the codes in random order
  npert::Integer=4, # number of codes to perturb in each ILS iterations
  knn::Integer=1000,
  sr_method::String="SR_D",
  schedule::Integer=1, # schedule to use
  p::AbstractFloat=0.5, # Temperature decay parameter
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  # TODO maybe epose this? cpp is always faster
  cpp     = true

  # Train LSQ
  d, _ = size(Xt)
  C, B, obj = train_sr(Xt, m, h, R, B, C, niter, ilsiter, icmiter, randord, npert, sr_method, schbedule, p, cpp, V)
  norms_B, norms_C = get_norms_codebook(B, C)

  # === Encode the base set ===
  B_base = convert(Matrix{Int16}, rand(1:h, m, size(Xb,2)))
  B_base = encoding_icm(Xb, B_base, C, ilsiter, icmiter, randord, npert, cpp, V)
  base_error = qerror(Xb, B_base, C)
  if V; @printf("Error in base is %e\n", base_error); end

  # Compute and quantize the database norms
  B_base_norms, db_norms_X = quantize_norms(B_base, C, norms_C)
  db_norms     = vec( norms_C[ B_base_norms ] )

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms, Matrix{Float32}(1.0*I, d, d), knn)
  if V; println("done"); end

  rec = eval_recall(gt, idx, knn)
end

"Runs an lsq experiment/demo"
function experiment_sr(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # OPQ initialization
  C, B, R, _ = train_opq(Xt, m, h, niter, "natural", V)

  # ChainQ (second initialization)
  C, B, R, train_error = train_chainq(Xt, m, h, R, B, C, niter, V)

  # Actual experiment
  experiment_sr(Xt, B, C, R, Xb, Xq, gt, m, h, niter, knn, V)
end

function experiment_sr_cuda(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  R::Matrix{T}, # rotation
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer=8, # number of ILS iterations
  icmiter::Integer=4, # number of ICM iterations
  randord::Bool=true, # whether to explore the codes in random order
  npert::Integer=4, # number of codes to perturb in each ILS iterations
  knn::Integer=1000,
  nsplits_train::Integer=1,
  nsplits_base::Integer=1, # Number of splits for training data so GPU does not run out of memory
  sr_method::String="SR_D",
  schedule::Integer=1, # schedule to use
  p::AbstractFloat=0.5, # Temperature decay parameter
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  if V
    println()
    println("**********************************************************************************************");
    println("Running LSQ++ ($(sr_method)) with $m codebooks, $npert perturbations, $icmiter icm iterations and random order = $randord");
    println("**********************************************************************************************");
  end

  # Train LSQ
  d, _ = size(Xt)
  C, B, train_error = train_sr_cuda(Xt, m, h, R, B, C, niter, ilsiter, icmiter, randord, npert, sr_method, schedule, p, nsplits_train, V)
  norms_B, norms_C = get_norms_codebook(B, C)

  # === Encode the base set ===
  B_base = convert(Matrix{Int16}, rand(1:h, m, size(Xb,2)))

  # ilsiters = [16, 32, 64, 128, 256]
  ilsiters = [ilsiter * 4]
  Bs_base, _ = encode_icm_cuda(Xb, B_base, C, ilsiters, icmiter, npert, randord, nsplits_base, V)

  # for (idx, ilsiter) in enumerate(ilsiters)
  B_base = Bs_base[end]
  base_error = qerror(Xb, B_base, C)
  if V; @printf("Error in base is %e\n", base_error); end

  # Compute and quantize the database norms
  B_base_norms, db_norms_X = quantize_norms(B_base, C, norms_C)
  db_norms = vec(norms_C[ B_base_norms ])

  if V; print("Querying m=$m ... "); end
  # @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms_X, eye(Float32, d), knn)
  @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms, Matrix{Float32}(1.0*I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  # end
  return C, B, R, train_error, B_base, recall
end

function experiment_sr_cuda_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  R::Matrix{T}, # rotation
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer=8, # number of ILS iterations
  icmiter::Integer=4, # number of ICM iterations
  randord::Bool=true, # whether to explore the codes in random order
  npert::Integer=4, # number of codes to perturb in each ILS iterations
  knn::Integer=1000,
  nsplits_train::Integer=1,
  sr_method::String="SR_D",
  schedule::Integer=1, # schedule to use
  p::AbstractFloat=0.5, # Temperature decay parameter
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  # Train LSQ
  d, _ = size(Xt)
  if V; print("Running CUDA $(sr_method) training... "); end
  C, B, train_error = train_sr_cuda(Xt, m, h, R, B, C, niter, ilsiter, icmiter, randord, npert, sr_method, schedule, p, nsplits_train, V)
  if V; @printf("done\n"); end

  norms_B, norms_C = get_norms_codebook(B, C)
  db_norms         = vec(norms_C[norms_B])

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B, Xq, C, db_norms, Matrix{Float32}(1.0*I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, R, train_error, recall
end

"Runs an lsq experiment/demo"
function experiment_sr_cuda(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  nsplits_train::Integer=1,
  nsplits_base::Integer=1,
  sr_method::String="SR_D",
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # OPQ initialization
  if V; @printf("Running OPQ initialization... "); end
  niter_opq = 25
  C, B, R, _ = train_opq(Xt, m, h, niter_opq, "natural", V)
  if V; @printf("done\n"); end

  # ChainQ (second initialization)
  # C, B, R, train_error = train_chainq(Xt, m, h, R, B, C, niter, V)

  # Actual experiment
  experiment_sr_cuda(Xt, B, C, R, Xb, Xq, gt, m, h, niter, knn, nsplits_train, nsplits_base, sr_method, V)
end

"Runs an lsq experiment/demo"
function experiment_sr_cuda_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  nsplits_train::Integer=1,
  sr_method::String="SR_D",
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # OPQ initialization
  if V; @printf("Running OPQ initialization... "); end
  niter_opq = 25
  C, B, R, opq_error = train_opq(Xt, m, h, niter_opq, "natural", V)
  if V; @printf("done\n"); end

  # ChainQ (second initialization)
  if V; @printf("Running ChainQ initialization... "); end
  niter_chainq = 25
  C, B, R, chainq_error = train_chainq(Xt, m, h, R, B, C, niter_chainq, V)
  if V; @printf("done\n"); end

  # Actual experiment
  experiment_sr_cuda_query_base(Xt, B, C, R, Xq, gt, m, h, niter, knn, nsplits_train, sr_method, V), opq_error
end


text:
module FDFD

using GeometryPrimitives, StaticArrays, AxisArrays, IntervalSets
using Logging, SparseArrays, LinearAlgebra

include("./types.jl");
include("./grid.jl");
include("./data.jl");
include("./device.jl");
include("./pml.jl");
include("./flux.jl");
include("./solver/solver.jl");
include("./solver/driven.jl");
include("./solver/eigen.jl");
include("./solver/modulation.jl");
include("./solver/nonlinear.jl");

export unwrap

iscallable(f) = !isempty(methods(f))

function unwrap(v, inplace=false)
  # currently assuming an array
  unwrapped = inplace ? v : copy(v)
  for i in 2:length(v)
    while unwrapped[i] - unwrapped[i-1] >= pi
      unwrapped[i] -= 2pi
    end
    while unwrapped[i] - unwrapped[i-1] <= -pi
      unwrapped[i] += 2pi
    end
  end
  return unwrapped
end

end

ext: jl


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
SUITE["Spacecraft"] = BenchmarkGroup()

include("SpacecraftRendezvous.jl")

SRNA01, options = spacecraft(abort_time=-1.)
SRA01, _ = spacecraft(abort_time=120.)
SRA02, _ = spacecraft(abort_time=[120., 125.])
SRA03, _ = spacecraft(abort_time=[120., 145.])
SRA04, _ = spacecraft(abort_time=240.)
SRA05, _ = spacecraft(abort_time=[235., 240.])
SRA06, _ = spacecraft(abort_time=[230., 240.])
SRA07, _ = spacecraft(abort_time=[50., 150.])
SRA08, _ = spacecraft(abort_time=[0., 240.])
SRU01, _ = spacecraft(abort_time=260.)
SRU02, _ = spacecraft(abort_time=[0., 260.])

# ==============================================================================
# Decomposition-based approach
# ==============================================================================

# algorithm-specific options
options[:mode] = "check"

𝑂_dense = Options(:partition => [1:5], :δ => 0.04)
𝑂_discrete = merge(𝑂_dense, Options(:discretization => "nobloating", :δ => 0.1))
opC_dense = BFFPSV18(𝑂_dense)
opC_discrete = BFFPSV18(𝑂_discrete)
opD = LazyDiscretePost(:lazy_R⋂I => true, :lazy_R⋂G => true)

# single run to verify that specification holds
res = solve(SRNA01, options, opC_dense, opD)
@assert res.satisfied
res = solve(SRNA01, options, opC_discrete, opD)
@assert res.satisfied
res = solve(SRA01, options, opC_dense, opD)
@assert res.satisfied
res = solve(SRA01, options, opC_discrete, opD)
@assert res.satisfied
# res = solve(SRA02, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA02, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA03, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA03, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA04, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA04, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA05, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA05, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA06, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA06, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA07, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA07, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA08, options, opC_dense, opD)
# @assert res.satisfied  # TODO NEW SETTING
# res = solve(SRA08, options, opC_discrete, opD)
# @assert res.satisfied  # TODO NEW SETTING
res = solve(SRU01, options, opC_dense, opD)
@assert !res.satisfied
res = solve(SRU01, options, opC_discrete, opD)
@assert !res.satisfied
res = solve(SRU02, options, opC_dense, opD)
@assert !res.satisfied
res = solve(SRU02, options, opC_discrete, opD)
@assert !res.satisfied

# benchmark
SUITE["Spacecraft"]["SRNA01-SR02", "dense"] =
    @benchmarkable solve($SRNA01, $options, $opC_dense, $opD)
SUITE["Spacecraft"]["SRNA01-SR02", "discrete"] =
    @benchmarkable solve($SRNA01, $options, $opC_discrete, $opD)
SUITE["Spacecraft"]["SRA01-SR02", "dense"] =
    @benchmarkable solve($SRA01, $options, $opC_dense, $opD)
SUITE["Spacecraft"]["SRA01-SR02", "discrete"] =
    @benchmarkable solve($SRA01, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA02-SR02", "dense"] =
#     @benchmarkable solve($SRA02, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA02-SR02", "discrete"] =
#     @benchmarkable solve($SRA02, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA03-SR02", "dense"] =
#     @benchmarkable solve($SRA03, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA03-SR02", "discrete"] =
#     @benchmarkable solve($SRA03, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA04-SR02", "dense"] =
#     @benchmarkable solve($SRA04, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA04-SR02", "discrete"] =
#     @benchmarkable solve($SRA04, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA05-SR02", "dense"] =
#     @benchmarkable solve($SRA05, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA05-SR02", "discrete"] =
#     @benchmarkable solve($SRA05, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA06-SR02", "dense"] =
#     @benchmarkable solve($SRA06, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA06-SR02", "discrete"] =
#     @benchmarkable solve($SRA06, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA07-SR02", "dense"] =
#     @benchmarkable solve($SRA07, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA07-SR02", "discrete"] =
#     @benchmarkable solve($SRA07, $options, $opC_discrete, $opD)
# SUITE["Spacecraft"]["SRA08-SR02", "dense"] =
#     @benchmarkable solve($SRA08, $options, $opC_dense, $opD)
# SUITE["Spacecraft"]["SRA08-SR02", "discrete"] =
#     @benchmarkable solve($SRA08, $options, $opC_discrete, $opD)
SUITE["Spacecraft"]["SRU01-SR02", "dense"] =
    @benchmarkable solve($SRU01, $options, $opC_dense, $opD)
SUITE["Spacecraft"]["SRU01-SR02", "discrete"] =
    @benchmarkable solve($SRU01, $options, $opC_discrete, $opD)
SUITE["Spacecraft"]["SRU02-SR02", "dense"] =
    @benchmarkable solve($SRU02, $options, $opC_dense, $opD)
SUITE["Spacecraft"]["SRU02-SR02", "discrete"] =
    @benchmarkable solve($SRU02, $options, $opC_discrete, $opD)

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=false)

# return the sample with the smallest time value in each test
println("minimum: ", minimum(results))

# return the median for each test
println("median: ", median(results))

# ==============================================================================
# Create plots
# ==============================================================================

options[:mode] = "reach"
options[:plot_vars] = [1, 2]
options[:project_reachset] = true
res = solve(SRNA01, options, opC_dense, opD)
plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"s_x\raisebox{-1mm}{\textcolor{white}{.}}",
     ylab=L"s_y\raisebox{-2mm}{\textcolor{white}{\rule{1mm}{4mm}}}",
     xtick=[-1000., -500., 0.], ytick=[-400., -300., -200., -100., 0., 100.],
     xlims=(-1000., 0.), ylims=(-500., 100.),
     bottom_margin=10mm, left_margin=10mm, top_margin=3mm,
     size=(1000, 1000))
savefig(@relpath "SRNA01_SR02.png")

res = solve(SRA01, options, opC_dense, opD)
plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"s_x\raisebox{-1mm}{\textcolor{white}{.}}",
     ylab=L"s_y\raisebox{-2mm}{\textcolor{white}{\rule{1mm}{4mm}}}",
     xtick=[-1000., -500., 0.], ytick=[-400., -300., -200., -100., 0., 100.],
     xlims=(-1000., 400.), ylims=(-500., 100.),
     bottom_margin=10mm, left_margin=10mm, top_margin=3mm,
     size=(1000, 1000))
savefig(@relpath "SRA01_SR02.png")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
import Mercury as Hg
import ZMQ
using Sockets
using Test

if !isdefined(@__MODULE__, :TestMsg)
    include("jlout/test_msg_pb.jl")
end

@testset "Serial Relay Tests" begin
    serial_port_name = "/dev/tty.usbmodem14201"
    baudrate = 57600
    ctx = ZMQ.Context()
    addr = ip"127.0.0.1"
    port_in = "5555"
    port_out = "5556"

    serial_relay = Hg.launch_relay(
        serial_port_name,
        baudrate,
        Hg.tcpstring(addr, port_in),
        Hg.tcpstring(addr, port_out),
    )

    do_publish = Threads.Atomic{Bool}(true)
    function pub_message(pub, msg)
        do_publish
        while (do_publish[])
            Hg.publish(pub, msg)
            sleep(0.5)
        end
    end

    pub = Hg.ZmqPublisher(ctx, addr, port_in, name = "TestPub")
    # Subscriber to subscribe directly to pub
    sub1 = Hg.ZmqSubscriber(ctx, addr, port_in, name = "TestSub")
    # Subscriber to subscribe directly to the serial_relay pub.
    # These should produce identical results
    sub2 = Hg.ZmqSubscriber(ctx, addr, port_out, name = "TestSub")

    # Test the serial relay with pure byte arrays
    do_publish[] = true
    bytes_out = Vector{UInt8}("Hello World!")
    pub_task = @async pub_message(pub, bytes_out)
    sleep(2.0)

    # Input arrays
    bytes_in1 = zeros(UInt8, sizeof(bytes_out))
    bytes_in2 = zeros(UInt8, sizeof(bytes_out))

    recieved1 = false
    recieved2 = false
    for i = 1:100
        if !recieved1
            recieved1 = Hg.receive(sub1, bytes_in1)
        elseif !recieved2
            recieved2 = Hg.receive(sub2, bytes_in2)
        else
            break
        end
        sleep(0.1)
    end
    @test recieved1
    @test recieved2
    @test all(bytes_out .== bytes_in1 .== bytes_in2)
    do_publish[] = false
    wait(pub_task)

    # Test the code with Protobuf messages
    msg_out = TestMsg(x = 10, y = 11, z = 12)
    do_publish[] = true
    pub_task = @async pub_message(pub, msg_out)
    sleep(2.0)

    msg_in1 = TestMsg(x = 1, y = 1, z = 1)
    msg_in2 = TestMsg(x = 1, y = 1, z = 1)

    recieved1 = false
    recieved2 = false
    for i = 1:100
        if !recieved1
            recieved1 = Hg.receive(sub1, msg_in1)
        elseif !recieved2
            recieved2 = Hg.receive(sub2, msg_in2)
        else
            break
        end
        sleep(0.1)
    end
    @test recieved1
    @test recieved2
    @test msg_out.x == msg_in1.x == msg_in2.x
    @test msg_out.y == msg_in1.y == msg_in2.y
    @test msg_out.z == msg_in1.z == msg_in2.z

    do_publish[] = false
    wait(pub_task)
    close(sub1)
    close(sub2)
    close(pub)
    close(serial_relay)
end


text:
include("HamiltonODE/HamiltonSystem.jl")
include("HamiltonODE/HamiltonSteps.jl")
include("HamiltonODE/HamiltonIntegrator.jl")
include("HamiltonODE/HamiltonPlot.jl")
include("HamiltonODE/HamiltonInterpolate.jl")

println("============================")
println("Finished loading HamiltonODE")

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
module Inputs

export AbstractInput

"""
    AbstractInput{T<:Real, N, L}

An abstract type for iterators producing SNN input data, where
`T` is the type of input produced, `N` is the number of neurons
to make input for, and `L` is the length of the input sequence.
As such, an AbstractInput produces a length-`N` `Vector{T}` each
iteration.
"""
abstract type AbstractInput{T<:Real, N, L} end

function Base.length(iter::AbstractInput{T, N, L}) where {T, N, L} return L end
function Base.size(iter::AbstractInput{T, N, L}) where {T, N, L} return (N, L) end
function Base.eltype(t::Type{AbstractInput{T, N, L}}) where {T, N, L} return T end

include("PoissonInputs.jl")
include("ConstantPoissonInputs.jl")
include("VariablePoissonInputs.jl")
include("SummedInputs.jl")
end


text:
using MarketData

facts("collapse operations") do

  context("collapse squishes correctly") do
    @fact collapse(cl, last).values[1]                  => 99.50
    @fact collapse(cl, last).timestamp[1]               => Date(2000,1,7)
    @fact collapse(cl, last, period=month).values[1]    => 103.75
    @fact collapse(cl, last, period=month).timestamp[1] => Date(2000,1,31)
  end
end

facts("merge works correctly") do

   context("takes colnames kwarg correctly") do
     @fact merge(cl,ohlc["High", "Low"], colnames=["a","b","c"]).colnames[1] => "a"
     @fact merge(cl,ohlc["High", "Low"], colnames=["a","b","c"]).colnames[2] => "b"
     @fact merge(cl,ohlc["High", "Low"], colnames=["a","b","c"]).colnames[3] => "c"
     @fact merge(cl,op, colnames=["a","b"]).colnames[1]                      => "a"
     @fact merge(cl,op, colnames=["a","b"]).colnames[2]                      => "b"
     @fact merge(cl,op, colnames=["a"]).colnames[1]                          => "Close"
     @fact merge(cl,op, colnames=["a"]).colnames[2]                          => "Open"
     @fact_throws merge(cl,op, colnames=["a","b","c"])
     @fact_throws merge(cl,op, colnames=["a","b","c"])
  end

  context("returns correct alignment with Dates and values") do
    @fact merge(cl,op).values[2,1] => cl.values[2,1]
    @fact merge(cl,op).values[2,2] => op.values[2,1]
  end

  context("aligns with disparate sized objects") do
     @fact merge(cl, op[2:5]).values[1,1]  => cl.values[2,1]
     @fact merge(cl, op[2:5]).values[1,2]  => op.values[2,1]
     @fact merge(cl, op[2:5]).timestamp[1] => Date(2000,1,4)
     @fact length(merge(cl, op[2:5]))      => 4
  end
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HairpinIncr.mso								$
#*------------------------------------------------------------------
type Results_Hairpin
	Results_Hairpin()=begin
		new(
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle Pressure Drop",
				:Default=>0.01,
				:Lower=>1e-10,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle Pressure Drop",
				:Default=>0.01,
				:Lower=>1e-10,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Total Pressure Drop",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle Velocity",
				:Default=>1,
				:Upper=>1e5
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle Velocity",
				:Default=>1,
				:Upper=>1e5
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle rho-V^2",
				:Default=>1,
				:Upper=>1e6,
				:Unit=>"kg/s^2/m"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle rho-V^2",
				:Default=>1,
				:Upper=>1e6,
				:Unit=>"kg/s^2/m"
			)),
			heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Average Film Coefficient",
				:Default=>1,
				:Lower=>1e-12,
				:Upper=>1e6,
				:DisplayUnit=>"W/m^2/K"
			)),
			[:Pdnozzle_in,:Pdnozzle_out,:Pdrop,:Vnozzle_in,:Vnozzle_out,:RVsquare_in,:RVsquare_out,:hcoeff,]
		)
	end
	Pdnozzle_in::press_delta
	Pdnozzle_out::press_delta
	Pdrop::press_delta
	Vnozzle_in::velocity
	Vnozzle_out::velocity
	RVsquare_in::positive
	RVsquare_out::positive
	hcoeff::heat_trans_coeff
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Results_Hairpin
function atributes(in::Results_Hairpin,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Results_Hairpin(_::Dict{Symbol,Any})=begin
	newModel=Results_Hairpin()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using LinearAlgebra

function preprocess(text)
    text = lowercase(text)
    text = replace(text, "." => " .")
    words = split(text, " ")

    word_to_id = Dict()
    id_to_word = Dict()

    for word in words
        if false == get(word_to_id, word, false)
            new_id = length(word_to_id) + 1     # 単語IDは配列のインデックス番号として使用するため0を含まないようにする
            push!(word_to_id, word => new_id)   # word_to_id = Dict(word_to_id..., word=>new_id) より高速(アロケーションは大きい)
            push!(id_to_word, new_id => word)
        end
    end

    corpus = [word_to_id[w] for w in words]
    return corpus, word_to_id, id_to_word
end

function create_co_matrix(corpus, vocab_size; windows_size = 1)
    """
    共起行列の作成

    パラメータ:
    corpus コーパス（単語IDのリスト）
    vocab_size 語彙の大きさ
    windows_size 窓の大きさ（これが1の場合，単語の左右1単語がコンテキスト）

    返り値:
    co_matrix 共起行列
    """

    corpus_size = length(corpus)
    co_matrix = zeros((vocab_size, vocab_size))

    for (idx, word_id) in enumerate(corpus)
        for i = 1:windows_size
            left_idx = idx - i
            right_idx = idx + i

            if left_idx >= 1
                left_word_id = corpus[left_idx]
                co_matrix[word_id, left_word_id] += 1
            end

            if right_idx <= corpus_size
                right_word_id = corpus[right_idx]
                co_matrix[word_id, right_word_id] += 1
            end
        end
    end

    return co_matrix
end

function cos_similarity(x, y; ϵ = 1e-8)
    """
    コサイン類似度の算出
    パラメータ:
    x ベクトル
    y ベクトル
    ϵ ゼロ除算の防止のための定数

    返り値:
    内積
    """
    nx = x ./ (sqrt(sum(x .^ 2)) .+ ϵ)
    ny = y ./ (sqrt(sum(y .^ 2)) .+ ϵ)
    return dot(nx, ny)
end

function most_similar(query, word_to_id, id_to_word, word_matrix; top = 5)
    """
    類似単語の検索
    パラメータ:
    query クエリ（テキスト）
    word_to_id 単語から単語IDへのディクショナリ
    id_to_word 単語IDから単語へのディクショナリ
    word_matrix 単語ベクトルをまとめた行列。各行に対応する単語のベクトルが格納されていることを想定
    top 上位何位まで表示するか
    """

    # 1.クエリを取り出す
    if false == get(word_to_id, query, false)
        println("$(query) is not found")
        return
    end

    println("[query] $(query)")
    query_id = word_to_id[query]
    query_vec = word_matrix[query_id, :]

    # 2.コサイン類似度の算出
    vocab_size = length(id_to_word)
    similarity = zeros(vocab_size)
    for i = 1:vocab_size
        similarity[i] = cos_similarity(word_matrix[i, :], query_vec)
    end

    # 3.コサイン類似度の結果から，その値を降順で出力
    count = 0
    for i in sortperm(-1 .* similarity)
        if id_to_word[i] == query
            continue
        end

        println(" $(id_to_word[i]), $(similarity[i])")

        count += 1
        if count >= top
            return
        end
    end
end

function ppmi(C; verbose = false, ϵ = 1e-8)
    """
    正の相互情報量（PPMI）の作成
    パラメータ:
    C 共起行列
    verbose 進捗を表示するか
    ϵ ゼロ除算回避のための定数

    返り値:
    M PPMI行列
    """
    M = zero(C)
    N = sum(C)
    S = sum(C, dims = 1)
    total = size(C, 1) .* size(C, 2)
    cnt = 0

    for i = 1:size(C, 1)
        for j = 1:size(C, 2)
            pmi = log2(C[i, j] .* N ./ (S[j] .* S[i]) + ϵ)
            M[i, j] = max(0, pmi...)

            if verbose == true
                cnt += 1
                if cnt % (total ÷ 100) == 0
                    println("$(round((100 .* cnt / total); digits=1))% done")
                end
            end
        end
    end
    return M
end

function create_contexts_target(corpus; window_size = 1)
    """
    コンテキストとターゲットの作成

    パラメータ:
    corpus コーパス（単語IDのリスト）
    windows_size 窓の大きさ（これが1の場合，単語の左右1単語がコンテキスト）

    返り値:
    contexts コーパス中の全単語IDのコンテキスト
    target ターゲットとする単語ID
    """
    target = corpus[window_size+1:end-window_size]
    cs = []

    for idx = (window_size+1):(length(corpus)-window_size) # window_size=1の場合，開始のインデックス番号は2
        # 各単語のコンテキストを取り出し，csに追加
        push!(cs, corpus[filter(!iszero, -window_size:window_size).+idx]')
    end
    # 全単語のコンテキストを作成
    contexts = vcat(cs...)
    return contexts, target
end

function convert_one_hot(corpus, vocab_size)
    """
    one-hot表現への変換

    パラメータ:
    corpus 単語IDのリスト（1次元もしくは2次元のNumPy配列）
    vocab_size 語彙の大きさ

    返り値
    one-hot表現（1次元もしくは2次元のNumPy配列）
    """
    N = size(corpus, 1)

    if ndims(corpus) == 1
        one_hot = zeros(Int32, N, vocab_size)
        for (idx, word_id) in enumerate(corpus)
            one_hot[idx, word_id] = 1
        end
    elseif ndims(corpus) == 2
        C = size(corpus, 2)
        one_hot = zeros(Int32, N, C, vocab_size)
        for i = 1:N
            for j = 1:C
                one_hot[i, j, corpus[i, j]] = 1
            end
        end
    end
    return one_hot
end


text:
if Sys.isapple()
    url = "https://github.com/sglyon/MINPACK.jl/releases/download/v0.0.0/libcminpack.dylib"
    download(url, joinpath(dirname(@__FILE__), "libcminpack.dylib"))
elseif Sys.iswindows()
    if Sys.WORD_SIZE == 64
        dll_url = "https://github.com/sglyon/MINPACK.jl/releases/download/v0.0.0/libcminpack.dll"
    else
        dll_url = "https://github.com/sglyon/MINPACK.jl/releases/download/v0.0.0/libcminpack32.dll"
    end
    download(dll_url, joinpath(dirname(@__FILE__), "libcminpack.dll"))
elseif Sys.islinux()
    if Sys.WORD_SIZE == 64
        so_url = "https://github.com/sglyon/MINPACK.jl/releases/download/v0.0.0/libcminpack.so"
    else
        so_url = "https://github.com/sglyon/MINPACK.jl/releases/download/v0.0.0/libcminpack32.so"
    end
    download(so_url, joinpath(dirname(@__FILE__), "libcminpack.so"))
end


text:
@testset "kmers" begin
    @testset "kmer_count" begin
        input = "AAAAAC"
        bin = kmer_count(input, 5)
        @test bin[1] == 1
        @test bin[2] == 1
    end

    @testset "corrected_kmer_dist" begin
        input1 = [5.0, 5.0, 3.0, 2.0]
        input2 = [5.0, 6.0, 6.0, 2.0]
        k = 2
        ret = corrected_kmer_dist(input1, input2; k = 2)

        @test ret == sqeuclidean(input1, input2)/ (k*(sum(input1) + sum(input2)))
    end

end


text:
include("grid.jl")
include("board.jl")
include("game.jl")
include("player.jl")
include("round.jl")
include("state.jl")
include("clock.jl")
include("bag.jl")
include("log.jl")
include("piece.jl")
include("block.jl")
include("entry.jl")


text:
# Recurrence quantifaction analysis measures

###########################################################################################
# 0. Recurrence rate based
###########################################################################################
"""
    recurrencerate(R[; theiler])

Calculate the recurrence rate of the recurrence matrix `R`.

## Description

The recurrence rate is calculated as:

```math
RR = \\frac{1}{S} \\sum R
```

where ``S`` is the size of `R` or the region of `R` with potential recurrent points.
There is not a unique definition of that denominator, which is defined as the
full size of the matrix in many sources (e.g. [1]), whereas
in others it is adjusted to remove the points of the LOI when they are
excluded from the count [2,3].

For matrices of type `RecurrenceMatrix` or `JointRecurrenceMatrix`, where the
points around the central diagonal are usually excluded, the denominator is
adjusted to the size of the matrix outside the Theiler window
(by default equal to the LOI, and adjustable with the keyword argument `theiler`;
see [`rqa`](@ref) for details). For matrices of type `CrossRecurrenceMatrix`,
where normally all points are analyzed, the denominator is always the full
size of the matrix, regardless of the Theiler window that might be defined
(none by default).

*Hint*: to reproduce the calculations done following the formulas that use
the full size of the matrix in the denominator, use
`CrossRecurrenceMatrix(s,s,ε)` to define the recurrence matrix, instead of
`RecurrenceMatrix(s,ε)`, setting `theiler=1` (or `theiler=n` in general) to
explicitly exclude the LOI or other diagonals around it.

## References

[1] : N. Marwan *et al.*, "Recurrence plots for the analysis of complex systems",
*Phys. Reports 438*(5-6), 237-329 (2007).
[DOI:10.1016/j.physrep.2006.11.001](https://doi.org/10.1016/j.physrep.2006.11.001)

[2] : C.L. Webber & J.P. Zbilut, "Recurrence Quantification Analysis of Nonlinear
Dynamical Systems", in: Riley MA & Van Orden GC, Tutorials in Contemporary
Nonlinear Methods for the Behavioral Sciences, 26-94 (2005).
URL: https://www.nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.pdf

[3] : N. Marwan & C.L. Webber, "Mathematical and computational foundations of
recurrence quantifications", in: Webber, C.L. & N. Marwan (eds.), *Recurrence
Quantification Analysis. Theory and Best Practices*, Springer, pp. 3-43 (2015).
"""
function recurrencerate(R::ARM; theiler::Integer=deftheiler(R), kwargs...)::Float64
    (theiler < 0) && throw(ErrorException(
        "Theiler window length must be greater than or equal to 0"))
    if theiler == 0
        return nnz(R)/length(R)
    end
    diags_remove = -(theiler-1):(theiler-1)
    theiler_nz = 0
    for d in diags_remove
        theiler_nz += nnz(diag(R,d))
    end
    return (nnz(R) - theiler_nz)/_rrdenominator(R; theiler=theiler, kwargs...)
end

# Calculate the denominator for the recurrence rate
_rrdenominator(R::ARM; theiler=0, kwargs...) = length(R)

function _rrdenominator(R::M; theiler=0, kwargs...) where
    M<:Union{RecurrenceMatrix,JointRecurrenceMatrix}

    (theiler == 0) && (return length(R))
    k = size(R,1) - theiler
    return k*(k+1)
end

###########################################################################################
# 0. Histograms
###########################################################################################
macro histogram_params(keyword, description, hist_fun)
    combined_descriptions = Dict(:average => "average of the $(description)s",
                                 :max     => "longest $(description)",
                                 :entropy => "Shannon entropy of the $(description)s")
    function_bodies = Dict(
        :average => quote
                (hist==[0]) && return 0.0
                points = (1:length(hist)) .* hist
                return sum(points)/sum(hist)
            end,
        :max => quote
                (hist==[0]) && return 0
                return length(hist)
            end,
        :entropy => quote
                (hist==[0]) && return NaN
                prob_bins = hist ./ sum(hist)
                prob_bins = prob_bins[findall(!iszero, prob_bins)]
                return -sum(prob_bins .* log.(prob_bins))
            end
        )

    ret = quote end
    for (name, param) in combined_descriptions
        fname = Symbol("$(keyword)_$(name)")
        _fname = Symbol("_$(keyword)_$(name)")
        fbody = function_bodies[name]
        doc = """
                $fname(R[; lmin=2, theiler])

            Calculate the $(param) contained in the recurrence matrix `R`,
            ruling out the lines shorter than `lmin` (2 by default) and all the
            points inside the Theiler window (see [`rqa`](@ref) for the
            default values and usage of the keyword argument `theiler`).
            """
        push!(ret.args, quote
            @doc $doc ->
            $fname(R::ARM; kwargs...) = $_fname($hist_fun(R; kwargs...))

            function $_fname(hist::Vector{<:Integer})
                $fbody
            end
        end)
    end
    return esc(ret)
end

###########################################################################################
# 1. Based on diagonal lines
###########################################################################################
@histogram_params dl "diagonal line" diagonalhistogram

"""
    determinism(R[; lmin=2, theiler])

Calculate the determinism of the recurrence matrix `R`:

## Description

The determinism is calculated as:

```math
DET = \\frac{\\sum_{l=lmin}{l P(l)}}{\\sum_{l=1}{l P(l)}} =
\\frac{\\sum_{l=lmin}{l P(l)}}{\\sum R}
```

where ``l`` stands for the lengths of diagonal lines in the matrix, and ``P(l)``
is the number of lines of length equal to ``l``.

`lmin` is set to 2 by default, and this calculation rules out all the
points inside the Theiler window (see [`rqa`](@ref) for the
default values and usage of the keyword argument `theiler`).
"""
function determinism(R::ARM; kwargs...)
    npoints = recurrencerate(R; kwargs...)*_rrdenominator(R; kwargs...)
    return _determinism(diagonalhistogram(R; kwargs...), npoints)
end

function _determinism(diag_hist::Vector{<:Integer}, npoints)::Float64
    (diag_hist==[0]) && return 0.0
    diag_points = (1:length(diag_hist)) .* diag_hist
    return sum(diag_points)/npoints
end

"""
    divergence(R[; theiler])

Calculate the divergence of the recurrence matrix `R`
(actually the inverse of [`dl_max`](@ref)).
"""
divergence(R::ARM; kwargs...) = ( 1.0/dl_max(R; kwargs...) )

"""
    trend(R[; border=10, theiler])

Calculate the trend of recurrences in the recurrence matrix `R`.

## Description

The trend is the slope of the linear regression that relates the density of
recurrent points in the diagonals parallel to the LOI and the distance between
those diagonals and the LOI. It quantifies the degree of system stationarity,
such that in recurrence plots where points "fade away" from the central diagonal,
the trend will have a negative value.

It is calculated as:

```math
TREND = 10^3\\frac{\\sum_{d=\\tau}^{\\tilde{N}}\\delta[d]\\left(RR[d]-\\langle RR[d]\\rangle\\right)}{\\sum_{d=\\tau}^{\\tilde{N}}\\delta[d]^2}
```

where ``RR[d]`` is the local recurrence rate of the diagonal ``d``,
``\\delta[d]`` is a balanced measure of the distance between that diagonal and the LOI,
``\\tau`` is the Theiler window (number of central diagonals that are excluded), and
``\\tilde{N}`` is the number of the outmost diagonal that is included.

This parameter is expressed in units of variation recurrence rate every
1000 data points, hence the factor ``10^3`` in the formula [1].

The 10 outermost diagonals (counting from the corners of the matrix)
are excluded by default to avoid "border effects". Use the keyword argument
`border` to define a different number of excluded lines, and `theiler`
to define the size of the Theiler window (see [`rqa`](@ref) for details).

*Note*: In rectangular cross-recurrence plots (i.e. when the time series that
originate them are not of the same length), the limits of the formula for TREND
are not clearly defined. For the sake of consistency, this function limits the
calculations to the biggest square matrix that contains the LOI.

## References
[1] C.L. Webber & J.P. Zbilut, "Recurrence Quantification Analysis of Nonlinear
Dynamical Systems", in: Riley MA & Van Orden GC, *Tutorials in Contemporary
Nonlinear Methods for the Behavioral Sciences*, 2005, 26-94.
https://www.nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.pdf
"""
trend(R::ARM; theiler=deftheiler(R), kwargs...) =
    _trend(tau_recurrence(R); theiler=theiler, kwargs...)

function tau_recurrence(R::ARM)
    n = minimum(size(R))
    rv = rowvals(R)
    rr_τ = zeros(n)
    @inbounds for col=1:n
        for i in nzrange(R, col)
            if (r=rv[i]) ≤ n
                d = abs(r-col)
                if d==0
                    rr_τ[1] += 1.0/(n-d)
                else
                    rr_τ[d+1] += 0.5/(n-d)
                end
            end
        end
    end
    return rr_τ
end

function _trend(rr_τ::Vector; theiler=1, border=10, kwargs...)::Float64
    nmax = length(rr_τ)
    a = 1+theiler
    b = nmax-border
    numerator = denominator = 0.0
    mean_rr = mean(@view rr_τ[a:b])
    for d = a:b
        δ = d - b/2
        numerator += δ*(rr_τ[d] - mean_rr)
        denominator += δ*δ
    end
    return 1000.0*numerator/denominator
end

# Number of l-length sequences, based on diagonals
# countsequences(R::ARM; kwargs...) = countsequences(diagonalhistogram(R; kwargs...))
#
# function countsequences(diag_hist::Vector; lmin=2, kwargs...)
#     overlap = (1:length(diag_hist))' .- (lmin+1)
#     overlap[overlap .< 0] = 0
#     overlap * diag_hist
# end

###########################################################################################
# 2. Based on vertical lines
###########################################################################################

@histogram_params vl "vertical line" vl_histogram

"""
    laminarity(R[; lmin=2, theiler])

Calculate the laminarity of the recurrence matrix `R`.

## Description

The laminarity is calculated as:

```math
LAM = \\frac{\\sum_{v=lmin}{v P(l)}}{\\sum_{v=1}{v P(v)}} =
\\frac{\\sum_{v=lmin}{v P(l)}}{\\sum R}
```

where ``v`` stands for the lengths of vertical lines in the matrix, and ``P(v)``
is the number of lines of length equal to ``v``.

`lmin` is set to 2 by default, and this calculation rules out all the
points inside the Theiler window (see [`rqa`](@ref) for the
default values and usage of the keyword argument `theiler`).
"""
function laminarity(R::ARM; kwargs...)
    npoints = recurrencerate(R)*_rrdenominator(R; kwargs...)
    return _laminarity(verticalhistograms(R; kwargs...)[1], npoints)
end

_laminarity(vert_hist::Vector{<:Integer}, npoints) = _determinism(vert_hist, npoints)

"""
    trappingtime(R[; lmin=2, theiler])

Calculate the trapping time of the recurrence matrix `R`, ruling out the
lines shorter than `lmin` (2 by default) and all the
points inside the Theiler window (see [`rqa`](@ref) for the
default values and usage of the keyword argument `theiler`).

The trapping time is the average of the vertical line structures and thus equal
to [`vl_average`](@ref).
"""
trappingtime(R::ARM; kwargs...) = vl_average(R; kwargs...)

###########################################################################################
# 3. Based on recurrence times
###########################################################################################

@histogram_params rt "recurrence time" rt_histogram

"""
    meanrecurrencetime(R[; lmin=2, theiler])

Calculate the mean recurrence time of the recurrence matrix `R`, ruling out the
lines shorter than `lmin` (2 by default) and all the
points inside the Theiler window (see [`rqa`](@ref) for the
default values and usage of the keyword argument `theiler`).

Equivalent to [`rt_average`](@ref).
"""
meanrecurrencetime(R::ARM; kwargs...) = rt_average(R; kwargs...)

"""
    nmprt(R[; lmin=2, theiler])

Calculate the number of the most probable recurrence time (NMPRT), ruling out the
lines shorter than `lmin` (2 by default) and all the
points inside the Theiler window (see [`rqa`](@ref) for the
default values and usage of the keyword argument `theiler`).

This number indicates how many times the system has recurred using the recurrence
time that appears most frequently, i.e it is the maximum value of the histogram
of recurrence times [1].

## References

[1] : E.J. Ngamga *et al.* "Recurrence analysis of strange nonchaotic dynamics",
*Physical Review E*, 75(3), 036222(1-8) (2007)
[DOI:10.1103/physreve.75.036222](https://doi.org/10.1103/physreve.75.036222)

"""
nmprt(R::ARM; kwargs) = maximum(verticalhistograms(R; kwargs...)[2])

###########################################################################################
# 4. All in one
###########################################################################################

# Transient type
struct RQA
    data::Dict
end

function Base.getproperty(result::RQA, name::Symbol)
    if name === :data
        return getfield(result, :data)
    end
    @warn "x.$name is deprecated for results of `rqa`; use x[:$name] instead"
    getindex(result.data, name)
end

for f in (:getindex, :get, :length, :keys, :values, :collect, :iterate)
    @eval Base.$f(result::RQA, args...; kwargs...) = $f(result.data, args...; kwargs...)
end

Dict(rqa::RQA) = rqa.data

function Base.show(io::IO, mime::MIME"text/plain", result::RQA)
    print(io, "RQA parameters in ")
    show(io, mime, result.data)
end

"""
    rqa(R; kwargs...)

Calculate all RQA parameters of a recurrence matrix `R`. See the functions
referred to below for the definition of
the different parameters and the default values of the arguments.
Using this function is much more efficient than calling all individual functions
one by one.

## Return
The returned value contains the following entries,
which can be retrieved as from a dictionary (e.g. `results[:RR]`, etc.):

* `:RR`: recurrence rate (see [`recurrencerate`](@ref))
* `:DET`: determinsm (see [`determinism`](@ref))
* `:L`: average length of diagonal structures (see [`dl_average`](@ref))
* `:Lmax`: maximum length of diagonal structures (see [`dl_max`](@ref))
* `:DIV`: divergence (see [`divergence`](@ref))
* `:ENTR`: entropy of diagonal structures (see [`dl_entropy`](@ref))
* `:TREND`: trend of recurrences (see [`trend`](@ref))
* `:LAM`: laminarity (see [`laminarity`](@ref))
* `:TT`: trapping time (see [`trappingtime`](@ref))
* `:Vmax`: maximum length of vertical structures (see [`vl_max`](@ref))
* `:VENTR`: entropy of vertical structures (see [`vl_entropy`](@ref))
* `:MRT`: mean recurrence time (see [`meanrecurrencetime`](@ref))
* `:RTE` recurrence time entropy (see [`rt_entropy`](@ref))
* `:NMPRT`: number of the most probable recurrence time (see [`nmprt`](@ref))

All the parameters returned by `rqa` are `Float64` numbers,
even for parameters like `:Lmax`, `:Vmax` or `:NMPRT` which are integer values.
In the case of empty histograms (e.g. no existing vertical lines
less than the keyword `lminvert`) the average and maximum values
(`:L`, `:Lmax`, `:TT`, `:Vmax`, `:MRT`)
are returned as `0.0` but their respective entropies (`:ENTR`, `:VENTR`, `:RTE`)
are returned as `NaN`.

## Keyword Arguments

Standard keyword arguments are the ones accepted by the functions listed below,
i.e. `theiler`, `lmin`, and `border`:

* `theiler` is used to define a "Theiler window" around the central diagonal or
  "line of identity" (LOI): a region of points that are excluded in the calculation
  of RQA parameters, in order to rule out self-recurrences and apparent recurrences
  for smooth or high resolution data. The LOI is excluded by default for matrices
  of the types `RecurrenceMatrix` or `JointRecurrenceMatrix`, but it is included
  for matrices of the type `CrossRecurrenceMatrix`. `theiler=0` means that the
  whole matrix is scanned for lines. `theiler=1` means that the LOI is excluded.
  In general, `theiler=n` means that the `n` central diagonals are excluded
  (at both sides of the LOI, i.e. actually `2n-1` diagonals are excluded).

* `lmin` is used to define the minimum line length in the parameters that
  describe the distributions of diagonal or vertical lines (it is set as 2 by
  default).

* `border` is used to avoid border effects in the calculation of `:TREND`
  (cf. [`trend`](@ref)).

In addition `theilerdiag`, `lmindiag` may be used to
declare specific values that override the values of `theiler` and `lmin` in the
calculation of parameters related to diagonal structures. Likewise, `theilervert` and
`lminvert` can be used for the calculation of parameters related to vertical
structures.

The keyword argument `onlydiagonal` (`false` by default) can be set to `true`
in order to restrict the analysis to the recurrence rate and the parameters related
to diagonal structures (`:RR`, `:DET`, `:L`, `:Lmax`, `:DIV` and `:ENTR`), which makes
this function slightly faster.

## Transitional note on the returned type

In older versions, the `rqa` function returned a `NamedTuple`,
and in future versions it is planned to return a `Dict` instead.
In both cases, the results can be indexed with square brackets and `Symbol` keys,
as `result[:RR]`, `result[:DET]`, etc.
However, named tuples can also be indexed with "dot syntax", e.g. `result.RR`,
whereas this will not be possible with dictionaries, and there are other
differences in the indexing and iteration of those two types.

In order to facilitate the transition between versions, this function currently
returns a `RQA` object that essentially works as a dictionary, but can
also be indexed with the dot syntax (logging a deprecation warning).
The returned type can also be specified as a first argument of `rqa`
in order to replicate the output of different versions:

* `rqa(NamedTuple, R...)` to obtain the output of the older version (as in 1.3).
* `rqa(Dict, R...)` to obtain the output of the planned future version.
* `rqa(RQA, R...)` to obtain the default current output (same as `rqa(R...)`)

"""
rqa(R; kwargs...) = rqa(RQA, R; kwargs...)

function rqa(::Type{RQA}, R; kwargs...)
    rqa_dict = rqa(Dict, R; kwargs...)
    RQA(rqa_dict)
end

function rqa(::Type{Dict}, R; onlydiagonal=false, kwargs...)
    # Parse arguments for diagonal and vertical structures
    kw_d = Dict(kwargs)
    haskey(kw_d, :theilerdiag) && (kw_d[:theiler] = kw_d[:theilerdiag])
    haskey(kw_d, :lmindiag) && (kw_d[:lmin] = kw_d[:lmindiag])
    dhist = diagonalhistogram(R; kw_d...)
    rr_d = recurrencerate(R; kw_d...)
    if onlydiagonal
        return Dict{Symbol, Float64}(
            :RR    => recurrencerate(R; kwargs...),
            :DET   => _determinism(dhist, rr_d*_rrdenominator(R; kw_d...)),
            :L     => _dl_average(dhist),
            :Lmax  => _dl_max(dhist),
            :DIV   => 1.0/_dl_max(dhist),
            :ENTR  => _dl_entropy(dhist)
        )
   else
        kw_v = Dict(kwargs)
        haskey(kw_v, :theilervert) && (kw_v[:theiler] = kw_v[:theilervert])
        haskey(kw_v, :lminvert) && (kw_v[:lmin] = kw_v[:lminvert])
        vhist, rthist = verticalhistograms(R; kw_v...)
        rr_v = recurrencerate(R; kw_v...)
        rqa_dict = Dict{Symbol, Float64}(
            :RR    => rr_d,
            :DET   => _determinism(dhist, rr_d*_rrdenominator(R; kw_v...)),
            :L     => _dl_average(dhist),
            :Lmax  => _dl_max(dhist),
            :DIV   => 1.0/_dl_max(dhist),
            :ENTR  => _dl_entropy(dhist),
            :TREND => trend(R; kw_d...),
            :LAM   => _laminarity(vhist, rr_v*_rrdenominator(R; kw_v...)),
            :TT    => _vl_average(vhist),
            :Vmax  => _vl_max(vhist),
            :VENTR => _vl_entropy(vhist),
            :MRT   => _rt_average(rthist),
            :RTE   => _rt_entropy(rthist),
            :NMPRT => maximum(rthist)
        )
        return rqa_dict
    end
end


text:
try
    # 1. do some work related to reading writing files
    # 2. invoke an HTTP request to a remote web service
    # 3. create a status report in PDF and save in a network drive
catch ex
    if ex isa FileNotFoundError
        println("Having trouble with reading local file")
        exit(1)
    elseif ex isa HTTPRequestError
        println("Unable to communicate with web service")
        exit(2)
    elseif ex isa NetworkDriveNotReadyError
        println("All done, except that the report cannot be saved")
        exit(3)
    else
        println("An unknown error has occured, please report. Error=", ex)
        exit(255)
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.7464435003747418 + -0.7676428720221384*1/(1 + exp(-(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])))) + -0.5937285886485864*1/(1 + exp(-(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])))) + 0.70936536775501*1/(1 + exp(-(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])))) + 0.9647291374973435*1/(1 + exp(-(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])))) + 0.527810271934432*1/(1 + exp(-(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))))) + 1/(1 + exp(-(-0.35636357579153266 + 0.8665784119884141*1/(1 + exp(-(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])))) + 0.03175667478129274*1/(1 + exp(-(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])))) + -0.7799384758999897*1/(1 + exp(-(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])))) + 0.4104070349373172*1/(1 + exp(-(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])))) + -0.7416668360909755*1/(1 + exp(-(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))))) + 1/(1 + exp(-(0.5038932713674513 + 0.7223377417739489*1/(1 + exp(-(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])))) + -0.17110213406275054*1/(1 + exp(-(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])))) + 0.7726394930743403*1/(1 + exp(-(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])))) + 0.6118790711423934*1/(1 + exp(-(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])))) + 0.16250007841990932*1/(1 + exp(-(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))))) + 1/(1 + exp(-(-0.7190092552996492 + 0.6952008118484776*1/(1 + exp(-(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])))) + 0.975536624348428*1/(1 + exp(-(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])))) + 0.03367235210941821*1/(1 + exp(-(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])))) + 0.8059867473566515*1/(1 + exp(-(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])))) + 0.5874603423029647*1/(1 + exp(-(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))))) + 1/(1 + exp(-(-0.5715134646750415 + 0.6325466436773857*1/(1 + exp(-(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])))) + -0.2633742432221431*1/(1 + exp(-(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])))) + -0.3786600783719858*1/(1 + exp(-(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])))) + 0.5923927477083897*1/(1 + exp(-(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])))) + 0.9114166052727355*1/(1 + exp(-(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using MacroTools: @q
using GeneralizedGenerated
import PyCall
using MLStyle

import SymPy
using SymPy: Sym, symbols, free_symbols
import SpecialFunctions
using SpecialFunctions: logfactorial

"""As type encoding a PyObject is unsafe without some hard
works with reference counting, we simply use a Julia proxy
to imitate the `sympy` module, e.g.,
    `sympy.attr = _pysympy.attr`
"""
struct PySymPyModule end
GeneralizedGenerated.NGG.@implement GeneralizedGenerated.NGG.Typeable{PySymPyModule}
const sympy = PySymPyModule()
Base.getproperty(::PySymPyModule, s::Symbol) = Base.getproperty(_pysympy, s)

const symfuncs = Dict()

_pow(a,b) = Base.:^(float(a),b)

function __init__()
    stats = PyCall.pyimport_conda("sympy.stats", "sympy")
    SymPy.import_from(stats)
    global stats = stats
    global _pysympy = SymPy.sympy

    # for dist in [:Normal, :Cauchy, :Laplace, :Beta, :Uniform]
    #     @eval begin
    #         function Distributions.$dist(μ::Sym, σ::Sym)
    #             println("Evaluating ",$dist, "(",μ," :: Sym, ", σ, " :: Sym)")
    #             stats.$dist(:dist, μ,σ) |> SymPy.density
    #         end

    #         Distributions.$dist(μ,σ) = $dist(promote(μ,σ)...)
    #     end
    # end

    # https://discourse.julialang.org/t/pyobjects-as-keys/26521/2
    merge!(symfuncs, Dict(
        sympy.log => Base.log
      , sympy.Pow => _pow
      , sympy.Abs => Base.abs
      , sympy.Indexed => getindex
    ))

    @eval begin
    @gg M function codegen(_::Type{M}, _m::Model, _args, _data) where M <: TypeLevel{Module}
        f = _codegen(type2model(_m))
        Expr(:let,
            Expr(:(=), :M, from_type(M)),
            :($f(_args, _data)))
    end
end

end

struct SymDist
    logpdf :: Function
end

logpdf(d::SymDist, x) = d.logpdf(sym(x))

# julia> logpdf(SymPy.density(Soss.stats.Poisson(:Poisson,sym(:λ))), sym(:x))
# x⋅log(λ) - λ - log(x!)

SpecialFunctions.logfactorial(x::Sym) = sympy.loggamma(x+1)

Distributions.Poisson(λ::Sym) = SymDist(x -> x * log(λ) - λ - logfactorial(x))
Distributions.Bernoulli(p::Sym) = SymDist(y -> y * log(p) + (1-y) * log(1-p))

for dist in [:Bernoulli, :Poisson]
    @eval begin
        logpdf(d::$dist, x::Sym) = logpdf($dist(sym.(Distributions.params(d))...), x)
    end
end

# "Half" distributions
for dist in [:Normal, :Cauchy]
    let half = Symbol(:Half, dist)
        @eval begin
            logpdf(d::$half, x::Sym) = 2 * logpdf($dist(0, sym.(d.σ)), x)
        end
    end
end

for dist in [:Normal, :Cauchy, :Laplace, :Beta, :Uniform]
    @eval begin
        function Distributions.$dist(μ::Sym, σ::Sym)
            stats.$dist(:dist, μ,σ) |> SymPy.density
        end

        Distributions.$dist(μ,σ) = $dist(promote(μ,σ)...)
    end
end

export sym
sym(s::Symbol) = sympy.symbols(s, real=true)
sym(s) = Base.convert(Sym, s)
function sym(expr::Expr)
    @match expr begin
        Expr(:call, f, args...) => :($f($(map(sym,args)...)))
        :($x[$j]) => begin
            j = symbols(:j, cls=sympy.Idx)
            x = sympy.IndexedBase(x, real=true)
            return x[j]
        end
        _ => begin
                 @show expr
                 error("sym: Argument type not implemented")
             end
    end
end

export expandSums
function expandSums(s::Sym)
    s.args==() && return s
    if s.func==sympy.Sum
        s2 = expandSum(s.args...)
        # Maybe it didn't simplify and we're done
        s2 == s && return s2
        # Or not, and we recurse
        s = s2
    end
    newargs = [expandSums(t) for t in s.args]
    result = s.func(newargs...)
    # @show result - s
    # @assert sympy.simplify(result - s) == sym(0)
    return result
end

function expandSum(s::Sym, limits::Sym...)

    func = s.func
    args = s.args
    args == () && return maybeSum(s,limits...)

    func == sympy.Add && return func([expandSum(t, limits...) for t in args]...)
    func == sympy.Mul && return expandMulSum(args, limits...)

    return maybeSum(s, limits...)

end

# Expand a sum of a `Mul`
# Currently a greedy algorithm, may need to optimize later
function expandMulSum(factors::NTuple{N,Sym}, limits::Sym...) where {N}
    limits == () && return prod(factors)

    p = sympy.expand_mul(prod(factors))
    p.func == sympy.Mul || return expandSum(p, limits...)
    factors = p.args
    for fac in factors
        for lim in limits
            (ix, ixlo, ixhi) = lim.args
            if ix ∉ atoms(fac)
                inSummand = prod(allbut(factors, fac))
                inSum = expandSum(inSummand, lim)
                outLims = allbut(limits, lim)
                return expandSum(fac*inSum, outLims...)
            end
        end
    end
    return maybeSum(prod(factors), limits...)
end

function atoms(s::Sym)
    s.func == sympy.Symbol && return [s]
    s.func == sympy.Idx && return [s]
    isempty(free_symbols(s)) && return []

    s.func == sympy.Sum && begin
        summand = s.args[1]
        ixs = [j.args[1] for j in s.args[2:end]]
        bounds = union([j.args[2:3] for j in s.args[2:end]]...)
        return setdiff(union(atoms(summand), atoms.(bounds)...), ixs)
    end
    result = union(map(atoms, s.args)...)
    return result
end

function allbut(tup, x)
    result = filter(collect(tup)) do v
        v ≠ x
    end
    tuple(result...)
end

# Force computation of sums that don't involve the index
function maybeSum(t::Sym, limits::Sym...)
    length(limits) > 0 || return t

    for lim in limits
        (ix, ixlo, ixhi) = lim.args
        ix ∈ atoms(t) || begin
            return maybeSum(t * (ixhi - ixlo + 1), allbut(limits, lim)...)
        end
    end

    return sympy.Sum(t, limits...)
end

# # integrate(exp(ℓ), (sym(:μ), -oo, oo), (sym(:logσ),-oo,oo))
export marginal
function marginal(ℓ,v)
    f = ℓ.func
    f == sympy.Add || return ℓ
    newargs = filter(t -> sym(v) in atoms(t), collect(ℓ.args))
    foldl(+,newargs)
end

export score
score(ℓ::Sym, v) = diff(ℓ, sym(v))

score(m::Model, v) = score(symlogpdf(m), v)

marginal(m::Model, v) = marginal(m |> symlogpdf, v)

symvar(st) = :($sympy.IndexedBase($(st.x)))

function symvar(st::Sample)
    st.rhs.args[1] ∈ [:For, :iid] && return :($sympy.IndexedBase($(st.x)))
    return :($sym($(st.x)))
end

export sourceSymlogpdf
function sourceSymlogpdf()
    function(_m::Model)
        function proc(_m, st :: Assign)
            x = st.x
            xname = QuoteNode(x)
            return :($x = $sympy.IndexedBase($xname))
        end

        function proc(_m, st :: Sample)
            s = :(_ℓ += symlogpdf($(st.rhs), $(symvar(st))))
            end
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing

        function wrap(kernel)
            q = @q begin
                _ℓ = 0.0
            end

            for x in arguments(_m)
                xname = QuoteNode(x)
                xsym = :($sympy.IndexedBase($xname))
                push!(q.args, :($x = $xsym))
            end

            for st in map(v -> findStatement(_m,v), toposort(_m))

                typeof(st) == Sample || continue
                x = st.x
                xname = QuoteNode(x)
                rhs = st.rhs
                xsym = ifelse(rhs.args[1] ∈ [:For, :iid]
                    , :($sympy.IndexedBase($xname))
                    , :($sym($xname))
                )
                push!(q.args, :($x = $xsym))
            end

            @q begin
                $q
                $kernel
                return _ℓ
            end
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end

For(f, θ::Sym) = For(f, (θ,))

function For(f::F, θ::NTuple{N,Sym}) where {F, N}
    For{F,NTuple{N,Sym},Sym,Sym}(f,θ)
end

# logpdf(d::For{F,N,T}, x::Array{Symbol, N}) where {F,N,T} = logpdf(d,sym.(x))

export symlogpdf
function symlogpdf(d::For{F,T,D,X}, x::Sym) where {F, N, J <: Union{Sym,Integer}, T <: NTuple{N,J}, D,  X}
    js = symbols.(Symbol.(:_j,1:N), cls=sympy.Idx)
    x = sympy.IndexedBase(x)
    result = symlogpdf(d.f(js...), x[js...])

    for k in N:-1:1
        result = sympy.Sum(result, (js[k], 1, d.θ[k]))
    end
    result
end

symlogpdf(d::iid, x::Sym) = symlogpdf(For(j -> d.dist, d.size), x)

symlogpdf(d::For{F,T,D,X}, x::Symbol) where {F,T,D,X} = symlogpdf(d,sym(x))

symlogpdf(d::Normal, x::Sym) = symlogpdf(Normal(sym(d.μ),sym(d.σ)), x)

symlogpdf(d::Cauchy, x::Sym) = symlogpdf(Cauchy(sym(d.μ),sym(d.σ)), x)

symlogpdf(d::Beta, x::Sym) = symlogpdf(Beta(sym(d.α),sym(d.β)), x)

symlogpdf(d::Poisson, x::Sym) = symlogpdf(Poisson(sym(d.λ)), x)

logpdf(d::Sym, x::Sym) = symlogpdf(d,x)

function symlogpdf(d::Sym, x::Sym)
    d.func
    result = d.pdf(x) |> log
    sympy.expand_log(result,force=true)
end

symlogpdf(d,x::Sym) = logpdf(d,x)

function symlogpdf(d::JointDistribution, simplify=true)
    symlogpdf(d.model)
end

function symlogpdf(m::Model, simplify=true)
    s = _symlogpdf(getmoduletypencoding(m), m)
    simplify && return foldConstants(expandSums(s))
    return s
end

@gg M function _symlogpdf(_::Type{M}, _m::Model) where M <: TypeLevel{Module}
    Expr(:let,
        Expr(:(=), :M, from_type(M)),
        type2model(_m) |> canonical |> sourceSymlogpdf())
end

export tolatex
function tolatex(ℓ::SymPy.Sym)
    r = r"Idx\\left\(_j(?<num>\d+)\\right\)"
    s = s"j_\g<num>"
    x = Base.replace(sympy.latex(ℓ), r => s)
end

export foldConstants
function foldConstants(s::Sym)
    s.func ==sympy.Integer && return s
    isempty(free_symbols(s)) && return Float64(SymPy.N(s))
    s.func == sympy.Sum && return sympy.Sum(foldConstants(s.args[1]), s.args[2:end]...)
    s.func == sympy.Indexed && return s
    s.func == sympy.IndexedBase && return s
    s.func == sympy.Symbol && return s
    newargs = foldConstants.(s.args)
    return s.func(newargs...)
end


text:
module NonconvexJuniper

export JuniperIpoptAlg, JuniperIpoptOptions, IpoptAlg, IpoptOptions

import JuMP, Ipopt, MathOptInterface
const MOI = MathOptInterface
using Reexport, Parameters, SparseArrays, Zygote
@reexport using NonconvexCore, NonconvexIpopt
using NonconvexCore: @params, VecModel, AbstractResult, JuMPProblem
using NonconvexCore: AbstractOptimizer, get_jump_problem, nvalues
import NonconvexCore: optimize!
using Juniper

include("juniper.jl")

end


text:
"""
    log_sum_exp(z::Array)

A stable implementation f(x) = log ∘ sum ∘ exp (x).
Since exponentiation can lead to very large numbers.
"""
log_sum_exp(z) = log_sum_exp(z, maximum(z, dims = 1))
log_sum_exp(z, m) = log.(sum(exp.(z .- m), dims = 1)) .+ m


text:
function -(S1::Set{T},S2::Set{T}) where T
    setdiff(S1,S2)
end

function ×(S1::Set{T},S2::Set{T}) where T
    S=Set()
    for s1 in S1
        for s2 in S2
            push!(S,(s1,s2))
        end
    end
    S
end

function ⊗(S1::Set{T},S2::Set{T}) where T
    S=Set()
    for s1 in S1
        for s2 in S2
            if s1<s2
                push!(S,(s1,s2))
            elseif s2<s1
                push!(S,(s2,s1))
            end
        end
    end
    S
end


text:
using Ramnet.Encoders
using Ramnet

using Random
using StaticArrays

using ..Approximators:QDiscriminator
using ..Buffers: MultiStepDynamicBuffer, add!
using ..Environments

mutable struct MCDiscriminatorAgent{AS,OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},A <: Real,E <: AbstractEncoder{T}} <: AbstractAgent{O,A}
    actions::SVector{AS,A}
    discount::Float64
    ϵ_start::Float64
    ϵ_end::Float64
    ϵ_decay::Int
    update_count::Int
    Q̂::QDiscriminator{AS,OS,T,O,A,E}
    buffer::MultiStepDynamicBuffer{O,A}
    observation::Union{Nothing,O}
    done::Bool
    action::Union{Nothing,A}
    rng::MersenneTwister
end

# TODO: Validate actions and encoder and ϵ
function MCDiscriminatorAgent(::Type{O}, actions::StaticArray{Tuple{AS},A,1}, n, forgetting_factor, discount, ϵ_start, ϵ_end, ϵ_decay, encoder::E; seed::Union{Nothing,UInt}=nothing) where {AS,OS,T,O <: StaticArray{Tuple{OS},T,1},A,E <: AbstractEncoder{T}}
    !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
    rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

    n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))

    Q̂ = QDiscriminator(O, actions, n, forgetting_factor, encoder, seed)
    buffer = MultiStepDynamicBuffer{O,A}()

    MCDiscriminatorAgent{AS,OS,T,O,A,E}(actions, discount, ϵ_start, ϵ_end, ϵ_decay, 0, Q̂, buffer, nothing, false, nothing, rng)
end

# TODO: Return type can only be determined at runtime. Make AbstractEnvironment parametric
function MCDiscriminatorAgent(
    env, encoder::E; tuple_size, forgetting_factor::Float64=1.0, discount::Float64=1.0, epsilon_start::Float64=0.01, epsilon_end::Float64=0.0, epsilon_decay::Int=1000, seed::Union{Nothing,UInt}=nothing) where {T <: Real,E <: AbstractEncoder{T}}
    # MCDiscriminatorAgent{observation_type(env),action_type(env),T,E}(action_set(env), tuple_size, observation_length(env), forgetting_factor, discount, epsilon, encoder; seed)
    MCDiscriminatorAgent(observation_type(env), action_set(env), tuple_size, forgetting_factor, discount, epsilon_start, epsilon_end, epsilon_decay, encoder; seed)
end

# function MCDiscriminatorAgent(env, encoder::E, n, forgetting_factor=1.0, discount=1.0, ϵ=0.01, seed=nothing) where {T <: Real,E <: AbstractEncoder{T}}
#     MCDiscriminatorAgent(observation_type(env), action_set(env), n, forgetting_factor, discount, ϵ, encoder; seed)
# end

function observe!(agent::MCDiscriminatorAgent{AS,OS,T,O,A,E}, observation::O) where {AS,OS,T,O,A,E} # {O,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    agent.observation = observation
    agent.done = false

    nothing
end

# TODO: All observe! methods are the same for all agents. Generalize.
# TODO: This method does not need to take in the action
function observe!(agent::MCDiscriminatorAgent{AS,OS,T,O,A,E}, action::A, reward::Float64, observation::O, done::Bool) where {AS,OS,T,O,A,E} # {O,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    add!(agent.buffer, agent.observation, agent.action, reward)

    agent.observation = observation
agent.done = done

    nothing
end

function epsilon(agent::MCDiscriminatorAgent)
    α = agent.update_count / agent.ϵ_decay
    (1 - α) * agent.ϵ_start + α * agent.ϵ_end
end

function select_action!(agent::MCDiscriminatorAgent{AS,OS,T,O,A,E}, observation::O) where {AS,OS,T,O,A,E} # {O,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    if rand(agent.rng) < epsilon(agent)
        agent.action = rand(agent.rng, agent.actions)
    else
        agent.action = select_action(agent.Q̂, observation)
    end

    return agent.action
end

function update!(agent::MCDiscriminatorAgent{AS,OS,T,O,A,E}) where {AS,OS,T,O,A,E} # {O,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    if agent.done
        G = 0.0
        for transition in Iterators.reverse(agent.buffer)
            G = transition.reward + agent.discount * G

            update!(agent.Q̂, transition.observation, transition.action, G)
        end

        agent.update_count = agent.update_count == agent.ϵ_decay ? agent.ϵ_decay : agent.update_count + 1
        reset!(agent.buffer)
    end

    nothing
end

function reset!(agent::MCDiscriminatorAgent{AS,OS,T,O,A,E}; seed::Union{Nothing,UInt}=nothing) where {AS,OS,T,O,A,E} # {O,A <: Real,S,T <: Real,E <: AbstractEncoder{T}}
    if !isnothing(seed)
        if seed ≥ 0
        seed!(agent.rng, seed)
        else
            throw(DomainError(seed, "Seed must be non-negative"))
        end
    end

    reset!(agent.Q̂; seed)
    reset!(agent.buffer)

    nothing
end

# function q_values!(agent::MCDiscriminatorAgent{O,A,S,T,E}, observation::O, dest::AbstractVector{Float64}) where {O,A <: Real,S,T <: Real,E <: AbstractEncoder{T}}
#     q_max = typemin(Float64)
#     q = 0.0
#     for (i, action) in Iterators.zip(eachindex(dest), agent.actions)
#         q = predict(agent.Q̂[action], observation)
#         dest[i] = q
#         q_max = q > q_max ? q : q_max
#     end

#     return q_max
# end


text:
using Pkg
using Random
using Test

tests = []

if length(ARGS) > 0
    tests = ARGS
end

Random.seed!(42)    # set a prng seed
tolerance = 1e-8    # set error tolerance

for t in tests
    include(t)      # include the different test suites
end


text:
include("SortSparse.jl")
include("KoLesky.jl")
using LinearAlgebra
using Distributions
using Distances
using Random
using Profile
using JLD
using MAT
Random.seed!(123)

NTrain = 8000
NTest = 100
N = NTest + NTrain
d =  2
y = randn(N)
x = randn(d, N)

setSplit = randperm(N)
trainInds = setSplit[1:NTrain]
testInds = setSplit[(NTrain + 1) : end]

xTrain = x[ :, trainInds ]
yTrain = y[ trainInds ]
xTest = x[:, testInds ]
yTest = y[ testInds ]

ρ = 3.5
λ = 1.0

l = 0.1
cov = r -> exp(-r/l)

@time colptr, rowval, P, revP, distances = sortSparseRev(xTrain,ρ, 1)
xTrain = xTrain[:, P]
yTrain = yTrain[P]
x = hcat(xTrain, xTest)

skeletons = construct_skeletons(colptr, rowval, distances, λ)
colptr = 0
length(skeletons)
rowval = 0
GC.gc()

#Comparison to full kodel
LinearAlgebra.BLAS.set_num_threads(8)
@time begin
  Θ = cov.(pairwise(Euclidean(), x, dims=2))
  μFull = Θ[(NTrain + 1) : end, 1 : NTrain] * ( Θ[1 : NTrain, 1 : NTrain] \ yTrain)
end

Θ = 0
GC.gc()
@show mean(abs.(yTest - μFull))
@show mean(abs.(yTest - μFull)) / std(yTest)

@time μ, σ = predict(skeletons, (NTrain+1): N, yTrain, x, cov)

mean(abs.(yTest - μ)) / std(yTest)
mean(abs.(yTest - μ))


text:
export intersect, renderscene!

function intersect(ray::Ray, sphere::Sphere)
    oc = (ray.origin - sphere.center)

    a = ray.dir ⋅ ray.dir
    bhalf = oc ⋅ ray.dir
    c = oc ⋅ oc - sphere.radius^2
    Δ = bhalf^2 - a * c

    if Δ > 0
        sqrtΔ = sqrt(Δ)

        t1 = (-bhalf - sqrtΔ) / a
        if t1 > 0
            return t1
        end

        t2 = (-bhalf + sqrtΔ) / a
        if t2 > 0
            return t2
        end
    end

    return zero(Δ)
end

function renderscene!(camera, shapes, image; background=RGB{Float64}(0, 0, 0))
    height, width = size(image)
    for y in 1:height
        for x in 1:width
            ray = cast_ray(camera, x, y, width, height)

            closest_dist = zero(Float64)
            closest_shape = 0
            closest_loc_ray = ray
            firsthit = true
            for i in eachindex(shapes)
                loc_ray = transform(shapes[i]) * ray
                cur_dist = intersect(loc_ray, shapes[i])
                if cur_dist > 0 && (firsthit || cur_dist < closest_dist)
                    closest_loc_ray = loc_ray
                    closest_dist = cur_dist
                    closest_shape = i
                    firsthit = false
                end
            end

            image[height - y + 1, x] = if closest_shape == 0
                background
            else
                getcolor(shapes[closest_shape], ray_at_dist(closest_loc_ray, closest_dist))
            end
        end
    end
end


text:
# Original file: https://github.com/JuliaGraphics/Luxor.jl/blob/master/src/play.jl
# Replaced whileloop with forloop in play macro to limit number of displayed frames to n

using MiniFB

include("key.jl")

current_key = Key()

function keyboard(window, key, mod, isPressed)::Cvoid

    if key == MiniFB.KB_KEY_LEFT
        if Bool(isPressed)
            current_key.left = true
        else
            current_key.left = false
        end
    elseif key == MiniFB.KB_KEY_RIGHT
        if Bool(isPressed)
            current_key.right = true
        else
            current_key.right = false
        end
    elseif key == MiniFB.KB_KEY_DOWN
        if Bool(isPressed)
            current_key.down = true
        else
            current_key.down = false
        end
    elseif key == MiniFB.KB_KEY_UP
        if Bool(isPressed)
            current_key.up = true
        else
            current_key.up = false
        end
    end
    return nothing
end

macro play(w, h, n, body)
    quote
        window = mfb_open_ex("Visualization", $(esc(w)), $(esc(h)), MiniFB.WF_RESIZABLE)

        mfb_set_keyboard_callback(window, keyboard)

        buffer = zeros(UInt32, $(esc(w)), $(esc(h)))
        for i = 1:$(esc(n))
            Drawing($(esc(w)), $(esc(h)), :image)
            $(esc(body))
            m = permutedims(image_as_matrix!(buffer), (2, 1))
            finish()
            state = mfb_update(window, m)
            if state != MiniFB.STATE_OK
                break
            end
        end
        mfb_close(window)
    end
end


text:
function countingsort(x::AbstractArray{<:Int})::AbstractArray{<:Int}
    l = length(x)
    n = maximum(x)+1
    count = zeros(Int, n)
    aux = copy(x)
    for i ∈ 1:l
        @inbounds count[x[i]+1]+=1
    end
    for i ∈ 2:n
        @inbounds count[i] += count[i-1]
    end
    for i ∈ l:-1:1
        @inbounds aux[count[x[i]+1]] = x[i]
        @inbounds count[x[i]+1]-=1
    end
    return aux
end