text:
# This file is a part of GERDADeepLearning.jl, licensed under the MIT License (MIT).

using MXNet, Formatting

import Base: getindex

function fc_layer(name, X, n_hidden, act_type, pdropout)
  fc = mx.FullyConnected(X, num_hidden=Int(n_hidden), name="$(name)")
  act = mx.Activation(fc, act_type=act_type, name="$(name)_act")
  return mx.Dropout(act, p=pdropout, name="$(name)_dropout")
end

function conv_layer(name, X, num_filter, filter_length, act_type, pool_size, pool_type, dropout)
  conv = mx.Convolution(X, kernel=(1, filter_length), pad=(0, 4), num_filter=num_filter, name="$name") # 2D->3D
  act = mx.Activation(conv, act_type=act_type, name="$(name)_act")
  pool = mx.Pooling(act, kernel=(1, pool_size), stride=(1, pool_size), pool_type=pool_type, name="$(name)_pool") # 256->64
  dropout = mx.Dropout(pool, p=dropout, name="$(name)_dropout")
  return dropout
end

function deconv_layer(name, X, num_filter, filter_length, act_type, pool_size, dropout)
  pad = Int( floor(filter_length/2) - pool_size/2 )
  if pad < 0
    throw(ArgumentError("upsampling not possible because padding is negative ($pad). Increase filter length or decrease pool size."))
  end
  X = mx.Deconvolution(X, kernel=(1, filter_length-1), stride=(1, pool_size), pad=(0, pad), num_filter=num_filter, name="$name")
  if act_type != nothing
    X = mx.Activation(X, act_type=act_type, name="$(name)_act")
  end
  X = mx.Dropout(X, p=dropout, name="$(name)_dropout")
  return X
end

type PlotCallback <: mx.AbstractEpochCallback
    graphs::Dict{Symbol, Array{Float32,1}}

    PlotCallback() = new(Dict{Symbol, Array{Float32,1}}())
end

function (cb::PlotCallback){T<:Real}(model :: Any, state :: mx.OptimizationState, metric :: Vector{Tuple{Base.Symbol, T}})
  for index in 1:length(metric)
    key = metric[index][1]
    if !haskey(cb.graphs, key)
      cb.graphs[key] = Float32[]
    end

    push!(cb.graphs[key], metric[index][2])
  end
end

function calculate_parameters(model, filepath)
  file = open(filepath, "w")
  total_parameter_count = 0
  for param in model.arg_params
    line = "Parameter $(param[1]) has shape $(size(param[2])) = $(length(param[2]))"
    write(file, line*"\n")
    total_parameter_count += length(param[2])
  end
  lastline = "Total parameter count: $total_parameter_count"
  write(file, lastline*"\n")
  close(file)
 end

 function get_total_parameter_count(model)
   return sum([length(param[2]) for param in model.arg_params])
 end

 function to_context(gpus::Vector{Int})
   if length(gpus) == 0
     return mx.cpu()
   else
     return [mx.gpu(i) for i in gpus]
   end
 end

# function exists_device(ctx::mx.Context)
#  try
#    mx.ones(Float32, (1,1), ctx)
#    return true
#  catch err
#    return false
#  end
# end
#
#
#  best_device_cache = nothing
#
# export best_device
#  function best_device()
#    global best_device_cache
#    if best_device_cache == nothing
#      gpus = list_xpus(mx.gpu)
#      if length(gpus) > 0
#        best_device_cache = gpus
#        info("$(length(gpus)) GPUs found.")
#      else
#        best_device_cache = mx.cpu()
#        info("No GPUs available, fallback to single CPU.")
#      end
#    end
#    return best_device_cache
#  end
#
# export use_gpu
# function use_gpu(id)
#   global best_device_cache
#   if exists_device(mx.gpu(id))
#     best_device_cache = mx.gpu(id)
#   else
#     info("Could not access GPU $id, fallback to CPU")
#     best_device_cache = mx.cpu()
#   end
# end
#
# export list_xpus
# function list_xpus(xpu=mx.gpu)
#   result = mx.Context[]
#   while exists_device(xpu(length(result))) && length(result) < 8
#     push!(result, xpu(length(result)))
#   end
#   return result
# end

 type NetworkInfo
   name::String
   dir::AbstractString
   config::Dict
   context
   model
   epoch::Integer # the current state of the model, initialized to 0.
   training_curve::Vector{Float64} # MSE, created during training
   xval_curve::Vector{Float64} # MSE, created on demand

   NetworkInfo(name::String, dir::AbstractString, config::Dict, context) =
      new(name, dir, config, context, nothing, 0, Float64[], Float64[])
 end

Base.getindex(n::NetworkInfo, key::AbstractString) = n.config[key]
Base.setindex!(n::NetworkInfo, value, key::AbstractString) = n.config[key] = value

function save_compatible_heckpoint(sym :: mx.SymbolicNode, arg_params :: Dict{Base.Symbol, mx.NDArray}, aux_params :: Dict{Base.Symbol, mx.NDArray}, prefix :: AbstractString, epoch :: Int)
  if epoch <= 1
    mx.save("$prefix-symbol.json", sym)
  end
  save_dict = merge(Dict{Base.Symbol, mx.NDArray}(map((x) -> Symbol("arg:$(x[1])") => x[2], arg_params)),
                    Dict{Base.Symbol, mx.NDArray}(map((x) -> Symbol("aux:$(x[1])") => x[2], aux_params)))
  save_filename = format("{1}-{2:04d}.params", prefix, epoch)
  mx.save(save_filename, save_dict)
end

function mx_create_kvstore(kv_type :: Base.Symbol, num_device :: Int, arg_params :: Dict{Base.Symbol,mx.NDArray})
  if num_device == 1 && !ismatch(r"dist", string(kv_type))
    kv = nothing
  else
    if kv_type == :local
      max_size = maximum([prod(size(param)) for (k,param) in arg_params])
      if max_size < 1024 * 1024 * 16
        kv_type = :local_update_cpu
      else
        kv_type = :local_allreduce_cpu
      end
      info("Auto-select kvstore type = $kv_type")
    end
    kv = mx.KVStore(kv_type)
  end

  update_on_kvstore = true
  if isa(kv, Void) || ismatch(r"local_allreduce", string(mx.get_type(kv)))
    update_on_kvstore = false
  end

  return (kv, update_on_kvstore)
end

function train(n::NetworkInfo,
      train_provider, eval_provider;
      verbosity=2)
  learning_rate = n["learning_rate"]
  epochs = n["epochs"]
  optimizer_name = n["optimizer"]

  if epochs <= n.epoch
    verbosity >= 2 && info("$(n.name): Target epoch already reached.")
    return
  end

  training_curve = PlotCallback()
  eval_curve = Float64[]

  metric = mx.MSE()

  if optimizer_name == "SGD"
    optimizer = mx.SGD(lr=learning_rate, momentum=n["momentum"])
  elseif optimizer_name == "ADAM"
    optimizer = mx.ADAM(lr=learning_rate)
  elseif optimizer_name == "None"
    optimizer = mx.SGD(lr=1e-9)
  end
  # optimizer = mx.SGD(lr=learning_rate, momentum=0.9)

  # Manual initialization
  if !isdefined(n.model, :arg_params)
    mx.init_model(n.model, mx.UniformInitializer(0.1); overwrite=false, [mx.provide_data(train_provider)..., mx.provide_label(train_provider)...]...)
  end

  verbosity >= 2 && info("$(train_provider.sample_count) data points, $(get_total_parameter_count(n.model)) learnable parameters, dropout=$(n["dropout"]), batch size $(train_provider.batch_size), using $optimizer_name with learning rate $learning_rate\u1b[K")
  verbosity >= 2 && info("Starting training on $(n.model.ctx) (from $(n.epoch+1) to $epochs)... \u1b[K")
  verbosity >= 3 && info("Untrained MSE: $(eval(n.model, eval_provider, mx.MSE())[1][2])\u1b[K")
  flush(STDOUT)

  # for arg_param in n.model.arg_params
  #   print("$(arg_param[1]): ")
  #   println(copy(arg_param[2])[1:min(8, length(arg_param[2]))])
  # end
  # kvstore, update_on_kvstore = mx_create_kvstore(:device, length(n.model.ctx), n.model.arg_params)

  # Train each step manually
  for epoch in (n.epoch+1) : epochs
    mx.fit(n.model, optimizer, train_provider,
           n_epoch=1,
           eval_metric=metric,
           callbacks=[training_curve],
           verbosity=0,
           kvstore=:device)
    save_compatible_heckpoint(n.model.arch, n.model.arg_params, n.model.aux_params, joinpath(n.dir,n.name), epoch)

    eval_mse = eval(n.model, eval_provider, mx.MSE())[1][2]
    push!(eval_curve, eval_mse)
    verbosity >= 3 && info("Epoch $epoch / $epochs: MSE = $eval_mse.")
    verbosity >= 3 && flush(STDOUT)
  end

  # mx.fit(n.model, optimizer, train_provider, n_epoch=epochs, eval_metric=metric, callbacks=[training_curve], kvstore=:device) # TODO

  n.epoch = epochs

  calculate_parameters(n.model, joinpath(n.dir,n.name*"-parameters.txt"))

  append!(n.training_curve, training_curve.graphs[:MSE])
  append!(n.xval_curve, eval_curve)
  writedlm(joinpath(n.dir,n.name*"-MSE-train.txt"), n.training_curve)
  writedlm(joinpath(n.dir,n.name*"-MSE-xval.txt"), n.xval_curve)
end

export build
function build(n::NetworkInfo, method::Symbol,
    train_provider, eval_provider, build_function;
    verbosity=2
  )
  target_epoch = n["epochs"]
  slim = n["slim"]

  if(slim > 0 && method != :load)
    verbosity >= 2 && info("$(n.name): slim $(train_provider.sample_count) -> $slim")
    train_provider = slim_provider(train_provider, slim)
    if eval_provider != nothing
      eval_provider = slim_provider(eval_provider, slim)
    end
  end

  if method == :train
    loss, net = build_function(n.config, size(train_provider.data_arrays[1],1))
    n.model = mx.FeedForward(loss, context=n.context)
    train(n, train_provider, eval_provider; verbosity=verbosity)
    load_network(n, target_epoch)
  elseif method == :load
    load_network(n, target_epoch)
  elseif method == :refine
    load_network(n, -1; pick_best=false)
    train(n, train_provider, eval_provider; verbosity=verbosity)
    load_network(n, target_epoch)
  else
    throw(ArgumentError("method must be train, load or refine. got $method"))
  end
end

function slim_provider(p::mx.ArrayDataProvider, slim)
  slim = min(slim, p.sample_count)
  return mx.ArrayDataProvider(p.data_names[1] => slim_array(p.data_arrays[1], slim),
      p.label_names[1] => slim_array(p.label_arrays[1], slim); batch_size=p.batch_size)
end

function slim_array(array, slim)
  if(length(size(array)) == 1)
    return array[1:slim]
  else
    return array[:,1:slim]
  end
end

function padded_array_provider(key, data::Matrix{Float32}, batch_size)
  if size(data, 2) < batch_size
    plot_waveforms_padded = zeros(size(data,1), batch_size)
    plot_waveforms_padded[:,1:size(data, 2)] = data
  else
    plot_waveforms_padded = data
  end
  return mx.ArrayDataProvider(key => plot_waveforms_padded, batch_size=batch_size)
end

function eval(model, provider::mx.ArrayDataProvider, metric::mx.AbstractEvalMetric)
  prediction = mx.predict(model, provider; verbosity=0)
  data = provider.label_arrays[1]
  mx.reset!(metric)

  data_nd = mx.NDArray(data)
  prediction_nd = mx.NDArray(prediction)

  mx.update!(metric, [data_nd], [prediction_nd])
  return mx.get(metric)
end

export load_network
 function load_network(n::NetworkInfo, max_epoch; output_name="softmax", delete_unneeded_arguments=true, pick_best=true)
   if max_epoch < 0
     max_epoch = last_epoch(n.dir, n.name)
   end

  n.training_curve = readdlm(joinpath(n.dir, "$(n.name)-MSE-train.txt"))[:,1]
  n.xval_curve = readdlm(joinpath(n.dir, "$(n.name)-MSE-xval.txt"))[:,1]

  if pick_best
    epoch = findmin(n.xval_curve)[2]
    info("$(n.name): best epoch is $epoch.")
  else
    epoch = max_epoch
  end

   load_network_checkpoint(n, epoch; output_name=output_name, delete_unneeded_arguments=delete_unneeded_arguments)

   return n
 end

 function load_network_checkpoint(n::NetworkInfo, epoch; output_name="softmax", delete_unneeded_arguments=true)
   sym, arg_params, aux_params = mx.load_checkpoint(joinpath(n.dir, n.name), epoch)
   n.model = subnetwork(sym, arg_params, aux_params, output_name, delete_unneeded_arguments, n.context)
   n.epoch = epoch
 end

 function subnetwork(sym, arg_params, aux_params, output_name, delete_unneeded_arguments, context)
   all_layers = mx.get_internals(sym)
   loss = all_layers[output_name*"_output"]

   model = mx.FeedForward(loss, context=context)
   model.arg_params = copy(arg_params)
   model.aux_params = copy(aux_params)

   if delete_unneeded_arguments
     needed_args = mx.list_arguments(loss)
     for (name, array) in model.arg_params
       if !(name in needed_args)
         delete!(model.arg_params, name)
       end
     end
   end

   return model
 end

 function subnetwork(network::mx.FeedForward, subnetwork::mx.FeedForward)
   subnetwork.arg_params = copy(network.arg_params)
   subnetwork.aux_params = copy(network.aux_params)

   needed_args = mx.list_arguments(subnetwork.arch)
   for (name, array) in subnetwork.arg_params
     if !(name in needed_args)
       delete!(subnetwork.arg_params, name)
     end
   end

   return subnetwork
 end

 function last_epoch(dir, prefix; start=1)
   i = start
   if !isfile("$dir/$prefix-$(lpad(i,4,0)).params")
     throw(ArgumentError("Model not found: $dir"))
   end

   while isfile("$dir/$prefix-$(lpad(i,4,0)).params")
     i += 1
   end
   return i-1
 end

 function exists_network(dir::AbstractString, name::AbstractString)
   return isfile("$dir/$name-symbol.json")
 end

function decide_best_action(n::NetworkInfo)
  if !exists_network(n.dir, n.name) return :train end
  if n["epochs"] > last_epoch(n.dir, n.name) return :refine end
  return :load
end


text:
include("../benchmark.jl")

using JuMP
using Ipopt

function formulation1(lambda,w,mu,mu1,A,B,T, lb_var)
    m = Model(solver = IpoptSolver())
    @variable(m, c[i in A,j in B] >= lb_var)# for (i,j) in T)
    @variable(m, y[i in A,j in B] >= lb_var)# for (i,j) in T )

    @NLobjective(m, Min, -sum( lambda[i,j] * (log(c[i,j]) - (y[i,j]/w[i])^mu1[j] / mu1[j]) for (i,j) = T) )

    for (i,j) in T
      for (p,q) in T
        if i != p || j != q
          @NLconstraint(m,(log(c[i,j]) - (y[i,j]/w[i])^mu1[j] / mu1[j]) - (log(c[p,q]) - (y[p,q]/w[i])^mu1[j] / mu1[j]) >= 0.0)
        end
      end
    end

    @NLconstraint(m, sum(lambda[i,j] * (y[i,j] - c[i,j]) for (i,j) in T)  >= 0.0);
    return m
end

na = 21
nb = 5

A = 1:na;
B = 1:nb;
T = [(i,j) for i in A for j in B];

lambda = ones(na,nb);
wmin = 1;
#wmax = wmin + na - 1;
wmax = 5
w = [wmin + ((wmax-wmin)/(na-1))*(i-1) for i = 1:na]
mu = [1.0, 2.0, 3.0, 5.0, 8.0]
#mu = [1.0, 1.5, 2.0, 2.5, 3.0]

gamma = 0.0

mu1 = mu + 1.0

#lb = 0.0
lb_var = 0.1

srand(1)
m = formulation1(lambda,w,mu,mu1,A,B,T, lb_var)
iter = one_phase_solve(m);

#=
ORG_STDOUT = STDOUT
file = open("test.txt", "w")
redirect_stdout(file)

srand(1)
m = create_tax_problem(10)
one_phase_solve(m)

STDOUT = ORG_STDOUT

close(file)=#

ext: jl


text:
# User interaction with N Queens Problem.

# (c) E. Adrian Garro S. Costa Rica Institute of Technology.

include("nqueens.jl")

type UI
	option::String
	problem::NQueens
end

"""
`input(prompt::AbstractString="")`

Read a string from STDIN. The trailing newline is stripped.

The prompt string, if given, is printed to standard output without a
trailing newline before reading input.
"""
function input(prompt::AbstractString="")
	print(prompt)
	return chomp(readline())
end

"""
`int_input(prompt::AbstractString="")`

Use input function to get integers from user.
"""
function int_input(string_num::AbstractString)
	possible_int = input(string_num)
	parse(Int, possible_int)
end

function read_option(app::UI)
	println("")
	println("--------(N Queens Problem)---------")
	println("    	1) Input problem.")
	println("    	2) See solutions.")
	println("    	0) End.")
	println("----------------------------------")
	println("")
	app.option = input("Select an option to continue: ")
end

function input_problem(app::UI)
	try
		queens_quant = int_input(
			"Please write the number of queens: "
		)
		try
			set_queens_quant!(app.problem, queens_quant)
			solve!(app.problem)
			println("Problem solve! Please see the solutions :)")
		catch e
			println(e)
		end
	catch e
		println(e)
	end
end

function see_solutions(app::UI)
	solutions_num = app.problem.solutions_num
	if solutions_num != 0
		try
			println("The number of solutions for this problem are: ", solutions_num)
			solution_num = int_input(
				"Please write the number of the solution you want to see: "
			)
			try
				print_solution(app.problem, solution_num)
			catch e
				println(e)
			end
		catch e
			println(e)
		end
	else
		println("See solutions is not available.")
	end
end

function menu(app::UI)
	if app.option == ""
		true
	elseif app.option == "1"
		input_problem(app::UI)
		true
	elseif app.option == "2"
		see_solutions(app::UI)
		true
	elseif app.option == "0"
		println("Thank you for your time. ;)")
		false
	else
		println("Pick a valid option from menu!")
		true
	end
end

function main()
	app = UI("", NQueens(0, [], 0))
	while menu(app::UI)
		read_option(app::UI)
	end
end

main()

max_stars_repo_name: adriangarro/nqueens

max_issues_repo_name: adriangarro/nqueens

max_forks_repo_name: adriangarro/nqueens


text:
using Flux, Gym, Printf, Zygote
using Flux.Tracker: data
using Flux.Optimise: update!
using Statistics: mean
using DataStructures: CircularBuffer
using Distributions: sample

#using CuArrays

#Load game environment

env = make("Pendulum-v0")
reset!(env)

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(state(env))
ACTION_SIZE = 1#length(env.actions)
ACTION_BOUND = env._env.action_space.high[1]
MAX_EP = 50_000
MAX_EP_LENGTH = 200

BATCH_SIZE = 64
MEM_SIZE = 100_000
MIN_EXP_SIZE = 50_000

γ = 99f-2     # discount rate

τ = 1f-3 # for running average while updating target networks
η_act = 1f-4   # Learning rate
η_crit = 1f-3
L2_DECAY = 1f-2

# Ornstein-Uhlenbeck Noise params
μ = 0f0
θ = 15f-2
σ = 2f-1

# --------------------------------- Memory ------------------------------------

memory = CircularBuffer{Any}(MEM_SIZE)

function getData(batch_size = BATCH_SIZE)
  # Getting data in shape
  minibatch = sample(memory, batch_size)
  x = hcat(minibatch...)

  s      =   hcat(x[1, :]...) |> gpu
  a      =   hcat(x[2, :]...) |> gpu
  r      =   hcat(x[3, :]...) |> gpu
  s′     =   hcat(x[4, :]...) |> gpu
  s_mask = .!hcat(x[5, :]...) |> gpu

  return s, a, r, s′, s_mask
end

# -------------------------------- Action Noise --------------------------------

struct OUNoise
  μ
  θ
  σ
  X
end

ou = OUNoise(μ, θ, σ, zeros(Float32, ACTION_SIZE) |> gpu)

function sample_noise(ou::OUNoise)
  dx     = ou.θ * (ou.μ .- ou.X)
  dx   .+= ou.σ * randn(Float32, length(ou.X)) |> gpu
  ou.X .+= dx
end

# Noise scale
τ_ = 25
ϵ  = exp(-1f0 / τ_)
noise_scale = 1f0 / ACTION_BOUND

# ----------------------------- Model Architecture -----------------------------

w_init(dims...) = 6f-3rand(Float32, dims...) .- 3f-3

actor = Chain(Dense(STATE_SIZE, 400, relu),
	      	  Dense(400, 300, relu),
              Dense(300, ACTION_SIZE, tanh, initW=w_init),
              x -> x * ACTION_BOUND) |> gpu
actor_target = deepcopy(actor)

# Critic model
struct crit
  state_crit
  act_crit
  sa_crit
end

Flux.@treelike crit

function (c::crit)(state, action)
  s = c.state_crit(state)
  a = c.act_crit(action)
  c.sa_crit(relu.(s .+ a))
end

Base.deepcopy(c::crit) = crit(deepcopy(c.state_crit),
                              deepcopy(c.act_crit),
			      			  deepcopy(c.sa_crit))

critic = crit(Chain(Dense(STATE_SIZE, 400, relu), Dense(400, 300)) |> gpu,
              		Dense(ACTION_SIZE, 300) |> gpu,
	      			Dense(300, 1, initW=w_init) |> gpu)
critic_target = deepcopy(critic)

# ---------------------- Param Update Functions --------------------------------

function update_target!(target, model; τ = 1f0)
  for (p_t, p_m) in zip(params(target), params(model))
    p_t.data .= (1f0 - τ) * p_t.data .+ τ * p_m.data
  end
end

function update_model!(model, opt, loss, inp...)
  grads = gradient(()->loss(inp...), params(model))
  update!(opt, params(model), grads)
end

# ---------------------------------- Training ----------------------------------
## Losses
function L2_loss(model)
  l2_loss = sum(map(p->sum(p.^2), params(model)))
  return L2_DECAY * l2_loss
end

loss_crit(y, s, a) = Flux.mse(critic(s, a), y) #+ L2_loss(critic)

function loss_act(s)
  actions = actor(s)
  crit_out = critic(s, actions)
  return -sum(crit_out)
end

## Optimizers
opt_crit = ADAM(η_crit)
opt_act  = ADAM(η_act)

function replay()
  s, a, r, s′, s_mask = getData()

  a′ = actor_target(s′)
  v′ = critic_target(s′, a′)
  y = data(r .+ γ * v′ .* s_mask)	# set v′ to 0 where s_ is terminal state

  update_model!(critic, opt_crit, loss_crit, y, s, a)
  update_model!(actor, opt_act, loss_act, s)

  # Update Target models
  update_target!(actor_target, actor; τ = τ)
  update_target!(critic_target, critic; τ = τ)
end

# ---------------------------- Helper Functions --------------------------------

# Stores tuple of state, action, reward, next_state, and done
remember(state, action, reward, next_state, done) =
  push!(memory, [data.((state, action, reward, next_state))..., done])

# Choose action according to policy PendulumPolicy
function action(state, train=true)
  state = reshape(state, size(state)..., 1)
  act_pred = actor(state |> gpu)
  if train
    act_pred = act_pred .+ noise_scale * sample_noise(ou)
  end
  clamp.(act_pred, -ACTION_BOUND, ACTION_BOUND) # returns action
end

function episode!(env::EnvWrapper)
  reset!(env)
  for ep=1:MAX_EP_LENGTH
	s = state(env)
    a = action(s, trainable(env))
    s′, r, done, _ = step!(env, a)
    if trainable(env)
      remember(s, a, r, s′, done)
      replay()
    end
  end
  env.total_reward
end

# -------------------------------- Testing -------------------------------------

# Returns average score over 100 episodes

function test(env::EnvWrapper)
  score_mean = 0f0
  testmode!(env)
  for e=1:100
    total_reward = episode!(env)
    score_mean += total_reward / 100
  end
  testmode!(env, false)
  return score_mean
end

# ------------------------------ Training --------------------------------------

# Populate memory with random actions

s = reset!(env)
for e=1:MIN_EXP_SIZE
  global s
  a = 2rand(Float32) * ACTION_BOUND - ACTION_BOUND
  s′, r, done, _ = step!(env, a)
  remember(s, a, r, s′, done)
  s = s′
end

for e=1:MAX_EP
  total_reward = episode!(env)
  total_reward = @sprintf "%9.3f" total_reward
  print("Episode: $e | Score: $total_reward | ")
  score_mean = test(env)
  score_mean = @sprintf "%9.3f" score_mean
  println("Mean score over 100 test episodes: $score_mean")
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
beucher_gradient_ref(img, dims, r) = beucher_gradient_ref(img, strel_box(img, dims; r))
beucher_gradient_ref(img, se) = float.(dilate(img, se)) .- float.(erode(img, se))
external_gradient_ref(img, dims, r) = external_gradient_ref(img, strel_box(img, dims; r))
external_gradient_ref(img, se) = float.(dilate(img, se)) .- img
internal_gradient_ref(img, dims, r) = internal_gradient_ref(img, strel_box(img, dims; r))
internal_gradient_ref(img, se) = img .- float.(erode(img, se))

@testset "mgradient" begin
    test_ranges = [Bool, 1:10, Float64, Gray{N0f8}, Gray{Float64}]
    modes = [
        (:beucher, beucher_gradient_ref),
        (:external, external_gradient_ref),
        (:internal, internal_gradient_ref),
    ]
    for (mode, ref_fun) in modes
        for T in test_ranges
            for N in (1, 2, 3)
                sz = ntuple(_ -> 32, N)
                img = rand(T, sz...)

                out = mgradient(img; mode)
                @test out == mgradient(img, strel_box(img, ntuple(identity, N); r=1); mode)
                @test out == ref_fun(img, ntuple(identity, N), 1)

                @test mgradient(img; dims=(1,), r=2, mode) == ref_fun(img, (1,), 2)

                se = rand_se_mask(3, N; symmetric=true)
                @test mgradient(img, se; mode) == ref_fun(img, se)
            end
        end
    end

    img = rand(1:5, 7, 7)
    out = similar(img)
    mgradient!(out, img, similar(img))
    @test out == mgradient(img)

    img = rand(RGB, 7, 7)
    @test_throws ArgumentError("color image is not supported") mgradient(img)

    img = rand(1:10, 7, 7)
    se = Bool[1 0 0; 0 1 0; 0 0 0]
    msg = "structuring element must be symmetric with respect to its center"
    @test_throws ArgumentError(msg) mlaplacian(img, se)
end


text:
# Tracing charged particle in the magnetic bottle.
# This example shows how to trace protons in a stationary magnetic field that
# corresponds to the magnetic bottle.
# Reference: https://en.wikipedia.org/wiki/Magnetic_mirror#Magnetic_bottles
#
# Hongyang Zhou, hyzhou@umich.edu

using TestParticle
using TestParticle: getB_bottle
using OrdinaryDiffEq
using StaticArrays
using PyPlot
using Statistics: mean
using Printf

## Obtain field

# Magnetic bottle parameters in SI units
const I1 = 80. # current in the solenoid
const N1 = 1500 # number of windings
const I2 = 100. # current in the central solenoid
const N2 = 8000 # number of windings
const distance = 10. # distance between solenoids
const a = 1.0 # radius of each coil
const b = 4.0 # radius of central coil

function getB(xu)
   SVector{3}(getB_Bottle(xu[1], xu[2], xu[3], distance, a, b, I1*N1, I2*N2))
end

function getE(xu)
   SA[0.0, 0.0, 0.0]
end

## Initialize particles
m = TestParticle.mₑ
q = TestParticle.qₑ
c = TestParticle.c

# initial velocity, [m/s]
v₀ = [0.8, 0.48, 0.3595] .* c
# initial position, [m]
r₀ = [0.0, 0.0, -0.8]
stateinit = [r₀..., v₀...]

# Theoretically we can take advantage of the fact that magnetic field does not
# accelerate particles, so that γ remains constant. However, we are not doing
# that here since it is not generally true in the EM field.

param = prepare(getE, getB; species=Electron)
tspan = (0.0, 3e-7)

prob = ODEProblem(trace_analytic_relativistic!, stateinit, tspan, param)

@printf "Speed = %6.4f %s\n" √(v₀[1]^2+v₀[2]^2+v₀[3]^2)/c*100 "% speed of light"
@printf "Energy = %6.4f MeV\n" (1/√(1-(v₀[1]/c)^2-(v₀[2]/c)^2-(v₀[3]/c)^2)-1)*m*c^2/abs(q)/1e6

# Default Tsit5() and many solvers does not work in this case!
sol = solve(prob, AB3(); dt=1.5e-12, save_idxs=[1,2,3])

## Visualization

using3D()
fig = plt.figure(figsize=(10,6))
ax = fig.gca(projection="3d")

ax.plot(sol[1,:], sol[2,:], sol[3,:], label="electron")

# Plot coils
θ = range(0, 2π, length=100)
x = a.*cos.(θ)
y = a.*sin.(θ)
z = fill(distance/2, size(x))
ax.plot(x, y, z, color="r", label="coil")

z = fill(-distance/2, size(x))
ax.plot(x, y, z, color="r")

x = b.*cos.(θ)
y = b.*sin.(θ)
z = fill(0.0, size(x))
ax.plot(x, y, z, color="r")

ax.legend()
title("Electron trajectory in the magnetic bottle")
xlabel("x [m]")
ylabel("y [m]")
zlabel("z [m]")

xMin, xMax = -5, 5
yMin, yMax = -5, 5
zMin, zMax = -5.2, 5.2

xlim(xMin, xMax)
ylim(yMin, yMax)
zlim(zMin, zMax)

TestParticle.set_axes_equal(ax)

#savefig("test.png", bbox_inches="tight")


text:
const STILL = 0
const NORTH = 1
const EAST = 2
const SOUTH = 3
const WEST = 4

const DIRECTIONS = [a for a in 0:4]
const CARDINALS = [a for a in 1:4]

const ATTACK = 0
const STOP_ATTACK = 1

type Location
  x :: Int64
  y :: Int64
  Location(x::Int64, y::Int64) = new(x, y)
  Location() = Location(0, 0)
end

type Site
  owner :: Int64
  strength :: Int64
  production :: Int64
  Site(owner::Int64, strength::Int64, production::Int64) = new(owner, strength, production)
  Site() = Site(0, 0, 0)
end

type Move
  loc :: Location
  direction :: Int64
  Move(loc::Location, direction::Int64) = new(loc, direction)
  Move() = Move(Location(), 0)
end

type GameMap
  width :: Int64
  height :: Int64
  contents :: Vector{Vector{Site}}
  function GameMap(width::Int64, height::Int64)
    contents = Vector{Vector{Site}}()
    for y in 0:height-1
      row = Vector{Site}()
      for x in 0:width-1
        push!(row, Site(0, 0, 0))
      end
      push!(contents, row)
    end
    new(width, height, contents)
  end
  GameMap() = GameMap(0, 0)
end

isBounds(gm::GameMap, l::Location) = l.x >= 0 && l.x < gm.width && l.y >= 0 && l.y < gm.height

function getDistance(gm::GameMap, l1::Location, l2::Location)
  dx = abs(l1.x - l2.x)
  dy = abs(l1.y - l2.y)

  if dx > gm.width / 2
    dx = gm.width - dx
  end

  if dy > gm.height / 2
    dy = gm.height - dy
  end

  dx + dy
end

function getAngle(gm::GameMap, l1::Location, l2::Location)
  dx = l2.x - l1.x
  dy = l2.y - l1.y

  if dx > gm.width - dx
    dx -= gm.width
  elseif -dx > gm.width + dx
    dx += gm.width
  end

  if dy > gm.height - dy
    dy -= gm.height
  elseif -dy > gm.height + dy
    dy += gm.height
  end

  atan2(dy, dx)
end

function getLocation(gm::GameMap, loc::Location, direction::Int64)
  l = deepcopy(loc)
  if direction != STILL
    if direction == NORTH
      if l.y == 0
        l.y = gm.height - 1
      else
        l.y = l.y - 1
      end
    elseif direction == EAST
      if l.x == gm.width - 1
        l.x = 0
      else
        l.x = l.x + 1
      end
    elseif direction == SOUTH
      if l.y == gm.height - 1
        l.y = 0
      else
        l.y = l.y + 1
      end
    elseif direction == WEST
      if l.x == 0
        l.x = gm.width - 1
      else
        l.x = l.x - 1
      end
    end
  end
  l
end

function getSite(gm::GameMap, l::Location, direction::Int64=STILL)
  l = getLocation(gm, l, direction)
  gm.contents[l.y + 1][l.x + 1]
end

max_stars_repo_name: yangle/HaliteIO

max_issues_repo_name: yangle/HaliteIO

max_forks_repo_name: yangle/HaliteIO


text:
module Wynn
    using SymPy

    # data structure to store the epsilon table
    # ϵ_ij indexed by ints i & j, stored within dict.
    struct EpsilonTable{T}
        series::T
        terms::Vector{T}
        etable::Dict{Tuple{Int,Int},T}
    end

    function EpsilonTable(terms::Vector{T}; simplified::Bool = true) where T<:Union{Float64,Sym}
        # maximum occuring i & j index
        max_ind = length(terms)

        # setting base case j = -1
        etable = Dict((i, -1) => eltype(terms)(0) for i in 0:max_ind)

        # setting base case j = 0, i = 0:max_ind-1
        merge!(etable, Dict((i-1, 0) => sum(terms[1:i]) for i in 1:max_ind))

        # setting base case i = -j-1, j even
        merge!(etable, Dict((-j-1, 2j) => eltype(terms)(0) for j in 0:max_ind))

        # setting base case for odd j
        merge!(etable, Dict((-j-1, 2j-1) => eltype(terms)(0) for j in 0:max_ind))

        # recursive calculations, j>=1, i = floor(Int,-j/2):(max_ind-j-1)
        for j in 1:2*(max_ind-1), i in floor(Int,-j/2):(max_ind-j-1)#j-1
            ϵ_ij = etable[i+1, j-2] + 1 / (etable[i+1, j-1] - etable[i, j-1])
            eltype(terms) <: Sym && simplified ? push!(etable, (i, j) =>
                simplify(ϵ_ij)) : push!(etable, (i, j) => ϵ_ij)
        end
        series = sum(terms)
        EpsilonTable(series, terms, etable)
    end

    EpsilonTable(terms::Vector{Int}; simplified::Bool = true) = EpsilonTable(convert(Vector{Float64}, terms); simplified = simplified)

    export EpsilonTable
end


text:
# Functionality for the promotion of dictionaries.

export promote_domaintype,
    promote_coefficienttype

# Convert the domaintype of the dictionary to `S`
convert(::Type{Dictionary{S}}, dict::Dictionary{S}) where {S} = dict
convert(::Type{Dictionary{S}}, dict::Dictionary) where {S} = similar(dict, S)

promote(d1::Dictionary, d2::Dictionary) =
    promote_coefficienttype(promote_domaintype(d1, d2)...)

## Domain type

"Promote the domain type of the dictionaries to a common type."
promote_domaintype(dict::Dictionary) = (dict,)

function promote_domaintype(d1::Dictionary, d2::Dictionary, d3::Dictionary...)
    d1b, d2b = promote_domaintype(d1, d2)
    (d1b, promote_domaintype(d2b, d3...))
end

promote_domaintype(d1::Dictionary, d2::Dictionary) =
    _promote_domaintype(d1, d2, promote_type(domaintype(d1), domaintype(d2)))
_promote_domaintype(d1, d2, ::Type{S}) where {S} =
    convert(Dictionary{S}, d1), convert(Dictionary{S}, d2)

## Coefficient type

"Promote the coefficient type of the dictionaries to a common type."
promote_coefficienttype(dict::Dictionary) = (dict,)

function promote_coefficienttype(d1::Dictionary, d2::Dictionary, d3::Dictionary...)
    d1b, d2b = promote_coefficienttype(d1, d2)
    (d1b, promote_coefficienttype(d2b, d3...))
end

promote_coefficienttype(d1::Dictionary, d2::Dictionary) =
    _promote_coefficienttype(d1, d2, coefficienttype(d1), coefficienttype(d2))

# Types are the same
_promote_coefficienttype(d1, d2, ::Type{T}, ::Type{T}) where {T} = (d1,d2)
# Types differ: compute a joined supertype
_promote_coefficienttype(d1, d2, ::Type{S}, ::Type{T}) where {S,T} =
    _promote_coefficienttype(d1, d2, S, T, promote_type(S,T))

_promote_coefficienttype(d1, d2, ::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U} =
    ensure_coefficienttype(U, d1), ensure_coefficienttype(U, d2)

"Ensure that the dictionary (or dictionaries) have coefficient type `T` or wider."
ensure_coefficienttype(::Type{C}, dict::Dictionary) where {C} =
    _ensure_coefficienttype(C, ensure_prectype(prectype(C),dict))
_ensure_coefficienttype(::Type{C}, dict::Dictionary) where {C} =
    _ensure_coefficienttype(dict, C, coefficienttype(dict))
_ensure_coefficienttype(dict::Dictionary, ::Type{T}, ::Type{T}) where {T} = dict
# No fallback is provided, in order to generate an error when a coefficient type
# is requested that can not be achieved. Complexification of a dictionary is
# achieved by ComplexifiedDict.

ensure_coefficienttype(::Type{T}, dicts::Dictionary...) where {T} =
    ensure_coefficienttype.(T, dicts)

## Precision type

"Ensure that the dictionary has precision type `T` or a wider type."
ensure_prectype(::Type{T}, dict::Dictionary) where {T} =
    _ensure_prectype(dict, T, prectype(dict))
_ensure_prectype(dict::Dictionary, ::Type{T}, ::Type{T}) where {T} = dict
_ensure_prectype(dict::Dictionary, ::Type{S}, ::Type{T}) where {S,T} =
    _ensure_prectype(dict, S, T, promote_type(S,T))

_ensure_prectype(dict::Dictionary, S, T, ::Type{U}) where {U} =
    similar(dict, widen_prectype(domaintype(dict), U))

ensure_prectype(::Type{T}, dicts::Dictionary...) where {T} =
    ensure_prectype.(T, dicts)

"Widen the type such that its precision type matches `U`."
widen_prectype(::Type{Complex{T}}, ::Type{U}) where {T,U} = Complex{U}
widen_prectype(::Type{<:Number}, ::Type{U}) where {U<:Number} = U
widen_prectype(::Type{SVector{T,N}}, ::Type{U}) where {T,N,U} = SVector{widen_prectype(T,U),N}
widen_prectype(::Type{Array{T,N}}, ::Type{U}) where {T,N,U} = Array{widen_prectype(T,U),N}
widen_prectype(::Type{NTuple{N,T}}, ::Type{U}) where {T,N,U} = NTuple{N,widen_prectype(T,U)}
widen_prectype(::Type{Tuple{A}}, ::Type{U}) where {A,U} = Tuple{widen_prectype(A,U)}
widen_prectype(::Type{Tuple{A,B}}, ::Type{U}) where {A,B,U} =
    Tuple{widen_prectype(A,U),widen_prectype(B,U)}
widen_prectype(::Type{Tuple{A,B,C}}, ::Type{U}) where {A,B,C,U} =
    Tuple{widen_prectype(A,U),widen_prectype(B,U),widen_prectype(C,U)}


text:
@testset "Brownian Bridge" begin

using DiffEqNoiseProcess, DiffEqBase, Test, Random, DiffEqBase.EnsembleAnalysis

Random.seed!(100)
W = BrownianBridge(0.0,1.0,0.0,1.0,0.0,0.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
@time sol = solve(ensemble_prob,dt=0.1,trajectories=100000)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 2:10
  q = qs[i]
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end
@test ≈(timestep_mean(sol,1)[1],0.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,1)[2],0.0,atol=1e-16)
@test ≈(timestep_mean(sol,11)[1],1.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,11)[2],0.0,atol=1e-16)

μ = 1.2
σ = 2.2
W = GeometricBrownianBridge(μ,σ,0.0,1.0,0.0,1.0,0.0,0.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
@time sol = solve(ensemble_prob,dt=0.1,trajectories=100)

Random.seed!(100)
r = 100 # should be independent of the rate, so make it crazy
rate(u,p,t) = r
W = CompoundPoissonBridge(rate,0.0,1.0,0.0,1.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
@time sol = solve(ensemble_prob,dt=0.1,trajectories=100000)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 2:10
  q = qs[i]
  # Mean and variance of binomial matches that of the Brownian bridge!
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end
@test ≈(timestep_mean(sol,1)[1],0.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,1)[2],0.0,atol=1e-16)
@test ≈(timestep_mean(sol,11)[1],1.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,11)[2],0.0,atol=1e-16)

# check VBT distributional properties

W = VirtualBrownianTree(0.0,0.0; Wend=1.0,tree_depth=3)
prob = NoiseProblem(W,(0.0,1.0))
function prob_func(prob,i,repeat)
  # to re-instantiate PRNG
  Wtmp = VirtualBrownianTree(0.0,0.0; Wend=1.0,tree_depth=3)
  remake(prob, noise=Wtmp)
end
ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
@time sol = solve(ensemble_prob,dt=0.125,trajectories=100000)

# Spot check the mean and the variance
qs = 0:0.125:1
for i in 2:8
  q = qs[i]
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end
@test ≈(timestep_mean(sol,1)[1],0.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,1)[2],0.0,atol=1e-16)
@test ≈(timestep_mean(sol,Int(2^(W.tree_depth)+1))[1],W.W[end],atol=1e-16)
@test ≈(timestep_meanvar(sol,Int(2^(W.tree_depth)+1))[2],0.0,atol=1e-16)

end


text:
struct NNParamKolmogorov{C,O,S,E} <: NeuralPDEAlgorithm
    chain::C
    opt::O
    sdealg::S
    ensemblealg::E
end
NNParamKolmogorov(chain  ; opt=Flux.ADAM(0.1) , sdealg = EM() , ensemblealg = EnsembleThreads()) = NNParamKolmogorov(chain , opt , sdealg , ensemblealg)

function DiffEqBase.solve(
    prob::ParamKolmogorovPDEProblem,
    alg::NNParamKolmogorov;
    abstol = 1f-6,
    verbose = false,
    maxiters = 300,
    trajectories = 1000,
    save_everystep = false,
    use_gpu = false,
    dγ = 0.01,
    dt,
    dx,
    kwargs...
    )

    tspan = prob.tspan
    g = prob.g
    f = prob.f
    noise_rate_prototype = prob.noise_rate_prototype
    Y_domain = prob.Y_domain
    xspan = prob.xspan
    d = prob.d
    phi = prob.phi
    ts = tspan[1]:dt:tspan[2]
    xs = xspan[1]:dx:xspan[2]
    γs_sigma = Y_domain.sigma[1]:dγ:Y_domain.sigma[2]
    γs_mu = Y_domain.mu[1]:dγ:Y_domain.mu[2]
    γs_phi = Y_domain.phi[1]:dγ:Y_domain.phi[2]
    γ_mu_prototype = prob.Y_mu_prototype
    γ_sigma_prototype = prob.Y_sigma_prototype
    γ_phi_prototype = prob.Y_phi_prototype

    #hidden layer
    chain  = alg.chain
    opt    = alg.opt
    sdealg = alg.sdealg
    ensemblealg = alg.ensemblealg
    ps     = Flux.params(chain)
    t = rand(ts , 1 , trajectories)
    x = rand(xs , d , 1 , trajectories)

    #Sample all the parameters from uniform distributions
    total_dims = d + 1
    if isnothing(γ_sigma_prototype)
      γ_sigma = nothing
    else
      γ_sigma = rand(γs_sigma , size(γ_sigma_prototype)[1] , size(γ_sigma_prototype)[2] , size(γ_sigma_prototype)[3] , trajectories)
      total_dims += d^2*(size(γ_sigma_prototype))[3]
    end
    if isnothing(γ_mu_prototype)
      γ_mu_1 = nothing
      γ_mu_2 = nothing
    else
      if !isnothing(γ_mu_prototype[1])
        γ_mu_1 = rand(γs_mu , size(γ_mu_prototype[1])[1], size(γ_mu_prototype[1])[2], 1 , trajectories)
        total_dims += d^2
      elseif !isnothing(γ_mu_prototype[1])
        γ_mu_2 = rand(γs_mu , size(γ_mu_prototype[2])[1], size(γ_mu_prototype[2])[2] , 1 ,trajectories)
        total_dims += d
      end
    end
    if isnothing(γ_phi_prototype)
      γ_phi = nothing
    else
      γ_phi = rand(γs_phi , size(γ_phi_prototype)[1] , 1 , trajectories)
      total_dims += size(γ_phi_prototype)[1]
    end

    #Declare wrapper functions to mu and sigma
    function sigma_(dx , x , γ_sigma , i)
      if isnothing(γ_sigma)
        return g(x , γ_sigma)
      else
        return g(x , γ_sigma[: , : , : , i])
      end
    end

    function mu_(dx, x , γ_mu_1, γ_mu_2 , i)
      if isnothing(γ_mu_prototype)
        return f(x , γ_mu_1 , γ_mu_2)
      else
        if isnothing(γ_mu_1) && isnothing(γ_mu_2)
          return f(x , γ_mu_1 , γ_mu_2)
        elseif isnothing(γ_mu_1) && !isnothing(γ_mu_2)
          return f(x , γ_mu_1 , γ_mu_2[: , : , : , i])
        elseif !isnothing(γ_mu_1) && isnothing(γ_mu_2)
          return f(x , γ_mu_1[: , : , : , i] , γ_mu_2)
        else
          return f(x , γ_mu_1[: , : , : , i] , γ_mu_2[: , : , : , i])
        end
      end
    end

    #Preparing training data
    X  = zeros(total_dims ,trajectories)
    for i in 1:length(t)
      K = vcat(t[i] , x[: , : , i])
      if !isnothing(γ_sigma_prototype)
        K = vcat(K , reshape(γ_sigma[: , : , : , i] , d^2*(size(γ_sigma_prototype))[3] , 1))
      end
      if !isnothing(γ_mu_prototype)
        K = vcat(K , reshape(γ_mu_1[: , : , : , i] , d^2 , 1) , reshape(γ_mu_2[: , : , : , i]  , d , 1))
      end
      if !isnothing(γ_phi_prototype)
        K = vcat(K , reshape(γ_phi[: , :  , i] , size(γ_phi_prototype)[1] , 1))
      end
      X[: , i] = K
    end
    print("Data Prepared to train the model")

    function sigma(dx , x , p , t)
      sigma_(dx , x , γ_sigma , 1)
    end
    function mu(dx, x , p , t)
      mu_(dx, x , γ_mu_1, γ_mu_2 , 1)
    end

    sdeproblem = SDEProblem(mu,sigma,x[: , : , 1], (0.00 , t[1]))
    function prob_func(prob,i,repeat)
      sigma(dx , x , p , t) = sigma_(dx , x , γ_sigma , i)
      mu(dx , x , p , t) = mu_(dx , x , γ_mu_1, γ_mu_2 , i)
      SDEProblem(prob.f , prob.g , x[: , : , i] , (0.00 , t[i]) ,noise_rate_prototype = prob.noise_rate_prototype)
    end
    output_func(sol,i) = (sol[end],false)
    ensembleprob = EnsembleProblem(sdeproblem , prob_func = prob_func , output_func = output_func)
    sim = solve(ensembleprob, sdealg, ensemblealg , dt=dt, trajectories=trajectories,adaptive=false)
    x_sde  = reshape([],d,0)
    # sol = solve(sdeproblem, sdealg ,dt=0.01 , save_everystep=false , kwargs...)
    # x_sde = sol[end]

    x_sde  = zeros(d,trajectories)
    for i in 1:length(sim.u)
        x_sde[: , i] = sim.u[i]
    end

    y = phi(x_sde , γ_phi)

    #Check the dimensions for the target
    if size(y)[1] != 1 && size(y)[1] != trajectories
      errMsg = "The output from phi function has dimensions "*string(size(y))*" while it should be (1, trajectories)"
      throw(DimensionMismatch(errMsg))
    end

    if use_gpu == true
        y = y |> gpu
        X = X |> gpu
    end
    data   = Iterators.repeated((X , y), maxiters)
    if use_gpu == true
        data = data |> gpu
    end

    #MSE Loss Function
    loss(x , y) =Flux.mse(chain(x) , y)

    callback = function ()
        l = loss(X, y)
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end

    Flux.train!(loss, ps, data, opt; cb = callback)

 end #solve


text:
#==========================================================================================#

# PARSE WEIGHTS

function parse_weights(nonmissing::BitVector, w::AnalyticWeights)
    n  = sum(nonmissing)
    v  = float(w[nonmissing])
    s  = n / sum(v)
    v .= v .* s
    AnalyticWeights(v, n)
end

function parse_weights(nonmissing::BitVector, w::FrequencyWeights)
    v = w[nonmissing]
    s = Int(sum(v))
    FrequencyWeights(float(v), s)
end

function parse_weights(nonmissing::BitVector, w::ProbabilityWeights)
    n  = sum(nonmissing)
    v  = float(w[nonmissing])
    s  = n / sum(v)
    v .= v .* s
    ProbabilityWeights(v, n)
end

parse_weights(nonmissing::BitVector, w::UnitWeights) = UnitWeights(Float64, sum(nonmissing))

function parse_weights(nonmissing::BitVector, w::Weights)
    n  = sum(nonmissing)
    v  = float(w[nonmissing])
    s  = n / sum(v)
    v .= v .* s
    Weights(v, n)
end

#==========================================================================================#

# INTERSECTION BETWEEN A FRAME AND CORRELATION STRUCTURE

parse_corr!(nonmissing::BitVector, corr::Homoscedastic)   = copy(corr)
parse_corr!(nonmissing::BitVector, corr::Heteroscedastic) = copy(corr)

function parse_corr!(nonmissing::BitVector, corr::Clustered)

    (nonmissing == corr.nonmissing) && (return copy(corr))

    nonmissing .*= corr.nonmissing
    touse        = nonmissing[corr.nonmissing]
    new_ic       = corr.ic[touse]
    iter         = sort(unique(corr.ic))
    new_iter     = sort(unique(new_ic))
    new_nc       = length(new_iter)
    new_mat      = corr.mat[findall((in)(new_iter), iter), touse]

    return Clustered(corr.adj, nonmissing, new_mat, new_ic, new_nc)
end

function parse_corr!(nonmissing::BitVector, corr::CrossCorrelated)

    (nonmissing == corr.nonmissing) && (return copy(corr))

    nonmissing .*= corr.nonmissing
    touse        = nonmissing[corr.nonmissing]
    new_mat      = Symmetric(corr.mat[touse, touse])

    return CrossCorrelated(corr.adj, nonmissing, SuiteSparse.CHOLMOD.Sparse(new_mat))
end

#==========================================================================================#

# RETRIEVAL

function getvector(obj::Microdata, x::Symbol)
    (length(obj.mapping[x]) == 1) || throw(ArgumentError(string(x) * " is not a vector"))
    return view(obj.matrix, :, obj.mapping[x]...)
end

function getmatrix(obj::Microdata, args...)
    x = mapreduce(i -> obj.mapping[i], vcat, args, init = Vector{Int64}())
    return view(obj.matrix, :, x)
end

function getnames(obj::Microdata, args...)
    n = coefnames(eterms(obj.model))
    x = mapreduce(i -> obj.mapping[i], vcat, args, init = Vector{Int64}())
    return n[x]
end

getcorr(obj::Microdata)       = obj.corr
getnonmissing(obj::Microdata) = obj.nonmissing
getweights(obj::Microdata)    = obj.weights


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Printf

# this macro tests for exceptions thrown at macro expansion
macro test_me(ty, ex)
    return quote
        @test_throws $(esc(ty)) try
            $(esc(ex))
        catch err
            @test err isa LoadError
            @test err.file === $(string(__source__.file))
            @test err.line === $(__source__.line)
            rethrow(err.error)
        end
    end
end

# printf
# int
@test (@sprintf "%d" typemax(Int64)) == "9223372036854775807"
@test (@sprintf "%a" typemax(Int64)) == "0x7.fffffffffffffffp+60"
@test (@sprintf "%A" typemax(Int64)) == "0X7.FFFFFFFFFFFFFFFP+60"

#printing an int value
for (fmt, val) in (("%i", "42"),
                   ("%u", "42"),
                   ("Test: %i", "Test: 42"),
                   ("%#x", "0x2a"),
                   ("%#o", "052"),
                   ("%x", "2a"),
                   ("%X", "2A"),
                   ("% i", " 42"),
                   ("%+i", "+42"),
                   ("%4i", "  42"),
                   ("%-4i", "42  "),
                   ("%a", "0x2.ap+4"),
                   ("%A", "0X2.AP+4"),
                   ("%20a","            0x2.ap+4"),
                   ("%-20a","0x2.ap+4            "),
                   ("%f", "42.000000"),
                   ("%g", "42")),
     num in (UInt16(42), UInt32(42), UInt64(42), UInt128(42),
              Int16(42), Int32(42), Int64(42), Int128(42), big"42")
            #big"42" causes stack overflow on %a ; gh #14409
    num isa BigInt && fmt in ["%a", "%#o", "%g"] && continue
    @test @eval(@sprintf($fmt, $num) == $val)
end

# pointers
if Sys.WORD_SIZE == 64
    @test (@sprintf "%20p" 0) == "  0x0000000000000000"
    @test (@sprintf "%-20p" 0) == "0x0000000000000000  "
elseif Sys.WORD_SIZE == 32
    @test (@sprintf "%20p" 0) == "          0x00000000"
    @test (@sprintf "%-20p" 0) == "0x00000000          "
else
    @test false
end

# float / BigFloat
for (fmt, val) in (("%7.2f", "   1.23"),
                   ("%-7.2f", "1.23   "),
                   ("%07.2f", "0001.23"),
                   ("%.0f", "1"),
                   ("%#.0f", "1."),
                   ("%.4e", "1.2345e+00"),
                   ("%.4E", "1.2345E+00"),
                   ("%.2a", "0x1.3cp+0"),
                   ("%.2A", "0X1.3CP+0")),
      num in (1.2345, big"1.2345")
    @test @eval(@sprintf($fmt, $num) == $val)
end

# numeric spacing and various flag tests
function _test_flags(val, vflag::AbstractString, fmt::AbstractString, res::AbstractString, prefix::AbstractString)
    vflag = string("%", vflag)
    space_fmt = string(length(res) + length(prefix) + 3, fmt)
    fsign = string((val < 0 ? "-" : "+"), prefix)
    nsign = string((val < 0 ? "-" : " "), prefix)
    osign = val < 0 ? string("-", prefix) : string(prefix, "0")
    esign = string(val < 0 ? "-" : "", prefix)
    esignend = val < 0 ? "" : " "

    for (flag::AbstractString, ans::AbstractString) in (
            ("", string("  ", nsign, res)),
            ("+", string("  ", fsign, res)),
            (" ", string("  ", nsign, res)),
            ("0", string(osign, "00", res)),
            ("-", string(esign, res, "  ", esignend)),
            ("0+", string(fsign, "00", res)),
            ("0 ", string(nsign, "00", res)),
            ("-+", string(fsign, res, "  ")),
            ("- ", string(nsign, res, "  ")),
        )
        fmt_string = string(vflag, flag, space_fmt)
        @test @eval(@sprintf($fmt_string, $val) == $ans)
    end
end

for i in (
        (42, "", "i", "42", ""),
        (42, "", "d", "42", ""),

        (42, "", "u", "42", ""),
        (42, "", "x", "2a", ""),
        (42, "", "X", "2A", ""),
        (42, "", "o", "52", ""),

        (42, "#", "x", "2a", "0x"),
        (42, "#", "X", "2A", "0X"),
        (42, "#", "o", "052", ""),

        (1.2345, "", ".2f", "1.23", ""),
        (1.2345, "", ".2e", "1.23e+00", ""),
        (1.2345, "", ".2E", "1.23E+00", ""),

        (1.2345, "#", ".0f", "1.", ""),
        (1.2345, "#", ".0e", "1.e+00", ""),
        (1.2345, "#", ".0E", "1.E+00", ""),

        (1.2345, "", ".2a", "1.3cp+0", "0x"),
        (1.2345, "", ".2A", "1.3CP+0", "0X"),
    )
    _test_flags(i...)
    _test_flags(-i[1], i[2:5]...)
end

# Inf / NaN handling
@test (@sprintf "%f" Inf) == "Inf"
@test (@sprintf "%f" NaN) == "NaN"
@test (@sprintf "%f" big"Inf") == "Inf"
@test (@sprintf "%f" big"NaN") == "NaN"

# scientific notation
@test (@sprintf "%.0e" 3e142) == "3e+142"
@test (@sprintf "%#.0e" 3e142) == "3.e+142"
@test (@sprintf "%.0e" big"3e142") == "3e+142"
@test (@sprintf "%#.0e" big"3e142") == "3.e+142"

@test (@sprintf "%.0e" big"3e1042") == "3e+1042"

@test (@sprintf "%e" 3e42) == "3.000000e+42"
@test (@sprintf "%E" 3e42) == "3.000000E+42"
@test (@sprintf "%e" 3e-42) == "3.000000e-42"
@test (@sprintf "%E" 3e-42) == "3.000000E-42"
@test (@sprintf "%a" 3e4) == "0x1.d4cp+14"
@test (@sprintf "%A" 3e4) == "0X1.D4CP+14"
@test (@sprintf "%.4a" 3e-4) == "0x1.3a93p-12"
@test (@sprintf "%.4A" 3e-4) == "0X1.3A93P-12"

# %g
for (val, res) in ((12345678., "1.23457e+07"),
                   (1234567.8, "1.23457e+06"),
                   (123456.78, "123457"),
                   (12345.678, "12345.7"),
                   (12340000.0, "1.234e+07"))
    @test (@sprintf("%.6g", val) == res)
end
for (val, res) in ((big"12345678.", "1.23457e+07"),
                   (big"1234567.8", "1.23457e+06"),
                   (big"123456.78", "123457"),
                   (big"12345.678", "12345.7"))
    @test (@sprintf("%.6g", val) == res)
end
for (fmt, val) in (("%10.5g", "     123.4"),
                   ("%+10.5g", "    +123.4"),
                   ("% 10.5g","     123.4"),
                   ("%#10.5g", "    123.40"),
                   ("%-10.5g", "123.4     "),
                   ("%-+10.5g", "+123.4    "),
                   ("%010.5g", "00000123.4")),
      num in (123.4, big"123.4")
    @test @eval(@sprintf($fmt, $num) == $val)
end
@test( @sprintf( "%10.5g", -123.4 ) == "    -123.4")
@test( @sprintf( "%010.5g", -123.4 ) == "-0000123.4")
@test( @sprintf( "%.6g", 12340000.0 ) == "1.234e+07")
@test( @sprintf( "%#.6g", 12340000.0 ) == "1.23400e+07")
@test( @sprintf( "%10.5g", big"-123.4" ) == "    -123.4")
@test( @sprintf( "%010.5g", big"-123.4" ) == "-0000123.4")
@test( @sprintf( "%.6g", big"12340000.0" ) == "1.234e+07")
@test( @sprintf( "%#.6g", big"12340000.0") == "1.23400e+07")

# %g regression gh #14331
@test( @sprintf( "%.5g", 42) == "42")
@test( @sprintf( "%#.2g", 42) == "42.")
@test( @sprintf( "%#.5g", 42) == "42.000")

# hex float
@test (@sprintf "%a" 1.5) == "0x1.8p+0"
@test (@sprintf "%a" 1.5f0) == "0x1.8p+0"
@test (@sprintf "%a" big"1.5") == "0x1.8p+0"
@test (@sprintf "%#.0a" 1.5) == "0x2.p+0"
@test (@sprintf "%+30a" 1/3) == "         +0x1.5555555555555p-2"

# chars
@test (@sprintf "%c" 65) == "A"
@test (@sprintf "%c" 'A') == "A"
@test (@sprintf "%3c" 'A') == "  A"
@test (@sprintf "%-3c" 'A') == "A  "
@test (@sprintf "%c" 248) == "ø"
@test (@sprintf "%c" 'ø') == "ø"

# escape %
@test (@sprintf "%%") == "%"
@test (@sprintf "%%s") == "%s"
@test_me ArgumentError("invalid printf format string: \"%\"") @macroexpand(@sprintf "%") #" (fixes syntax highlighting)

# argument count
@test_me ArgumentError("@sprintf: wrong number of arguments (0) should be (1)") @macroexpand(@sprintf "%s")
@test_me ArgumentError("@sprintf: wrong number of arguments (2) should be (1)") @macroexpand(@sprintf "%s" "1" "2")

# no interpolation
@test_me ArgumentError("@sprintf: format must be a plain static string (no interpolation or prefix)") @macroexpand(@sprintf "$n")

# type width specifier parsing (ignored)
@test (@sprintf "%llf" 1.2) == "1.200000"
@test (@sprintf "%Lf" 1.2) == "1.200000"
@test (@sprintf "%hhu" 1) == "1"
@test (@sprintf "%hu" 1) == "1"
@test (@sprintf "%lu" 1) == "1"
@test (@sprintf "%llu" 1) == "1"
@test (@sprintf "%Lu" 1) == "1"
@test (@sprintf "%zu" 1) == "1"
@test (@sprintf "%ju" 1) == "1"
@test (@sprintf "%tu" 1) == "1"

# strings
@test (@sprintf "%s" "test") == "test"
@test (@sprintf "%s" "tést") == "tést"

@test (@sprintf "%8s" "test") == "    test"
@test (@sprintf "%-8s" "test") == "test    "

@test (@sprintf "%s" "tést") == "tést"

@test (@sprintf "%s" :test) == "test"
@test (@sprintf "%#s" :test) == ":test"
@test (@sprintf "%#8s" :test) == "   :test"
@test (@sprintf "%#-8s" :test) == ":test   "

@test (@sprintf "%8.3s" "test") == "     tes"
@test (@sprintf "%#8.3s" "test") == "     \"te"
@test (@sprintf "%-8.3s" "test") == "tes     "
@test (@sprintf "%#-8.3s" "test") == "\"te     "
@test (@sprintf "%.3s" "test") == "tes"
@test (@sprintf "%#.3s" "test") == "\"te"
@test (@sprintf "%-.3s" "test") == "tes"
@test (@sprintf "%#-.3s" "test") == "\"te"

# reasonably complex
@test (@sprintf "Test: %s%c%C%c%#-.0f." "t" 65 66 67 -42) == "Test: tABC-42.."

#test simple splatting
@test (@sprintf "%d%d" [1 2]...) == "12"

# invalid format specifiers, not "diouxXDOUeEfFgGaAcCsSpn"
for c in "bBhHIjJkKlLmMNPqQrRtTvVwWyYzZ"
    fmt_str = string("%", c)
    @test_me ArgumentError("@sprintf: first argument must be a format string") @macroexpand(@sprintf $fmt_str 1)
end

# combo
@test (@sprintf "%f %d %d %f" 1.0 [3 4]... 5) == "1.000000 3 4 5.000000"

# multi
@test (@sprintf "%s %f %9.5f %d %d %d %d%d%d%d" [1:6;]... [7,8,9,10]...) == "1 2.000000   3.00000 4 5 6 78910"

# comprehension
@test (@sprintf "%s %s %s %d %d %d %f %f %f" Any[10^x+y for x=1:3,y=1:3 ]...) == "11 101 1001 12 102 1002 13.000000 103.000000 1003.000000"

# @printf
@test_me ArgumentError("@printf: called with no arguments") @macroexpand(@printf)
@test_me ArgumentError("@printf: first or second argument must be a format string") @macroexpand(@printf 1)

# Check bug with trailing nul printing BigFloat
@test (@sprintf("%.330f", BigFloat(1)))[end] != '\0'

# Check utf8 strings #23880
@test (@sprintf("X%d", 2)) == "X2"
@test (@sprintf("\u00d0%d", 2)) == "\u00d02"
@test (@sprintf("\u0f00%d", 2)) == "\u0f002"
@test (@sprintf("\U0001ffff%d", 2)) == "\U0001ffff2"
@test (@sprintf("%dX%d", 1, 2)) == "1X2"
@test (@sprintf("%d\u00d0%d", 1, 2)) == "1\u00d02"
@test (@sprintf("%d\u0f00%d", 1, 2)) == "1\u0f002"
@test (@sprintf("%d\U0001ffff%d", 1, 2)) == "1\U0001ffff2"
@test (@sprintf("%d\u2203%d\u0203", 1, 2)) == "1\u22032\u0203"
@test_me ArgumentError @macroexpand(@sprintf("%y%d", 1, 2))
@test_me ArgumentError @macroexpand(@sprintf("%\u00d0%d", 1, 2))
@test_me ArgumentError @macroexpand(@sprintf("%\u0f00%d", 1, 2))
@test_me ArgumentError @macroexpand(@sprintf("%\U0001ffff%d", 1, 2))

# test at macro execution time
@test_throws ArgumentError("@sprintf: wrong number of arguments (2) should be (3)") (@sprintf "%d%d%d" 1:2...)

# issue #29662
@test (@sprintf "%12.3e" pi*1e100) == "  3.142e+100"

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: milesfrain/julia

max_forks_repo_name: milesfrain/julia


text:
module Cowsay

using TextWrap

struct Balloon
    thoughts::String
    left::String
    right::String
    upleft::String
    upright::String
    midleft::String
    midright::String
    downleft::String
    downright::String
end

struct CowFace
    eyes::String
    tongue::String
end

const DEFAULT_WRAPWIDTH = 30
const SAY_BALLOON = Balloon("\\", "<", ">", "/", "\\", "|", "|", "\\", "/")
const THINK_BALLOON = Balloon("o", "(", ")", "(", ")", "(", ")", "(", ")")
const DEFAULT_FACE = CowFace("oo", "  ")
const COW_TYPES = Dict{String,CowFace}(
    "default" => DEFAULT_FACE,
    "borg" => CowFace("==", "  "),
    "dead" => CowFace("xx", "U "),
    "greedy" => CowFace("\$\$", "  "),
    "paranoid" => CowFace("@@", "  "),
    "stoned" => CowFace("**", "U "),
    "tired" => CowFace("--", "  "),
    "wired" => CowFace("OO", "  "),
    "young" => CowFace("..", "  ")
    )

function cowsay(message::String=""; type::String="default", wrapwidth::Integer=DEFAULT_WRAPWIDTH, wordwrap::Bool=true)
    face = get(COW_TYPES, type, DEFAULT_FACE)

    the_cow = construct_cow(message, balloon=SAY_BALLOON, face=face, wrapwidth=wrapwidth, wordwrap=wordwrap)

    print(the_cow)
end

function cowthink(message::String=""; type::String="default", wrapwidth::Integer=DEFAULT_WRAPWIDTH, wordwrap::Bool=true)
    face = get(COW_TYPES, type, DEFAULT_FACE)

    the_cow = construct_cow(message, balloon=THINK_BALLOON, face=face, wrapwidth=wrapwidth, wordwrap=wordwrap)

    print(the_cow)
end

function construct_cow(message::String; balloon::Balloon=SAY_BALLOON, face::CowFace=DEFAULT_FACE, wrapwidth::Integer=DEFAULT_WRAPWIDTH, wordwrap::Bool=true)::String

    wrapped_message = wordwrap ? wrap(message, width=wrapwidth) : message
    lines = split(wrapped_message, "\n")
    textwidths = map(x->textwidth(x), lines)
    max_line_length = maximum(textwidths, dims=1)[1]
    padded_lines = map(x->rpad_textwidth(x, max_line_length), lines)

    message_width = max_line_length + 2
    message_height = length(padded_lines)

    cow = " $(repeat("_", message_width)) \n"

    for (index, line) in enumerate(padded_lines)
        if message_height == 1
            cow = cow * "$(balloon.left) $line $(balloon.right)\n"
        elseif index == 1
            cow = cow * "$(balloon.upleft) $line $(balloon.upright)\n"
        elseif index == message_height
            cow = cow * "$(balloon.downleft) $line $(balloon.downright)\n"
        else
            cow = cow * "$(balloon.midleft) $line $(balloon.midright)\n"
        end
    end

    cow = cow * " $(repeat("-", message_width)) \n"

    cow = cow * """
        $(balloon.thoughts)   ^__^
         $(balloon.thoughts)  ($(face.eyes))\\_______
            (__)\\       )\\/\\
             $(face.tongue) ||----w |
                ||     ||
"""

    return cow
end

function rpad_textwidth(text::AbstractString, maxwidth::Integer, padding::String=" ")::String

    padded_text = text

    while textwidth(padded_text) < maxwidth
        padded_text = padded_text * padding
    end

    return padded_text
end

end # module


text:
# count_sites_bitpar.jl
# =====================
#
# Counting sites in a bitparallel fashion.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

@inline bp_counter_type(::Type{<:Site}, ::Type{<:Alphabet}) = Int
@inline bp_start_counter(::Type{S}, ::Type{A}) where {S,A<:Alphabet} = zero(bp_counter_type(S, A))
@inline bp_update_counter(acc::Int, up::Int) = acc + up
@inline bp_correct_emptyspace(::Type{<:Site}, ::Type{<:Alphabet}) = false
@inline bp_emptyspace_correction(nempty::Int, count::Int) = count - nempty

@generated function bitpar_counter(::Type{S}, a::BioSequence{A}, b::BioSequence{A}) where {S<:Site,A<:NucAlphs}
    n = bitsof(A)
    n_elems = div(64, n)

    quote
        # Ensure that `a` is always the shorter of the two sequences.
        if length(a) > length(b)
            return bitpar_counter(S, b, a)
        end
        @assert length(a) ≤ length(b)

        nexta = bitindex(a, 1)
        stopa = bitindex(a, lastindex(a) + 1)
        nextb = bitindex(b, 1)
        stopb = bitindex(b, lastindex(b) + 1)
        counts = bp_start_counter(S, A)
#=
        println(A)
        println("Start indexes.")
        println("Index A: ", index(nexta), ", Offset A: ", offset(nexta))
        println("Index B: ", index(nextb), ", Offset B: ", offset(nextb))
        println("Stop Index A: ", index(stopa), ", Stop Offset A:", offset(stopa))
        println("Stop Index B: ", index(stopb), ", Stop Offset B:", offset(stopb))
=#

        # The first thing we need to sort out is to correctly align the head of
        # sequence / subsequence `a`s data is aligned such that the offset of
        # `nexta` is essentially reduced to 0.
        # With sequence / subsequence `a` aligned, from there, we only need to
        # worry about the alignment of sequence / subsequence `b` with respect
        # to `a`.
        if nexta < stopa && offset(nexta) != 0

            #println("Initial aligning, of a.data...")

            # Here we shift the first data chunks to the right so as the first
            # nucleotide of the seq/subseq is the first nibble / pair of bits.
            x = a.data[index(nexta)] >> offset(nexta)
            y = b.data[index(nextb)] >> offset(nextb)

            #println("x: ", hex(x), "\ny: ", hex(y))

            # Here it was assumed that there is something to go and get from
            # the next chunk of `b`, yet that may not be true.
            # We know that if this is not true of `b`, then it is certainly not
            # true of `a`.
            # We check if the end of the sequence is not contained in the same
            # integer like so: `64 - offset(nextb) < stopb - nextb`.
            #
            # This edge case was found and accounted for by Ben Ward @Ward9250.
            # Ask this maintainer for more information.

            #=
            println("64 - offset(nextb): ", 64 - offset(nextb))
            println("stopb - nextb: ", stopb - nextb)
            println("64 - offset(nextb) < stopb - nextb: ", 64 - offset(nextb) < stopb - nextb)
            =#
            if offset(nextb) > offset(nexta) && 64 - offset(nextb) < stopb - nextb
                y |= b.data[index(nextb) + 1] << (64 - offset(nextb))
                #println("Modified y: ", hex(y))
            end

            # Here we need to check something, we need to check if the
            # integer of `a` we are currently aligning contains the end of
            # seq/subseq `a`. Because if so it's something we need to take into
            # account of when we mask x and y.
            #
            # In other words if `64 - offset(nexta) > stopa - nexta, we know
            # seq or subseq a's data ends before the end of this data chunk,
            # and so the mask used needs to be defined to account for this:
            # `mask(stopa - nexta)`, otherwise the mask simply needs to be
            # `mask(64 - offset(nexta))`.
            #
            # This edge case was found and accounted for by Ben Ward @Ward9250.
            # Ask this maintainer for more information.

            #=
            println("stopa and nexta stuff:")
            println("64 - offset(nexta): ", 64 - offset(nexta))
            println("stopa - nexta: ", stopa - nexta)
            println("64 - offset(nexta) > stopa - nexta: ", 64 - offset(nexta) > stopa - nexta)
            println("stopa based mask: ", hex(mask(stopa - nexta)))
            println("64 - offset based mask: ", hex(mask(64 - offset(nexta))))
            =#

            k = ifelse(64 - offset(nexta) > stopa - nexta, stopa - nexta, 64 - offset(nexta))
            #println("k: ", k)
            m = mask(k)
            #=
            println("mask used: ", hex(m))
            println("masked x: ", hex(x & m))
            println("masked y: ", hex(y & m))
            =#
            counts = bp_update_counter(counts, bp_chunk_count(S, A, x & m, y & m))
            if bp_correct_emptyspace(S, A)
                #println("Correcting for emptyspace...")
                nempty = $n_elems - div(k, $n)
                #println("nempty: ", nempty)
                counts = bp_emptyspace_correction(nempty, counts)
                #println("counts: ", counts)
            end
            # Here we move our current position markers by k, meaning they move
            # to either, A). The next integer, or B). The end of the sequence if
            # it is in the current integer.
            nexta += k
            nextb += k
            #=
            println("Index A: ", index(nexta), ", Offset A: ", offset(nexta))
            println("Index B: ", index(nextb), ", Offset B: ", offset(nextb))
            println("Finished Aligning of a.data...")
            =#
        end
        #@assert offset(nexta) == 0

        if offset(nextb) == 0  # data are aligned with each other
            #println("Data are aligned...")
            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                #println("x: ", hex(x))
                #println("y: ", hex(y))
                counts = bp_update_counter(counts, bp_chunk_count(S, A, x, y))
                #println("counts: ", counts)
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                #println("Processing tail...")

                x = a.data[index(nexta)]
                #println("x: ", hex(x))
                y = b.data[index(nextb)]
                #println("y: ", hex(y))

                offs = stopa - nexta
                #println("offs: ", offs)
                m = mask(offs)
                #println("mask: ", hex(m))
                #println("masked x: ", hex(x & m))
                #println("masked y: ", hex(y & m))
                counts = bp_update_counter(counts, bp_chunk_count(S, A, x & m, y & m))
                #println("counts: ", counts)
                if bp_correct_emptyspace(S, A)
                    #println("Correcting for emptyspace...")
                    nempty = $n_elems - div(offs, $n)
                    #println("nempty: ", nempty)
                    counts = bp_emptyspace_correction(nempty, counts)
                    #println("counts: ", counts)
                end
            end
        elseif nexta < stopa
            #println("Data are unaligned...")

            #println("Index B: ", index(nextb), " Offset B: ", offset(nextb))
            y = b.data[index(nextb)]
            #println("y: ", hex(y))
            nextb += 64
            # Note that here, updating `nextb` by 64, increases the chunk index,
            # but the `offset(nextb)` will remain the same.

            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                z = b.data[index(nextb)]
                y = y >> offset(nextb) | z << (64 - offset(nextb))
                #println("x: ", hex(x))
                #println("z: ", hex(z))
                #println("y: ", hex(y))

                counts = bp_update_counter(counts, bp_chunk_count(S, A, x, y))
                #println("counts: ", counts)
                y = z
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                #println("Processing tail...")

                #println("Index A: ", index(nexta), " Offset A: ", offset(nexta))
                #println("Index B: ", index(nextb), " Offset B: ", offset(nextb))

                # Get the data chunk of `a`.
                x = a.data[index(nexta)]
                #println("x: ", hex(x))
                y = y >> offset(nextb)
                #println("y: ", hex(y))

                if 64 - offset(nextb) < stopa - nexta
                    #y |= a.data[index(nextb)] << (64 - offset(nextb))
                    y |= b.data[index(nextb)] << (64 - offset(nextb))
                end
                #println("modified y: ", hex(y))

                offs = stopa - nexta
                #println("offs: ", offs)
                m = mask(offs)
                #println("mask: ", hex(m))
                #println("masked x: ", hex(x & m))
                #println("masked y: ", hex(y & m))
                counts = bp_update_counter(counts, bp_chunk_count(S, A, x & m, y & m))
                #println("counts: ", counts)
                if bp_correct_emptyspace(S, A)
                    #println("Correcting for emptyspace...")
                    nempty = $n_elems - div(offs, $n)
                    #println("nempty: ", nempty)
                    counts = bp_emptyspace_correction(nempty, counts)
                    #println("counts: ", counts)
                end
            end
        end
        return counts
    end
end


text:
#
# Linear Algebra of order 2 Tensors
#
# Even though CuDenseTensor{_,2} is strided
# and passable to BLAS/LAPACK, it cannot
# be made <: StridedArray

function Base.:*(T1::Tensor{ElT1,2,StoreT1,IndsT1},
                 T2::Tensor{ElT2,2,StoreT2,IndsT2}) where
                                        {ElT1,StoreT1<:CuDense,IndsT1,
                                        ElT2,StoreT2<:CuDense,IndsT2}
  RM    = matrix(T1)*matrix(T2)
  indsR = IndsT1(ind(T1,1),ind(T2,2))
  pT    = promote_type(ElT1,ElT2)
  return Tensor(Dense(vec(RM)),indsR)
end
#= FIX ME
function LinearAlgebra.exp(T::CuDenseTensor{ElT,2}) where {ElT,IndsT}
  expTM = exp(matrix(T))
  return Tensor(Dense(vec(expTM)),inds(T))
end

function expHermitian(T::CuDenseTensor{ElT,2}) where {ElT,IndsT}
  # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric,
  # so extract the parent matrix
  expTM = parent(exp(Hermitian(matrix(T))))
  return Tensor(Dense(vec(expTM)),inds(T))
end
=#

# svd of an order-2 tensor
function LinearAlgebra.svd(T::CuDenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  fastSVD::Bool = get(kwargs,:fastSVD,false)
  aT = array(T)
  @timeit "CUSOLVER svd" begin
      MU,MS,MV = CUSOLVER.svd!(aT)
  end
  #conj!(MV)
  P = MS.^2
  truncerr, docut, P = truncate!(P;mindim=mindim,
              maxdim=maxdim,
              cutoff=cutoff,
              absoluteCutoff=absoluteCutoff,
              doRelCutoff=doRelCutoff)
  spec = Spectrum(P,truncerr)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:,1:dS]
    MS = MS[1:dS]
    MV = MV[:,1:dS]
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dS)
  v = eltype(IndsT)(dS)
  Uinds = IndsT((ind(T,1),u))
  Sinds = IndsT((u,v))
  Vinds = IndsT((ind(T,2),v))
  U = Tensor(Dense(vec(MU)),Uinds)
  Sdata      = CUDA.zeros(ElT, dS * dS)
  dsi        = diagind(reshape(Sdata, dS, dS), 0)
  Sdata[dsi] = MS
  S = Tensor(Dense(Sdata),Sinds)
  V = Tensor(Dense(vec(MV)),Vinds)
  return U,S,V,spec
end

function LinearAlgebra.eigen(T::Hermitian{ElT,<:CuDenseTensor{ElT,2,IndsT}};
                        kwargs...) where {ElT<:Union{Real,Complex},IndsT}
  ispossemidef::Bool = get(kwargs,:ispossemidef,false)
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  @timeit "CUSOLVER eigen" begin
      local DM, UM
      if ElT <: Complex
          DM, UM = CUSOLVER.heevd!('V', 'U', matrix(parent(T)))
      else
          DM, UM = CUSOLVER.syevd!('V', 'U', matrix(parent(T)))
      end
  end
  DM_ = reverse(DM)
  @timeit "truncate" begin
      truncerr, docut, DM = truncate!(DM_;maxdim=maxdim, cutoff=cutoff, absoluteCutoff=absoluteCutoff, doRelCutoff=doRelCutoff)
  end
  spec = Spectrum(DM,truncerr)
  dD = length(DM)
  dV = reverse(UM, dims=2)
  if dD < size(dV,2)
      dV = CuMatrix(dV[:,1:dD])
  end
  # Make the new indices to go onto U and V
  l = eltype(IndsT)(dD)
  r = eltype(IndsT)(dD)
  Vinds = IndsT((dag(ind(T, 2)), dag(r)))
  Dinds = IndsT((l, dag(r)))
  U = Tensor(Dense(vec(dV)),Vinds)
  D = Tensor(Diag(real.(DM)),Dinds)
  return D,U,spec
end

function LinearAlgebra.qr(T::CuDenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
  QM,RM = qr(matrix(T))
  # Make the new indices to go onto Q and R
  q,r = inds(T)
  q = dim(q) < dim(r) ? sim(q) : sim(r)
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  QM = CuMatrix(QM)
  Q = Tensor(Dense(vec(QM)),Qinds)
  R = Tensor(Dense(vec(RM)),Rinds)
  return Q,R
end

function polar(T::CuDenseTensor{ElT,2,IndsT}) where {ElT,IndsT}
  QM,RM = polar(matrix(T))
  dim = size(QM,2)
  # Make the new indices to go onto Q and R
  q = eltype(IndsT)(dim)
  # TODO: use push/pushfirst instead of a constructor
  # call here
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = Tensor(Dense(vec(QM)),Qinds)
  R = Tensor(Dense(vec(RM)),Rinds)
  return Q,R
end


text:
"""
$(TYPEDSIGNATURES)

Generate `ODESystem`, dependent variables, and parameters from an `ODEProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.ODEProblem; kwargs...)
    prob.f isa DiffEqBase.AbstractParameterizedFunction &&
        return prob.f.sys
    @parameters t

    p = prob.p
    has_p = !(p isa Union{DiffEqBase.NullParameters, Nothing})

    _vars = define_vars(prob.u0, t)

    vars = prob.u0 isa Number ? _vars : ArrayInterfaceCore.restructure(prob.u0, _vars)
    params = if has_p
        _params = define_params(p)
        p isa Number ? _params[1] :
        (p isa Tuple || p isa NamedTuple ? _params :
         ArrayInterfaceCore.restructure(p, _params))
    else
        []
    end

    var_set = Set(vars)

    D = Differential(t)
    mm = prob.f.mass_matrix

    if mm === I
        lhs = map(v -> D(v), vars)
    else
        lhs = map(mm * vars) do v
            if iszero(v)
                0
            elseif v in var_set
                D(v)
            else
                error("Non-permuation mass matrix is not supported.")
            end
        end
    end

    if DiffEqBase.isinplace(prob)
        rhs = ArrayInterfaceCore.restructure(prob.u0, similar(vars, Num))
        prob.f(rhs, vars, params, t)
    else
        rhs = prob.f(vars, params, t)
    end

    eqs = vcat([lhs[i] ~ rhs[i] for i in eachindex(prob.u0)]...)

    sts = vec(collect(vars))

    params = if params isa Number || (params isa Array && ndims(params) == 0)
        [params[1]]
    else
        vec(collect(params))
    end
    default_u0 = Dict(sts .=> vec(collect(prob.u0)))
    default_p = has_p ? Dict(params .=> vec(collect(prob.p))) : Dict()

    de = ODESystem(eqs, t, sts, params,
                   defaults = merge(default_u0, default_p);
                   name = gensym(:MTKizedODE),
                   kwargs...)

    de
end

_defvaridx(x, i, t) = variable(x, i, T = SymbolicUtils.FnType{Tuple, Real})
_defvar(x, t) = variable(x, T = SymbolicUtils.FnType{Tuple, Real})

function define_vars(u, t)
    _vars = [_defvaridx(:x, i, t)(t) for i in eachindex(u)]
end

function define_vars(u::Union{SLArray, LArray}, t)
    _vars = [_defvar(x, t)(t) for x in LabelledArrays.symnames(typeof(u))]
end

function define_vars(u::Tuple, t)
    _vars = tuple((_defvaridx(:x, i, t)(ModelingToolkit.value(t)) for i in eachindex(u))...)
end

function define_vars(u::NamedTuple, t)
    _vars = NamedTuple(x => _defvar(x, t)(ModelingToolkit.value(t)) for x in keys(u))
end

function define_params(p)
    [toparam(variable(:α, i)) for i in eachindex(p)]
end

function define_params(p::Union{SLArray, LArray})
    [toparam(variable(x)) for x in LabelledArrays.symnames(typeof(p))]
end

function define_params(p::Tuple)
    tuple((toparam(variable(:α, i)) for i in eachindex(p))...)
end

function define_params(p::NamedTuple)
    NamedTuple(x => toparam(variable(x)) for x in keys(p))
end

"""
$(TYPEDSIGNATURES)

Generate `SDESystem`, dependent variables, and parameters from an `SDEProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.SDEProblem; kwargs...)
    prob.f isa DiffEqBase.AbstractParameterizedFunction &&
        return (prob.f.sys, prob.f.sys.states, prob.f.sys.ps)
    @parameters t
    p = prob.p
    has_p = !(p isa Union{DiffEqBase.NullParameters, Nothing})

    _vars = define_vars(prob.u0, t)

    vars = prob.u0 isa Number ? _vars : ArrayInterfaceCore.restructure(prob.u0, _vars)
    params = if has_p
        _params = define_params(p)
        p isa Number ? _params[1] :
        (p isa Tuple || p isa NamedTuple ? _params :
         ArrayInterfaceCore.restructure(p, _params))
    else
        []
    end

    D = Differential(t)

    rhs = [D(var) for var in vars]

    if DiffEqBase.isinplace(prob)
        lhs = similar(vars, Any)

        prob.f(lhs, vars, params, t)

        if DiffEqBase.is_diagonal_noise(prob)
            neqs = similar(vars, Any)
            prob.g(neqs, vars, params, t)
        else
            neqs = similar(vars, Any, size(prob.noise_rate_prototype))
            prob.g(neqs, vars, params, t)
        end
    else
        lhs = prob.f(vars, params, t)
        if DiffEqBase.is_diagonal_noise(prob)
            neqs = prob.g(vars, params, t)
        else
            neqs = prob.g(vars, params, t)
        end
    end
    deqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)

    params = if ndims(params) == 0
        [params[1]]
    else
        Vector(vec(params))
    end

    de = SDESystem(deqs, neqs, t, Vector(vec(vars)), params;
                   name = gensym(:MTKizedSDE),
                   kwargs...)

    de
end

"""
$(TYPEDSIGNATURES)

Generate `OptimizationSystem`, dependent variables, and parameters from an `OptimizationProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.OptimizationProblem; kwargs...)
    if prob.p isa Tuple || prob.p isa NamedTuple
        p = [x for x in prob.p]
    else
        p = prob.p
    end

    vars = reshape([variable(:x, i) for i in eachindex(prob.u0)], size(prob.u0))
    params = p isa DiffEqBase.NullParameters ? [] :
             reshape([variable(:α, i) for i in eachindex(p)], size(Array(p)))

    eqs = prob.f(vars, params)
    de = OptimizationSystem(eqs, vec(vars), vec(params);
                            name = gensym(:MTKizedOpt),
                            kwargs...)
    de
end


text:
using MathOptInterface, Test

const CleverDicts = MathOptInterface.Utilities.CleverDicts

# Note: `MyKey` is just for testing. You wouldn't want to use it in practice
# because the key type of the dictionary isn't a concrete type.
struct MyKey{X} end
CleverDicts.key_to_index(::MyKey{X}) where {X} = X
CleverDicts.index_to_key(::Type{MyKey}, index::Int) = MyKey{index}()

@testset "CleverDict" begin
    @testset "MyKey type" begin
        d = CleverDicts.CleverDict{MyKey, String}()
        key = CleverDicts.add_item(d, "first")
        @test key == MyKey{1}()
        @test d[MyKey{1}()] == "first"
    end

    @testset "get/set" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        @test key == MathOptInterface.VariableIndex(1)
        sizehint!(d, 1)
        @test d[key] == "first"
        @test haskey(d, key) == true
        @test_throws KeyError d[MathOptInterface.VariableIndex(2)]
        delete!(d, key)
        sizehint!(d, 2)
        @test_throws KeyError d[key]
        # set index is valid now
        # @test_throws KeyError d[key] = "key"
        @test haskey(d, key) == false
        key2 = CleverDicts.add_item(d, "second")
        @test key2 == MathOptInterface.VariableIndex(2)
        @test d[key2] == "second"
        d[key2] = "third"
        @test d[key2] == "third"

        empty!(d)

        key = CleverDicts.add_item(d, "first")
        @test key == MathOptInterface.VariableIndex(1)
        @test d[key] == "first"
        d[key] = "zeroth"
        @test d[key] == "zeroth"
        @test haskey(d, key) == true
        @test_throws KeyError d[MathOptInterface.VariableIndex(2)]
        delete!(d, key)
        @test_throws KeyError d[key]
        # set index is valid now
        # @test_throws KeyError d[key] = "key"
        @test haskey(d, key) == false
        key2 = CleverDicts.add_item(d, "second")
        @test key2 == MathOptInterface.VariableIndex(2)
        @test d[key2] == "second"
    end

    @testset "LinearIndex" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        @test d[CleverDicts.LinearIndex(1)] == "first"
        key2 = CleverDicts.add_item(d, "second")
        @test d[CleverDicts.LinearIndex(2)] == "second"
        @test length(d) == 2
        delete!(d, key)
        @test d[CleverDicts.LinearIndex(1)] == "second"
        @test_throws KeyError d[CleverDicts.LinearIndex(2)]
        @test length(d) == 1
    end

    @testset "keys/values" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        key2 = CleverDicts.add_item(d, "second")
        @test collect(keys(d)) == [MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)]
        @test collect(values(d)) == ["first", "second"]
        delete!(d, key)
        key3 = CleverDicts.add_item(d, "third")
        @test collect(keys(d)) == [MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(3)]
        @test collect(values(d)) == ["second", "third"]
    end

    @testset "iterate" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        key2 = CleverDicts.add_item(d, "second")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
           push!(my_keys, k)
           push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)]
        @test my_values == ["first", "second"]
        delete!(d, key)
        key3 = CleverDicts.add_item(d, "third")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
           push!(my_keys, k)
           push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(3)]
        @test my_values == ["second", "third"]
    end

    @testset "iterate ii" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        key2 = CleverDicts.add_item(d, "second")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
            push!(my_keys, k)
            push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)]
        @test my_values == ["first", "second"]
        delete!(d, key)
        @test d[CleverDicts.LinearIndex(1)] == "second"
        key3 = CleverDicts.add_item(d, "third")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
            push!(my_keys, k)
            push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(3)]
        @test my_values == ["second", "third"]
    end

    @testset "iterate iii" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        y = 0
        for (k, v) in d
            y += 1
        end
        @test y == 0
    end

    @testset "haskey" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        @test !haskey(d, 1)
        k = CleverDicts.add_item(d, "a")
        @test haskey(d, k)
        j = CleverDicts.add_item(d, "b")
        @test haskey(d, j)
        delete!(d, k)
        @test !haskey(d, k)
        @test haskey(d, j)
    end

    @testset "haskey" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        @test isempty(d) == true
        k = CleverDicts.add_item(d, "a")
        @test isempty(d) == false
        delete!(d, k)
        @test isempty(d) == true
        j = CleverDicts.add_item(d, "b")
        @test isempty(d) == false
    end

    @testset "delete!" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        @test length(d) == 0
        @test delete!(d, MathOptInterface.VariableIndex(0)) == nothing
        k1 = CleverDicts.add_item(d, "a")
        k2 = CleverDicts.add_item(d, "b")
        d[CleverDicts.LinearIndex(2)] == "b"
        delete!(d, k1)
        d[CleverDicts.LinearIndex(1)] == "b"
        k3 = CleverDicts.add_item(d, "c")
        @test d[k3] == "c"
        @test d[CleverDicts.LinearIndex(1)] == "b"
        @test d[CleverDicts.LinearIndex(2)] == "c"
    end

    @testset "Dense Operations" begin
        mul2(x) = x * 2
        div2(x) = div(x, 2)
        d = CleverDicts.CleverDict{Int, Float64}(div2, mul2, 3)

        d[4] = 0.25
        @test !haskey(d, 2)
        @test haskey(d, 4)
        @test !haskey(d, 6)
        @test length(d) == 1
        @test collect(d) == [4 => 0.25]
        @test d[4] == 0.25

        d[2] = 1.5
        @test haskey(d, 2)
        @test haskey(d, 4)
        @test !haskey(d, 6)
        @test length(d) == 2
        @test collect(d) == [2 => 1.5, 4 => 0.25]
        @test d[2] == 1.5
        @test d[4] == 0.25

        d[6] = 0.75
        @test haskey(d, 2)
        @test haskey(d, 4)
        @test haskey(d, 6)
        @test length(d) == 3
        @test collect(d) == [2 => 1.5, 4 => 0.25, 6 => 0.75]
        @test d[2] == 1.5
        @test d[4] == 0.25
        @test d[6] == 0.75
    end
end


text:
using Test
using StringCases

import StringCases._defaultcase
import StringCases._decamelize

# -------------------
#  test plain string
# -------------------

test_string = "foo BAR bAZ"

StringCases.camelize("dwed_dewew")

@test _defaultcase(test_string) == "foo bar baz"

@test spacecase(test_string) == "foo bar baz"

@test pascalcase(test_string) == "FooBarBaz"

@test camelcase(test_string) == "fooBarBaz"

@test snakecase(test_string) == "foo_bar_baz"

@test kebabcase(test_string) == "foo-bar-baz"

# -------------
#  test delims
# -------------

test_string = "this_var-here"

@test _defaultcase(test_string) == "this var here"

@test spacecase(test_string) == "this var here"

@test pascalcase(test_string) == "ThisVarHere"

@test camelcase(test_string) == "thisVarHere"

@test snakecase(test_string) == "this_var_here"

@test kebabcase(test_string) == "this-var-here"

# -----------------
#  test underscore
# -----------------

test_string = "foo_bar_baz"

@test underscore("foo_bar_baz") == test_string

@test underscore("foo bar baz") == test_string

@test underscore("foo-bar-baz") == test_string

@test underscore("fooBarBaz") == test_string

@test underscore("FooBarBaz") == test_string

@test underscore("modelA") == "model_a"

@test underscore("ModelA") == "model_a"

@test underscore("aModel") == "a_model"

@test underscore("AModel") == "amodel"

# -----------------
#  test dasherize
# -----------------

test_string = "foo-bar-baz"

@test dasherize("foo_bar_baz") == test_string

@test dasherize("foo bar baz") == test_string

@test dasherize("foo-bar-baz") == test_string

@test dasherize("fooBarBaz") == test_string

@test dasherize("FooBarBaz") == test_string

@test dasherize("modelA") == "model-a"

@test dasherize("ModelA") == "model-a"

@test dasherize("aModel") == "a-model"

@test dasherize("AModel") == "amodel"

# -----------------
#  test decamelize
# -----------------

test_string = "foo_bar_baz"

@test _decamelize("foo_bar_baz") == test_string

@test _decamelize("foo bar baz") == "foo bar baz"

@test _decamelize("foo-bar-baz") == "foo-bar-baz"

@test _decamelize("fooBarBaz") == test_string

@test _decamelize("FooBarBaz") == test_string

@test _decamelize("modelA") == "model_a"

@test _decamelize("ModelA") == "model_a"

@test _decamelize("aModel") == "a_model"

@test _decamelize("AModel") == "amodel"

# -------------
#  test spaces
# -------------

test_string = "boo"

@test _defaultcase("boo") == test_string

@test _defaultcase("  boo") == test_string

@test _defaultcase("boo  ") == test_string

@test _defaultcase(" boo ") == test_string

# ============
#  tests from
#   ember.js
# ============

include("ember_js_test.jl")


text:
export DryModel, EquilMoist

#### Moisture component in atmosphere model
abstract type MoistureModel end

vars_state(::MoistureModel, FT) = @vars()
vars_gradient(::MoistureModel, FT) = @vars()
vars_diffusive(::MoistureModel, FT) = @vars()
vars_aux(::MoistureModel, FT) = @vars()

function atmos_nodal_update_aux!(::MoistureModel, m::AtmosModel, state::Vars,
                                 aux::Vars, t::Real)
end
function flux_moisture!(::MoistureModel, ::AtmosModel, flux::Grad, state::Vars, aux::Vars, t::Real)
end
function diffusive!(::MoistureModel, diffusive, ∇transform, state, aux, t)
end
function flux_diffusive!(::MoistureModel, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real, D_t)
end
function flux_nondiffusive!(::MoistureModel, flux::Grad, state::Vars, aux::Vars, t::Real)
end
function gradvariables!(::MoistureModel, transform::Vars, state::Vars, aux::Vars, t::Real)
end

@inline function internal_energy(moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars)
  MoistThermodynamics.internal_energy(state.ρ, state.ρe, state.ρu, gravitational_potential(orientation, aux))
end
@inline temperature(moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars) = air_temperature(thermo_state(moist, orientation, state, aux))
@inline pressure(moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars) = air_pressure(thermo_state(moist, orientation, state, aux))
@inline soundspeed(moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars) = soundspeed_air(thermo_state(moist, orientation, state, aux))

@inline function total_specific_enthalpy(moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars)
  phase = thermo_state(moist, orientation, state, aux)
  R_m = gas_constant_air(phase)
  T = air_temperature(phase)
  e_tot = state.ρe * (1/state.ρ)
  e_tot + R_m*T
end

"""
    DryModel

Assumes the moisture components is in the dry limit.
"""
struct DryModel <: MoistureModel
end

vars_aux(::DryModel,FT) = @vars(θ_v::FT)
@inline function atmos_nodal_update_aux!(moist::DryModel, atmos::AtmosModel,
                                         state::Vars, aux::Vars, t::Real)
  e_int = internal_energy(moist, atmos.orientation, state, aux)
  TS = PhaseDry(e_int, state.ρ)
  aux.moisture.θ_v = virtual_pottemp(TS)
  nothing
end

thermo_state(moist::DryModel, orientation::Orientation, state::Vars, aux::Vars) = PhaseDry(internal_energy(moist, orientation, state, aux), state.ρ)

"""
    EquilMoist

Assumes the moisture components are computed via thermodynamic equilibrium.
"""
Base.@kwdef struct EquilMoist <: MoistureModel
  maxiter::Int = 3
end
vars_state(::EquilMoist,FT) = @vars(ρq_tot::FT)
vars_gradient(::EquilMoist,FT) = @vars(q_tot::FT, h_tot::FT)
vars_diffusive(::EquilMoist,FT) = @vars(∇q_tot::SVector{3,FT})
vars_aux(::EquilMoist,FT) = @vars(temperature::FT, θ_v::FT, q_liq::FT)

@inline function atmos_nodal_update_aux!(moist::EquilMoist, atmos::AtmosModel,
                                         state::Vars, aux::Vars, t::Real)
  e_int = internal_energy(moist, atmos.orientation, state, aux)
  TS = PhaseEquil(e_int, state.ρ, state.moisture.ρq_tot/state.ρ, moist.maxiter)
  aux.moisture.temperature = air_temperature(TS)
  aux.moisture.θ_v = virtual_pottemp(TS)
  aux.moisture.q_liq = PhasePartition(TS).liq
  nothing
end

function thermo_state(moist::EquilMoist, orientation::Orientation, state::Vars, aux::Vars)
  e_int = internal_energy(moist, orientation, state, aux)
  FT = eltype(state)
  return PhaseEquil{FT}(e_int, state.ρ, state.moisture.ρq_tot/state.ρ, aux.moisture.temperature)
end

function gradvariables!(moist::EquilMoist, transform::Vars, state::Vars, aux::Vars, t::Real)
  ρinv = 1/state.ρ
  transform.moisture.q_tot = state.moisture.ρq_tot * ρinv
end

function diffusive!(moist::EquilMoist, diffusive::Vars, ∇transform::Grad, state::Vars, aux::Vars, t::Real)
  # diffusive flux of q_tot
  diffusive.moisture.∇q_tot = ∇transform.moisture.q_tot
end

function flux_moisture!(moist::EquilMoist, atmos::AtmosModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  ρ = state.ρ
  u = state.ρu / ρ
  flux.moisture.ρq_tot += u * state.moisture.ρq_tot
end

function flux_diffusive!(moist::EquilMoist, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real, D_t)
  d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
  flux_diffusive!(moist, flux, state, d_q_tot)
end
#TODO: Consider whether to not pass ρ and ρu (not state), foc BCs reasons
function flux_diffusive!(moist::EquilMoist, flux::Grad, state::Vars, d_q_tot)
  flux.ρ += d_q_tot * state.ρ
  flux.ρu += d_q_tot .* state.ρu'
  flux.moisture.ρq_tot += d_q_tot * state.ρ
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

Makie.plottype(::SimpleMesh) = Viz{<:Tuple{SimpleMesh}}

function Makie.plot!(plot::Viz{<:Tuple{SimpleMesh}})
  # retrieve mesh object
  mesh = plot[:object][]

  color      = plot[:color][]
  facetcolor = plot[:facetcolor][]
  showfacets = plot[:showfacets][]

  # relevant settings
  dim   = embeddim(mesh)
  nvert = nvertices(mesh)
  nelem = nelements(mesh)
  verts = vertices(mesh)
  topo  = topology(mesh)
  elems = elements(topo)

  # coordinates of vertices
  coords = coordinates.(verts)

  # fan triangulation (assume convexity)
  tris4elem = map(elems) do elem
    I = indices(elem)
    [[I[1], I[i], I[i+1]] for i in 2:length(I)-1]
  end

  # flatten vector of triangles
  tris = [tri for tris in tris4elem for tri in tris]

  # element vs. vertex coloring
  if color isa AbstractVector
    ncolor = length(color)
    if ncolor == nelem # element coloring
      # duplicate vertices and adjust
      # connectivities to avoid linear
      # interpolation of colors
      nt = 0
      elem4tri = Dict{Int,Int}()
      for e in 1:nelem
        Δs = tris4elem[e]
        for _ in 1:length(Δs)
          nt += 1
          elem4tri[nt] = e
        end
      end
      nv = 3nt
      tcoords = [coords[i] for tri in tris for i in tri]
      tconnec = [collect(I) for I in Iterators.partition(1:nv, 3)]
      tcolors = map(1:nv) do i
        t = ceil(Int, i/3)
        e = elem4tri[t]
        color[e]
      end
    elseif ncolor == nvert # vertex coloring
      # nothing needs to be done because
      # this is the default in Makie and
      # because the triangulation above
      # does not change the vertices in
      # the original polygonal mesh
      tcoords = coords
      tconnec = tris
      tcolors = color
    else
      throw(ArgumentError("Provided $ncolor colors but the mesh has
                           $nvert vertices and $nelem elements"))
    end
  else # single color
    # nothing needs to be done
    tcoords = coords
    tconnec = tris
    tcolors = color
  end

  # convert connectivities to matrix format
  tmatrix = reduce(hcat, tconnec) |> transpose

  # enable shading in 3D
  shading = dim == 3

  Makie.mesh!(plot, tcoords, tmatrix,
    color = tcolors,
    colormap = plot[:colormap],
    shading = shading,
  )

  if showfacets
    # use a sophisticated data structure
    # to extract the edges from the n-gons
    t = convert(HalfEdgeTopology, topo)
    ∂ = Boundary{1,0}(t)

    # append indices of incident vertices
    # interleaved with a sentinel index
    inds = Int[]
    for i in 1:nfacets(t)
      append!(inds, ∂(i))
      push!(inds, nvert+1)
    end

    # fill sentinel index with NaN coordinates
    push!(coords, Vec(ntuple(i->NaN, dim)))

    # extract incident vertices
    coords = coords[inds]

    # split coordinates to match signature
    xyz = [getindex.(coords, j) for j in 1:dim]

    Makie.lines!(plot, xyz...,
      color = facetcolor,
    )
  end
end


text:
function uptimes_bar(S::SeisData, fmt::String, use_name::Bool, nxt::Int64)
  xmi = 2^63-1
  xma = xmi+1
  fig = PyPlot.figure(figsize=[8.0, 6.0], dpi=150)
  ax = PyPlot.axes([0.20, 0.10, 0.72, 0.85])

  for i = 1:S.n
    if S.fs[i] == 0.0
      t = view(S.t[i], :, 2)
      x = rescaled(S.x[i].-mean(S.x[i]),i)
      plot(t, x, "ko", markersize = 2.0)
    else
      rect_x = Array{Float64,1}(undef, 0)
      rect_y = Array{Float64,1}(undef, 0)
      t = t_win(S.t[i], S.fs[i])
      for j = 1:size(t,1)
        ts = Float64(t[j,1])
        te = Float64(t[j,2])
        append!(rect_x, [ts, te, te, ts, ts])
        append!(rect_y, [i-0.48, i-0.48, i+0.48, i+0.48, i-0.48])
      end

      p = fill(rect_x, rect_y, linewidth=1.0, edgecolor="k")
    end
    xmi = min(xmi, first(t))
    xma = max(xma, last(t))
  end

  PyPlot.title("Channel Uptimes", fontweight="bold", fontsize=13.0, family="serif", color="black")

  # Y scaling and axis manipulation
  PyPlot.yticks(1:S.n, map((i) -> replace(i, " " => ""), use_name ? S.name : S.id))
  PyPlot.ylim(0.5, S.n+0.5)
  PyPlot.setp(gca().get_yticklabels(), fontsize=8.0, color="black", fontweight="bold", family="serif")

  # X scaling and axis manipulation
  xfmt(xmi, xma, fmt, nxt)
  PyPlot.setp(gca().get_xticklabels(), fontsize=10.0, color="black", fontweight="bold", family="serif")

  return fig
end

function uptimes_sum(S::SeisData, fmt::String, use_name::Bool, nxt::Int64)
  ntr = sum(S.fs.>0)
  W = Array{Int64, 1}(undef, 0)
  for i = 1:S.n
    S.fs[i] == 0.0 && continue
    w = div.(t_win(S.t[i], S.fs[i]), 1000000)
    for i = 1:size(w,1)
      append!(W, collect(w[i,1]:w[i,2]))
    end
  end
  sort!(W)
  t = collect(first(W):last(W))
  y = zeros(Int64, length(t))
  i = 0
  j = 1
  τ = first(t)
  while i < length(W)
    i = i + 1
    if W[i] == τ
      y[j] += 1
    else
      while W[i] != τ
        j += 1
        j > length(t) && break
        τ = t[j]
      end
      y[j] += 1
    end
  end
  fig = PyPlot.figure(figsize=[8.0, 6.0], dpi=150)
  step(t.*1000000, y, color="k", linewidth=2.0)

  # Y scaling and axis manipulation
  dy = ntr > 24 ? 5.0 : ntr > 19 ? 4.0 : ntr > 14 ? 3.0 : ntr > 9 ? 2.0 : 1.0
  PyPlot.ylim(0, ntr)
  PyPlot.yticks(ntr:-dy:0)
  PyPlot.ylabel("Active Channels", fontweight="bold", fontsize=12.0, family="serif", color="black")
  PyPlot.setp(gca().get_yticklabels(), fontsize=10.0, color="black", fontweight="bold", family="serif")

  # X scaling and axis manipulation
  xfmt(first(t)*1000000, last(t)*1000000, fmt, nxt)
  PyPlot.setp(gca().get_xticklabels(), fontsize=10.0, color="black", fontweight="bold", family="serif")
  return fig
end

"""
    uptimes(S[, summed=false])

Bar plot of uptimes for each channel in S.

If summed==true, plot uptimes for all channels in S that record timeseries data,
scaled so that y=1 corresponds to 100% of channels active. Non-timeseries
channels in S are not counted toward the cumulative total in a summed uptime plot.
"""
function uptimes(S::SeisData; summed::Bool=false, fmt::String="auto", use_name::Bool=false, nxt::Int64=5)
  if summed
    fig = uptimes_sum(S, fmt, use_name, nxt)
  else
    fig = uptimes_bar(S, fmt, use_name, nxt)
  end
  return fig
end


text:
    # Test for GaussLobatto()
    tol = 1e-14
    n = 42
    x = gausslobatto(n)
    w = x[2]; x = x[1]
    @test ( length(x) == n) && ( length(w) == n )
    @test abs(x[37] - 0.922259214258616) < tol
    @test abs(w[37] - 0.029306411216166) < tol
    @test ( x[1] == -1 && x[n] == 1 )