text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6e3ef066-d115-11eb-2338-013a707dfe8a
begin
	using PlutoUI
	using DataFrames
	using Statistics
	using RollingFunctions
	using HTTP
	using JSON
	using Chain
	using Tables
	using PerceptualColourMaps
	using Gadfly
	Gadfly.set_default_plot_size(16cm, 14cm)
end

# ╔═╡ 8798b46c-1b0b-4d34-93ef-1d5062f8a632
md"""
## Summary
This report compares the performance of a standard A2C baseline against a "GPT-in-the-middle" variant. The architectures are as follows:

#### A2C

- Convolution with
  - output-size $32$
  - kernel-shape $8\times 8$
  - stride $4\times 4$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $4\times 4$
  - stride $2\times 2$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $3\times 3$
  - stride $1\times 1$
- ReLU
- Linear layer with output-size $512$
- Parallel value head and action head

#### GPT
(differences in **bold**)

- Convolution with
  - output-size $32$
  - kernel-shape $8\times 8$
  - stride $4\times 4$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $4\times 4$
  - stride $2\times 2$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $3\times 3$
  - stride $1\times 1$
- ReLU
- **Linear layer with output-size $1600$**
- **XL-size GPT torso**
- Parallel value head and action head
"""

# ╔═╡ 2954c622-ec01-4b1f-8f44-a48aa5413a16
md"""
## Results
The following graphs compare the performance of these two architectures on the Atari game Breakout.

The baseline uses 8 seeds using the best parameters from a random hyperparameter search.

The GPT uses the top 4 hyperparameter sets from our hyperparameter search.
"""

# ╔═╡ 5d3bfc57-6a5a-41e5-8f63-6902e4958936
begin
	struct HTMLDocument
		embedded
	end
	function Base.show(io::IO, mime::MIME"text/html", doc::HTMLDocument)
		println(io, "<html>")
		show(io, mime, doc.embedded)
		println(io, "</html>")
	end
end;

# ╔═╡ 5cba55d9-c7fe-4918-88b3-abbabe08c56b
function smooth_curves(df::AbstractDataFrame, y::Symbol, group_col::Symbol, window_size::Int)::DataFrame
	bounded_window(d) = max(1, min(window_size, nrow(d) ÷ 5))
	trim_window(d) = first(d, nrow(d) - mod(nrow(d), bounded_window(d)))
	@chain df begin
		dropmissing!([y, group_col])
		groupby(group_col)
		pairs
		[@chain g begin
				DataFrame
				sort!(:step)
				trim_window
				DataFrame(
					group_col => first(k),
					:step => round.(Int, rollmedian(_.step, bounded_window(_))),
					:mean => rollmean(_[:, y],  bounded_window(_)),
					:std => rollstd(_[:, y],  bounded_window(_)),
				)
				transform(
					[:mean, :std] => ((x, y) -> x - y) => [:lower],
					[:mean, :std] => ((x, y) -> x + y) => [:upper]
				)
				leftjoin(g; on=:step, makeunique=true)
			end for (k, g) in _]
		vcat(_...)
	end
end;

# ╔═╡ f7306428-2a04-43da-822e-85bed96216d8
md"""
## Analysis
The Abbeel paper ([Pretrained Transformers as Universal Computation Engines](https://arxiv.org/abs/2103.05247)) does not present state of the art performance after convergence. Instead it presents faster learning. It achieves this result by learning a very small number of parameters (relative to the baselines), forcing the networks to take advantage of information stored in the pretrained GPT network.

In contrast, our GPT model learns just as many parameters as the baseline and therefore it is not forced to take advantage of the GPT network. Instead the GPT network simply adds noise to the learning process.
"""

# ╔═╡ 5d99bf99-8973-4c4d-b88a-1e478b1678fb
md"""
### Analysis of Action Hidden Size

We explored a wide range of parameters including action hidden sizes (by this we mean then additional perceptron layer in the action head). We found that in general, layers of size $256$ or greater seemed to performed poorly on the first 30M steps of Breakout as the following graph indicates.
"""

# ╔═╡ f35a94e5-7568-4b89-af94-6167aebd3dab
md"(Group runs with error bands?)"

# ╔═╡ 3b614ff9-663e-459e-9621-98edd8922fdc
@bind collapse_runs CheckBox()

# ╔═╡ b721ca18-e5cd-4762-9ff0-9d94294b6969
md"""
###### Lit review for Deep Learning tricks

- [STABILIZING TRANSFORMERS FOR REINFORCEMENT LEARNING](https://arxiv.org/pdf/1910.06764.pdf) offers two tricks:
  - Replacing residual layers with gates. The reason this helps is not explained.
  - Moving the layer norm before the inputs (instead of after). This is supposed to help the agent learn a Markov policy before learning a history-dependent policy. Not relevant to what we are doing.

- [Pretrained Transformers As Universal Computation Engines](https://arxiv.org/pdf/2103.05247.pdf) had some additional tricks that we should try:
  - Fine-tuning the layer norm.
  - breaking up the input into patches and feeding in sequentially.
- [WHAT MATTERS FOR ON-POLICY DEEP ACTORCRITIC METHODS? A LARGE-SCALE STUDY](https://arxiv.org/pdf/2006.05990.pdf) suggests some good general techniques to use with actor-critic algorithms:
  - Use PPO loss instead of generic actor-critic
  - Use separate networks for value and policy. Try using the GPT only for value or only for policy.
  - Try Tanh instead of ReLU
  - Try normalizing observations and GPT outputs
  - Try using GAE
  - Tune gamma
  - Switch to Adam
"""

# ╔═╡ 731ba6f9-fa51-4fa3-a2dd-7748c9ecbea3
md"""
###### To do:
- Deep learning tricks
- Try minimal parameter learning
- Try additional Atari games
- Debug memory issues (get to 8)
  - https://github.com/pytorch/pytorch/issues/3022
  - Vivek
  - Zeyu
  - Laura
- Look into performance issues
- time forward and backward pass and compare against huggingface implementation
"""

# ╔═╡ 939f57a0-cb34-471b-aa24-10e5a95f410d
md"""
###### Open questions:
- [Do Frameskip versions work?](http://rldl12.eecs.umich.edu:8081/#sweeps/49)
- [Does normalize observation work?](http://rldl12.eecs.umich.edu:8081/#sweeps/50)
- [What are the best GPT parameters for Pong?](http://rldl12.eecs.umich.edu:8081/#sweeps/52)
"""

# ╔═╡ 41bd06a2-e7bd-46c6-9249-6e69223b0e11
begin
	HASRUA_ENDPOINT_URL = "http://rldl12.eecs.umich.edu:1200/v1/graphql"
	HASURA_ADMIN_SECRET = "hNuxFSTJMk28GMiFDpZDbviKtelFQamcb20UJiYZSIZ4LYsyLCCwQFYA6Y8HRSUo"
end;

# ╔═╡ 94185c13-b6d2-4337-b5ce-336f5e128032
function gql_query(query:: String; variables:: Dict = nothing)
	r = HTTP.request(
		"POST",
		HASRUA_ENDPOINT_URL;
		verbose=3,
		body= JSON.json(Dict("query" => query, "variables" => variables)) ,
		headers=Dict("x-hasura-admin-secret" => HASURA_ADMIN_SECRET),
	)
	@chain r begin
		_.body
		String
		JSON.parse(_, null=missing)
		_["data"]
	end
end;

# ╔═╡ 8f417ce0-fa72-4b12-8d3b-8f88b2e5a2bf
function latest_runs(n::Integer)::Vector{AbstractDict}
	query = """
	  query latestRuns(\$limit: Int!){
		run(order_by: {id: desc}, limit: \$limit) {
		  id
		  metadata
		  sweep{
	        id
		  }
		}
	  }
	  """
	gql_query(query; variables=Dict("limit" => n))["run"]
end;

# ╔═╡ 48be0808-f22a-4b72-a078-b2f2552db327
function format_logs(run::AbstractDict)
	sweep_id = ismissing(run["sweep"]) ? nothing : run["sweep"]["id"]
	@chain run["run_logs"] begin
		map(x -> Dict((Symbol(replace(k, " " => "_")) => v for (k,v) in x["log"])), _)
		DataFrame
		insertcols!(:run_id => run["id"], :sweep_id => sweep_id)
	end
end;

# ╔═╡ f4d4785e-0df1-4882-8eb9-005ff2a50614
function format_run(run::AbstractDict)
	sweep_id = something(get(run, "sweep", nothing), Dict("id" => nothing))["id"]
	data = run["metadata"]
	push!(data, pop!(data, "parameters")...)
	items = [Symbol(replace(k, " " => "_")) => v for (k, v) in data]
	push!(items, :run_id => run["id"], :sweep_id => sweep_id)
	Dict(items)
end;

# ╔═╡ 0ca92cdc-102c-45e7-a5b0-17087b9ba634
function format_runs(runs::Vector)::DataFrame
	rows = map(format_run, runs)
	cols = map(keys, rows) |> Iterators.flatten |> Set
	for r in rows
		for c in setdiff(cols, keys(r))
			r[c] = nothing
		end
	end
	df = DataFrame(rows)
	leading = [:run_id, :sweep_id]
	cols = filter(x -> !(x in leading), propertynames(df))
	df[!, vcat(leading, cols)]
end;

# ╔═╡ 68e80154-126b-403c-a593-d896cdf96ef2
function run_data(ids::AbstractVector{Int}, runs::AbstractDataFrame)
	query = """
		query getRuns(\$ids: [Int!]) {
		  run(where: {id: {_in: \$ids}}) {
			id
			run_logs{
			  log
			}
			metadata
			charts{
			  spec
			}
			sweep {
				id
			}
		  }
		}
  	"""
	@chain gql_query(query; variables=Dict("ids" => ids)) begin
			_["run"]
			format_logs.(_)
			DataFrame.(_)
			vcat(_...;cols=:union)
			leftjoin(runs[!, Not(:sweep_id)]; on=:run_id)
	end
end;

# ╔═╡ 0f54629c-e045-47bf-90a5-36e44c05b8f0
function sweep_runs(sweep_ids::AbstractVector{Int})
	query = """
		query getSweepRuns(\$ids: [Int!]) {
		  sweep(where: {id: {_in: \$ids}}) {
			runs {
           		id
			    metadata
			    sweep{
				 id
			    }
        	}
		  }
		}
  	"""
	rows = @chain gql_query(query; variables=Dict("ids" => sweep_ids)) begin
		_["sweep"]
		map(values, _)
		Iterators.flatten
		Iterators.flatten
		collect
	end
	format_runs(rows)
end;

# ╔═╡ 18147da2-1945-42e4-b140-c1375f46962b
function sweep_data(sweep_runs::AbstractDataFrame; ids = nothing)::DataFrame
	s_run_ids = sweep_runs.run_id |> unique |> collect
	if !isnothing(ids)
		s_run_ids = filter(in(ids), s_run_ids)
	end
	s_df = [run_data(collect(i), sweep_runs) for i in Iterators.partition(s_run_ids, 10)]
	s_df = vcat(s_df...; cols=:union)
	dropmissing!(s_df, :step)
end;

# ╔═╡ 1f0c13a6-ae73-4aec-bee3-745610c23175
begin
	gpt_plot_group = collapse_runs ? :sweep_id : :run_id
	plot(
		# smooth_curves(
		# 	sweep_data(sweep_runs([30, 31])),
		# 	:episode_return,
		# 	gpt_plot_group,
		# 	parse(Int, window_size)
		# ),
		sweep_data(sweep_runs([30, 32])),
		x=:step, y=:episode_return,
		# ymin=:lower, ymax=:upper,
		group=gpt_plot_group,
		color=:sweep_id,
		Guide.xlabel("Step"),
		Guide.ylabel("Episode Return"),
		Geom.line,
		# Geom.ribbon,
		Scale.color_discrete(),
		Guide.colorkey(labels=["Yes", "No"], title="GPT?"),
		Guide.title("GPT vs. No GPT"),
		alpha=[0.5]
	) |> HTMLDocument
end

# ╔═╡ 0135af1a-d8c8-4027-af38-7f57bbfef282
gpt_runs = sweep_runs([26, 23, 22, 19, 17]);

# ╔═╡ 3ba73b7c-f5dc-4639-ab9c-c43cc992634f
gpt_df = sweep_data(gpt_runs);

# ╔═╡ bbed8889-d401-4167-9370-9927aabbc83b
# @bind window_size Select(string.(collapse_runs ? (100:100:500) : (1:10)), default=collapse_runs ? "200" : "5")
window_size = collapse_runs ? "500" : "11";

# ╔═╡ fe6f09ca-dc39-4d77-aae1-757c87ccd3e7
begin
	action_plot_group = collapse_runs ? :action_hidden_size : :run_id
	action_sizes = gpt_df.action_hidden_size |> unique |> sort
	action_colors(n) = cmap("L08"; N=n)
	plot(
		smooth_curves(
			gpt_df,
			:episode_return,
			action_plot_group,
			parse(Int, window_size)
		),
		x=:step, y=:mean,
		ymin=:lower, ymax=:upper,
		group=action_plot_group,
		color=:action_hidden_size,
		Guide.xlabel("Step"),
		Guide.ylabel("Episode Return"),
		Geom.line,
		Geom.ribbon,
		Scale.color_discrete(action_colors; levels=action_sizes),
		Guide.colorkey(title="Action Hidden Size"),
		Guide.title("Performance Across Action Sizes"),
		alpha=[0.5]
	) |> HTMLDocument
end

# ╔═╡ Cell order:
# ╟─8798b46c-1b0b-4d34-93ef-1d5062f8a632
# ╟─2954c622-ec01-4b1f-8f44-a48aa5413a16
# ╟─5d3bfc57-6a5a-41e5-8f63-6902e4958936
# ╟─5cba55d9-c7fe-4918-88b3-abbabe08c56b
# ╟─1f0c13a6-ae73-4aec-bee3-745610c23175
# ╟─f7306428-2a04-43da-822e-85bed96216d8
# ╟─5d99bf99-8973-4c4d-b88a-1e478b1678fb
# ╟─fe6f09ca-dc39-4d77-aae1-757c87ccd3e7
# ╟─f35a94e5-7568-4b89-af94-6167aebd3dab
# ╟─3b614ff9-663e-459e-9621-98edd8922fdc
# ╟─b721ca18-e5cd-4762-9ff0-9d94294b6969
# ╠═731ba6f9-fa51-4fa3-a2dd-7748c9ecbea3
# ╠═939f57a0-cb34-471b-aa24-10e5a95f410d
# ╟─6e3ef066-d115-11eb-2338-013a707dfe8a
# ╟─41bd06a2-e7bd-46c6-9249-6e69223b0e11
# ╟─94185c13-b6d2-4337-b5ce-336f5e128032
# ╟─8f417ce0-fa72-4b12-8d3b-8f88b2e5a2bf
# ╟─48be0808-f22a-4b72-a078-b2f2552db327
# ╟─f4d4785e-0df1-4882-8eb9-005ff2a50614
# ╟─0ca92cdc-102c-45e7-a5b0-17087b9ba634
# ╟─68e80154-126b-403c-a593-d896cdf96ef2
# ╟─0f54629c-e045-47bf-90a5-36e44c05b8f0
# ╟─18147da2-1945-42e4-b140-c1375f46962b
# ╟─0135af1a-d8c8-4027-af38-7f57bbfef282
# ╟─3ba73b7c-f5dc-4639-ab9c-c43cc992634f
# ╟─bbed8889-d401-4167-9370-9927aabbc83b

max_stars_repo_name: GPT-RL/generalization

max_issues_repo_name: GPT-RL/generalization

max_forks_repo_name: GPT-RL/generalization


text:
using JLD
using Printf

scratch_dir = ARGS[1]

for set in ["set1", "set2", "set3", "set4"]
    power = zeros(2, length(-0.75:0.05:0.75))

    i = 0
    for del in -0.75:0.05:0.75
        i += 1
        nrep = 0
        for batch in 1:10
            try
                file = jldopen("$(scratch_dir)/res_$(set)_$(del)_$(batch).jld", "r")
                T1 = read(file, "T1")
                T2 = read(file, "T2")
                close(file)

                power[1, i] += sum(T1)
                power[2, i] += sum(T2)

                nrep += 1
            catch
                @warn "could not open $(set) $(del) $(batch)"
            end
        end
	    power[1, i] = power[1, i] / (nrep * 100)
	    power[2, i] = power[2, i] / (nrep * 100)
    end

    @save "power_$(set).jld" power
end

ext: jl


text:
using Knet
using ArgParse
using JSON
using MAT
using JLD

include("util.jl")
include("model.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "CNN-HMM hybrid ASR system in Knet by Ilker Kesen."

    @add_arg_table s begin
        ("--features"; required=true; help="MFSC features file")
        ("--jsonfile"; required=true)
        ("--loadfile"; default=nothing)
        ("--savefile"; default=nothing)
        ("--hidden"; default=128; arg_type=Int64)
        ("--epochs"; default=100; arg_type=Int64)
        ("--batchsize"; default=200; arg_type=Int64)
        ("--atype"; default=(gpu()>=0?"KnetArray{Float32}":"Array{Float32}"))
        ("--optim"; default="Adam(;gclip=5.0)")
        ("--pdrop"; default=0.0; arg_type=Float64)
        ("--seed"; default=-1; arg_type=Int64)
        ("--nvalid"; default=200; arg_type=Int64)
        ("--generate"; action=:store_true)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:atype] = eval(parse(o[:atype]))
    sr = o[:seed] > 0 ? srand(o[:seed]) : srand()

    # load data
    jsondata = JSON.parsefile(abspath(o[:jsonfile]))

    # TODO: data loading
    p2i, i2p = build_vocab(jsondata)
    @time trn = make_data(o[:features], jsondata, "train")
    shuffle!(trn)
    val = []
    for k = 1:o[:nvalid]
        push!(val, pop!(trn))
    end
    trn = reduce(vcat, trn)
    val = reduce(vcat, val)

    # load model
    w = opts = bestacc = nothing
    if o[:loadfile] == nothing
        w = initweights(o[:atype], o[:hidden])
        bestacc = -Inf
    else
        w = load(o[:loadfile], "w")
        bestacc = load(o[:loadfile], "acc")
        w = map(wi->convert(o[:atype],wi), w)
    end
    opts = map(wi->eval(parse(o[:optim])), w)

    if o[:generate]
        tst = make_data(o[:features], jsondata, "test")
        tst = reduce(vcat, tst)
        @time lossval, acc = validate(w,val,p2i)
        println("lossval:$lossval,acc:$acc")
        return
    end

    bestacc = -Inf
    println("training has been started. ", now())
    for epoch = 1:o[:epochs]
        # shuffle train data
        shuffle!(trn)

        # one epoch training
        losstrn = 0
        @time for k = 1:o[:batchsize]:length(trn)
            samples = trn[k:min(length(trn),k+o[:batchsize]-1)]
            x = make_input(samples)
            ygold = make_output(samples,p2i)
            this_loss = train!(w,convert(o[:atype],x),ygold,opts;o=o)
            losstrn += this_loss
        end
        losstrn = losstrn / length(trn)

        # compute validation loss and accuracy
        @time lossval, acc = validate(w,val,p2i; o=o)

        # report and save
        print("(epoch:$epoch,losstrn:$losstrn,lossval:$lossval,acc:$acc) ")
        print(now())
        if acc > bestacc
            bestacc = acc
            save(o[:savefile], "w", map(wi->Array(wi),w), "acc", acc)
            print(" model saved!")
        end
        println()
    end
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/phoneme-convnet

max_issues_repo_name: ilkerkesen/phoneme-convnet

max_forks_repo_name: ilkerkesen/phoneme-convnet


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1c3ff5c4-90d5-11eb-1459-7d7d43335a80
using Yao, YaoPlots, YaoExtensions, PlutoUI, Flux, Plots

# ╔═╡ e9342104-91cf-11eb-01d7-ffd747975c7c
using LinearAlgebra

# ╔═╡ aa09eef8-8eac-11eb-00e6-6794d5c9e0b5
md"""
# Quantum Computing Basics

In case you are not familar with some basic concept of quantum computing or quantum algorithms, we will have a quick and short introduction for quantum computing.
"""

# ╔═╡ 22f46d6e-90d5-11eb-3094-e7015012995c
md"""

## Qubits

The computation states are often represented as qubits in quantum computing model. A qubit is a quantum state that has two basis states ``|0⟩`` and ``|1⟩``

## Quantum Registers

A quantum register is an abstraction of the qubit memory, it usualy contains many qubits.

## Quantum Gates

A quantum gate is a unitary operator that can be applied on the (quantum) register. we usually visualize it as a square block with two wires, e.g the following is a Pauli X gate.
"""

# ╔═╡ 920da102-90d5-11eb-0437-495807d8b345
YaoPlots.plot(X)

# ╔═╡ a5d5d56a-90d5-11eb-1e4a-8bf893d8534a
md"""
## Parameterized Quantum Gates

There are some quantum gates that can take classical parameters, e.g the rotation gates on X, Y, Z or other directions.
"""

# ╔═╡ be1f51c6-90d5-11eb-3db1-99708f112467
YaoPlots.plot(Rx(0.1))

# ╔═╡ 4f4a700e-90d6-11eb-3058-b1809082853a
md"""
## Controlled Gate

The gate operation can be controlled by another qubits, we usually represent them as the following diagram
"""

# ╔═╡ 6b2079e2-90d6-11eb-330c-0194f2906114
YaoPlots.plot(control(2, 1, 2=>Rx(0.1)))

# ╔═╡ 7854398c-90d6-11eb-2aad-b75727ab7b79
md"for controlled Pauli X gate, we also call it the CNOT gate, it can be plotted in a simpler way"

# ╔═╡ 9003bb5c-90d6-11eb-102e-f1954fe71ec6
YaoPlots.plot(cnot(2, 1, 2))

# ╔═╡ debe78d2-90d5-11eb-1526-87f514b030bf
md"""
## Quantum Circuits & Quantum Algorithms

Now with the concept of register and gates, we can define quantum circuits. The quantum circuits are composition of gates and measurements.

e.g the following is a quantum fourier transformation circuit, it is defined using hadamard gate ``H`` and control phase shift gate.
"""

# ╔═╡ 0a379552-90d6-11eb-1a07-130ebf3c1264
begin
	A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1))))
	B(n, k) = chain(n, j==k ? put(k=>H) : A(j, k) for j in k:n)
	qft(n) = chain(B(n, k) for k in 1:n)
	YaoPlots.plot(qft(3))
end

# ╔═╡ e4f1f268-91ca-11eb-231a-712cb8cf4a55
md"""
## Program with Yao

checkout the website: [yaoquantum.org](https://yaoquantum.org/)
"""

# ╔═╡ 0478efb2-91cb-11eb-034e-0b0c3be90bb7
html"""
<iframe src="https://yaoquantum.org/" width=100% height=1000></iframe>
"""

# ╔═╡ 2a546ef0-91cb-11eb-17fd-4b8c74878e6a
md"""
### the quantum blocks

Yao represents the quantum circuits as blocks

![](https://github.com/Roger-luo/PIML/blob/master/notebooks/assets/qblock.png?raw=true)
"""

# ╔═╡ 76b6bf64-91cb-11eb-3c39-576a66e70cf4
qft(3)

# ╔═╡ c75b0244-8eac-11eb-2528-b9f59f61634d
md"""
## Variational Quantum Algorithms

*in some context, this is also called quantum machine learning.*

Traditionally, a quantum algorithm usually requires a well designed circuits that provides an exponential speed up comparing to its classical counterparts, e.g the well known Shor algorithm.

However, designing and finding such algorithms can be extremely hard, and may not be possible implemented on near term quantum devices that is noisy.
"""

# ╔═╡ e1345ee8-8eac-11eb-3c22-0db185142833
md"""
## Differentiating Quantum Circuits: the forward mode

In literature, this is usually refered as the **Parameter Shift Rule**, it was first discovered in [Quantum Circuit Learning, arXiv: 1803.00745](https://arxiv.org/pdf/1803.00745.pdf)

we will now derive the rules, assume the unitary of our parameterized quantum circuit is ``U(\mathbf{θ}) = \prod_{i=1}^n U(θ_i)``, and for convenience, we will denote ``U_{j:k} = U_j⋯U_k``

then for an observable ``B``, the expectation value is ``⟨B(\mathbf{θ})⟩ = Tr(B U_{l:1}ρU_{l:1}^{†})``

now if we assume our unitary is generated by Pauli strings (which is not always true), we have ``U_j(θ) = exp(-iθ_j P_j/2)``.

The gradient of ``⟨B(\mathbf{θ})⟩`` can be calculated as following

```math
\frac{\partial ⟨B⟩}{\partial θ_j} = -\frac{i}{2} (B U_{l:j} [P_j, U_{j-1:1}ρU_{j-1:1}^{†}] U_{l:j}^{†})
```

and because we have the following relationship for arbitrary ``ρ``

```math
[P_j, ρ] = i[U_j(\frac{π}{2})ρ U_j^{†}(\frac{π}{2}) - U_j(-\frac{π}{2})ρ U_j(-\frac{π}{2})]
```

thus the gradient can be evaluated by

```math
\begin{aligned}
\frac{\partial B}{\partial θ_j} =& \frac{1}{2} (Tr[B U_{l:j+1}U(θ_j+\frac{π}{2})U_{j-1:1}ρU_{j-1:1}^{†}U_j^{†}(θ_j+\frac{π}{2})U_{l:j+1}^{†}] -\\
&Tr[B U_{l:j+1}U(θ_j-\frac{π}{2})U_{j-1:1}ρU_{j-1:1}^{†}U_j^{†}(θ_j-\frac{π}{2})U_{l:j+1}^{†}])
\end{aligned}
```

which means

```math
\frac{\partial B}{\partial θ_j} = \frac{1}{2} (⟨B(\theta_j + \frac{π}{2})⟩ - ⟨B(\theta_j - \frac{π}{2})⟩)
```
"""

# ╔═╡ f8270b02-8eac-11eb-2ab3-3323e01f4f62
md"""
## Differentiating Quantum Circuits: the reverse mode

is it possible to differentiate quantum circuits using reverse mode?
"""

# ╔═╡ c97fffa8-91d0-11eb-0a0c-092545e368f8
md"""
Only for simulation

![](https://github.com/Roger-luo/PIML/blob/master/notebooks/assets/yao-reverse-ad.png?raw=true)
"""

# ╔═╡ 3c919a50-8ead-11eb-058f-27eceb1deeaa
md"""
## Implementing A Simple Variational Quantum Eigensolver (VQE)

You have learned how variational Monte Carlo works, now let's:

1. use a quantum circuit as our ansatz
2. directly calculate the energy by measure the Hamiltonian
3. use gradient based optimization to optimize the energy
"""

# ╔═╡ 68e75bba-91ca-11eb-128f-1d6b0b7cbb78
md"first let's create a variational circuit using the `variational_circuit` function"

# ╔═╡ 4422a0fa-91ca-11eb-3bda-f37b2dab68b3
YaoPlots.plot(variational_circuit(5, 3))

# ╔═╡ 82f2d962-91ca-11eb-2e66-597a90047e87
md"then create a quantum register at state ``|00⋯00⟩``"

# ╔═╡ 895600fc-91ca-11eb-3d03-e5aedfa44576
r = zero_state(5)

# ╔═╡ 940657d8-91ca-11eb-33fd-63a019b4ebfc
md"we can evaluate this circuit on the emulator"

# ╔═╡ a0527026-91ca-11eb-0604-21a298350c8d
r |> variational_circuit(5, 3)

# ╔═╡ b8a71bf4-91ca-11eb-2f81-eb67cfeeada0
md"""
we can now create a 1-D heisenberg hamiltonian

```math
\sum_i \mathbf{S}_i \mathbf{S}_{i+1}
```
"""

# ╔═╡ da9ba108-91ca-11eb-28a4-755cf76ea605
h = heisenberg(5)

# ╔═╡ b3cc56a2-91cb-11eb-3736-87df27d04545
md"we can calculate the expectation value of this hamiltonian"

# ╔═╡ 848d146c-91cb-11eb-17cc-87a1aa3c0d76
expect(heisenberg(5), r)

# ╔═╡ bd3c4bca-91cb-11eb-29cf-657eb95db73d
@doc expect

# ╔═╡ 9fe3db9c-91cb-11eb-1888-cf7dd2aa3aa1
expect(heisenberg(5), r=>variational_circuit(5, 3))

# ╔═╡ a6c13932-91cb-11eb-3822-7da569e520e0
reg, ∇θ = expect'(heisenberg(5), zero_state(5)=>variational_circuit(5, 3))

# ╔═╡ cbacc1f8-91cb-11eb-316c-37e3da59855c
md"we can use the optimizers from the ML framework `Flux`"

# ╔═╡ dddcbf20-91cd-11eb-2d73-2daa1643074d
opt = ADAM()

# ╔═╡ 18d83a62-91d0-11eb-2d8d-2bc2e69785e8
md"depth=$(@bind depth Slider(1:10; show_value=true))"

# ╔═╡ 12aff314-91ce-11eb-36d3-772d32cb8f3e
circuit = variational_circuit(5, depth)

# ╔═╡ e172b6f0-91cd-11eb-35be-2d44f6f26a80
begin
	history = []
	θ = rand(nparameters(circuit))
	for _ in 1:1000
		reg, ∇θ = expect'(heisenberg(5), zero_state(5)=>dispatch!(circuit, θ))
		Flux.Optimise.update!(opt, θ, ∇θ)
		current = expect(heisenberg(5), zero_state(5)=>dispatch!(circuit, θ))
		push!(history, real(current))
	end
	Plots.plot(history)
end

# ╔═╡ ec0935e0-91cf-11eb-1837-334e0e92edee
minimum(eigvals(Matrix(heisenberg(5))))

# ╔═╡ ad0071da-8ead-11eb-31e8-b94b8bff7635
md"""
## Placement Ads: Google Summer of Code 2022

https://summerofcode.withgoogle.com/
"""

# ╔═╡ Cell order:
# ╠═1c3ff5c4-90d5-11eb-1459-7d7d43335a80
# ╟─aa09eef8-8eac-11eb-00e6-6794d5c9e0b5
# ╟─22f46d6e-90d5-11eb-3094-e7015012995c
# ╟─920da102-90d5-11eb-0437-495807d8b345
# ╟─a5d5d56a-90d5-11eb-1e4a-8bf893d8534a
# ╟─be1f51c6-90d5-11eb-3db1-99708f112467
# ╟─4f4a700e-90d6-11eb-3058-b1809082853a
# ╟─6b2079e2-90d6-11eb-330c-0194f2906114
# ╟─7854398c-90d6-11eb-2aad-b75727ab7b79
# ╟─9003bb5c-90d6-11eb-102e-f1954fe71ec6
# ╟─debe78d2-90d5-11eb-1526-87f514b030bf
# ╟─0a379552-90d6-11eb-1a07-130ebf3c1264
# ╟─e4f1f268-91ca-11eb-231a-712cb8cf4a55
# ╟─0478efb2-91cb-11eb-034e-0b0c3be90bb7
# ╟─2a546ef0-91cb-11eb-17fd-4b8c74878e6a
# ╠═76b6bf64-91cb-11eb-3c39-576a66e70cf4
# ╟─c75b0244-8eac-11eb-2528-b9f59f61634d
# ╟─e1345ee8-8eac-11eb-3c22-0db185142833
# ╟─f8270b02-8eac-11eb-2ab3-3323e01f4f62
# ╟─c97fffa8-91d0-11eb-0a0c-092545e368f8
# ╟─3c919a50-8ead-11eb-058f-27eceb1deeaa
# ╟─68e75bba-91ca-11eb-128f-1d6b0b7cbb78
# ╠═4422a0fa-91ca-11eb-3bda-f37b2dab68b3
# ╟─82f2d962-91ca-11eb-2e66-597a90047e87
# ╠═895600fc-91ca-11eb-3d03-e5aedfa44576
# ╟─940657d8-91ca-11eb-33fd-63a019b4ebfc
# ╠═a0527026-91ca-11eb-0604-21a298350c8d
# ╟─b8a71bf4-91ca-11eb-2f81-eb67cfeeada0
# ╠═da9ba108-91ca-11eb-28a4-755cf76ea605
# ╟─b3cc56a2-91cb-11eb-3736-87df27d04545
# ╠═848d146c-91cb-11eb-17cc-87a1aa3c0d76
# ╠═bd3c4bca-91cb-11eb-29cf-657eb95db73d
# ╠═9fe3db9c-91cb-11eb-1888-cf7dd2aa3aa1
# ╠═a6c13932-91cb-11eb-3822-7da569e520e0
# ╟─cbacc1f8-91cb-11eb-316c-37e3da59855c
# ╠═dddcbf20-91cd-11eb-2d73-2daa1643074d
# ╠═18d83a62-91d0-11eb-2d8d-2bc2e69785e8
# ╠═12aff314-91ce-11eb-36d3-772d32cb8f3e
# ╠═e172b6f0-91cd-11eb-35be-2d44f6f26a80
# ╠═e9342104-91cf-11eb-01d7-ffd747975c7c
# ╠═ec0935e0-91cf-11eb-1837-334e0e92edee
# ╟─ad0071da-8ead-11eb-31e8-b94b8bff7635

max_stars_repo_name: Roger-luo/PIML

max_issues_repo_name: Roger-luo/PIML

max_forks_repo_name: Roger-luo/PIML


text:
# This is the precompilation package
import Plots, CSV, StatsPlots, Markdown, DataFrames, Query, ColorSchemes, StatsBase, Random, Distributions

p = Plots.plot()
fname = "$(tempname()).png"
Plots.savefig(p, fname)

try
    CSV.File("a")
catch
end

ColorSchemes.acton.colors

DataFrames.DataFrame()

dist = Distributions.Gamma(2)
p = StatsPlots.scatter(dist, leg = false)

StatsBase.AnalyticWeights([0.2, 0.1, 0.3])

Markdown.md"a"

nothing

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
using LinearAlgebra
using Statistics
using PyCall

#I'm just used to PyPlot, feel free to change to Plots
using PyPlot

using NCFlex
using Optim
using LineSearches
using JuLIP
using ASE
using FiniteDifferences
using Roots

using StaticArrays
using SparseArrays
using Printf

using Parameters, Setfield
using BifurcationKit
const BK = BifurcationKit

@pyimport ase.units as ase_units
@pyimport ase.lattice.cubic as ase_lattice_cubic
@pyimport matscipy.fracture_mechanics.crack as crack
@pyimport matscipy.fracture_mechanics.clusters as clusters
@pyimport matscipy.elasticity as elasticity

include("tungsten_MB.jl")
include("NCF_JuLIP_setup.jl")

ncf_t = NCF_1(clust);
N1 = ncf_t.at["N1"]

x0 = copy(dofs(ncf_t.at));
Ubar0 = zeros(3*N1)

### example of a plot:
# plot_at_ncf(ncf_t,Ubar0,1.0,aa)
# plot(aa,0.0,"*")

####################
# Preparation Part 1:
#
# a clean NCF struct just in case:
ncf_t = NCF_1(clust)

# runs the CLE-only estimation by finding K for a fixed alpha
# each alpha is taken from αss, which is formed as
# αss = (a1,a2,length=ll)
# if we have no clue about K, best to set a1,a2 far apart and ll large
# but this can be very slow
# here we do not have to do this step as manually found the good K, so I just set ll=2
# can be commented out too
kss, αss = preparation_part_1(ncf_t,ll = 2,a1=-0.6,a2=-0.5)

####################
# Preparation Part 2:
#
# a clean NCF struct just in case:
ncf_t = NCF_1(clust);

# fix some alpha and hope for the best
# if this one does not work, change it and re-run
a_st = -0.2

# typically k_st is taken as a member of kss e.g.
#k_st = minimum(kss)

# here is a manually-found near perfect k_st:
k_st = 0.8877816860333725

#using k_st and a_st, we find the first static equilibrium using a hessian-free Optim algorithm:
xbar1_1 = preparation_part_2(ncf_t,k_st,a_st)
# to speed things up Optim algorithm has very large error tolerance, so best to 'top it up' with a quick Newton solver:
xbar1_2 = Newton_static_1(ncf_t,xbar1_1,k_st,a_st; show = true);
# copying just in case
xbar_c = copy(xbar1_2)

# so now we have first static equilibrium triplet (x_bar1_2,k_st,a_st)
# but not quite giving us f_alpha = 0.0:
@show ncf_t.g(xbar1_2,k_st,a_st)[end]

# we can "manually" make it a bit better by finding a second static equilibrium triplet is (xbar_b,kk,a_st):
kk = k_st + 0.000206
xbar_b, gg_b = find_k(ncf_t,xbar1_2,kk,a_st)

# this is now close enough to what we want that we can safely call a Roots.find_zero routine, which can be done as follows:
K_best, xbar_best = preparation_part_3(ncf_t,kk-0.001,kk+0.001,xbar_b,a_st)

# as a result we have our first flex equilibrium triplet (xbar_best,K_best,a_st):
@show ncf_t.g(xbar_best,K_best,a_st)[end];

#we can further run Newton_flex:
kbar_fin = K_best
xbar_fin = xbar_best
abar_fin = a_st
flex1 = Newton_flex_1(ncf_t,[xbar_fin;abar_fin],kbar_fin;show=true)

#and just in case find a second flex equilibrium a bit away to make kick-starting the continuation easier:
ee = 0.0001
flex2 = Newton_flex_1(ncf_t,flex1.+ee,kbar_fin+ee;show=true)
;

######################
# Continuation:

#clean NCF struct just in case:
ncf_t = NCF_1(clust)

#populate it with the first two flex solutions:
Us = []
push!(Us,flex1[1:end-1])
ncf_t.Us = Us
ncf_t.Ks = [kbar_fin]
ncf_t.αs = [flex1[end]]

push!(ncf_t.Us,flex2[1:end-1])
push!(ncf_t.Ks,ncf_t.Ks[1]+ee)
push!(ncf_t.αs,flex2[end])

# check that we have two data points:
@show ncf_t.Ks
@show ncf_t.αs

# and off we go:
simple_continuation_1(ncf_t, theta=0.5, maxSteps=10,
        dsmin = 0.002, dsmax = 0.002, ds= 0.002,
        tangentAlgo = BorderedPred(),dotPALC = (x, y) -> dot(x,y)/length(x),
        linsolver = GMRESKrylovKit(),linearAlgo = BorderingBLS())

#some comments on the function:
# this is a custom-made routine based on BifurcationKit.jl which stores every data point
# and manually halves the ds near fold points to make sure we capture them

# dsmin, dsmax - minimal/maximal ds the algorithm will take in
# ds - the one we start with
# the rest is quite standard

# it can be called on repeatedly and will just pick up where it left

#using PyPlot we can then create a handy figure:
figure(figsize=(10,4))

iii = length(ncf_t.Ks) - 2
subplot(121)
plot(ncf_t.Ks,ncf_t.αs,".-")
plot(ncf_t.Ks[iii],ncf_t.αs[iii],"r*")

subplot(122)
plot_at_ncf(ncf_t,ncf_t.Us[iii],ncf_t.Ks[iii],ncf_t.αs[iii])

plot(ncf_t.αs[iii],0.0,"*")


text:
# This file is not part of the automated test suite, as it requires Docker
# to be present on the system.

using Test, Joseki, JSON

@testset "Test sample dockerfile" begin
    @info "Starting server from included dockerfile"
    cd(dirname(@__FILE__))
    cd("../")
    run(`docker build -t joseki .`)
    run(`docker run -d --rm -p 8000:8000 --name joseki-test joseki`)
    sleep(5)

    res = HTTP.get("http://localhost:8000/pow/?x=2&y=3")
    @test res.status == 200
    @test HTTP.hasheader(res, "Access-Control-Allow-Methods", "GET,PUT,POST,DELETE,PATCH,OPTIONS")

    @info "Shutting down dockerfile server"
    run(`docker rm -f joseki-test`)
end


text:
using Test
using ConstraintSolver
using JSON

CS = ConstraintSolver

include("sudoku_fcts.jl")
include("small_special_tests.jl")
include("sudoku_tests.jl")
include("killer_sudoku_tests.jl")
include("graph_color_tests.jl")


text:
export zoo_cnf, zoo_dnf,
    zoo_cnf_file, zoo_dnf_file,
    FnfFormat

struct FnfFormat <: FileFormat end

##############################################
# Read CNF/DNF
##############################################

"""
    zoo_cnf(name)

Loads CNF file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_cnf(name) =
    read(zoo_cnf_file(name), LogicCircuit, FnfFormat())::Plain⋀Node

"""
    zoo_dnf(name)

Loads DNF file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_dnf(name) =
    read(zoo_dnf_file(name), LogicCircuit, FnfFormat())::Plain⋁Node

zoo_cnf_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/cnfs/$name"

zoo_dnf_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/dnfs/$name"

const fnf_grammar = raw"""
    start: header (_NL clause)+ _NL?

    header : "p" _WS (cnf | dnf) _WS INT _WS INT

    cnf : "cnf"
    dnf : "dnf"

    clause : (SIGNED_INT _WS)+ "0"

    TRAILING_COMMENT : "%" (/./|_NL)* /$/
    %ignore TRAILING_COMMENT

    %import common.INT
    %import common.SIGNED_INT
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const fnf_parser = Lark(fnf_grammar)

struct FnfParse <: JuiceTransformer
    lit_cache::Dict{Lit,PlainLiteralNode}
    FnfParse() = new(Dict{Lit,PlainLiteralNode}())
end

@rule header(t::FnfParse, x) =
    (x[1], Base.parse(Int,x[2]), Base.parse(Int,x[3]))

@rule cnf(t::FnfParse, _) = true
@rule dnf(t::FnfParse, _) = false

@rule start(t::FnfParse, x) = begin
    polarity, numvars, numclauses = x[1]
    clauses = x[2:end]
    @assert numclauses == length(clauses)
    circuit = if polarity
        Plain⋀Node(map(Plain⋁Node, clauses))
    else
        Plain⋁Node(map(Plain⋀Node, clauses))
    end
    @assert num_variables(circuit) == numvars
    circuit
end

@rule clause(t::FnfParse, x) = begin
    map(x) do l
        literal = Base.parse(Lit,l)
        get!(t.lit_cache, literal) do
            PlainLiteralNode(literal)
        end
    end
end

function Base.parse(::Type{PlainLogicCircuit}, str, ::FnfFormat)
    ast = Lerche.parse(fnf_parser, str)
    Lerche.transform(FnfParse(), ast)
end

"""
    Base.read(io::IO, ::Type{PlainLogicCircuit}, ::FnfFormat)

Read CNF/DNF from file.
"""
Base.read(io::IO, ::Type{PlainLogicCircuit}, ::FnfFormat) =
    parse(PlainLogicCircuit, read(io, String), FnfFormat())

##############################################
# Write CNF/DNF
##############################################

const FNF_FORMAT = """c this file was saved by LogicCircuits.jl
c literals are numbered beginning at 1
c positive/negative number indicates positive/negative literal
c file syntax:
c p cnf num-vars num-clauses
c {literal-in-clause}* 0
c"""

"""
    Base.write(io::IO, fnf::LogicCircuit, ::FnfFormat)

Write CNF/DNF to file.
"""
function Base.write(io::IO, fnf::LogicCircuit, ::FnfFormat)
    @assert isflat(fnf)
    println(io, FNF_FORMAT)
    polarity = iscnf(fnf) ? "cnf" : "dnf"
    println(io, "p $polarity $(num_variables(fnf)) $(num_children(fnf))")
    for clause in children(fnf)
        for lit_node in children(clause)
            print(io, "$(literal(lit_node)) ")
        end
        println(io, "0")
    end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
#
# This difference is important because for Orthorhombic cells it is
# possible to run over only half of the cells, and wrapping coordinates
# in Orthorhombic cells is slightly cheaper.
#
struct TriclinicCell end
struct OrthorhombicCell end

struct UnitCell{UnitCellType,N,T,M}
  matrix::SMatrix{N,N,T,M}
end

"""

$(TYPEDEF)

$(TYPEDFIELDS)

Structure that contains some data required to compute the linked cells. To
be initialized with the box size and cutoff.

## Examples

```julia-repl
julia> sides = [250,250,250];

julia> cutoff = 10;

julia> box = Box(sides,cutoff)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [250.0 0.0 0.0; 0.0 250.0 0.0; 0.0 0.0 250.0]
  cutoff: 10.0
  number of computing cells on each dimension: [27, 27, 27]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 19683

julia> box = Box([ 10  0  0
                    0 10  5
                    0  0 10 ], 1)
Box{TriclinicCell, 3, Float64, 9}
unit cell matrix: [10.0 0.0 0.0; 0.0 10.0 5.0; 0.0 0.0 10.0]
cutoff: 1.0
number of computing cells on each dimension: [12, 17, 12]
computing cell sizes: [1.0, 1.0, 1.0] (lcell: 1)
Total number of cells: 2448

```

"""
Base.@kwdef struct Box{UnitCellType,N,T,M}
  unit_cell::UnitCell{UnitCellType,N,T,M}
  lcell::Int
  nc::SVector{N,Int}
  cutoff::T
  cutoff_sq::T
  ranges::SVector{N,UnitRange{Int}}
  cell_size::SVector{N,T}
  unit_cell_max::SVector{N,T}
end

"""

```
Box(
  unit_cell_matrix::AbstractMatrix,
  cutoff,
  T::DataType,
  lcell::Int=1,
  UnitCellType=TriclinicCell
)
```

Construct box structure given the cell matrix of lattice vectors. This
constructor will always return a `TriclinicCell` box type, unless the
`UnitCellType` parameter is set manually to `OrthorhombicCell`

### Example
```julia
julia> unit_cell = [ 100   50    0
                       0  120    0
                       0    0  130 ];

julia> box = Box(unit_cell,10)
Box{TriclinicCell, 3, Float64, 9}
  unit cell matrix: [100.0 50.0 0.0; 0.0 120.0 0.0; 0.0 0.0 130.0]
  cutoff: 10.0
  number of computing cells on each dimension: [17, 14, 15]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 3570

```

"""
function Box(
  unit_cell_matrix::AbstractMatrix,
  cutoff,
  T::DataType,
  lcell::Int=1,
  UnitCellType=TriclinicCell
)

  s = size(unit_cell_matrix)
  unit_cell_matrix = SMatrix{s[1],s[2],Float64,s[1]*s[2]}(unit_cell_matrix)

  @assert lcell >= 1 "lcell must be greater or equal to 1"

  N = size(unit_cell_matrix)[1]
  @assert N == size(unit_cell_matrix)[2] "Unit cell matrix must be square."
  @assert check_unit_cell(unit_cell_matrix,cutoff) " Unit cell matrix does not satisfy required conditions."

  unit_cell = UnitCell{UnitCellType,N,T,N*N}(SMatrix{N,N,T,N*N}(unit_cell_matrix))
  unit_cell_max = sum(@view(unit_cell_matrix[:,i]) for i in 1:N)

  nc_min = @. floor.(Int,unit_cell_max/cutoff)
  cell_size_max = unit_cell_max ./ nc_min
  cell_size = cell_size_max/lcell

  nc = SVector{N,Int}(ceil.(Int,unit_cell_max ./ cell_size) .+ 2*lcell)

  #ranges = ranges_of_replicas(cell_size, lcell, nc, unit_cell_matrix)
  ranges = SVector{N,UnitRange{Int}}(ntuple(i->-1:1,N))
  return Box{UnitCellType,N,T,N*N}(
    unit_cell,
    lcell,
    nc,
    cutoff,
    cutoff^2,
    ranges,
    cell_size,
    unit_cell_max
  )
end
Box(
  unit_cell_matrix::AbstractMatrix,
  cutoff;
  T::DataType=Float64,
  lcell::Int=1,
  UnitCellType=TriclinicCell
) = Box(unit_cell_matrix,cutoff,T,lcell,UnitCellType)

function Base.show(io::IO,::MIME"text/plain",box::Box)
  println(typeof(box))
  println("  unit cell matrix: ", box.unit_cell.matrix)
  println("  cutoff: ", box.cutoff)
  println("  number of computing cells on each dimension: ",box.nc)
  println("  computing cell sizes: ", box.cell_size, " (lcell: ",box.lcell,")")
  print("  Total number of cells: ", prod(box.nc))
end

"""

```
Box(
  sides::AbstractVector,
  cutoff,
  T::DataType,
  lcell::Int=1,
  UnitCellType=OrthorhombicCell
)
```

For orthorhombic unit cells, `Box` can be initialized with a vector of the
length of each side.

### Example
```julia
julia> box = Box([120,150,100],10)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [120.0 0.0 0.0; 0.0 150.0 0.0; 0.0 0.0 100.0]
  cutoff: 10.0
  number of computing cells on each dimension: [14, 17, 12]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 2856

```

"""
function Box(
  sides::AbstractVector,
  cutoff,
  T::DataType,
  lcell::Int=1,
  UnitCellType=OrthorhombicCell
)
  N = length(sides)
  cart_idxs = CartesianIndices((1:N,1:N))
  # Build unit cell matrix from lengths
  unit_cell_matrix = SMatrix{N,N,T,N*N}(
    ntuple(N*N) do i
      c = cart_idxs[i]
      if c[1] == c[2]
        return sides[c[1]]
      else
        return zero(T)
      end
    end
  )
  return Box(
    unit_cell_matrix,
    cutoff,
    T,
    lcell,
    UnitCellType
  )
end
Box(
  sides::AbstractVector,
  cutoff;
  T::DataType=Float64,
  lcell::Int=1,
  UnitCellType=OrthorhombicCell
) = Box(sides,cutoff,T,lcell,UnitCellType)

"""

$(TYPEDEF)

$(TYPEDFIELDS)

Copies particle coordinates and associated index, to build contiguous particle lists
in memory when building the cell lists. This strategy duplicates the particle coordinates
data, but is probably worth the effort.

"""
struct ParticleWithIndex{N,T}
  index::Int
  index_original::Int
  coordinates::SVector{N,T}
  real::Bool
end
Base.zero(::Type{ParticleWithIndex{N,T}}) where {N,T} =
  ParticleWithIndex{N,T}(0,0,zeros(SVector{N,T}),false)

"""

$(TYPEDEF)

$(TYPEDFIELDS)

This structure contains the cell linear index and the information
about if this cell is in the border of the box (such that its
neighbouring cells need to be wrapped)

"""
struct Cell{N,T}
  icell::Int
  cartesian::CartesianIndex{N}
  center::SVector{N,T}
end
Base.zero(::Type{Cell{N,T}}) where {N,T} =
  Cell{N,T}(0,CartesianIndex{N}(ntuple(i->0,N)),zeros(SVector{N,T}))

"""

$(TYPEDEF)

$(TYPEDFIELDS)

Auxiliary structure to contain projected particles in large and
and dense systems.

"""
struct ProjectedParticle{N,T}
  index_original::Int
  xproj::T
  coordinates::SVector{N,T}
  real::Bool
end

"""

$(TYPEDEF)

$(TYPEDFIELDS)

Structure that contains the cell lists information.

"""
Base.@kwdef struct CellList{N,T}
  " *mutable* number of cells with real particles. "
  ncwp::Vector{Int}
  " *mutable* number of particles in the computing box "
  ncp::Vector{Int}
  " Auxiliary array to annotate if the cell contains real particles. "
  contains_real::Vector{Bool}
  " Indices of the unique cells with real particles. "
  cwp::Vector{Cell{N,T}}
  " First particle of cell. "
  fp::Vector{ParticleWithIndex{N,T}}
  " Next particle of cell "
  np::Vector{ParticleWithIndex{N,T}}
  " Number of particles of cell. "
  npcell::Vector{Int}
  " Auxiliar array to store projected particles. "
  projected_particles::Vector{Vector{ProjectedParticle{N,T}}}
end
function Base.show(io::IO,::MIME"text/plain",cl::CellList)
  println(typeof(cl))
  println("  $(cl.ncwp[1]) cells with real particles.")
  println("  $(cl.ncp[1]) particles in computing box, including images.")
end

"""

$(TYPEDEF)

$(TYPEDFIELDS)

Structure that will cointain the cell lists of two independent sets of
particles for cross-computation of interactions

"""
@with_kw struct CellListPair{V,N,T}
  small::V
  large::CellList{N,T}
  swap::Bool
end
function Base.show(io::IO,::MIME"text/plain",cl::CellListPair)
  print(typeof(cl),"\n")
  print("   $(length(cl.small)) particles in the smallest vector.\n")
  print("   $(cl.large.ncwp[1]) cells with particles.")
end

"""

```
CellList(
  x::AbstractVector{AbstractVector},
  box::Box{UnitCellType,N,T};
  parallel::Bool=true
) where {UnitCellType,N,T}
```

Function that will initialize a `CellList` structure from scracth, given a vector
or particle coordinates (a vector of vectors, typically of static vectors)
and a `Box`, which contain the size ofthe system, cutoff, etc.

### Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> cl = CellListMap.CellList(x,box)
CellList{3, Float64}
  1800 cells with real particles.
  161084 particles in computing box, including images.

```

"""
function CellList(
  x::AbstractVector{<:AbstractVector},
  box::Box{UnitCellType,N,T};
  parallel::Bool=true
) where {UnitCellType,N,T}
  number_of_cells = prod(box.nc)
  # number_of_particles is a lower bound, will be resized when necessary to incorporate particle images
  number_of_particles = ceil(Int,1.2*length(x))
  ncwp = [0]
  ncp = [0]
  cwp = Vector{Cell{N,T}}(undef,number_of_cells)
  contains_real = Vector{Bool}(undef,number_of_cells)
  fp = Vector{ParticleWithIndex{N,T}}(undef,number_of_cells)
  np = Vector{ParticleWithIndex{N,T}}(undef,number_of_particles)
  npcell = Vector{Int}(undef,number_of_cells)
  projected_particles = [ Vector{ProjectedParticle{N,T}}(undef,0) for i in 1:nthreads() ]

  cl = CellList{N,T}(
    ncwp,
    ncp,
    contains_real,
    cwp,
    fp,
    np,
    npcell,
    projected_particles
  )
  return UpdateCellList!(x,box,cl,parallel=parallel)
end

"""

```
CellList(
  x::AbstractVector{<:AbstractVector},
  y::AbstractVector{<:AbstractVector},
  box::Box{UnitCellType,N,T};
  parallel::Bool=true
) where {UnitCellType,N,T}
```

Function that will initialize a `CellListPair` structure from scracth, given two vectors
of particle coordinates and a `Box`, which contain the size ofthe system, cutoff, etc. The cell lists
will be constructed for the largest vector, and a reference to the smallest vector is annotated.

### Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> y = [ 250*rand(SVector{3,Float64}) for i in 1:10000 ];

julia> cl = CellList(x,y,box)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   1000 particles in the smallest vector.
   1767 cells with particles.

```

"""
function CellList(
  x::AbstractVector{SVector{N,T}},
  y::AbstractVector{SVector{N,T}},
  box::Box;
  parallel::Bool=true
) where {N,T}
  if length(x) <= length(y)
    y_cl = CellList(y,box,parallel=parallel)
    cl_pair = CellListPair(small=x,large=y_cl,swap=false)
  else
    x_cl = CellList(x,box,parallel=parallel)
    cl_pair = CellListPair(small=y,large=x_cl,swap=true)
  end
  return cl_pair
end

"""

```
particles_per_cell(cl::CellList,box::Box)
```

Returns the average number of particles per computing cell.

"""
particles_per_cell(cl::CellList,box::Box) = cl.ncp[1] / prod(box.nc)

"""

```
UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,cl:CellList{N,T},
    parallel=true
) where {N,T}
```

Function that will update a previously allocated `CellList` structure, given new
updated particle positions.

## Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> cl = CellList(x,box);

julia> box = Box([260,260,260],10);

julia> x = [ 260*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> cl = UpdateCellList!(x,box,cl); # update lists

```

"""
function UpdateCellList!(
  x::AbstractVector{<:AbstractVector},
  box::Box,
  cl::CellList{N,T};
  parallel::Bool=true
) where {N,T}
  @unpack contains_real, cwp, fp, np, npcell, projected_particles = cl

  number_of_cells = prod(box.nc)
  if number_of_cells > length(cwp)
    number_of_cells = ceil(Int,1.2*number_of_cells) # some margin in case of box size variations
    resize!(contains_real,number_of_cells)
    resize!(cwp,number_of_cells)
    resize!(fp,number_of_cells)
    resize!(npcell,number_of_cells)
  end

  cl.ncwp[1] = 0
  fill!(contains_real,false)
  fill!(cwp,zero(Cell{N,T}))
  fill!(fp,zero(ParticleWithIndex{N,T}))
  fill!(np,zero(ParticleWithIndex{N,T}))
  fill!(npcell,0)

  #
  # The following part cannot be *easily* paralelized, because
  # there is concurrency on the construction of the cell lists
  #

  #
  # Add virtual particles to edge cells
  #
  for (ip,particle) in pairs(x)
    p = SVector{N,T}(wrap_to_first(particle,box))
    cl = replicate_particle!(ip,p,box,cl)
  end
  #
  # Add true particles, such that the first particle of each cell is
  # always a true particle
  #
  for (ip,particle) in pairs(x)
    p = SVector{N,T}(wrap_to_first(particle,box))
    cl = add_particle_to_celllist!(ip,p,box,cl)
  end

  maximum_npcell = maximum(npcell)
  if maximum_npcell > length(projected_particles[1])
    for i in 1:nthreads()
      resize!(projected_particles[i],ceil(Int,1.2*maximum_npcell))
    end
  end

  return cl
end

"""

```
add_particle_to_celllist!(
  ip,
  x::SVector{N,T},
  box,
  cl::CellList{N,T};
  real_particle::Bool=true
) where {N,T}
```

Adds one particle to the cell lists, updating all necessary arrays.

"""
function add_particle_to_celllist!(
  ip,
  x::SVector{N,T},
  box,
  cl::CellList{N,T};
  real_particle::Bool=true
) where {N,T}
  @unpack contains_real, ncp, ncwp, cwp, fp, np, npcell = cl
  ncp[1] += 1
  icell_cartesian = particle_cell(x,box)
  icell = cell_linear_index(box.nc,icell_cartesian)
  #
  # Cells starting with real particles are annotated to be run over
  #
  if real_particle && (!contains_real[icell])
    contains_real[icell] = true
    ncwp[1] += 1
    cwp[ncwp[1]] = Cell{N,T}(
      icell,
      icell_cartesian,
      cell_center(icell_cartesian,box)
    )
  end
  if fp[icell].index == 0
    npcell[icell] = 1
  else
    npcell[icell] += 1
  end
  if ncp[1] > length(np)
    old_length = length(np)
    resize!(np,ceil(Int,1.2*old_length))
    for i in old_length+1:length(np)
      np[i] = zero(ParticleWithIndex{N,T})
    end
  end
  np[ncp[1]] = fp[icell]
  fp[icell] = ParticleWithIndex(ncp[1],ip,x,real_particle)
  return cl
end

"""

```
UpdateCellList!(
  x::AbstractVector{<:AbstractVector},
  y::AbstractVector{<:AbstractVector},
  box::Box{UnitCellType,N,T},cl:CellListPair,parallel=true
) where {UnitCellType,N,T}
```

Function that will update a previously allocated `CellListPair` structure, given new updated particle positions, for example.

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> y = [ 250*rand(SVector{3,Float64}) for i in 1:10000 ];

julia> cl = CellList(x,y,box);

julia> cl = UpdateCellList!(x,y,box,cl); # update lists

```

"""
function UpdateCellList!(
  x::AbstractVector{<:AbstractVector},
  y::AbstractVector{<:AbstractVector},
  box::Box{UnitCellType,N,T},cl_pair::CellListPair;
  parallel::Bool=true
) where {UnitCellType,N,T}

  if length(x) <= length(y)
    UpdateCellList!(y,box,cl_pair.large,parallel=parallel)
  else
    UpdateCellList!(x,box,cl_pair.large,parallel=parallel)
  end

  return cl_pair
end


text:
# Mon 31 Dec 2012 01:40:11 PM EST
#
# NeHe Tut 17 - 2D texture font
#
# Q - quit
# L - turn lights on/off
# F - change texture filter (linear, nearest, mipmap)
# PageUp/Down - move camera closer/further away
# Up/Down - increase/decrease x-rotation speed
# Left/Right - increase/decrease y-rotation speed

# TODO: This example runs, but it produces a very glitchy output.

# load necessary GLUT/OpenGL routines

global OpenGLver="1.0"
using OpenGL
using GLUT

### auxiliary functions

function glPrint(x::Integer, y::Integer, string::String, set::Integer)
    global base

    if set > 1
        set = 1
    end

    glBindTexture(GL_TEXTURE_2D, tex[1])
    glDisable(GL_DEPTH_TEST)

    glMatrixMode(GL_PROJECTION)
    glpushmatrix()

    glLoadIdentity()
    glOrtho(0, 640, 0, 480, -1, 1)

    glMatrixMode(GL_MODELVIEW)
    glPushMatrix()
    glLoadIdentity()

    glTranslate(x, y, 0)
    glListBase(uint32(base-32+(128*set)))
    glCallLists(strlen(string), GL_BYTE, string)

    glMatrixMode(GL_PROJECTION)
    glPopMatrix()

    glMatrixMode(GL_MODELVIEW)
    glPopMatrix()
end

### end of auxiliary functions

# initialize variables

global window
global base

global tex       = Array(Uint32,2) # generating 2 textures

T0                 = 0
Frames             = 0

cnt1               = 0
cnt2               = 0

width            = 640
height           = 480

# load textures from images

function LoadGLTextures()
    global tex

    imgFont, wFont, hFont = glimread(expanduser("~/.julia/SDL/Examples/NeHe/tut17/font.bmp"))

    imgBumps, wBumps, hBumps = glimread(expanduser("~/.julia/SDL/Examples/NeHe/tut17/bumps.bmp"))

    glGenTextures(2,tex)
    glBindTexture(GL_TEXTURE_2D,tex[1])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, 3, wFont, hFont, 0, GL_RGB, GL_UNSIGNED_BYTE, imgFont)

    glBindTexture(GL_TEXTURE_2D,tex[2])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, 3, wBumps, hBumps, 0, GL_RGB, GL_UNSIGNED_BYTE, imgBumps)
end

# function to init OpenGL context

function initGL(w::Integer,h::Integer)
    global base

    glViewport(0,0,w,h)
    LoadGLTextures()
    glClearColor(0.0, 0.0, 0.0, 0.0)
    glClearDepth(1.0)
    glDepthFunc(GL_LEQUAL)
    glEnable(GL_DEPTH_TEST)
    glShadeModel(GL_SMOOTH)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    # enable texture mapping & blending
    glEnable(GL_TEXTURE_2D)
    glBlendFunc(GL_SRC_ALPHA, GL_ONE)

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    # build the fonts

    base = glGenLists(256)
    glBindTexture(GL_TEXTURE_2D, tex[1])

    for loop = 1:256
        cx = (loop%16)/16
        cy = (loop/16)/16

        glNewList(uint32(base+(loop-1)), GL_COMPILE)
            glBegin(GL_QUADS)
                glTexCoord(cx, 1-cy-0.0625)
                glVertex(0, 0)
                glTexCoord(cx+0.0625, 1-cy-0.0625)
                glVertex(16, 0)
                glTexCoord(cx+0.0625, 1-cy)
                glVertex(16, 16)
                glTexCoord(cx, 1-cy)
                glVertex(0, 16)
            glEnd()
            glTranslate(10, 0, 0)
        glEndList()
    end
end

# prepare Julia equivalents of C callbacks that are typically used in GLUT code

function ReSizeGLScene(w::Int32,h::Int32)
    if h == 0
        h = 1
    end

    glViewPort(0,0,w,h)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    return nothing
end

_ReSizeGLScene = cfunction(ReSizeGLScene, Void, (Int32, Int32))

function DrawGLScene()
    global tex
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    glBindTexture(GL_TEXTURE_2D,tex[2])

    glTranslate(0.0,0.0,-5.0)
    glRotate(45.0, 0.0, 0.0, 1.0)
    glRotate(cnt1*30.0, 1.0, 1.0, 0.0)

    glDisable(GL_BLEND)
    glColor(1.0, 1.0, 1.0)

    glBegin(GL_QUADS)
        glTexCoord(0.0, 0.0)
        glVertex(-1.0, 1.0)
        glTexCoord(1.0, 0.0)
        glVertex(1.0, 1.0)
        glTexCoord(1.0, 1.0)
        glVertex(1.0, -1.0)
        glTexCoord(0.0, 1.0)
        glVertex(-1.0, -1.0)
    glEnd()

    glRotate(90.0, 1.0, 1.0, 0.0)
    glBegin(GL_QUADS)
        glTexCoord(0.0, 0.0)
        glVertex(-1.0, 1.0)
        glTexCoord(1.0, 0.0)
        glVertex(1.0, 1.0)
        glTexCoord(1.0, 1.0)
        glVertex(1.0, -1.0)
        glTexCoord(0.0, 1.0)
        glVertex(-1.0, -1.0)
    glEnd()

    glEnable(GL_BLEND)
    glLoadIdentity()

    glColor(1.0cos(cnt1), 1.0sin(cnt2), 1.0-0.5cos(cnt1+cnt2))
    glPrint(int(280+250cos(cnt1)), int(235+200sin(cnt2)), "NeHe", 0)
    glColor(1.0sin(cnt2), 1.0-0.5cos(cnt1+cnt2), 1.0cos(cnt1))
    glPrint(int(280+230cos(cnt2)), int(235+200sin(cnt1)), "OpenGL", 1)

    glColor(0.0, 0.0, 1.0)
    glPrint(int(240+200cos((cnt2+cnt1)/5)), 2, "JuliaLang", 0)
    glColor(1.0, 1.0, 1.0)
    glPrint(int(242+200cos((cnt2+cnt1)/5)), 2, "JuliaLang", 0)

    cnt1 +=0.01
    cnt2 +=0.0081

    glutSwapBuffers()

    return nothing
end

_DrawGLScene = cfunction(DrawGLScene, Void, ())

function keyPressed(key::Char,x::Int32,y::Int32)
    global base

    if key == int('q')
        glDeleteLists(base,256)
        glDeleteTextures(2,tex)
        glutDestroyWindow(window)
    end

    return nothing # keyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_keyPressed = cfunction(keyPressed, Void, (Char, Int32, Int32))

# run GLUT routines

glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH)
glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)

window = glutCreateWindow("NeHe Tut 17")

glutDisplayFunc(_DrawGLScene)
glutFullScreen()

glutIdleFunc(_DrawGLScene)
glutReshapeFunc(_ReSizeGLScene)
glutKeyboardFunc(_keyPressed)

initGL(width, height)

glutMainLoop()


text:
const sleeptime = 0.5

println("... testing multithreading  1, (sleeptime=$sleeptime) ...")
println("number of available threads: ", nthreads())

clk = PClock()
Δt = clk.Δt
sleep(sleeptime)

print(clk)
@test clk.id == 0
m = match(r"Clock 0, thrd 1 \(\+ (\d+) ac\)", repr(clk))
@test parse(Int, m.captures[1]) > 0
@test length(clk.ac) == nthreads()-1
@test clk.ac[1].thread == 2

println("... parallel clock identification ...")
c1 = pclock(clk, 1)
@test c1.thread == 2
@test c1.clock.id == 1
m = match(r"Active clock 1 on thrd (\d+)\:", repr(c1))
@test parse(Int, m.captures[1]) > 1

if DiscreteEvents._handle_exceptions[1]
    println("... remote error handling ...")
    put!(clk.ac[1].forth, DiscreteEvents.Clear())
    sleep(sleeptime)
    err = diagnose(clk, 1)
    @test err[1] isa ErrorException
    @test occursin(r"^step\!\(\:\:DiscreteEvents.ActiveClock.+ at threads\.jl", string(err[2][2]))
end

println("... testing channel and active clock ...")
a = 0
ev = DiscreteEvents.DiscreteEvent(fun(()->global a+=1),1.0,0.0)
put!(clk.ac[1].forth, DiscreteEvents.Register(ev))
sleep(sleeptime)
@test length(c1.clock.sc.events) == 1
cev = DiscreteEvents.DiscreteCond(fun(≥, fun(tau, c1), 5), fun(()->global a+=1))
put!(clk.ac[1].forth, DiscreteEvents.Register(cev))
sleep(sleeptime)
@test length(c1.clock.sc.cevents) == 1
b = 0
sp = DiscreteEvents.Sample(fun(()->global b+=1))
put!(clk.ac[1].forth, DiscreteEvents.Register(sp))
sleep(sleeptime)
@test length(c1.clock.sc.samples) == 1

println("... run parallel clock 1 (thread 2) ...")
# put!(clk.ac[1].forth, DiscreteEvents.Run(10.0))
# @test take!(clk.ac[1].back).t > 0
# sleep(sleeptime)
tns = 0
t1 = time_ns()
iter = 0
while (c1.clock.time+Δt) ≤ 10
    put!(clk.ac[1].forth, DiscreteEvents.Run(Δt, false))
    global tns += take!(clk.ac[1].back).t
    global iter += 1
end
t1 = time_ns() - t1
println("$iter ticks took $(Int(t1)*1e-9) s, clock time $(Int(tns)*1e-9) s")
@test c1.clock.time ≈ 10
@test c1.clock.scount == 1000
@test a == 2
@test b == 1000

println("... parallel clock 1 reset ...")
put!(clk.ac[1].forth, DiscreteEvents.Reset(true))
@test take!(clk.ac[1].back).x == 1
@test c1.clock.time == 0
@test c1.clock.scount == 0

println("... testing register(!) five cases ... ")
ev1 = DiscreteEvents.DiscreteEvent(fun(()->global a+=1),1.0,0.0)
ev2 = DiscreteEvents.DiscreteEvent(fun(()->global a+=1),2.0,0.0)
DiscreteEvents._register(clk, ev1, 0)              # 1. register ev1 to clk
@test DiscreteEvents._nextevent(clk) == ev1
DiscreteEvents._register(clk, ev1, 1)              # 2. register ev1 to 1st parallel clock
sleep(sleeptime)
@test DiscreteEvents._nextevent(c1.clock) == ev1
DiscreteEvents._register(c1, ev2, 1)               # 3. register ev2 directly to 1st parallel clock
sleep(sleeptime)
@test length(c1.clock.sc.events) == 2
DiscreteEvents._register(c1, ev2, 0)               # 4. register ev2 back to master
sleep(sleeptime)
@test length(clk.sc.events) == 2

if nthreads() > 2                           # This fails on CI (only 2 threads)
    DiscreteEvents._register(c1, ev2, 2)           # 5. register ev2 to another parallel clock
    c2 = pclock(clk, 2)
    sleep(sleeptime)
    @test DiscreteEvents._nextevent(c2.clock) == ev2
end

println("... testing API on multiple threads ...")

function clocks_ok(clk::Clock)
    ok = true
    for i in eachindex(clk.ac)
        while isready(clk.ac[i].back)
            token = take!(clk.ac[i].back)
            if token isa Error
                println("clock $i error: ", token.exc)
                ok = false
            end
        end
        if istaskfailed(clk.ac[i].ref[])
            println("clock $i error, stacktrace:")
            println(clk.ac[i].ref[])
            ok = false
        end
    end
    return ok
end

clock = Vector{Clock}()
push!(clock, clk)
for i in eachindex(clk.ac)
    push!(clock, pclock(clk, i).clock)
end
println("    got $(length(clock)) clocks")
@test sum(length(c.sc.events) for c in clock) ≥ 4
resetClock!(clk)                                 # test resetClock! on multiple clocks
@test clk.Δt == Δt                          # must still be the same
@test sum(length(c.sc.events) for c in clock) == 0
println("    resetClock!     ok")

a = zeros(Int, nthreads())
event!(clk, ()->a[1]+=1, 1)
sleep(sleeptime); @test clocks_ok(clk)
event!(clk, ()->a[2]+=1, 1, cid=1)
sleep(sleeptime); @test clocks_ok(clk)
event!(c1,  ()->a[2]+=1, 2)
sleep(sleeptime); @test clocks_ok(clk)
event!(c1,  ()->a[1]+=1, 2, cid=0)
sleep(sleeptime); @test clocks_ok(clk)
if nthreads() > 2
    event!(c1,  ()->a[3]+=1, 1, cid=2)
    sleep(sleeptime); @test clocks_ok(clk)
end
sleep(sleeptime)
@test length(clk.sc.events) == 2
@test length(c1.clock.sc.events) == 2
nthreads() > 2 && @test length(c2.clock.sc.events) == 1
error = false
@test !error
println("    distributed events ok")
run!(clk, 3)
@test a[1] == 2
@test a[2] == 2
nthreads() > 2 && @test a[3] == 1
println("    first parallel run ok")

println("... collapse! ...")
ac = clk.ac
collapse!(clk)
@test all(x->istaskdone(x.ref[]), ac)
@test all(x->x.ch.state == :closed, ac)


text:
## ----------------------------------------------------------------------------
# BUNDLE from the glucose-limited culture reported at
# Folsom, James Patrick, Albert E. Parker, and Ross P. Carlson.
# “Physiological and Proteomic Analysis of Escherichia Coli Iron-Limited
# Chemostat Growth.” Journal of Bacteriology 196, no. 15 (August 1, 2014):
# 2748–61. https://doi.org/10.1128/JB.01606-14.

## ----------------------------------------------------------------------------
# Sumpplementary table 2
TABLE_S2 = Dict()
function _populate_table_s2()
    empty!(TABLE_S2)
    TABLE_S2["D"] = Dict(
        "unit" => "1/h",
        "val" => [0.1, 0.2, 0.3, 0.4],
        "err" => [0.0, 0.0, 0.0, 0.0],
    )
    TABLE_S2["μ"] = deepcopy(TABLE_S2["D"])
    TABLE_S2["YX_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.426, 0.525, 0.458, 0.459],
        "err" => [0.032, 0.028, 0.034, 0.038],
    )
    TABLE_S2["Ypyr_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.0, 0.0, 0.0, 0.0],
        "err" => [0.0, 0.0, 0.0, 0.0],
    )
    TABLE_S2["Ysucc_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.004, 0.003, 0.0, 0.0],
        "err" => [0.001, 0.001, 0.0, 0.0],
    )
    TABLE_S2["Ylac_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.013, 0.003, 0.026, 0.023],
        "err" => [0.003, 0.001, 0.002, 0.003],
    )
    TABLE_S2["Yform_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.041, 0.005, 0.001, 0.002],
        "err" => [0.015, 0.001, 0.000, 0.002],
    )
    TABLE_S2["Yac_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.000, 0.000, 0.001, 0.004],
        "err" => [0.000, 0.000, 0.001, 0.002],
    )
    TABLE_S2["Yco2_glc"] = Dict(
        "unit" => "Cmmol/Cmmol glc",
        "val" => [0.515, 0.465, 0.514, 0.512],
        "err" => [0.051, 0.031, 0.037, 0.047],
    )

    TABLE_S2["uglc"] = Dict(
        "unit" => "Cmmol/gCDW-h",
        "val" => [9.1, 14.8, 25.5, 33.9],
        "err" => [0.68, 0.78, 1.87, 2.79],
    )
    TABLE_S2["uO2"] = Dict(
        "unit" => "Cmmol/gCDW-h",
        "val" => [4.74, 6.64, 12.68, 16.82],
        "err" => [0.826, 0.786, 1.85, 2.93],
    )
    TABLE_S2
end

## ----------------------------------------------------------------------------
#  M9 minima medium, all in mM
MEDIUM = Dict()
function _populate_medium()
    empty!(MEDIUM)
    MEDIUM["glc"] = 2.2 # mM
    MEDIUM
end

## ----------------------------------------------------------------------------
# Bundle, converting
BUNDLE = Dict()
function _populate_bundle()

    _populate_table_s2()
    _populate_medium()
    empty!(BUNDLE)

    # flxs
    BUNDLE["D"] = deepcopy(TABLE_S2["D"])
    BUNDLE["μ"] = deepcopy(TABLE_S2["μ"])
    BUNDLE["uGLC"] = Dict(
        # q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "glucose uptake rate",
        "unit" => "mmol/gCDW h",
        "val" => -TABLE_S2["uglc"]["val"] ./ 6,
        "err" => TABLE_S2["uglc"]["err"] ./ 6,
    )
    BUNDLE["uO2"] = Dict(
        # q[mmol/gCDW-h]
        "name" => "oxigen uptake rate",
        "unit" => "mmol/gCDW h",
        "val" => -TABLE_S2["uO2"]["val"],
        "err" => TABLE_S2["uO2"]["err"]
    )

    # yields
    # qglc[Cmmol/gCDW-h]
    qglc = TABLE_S2["uglc"]["val"]

    BUNDLE["uPYR"] = Dict(
        # Y[Cmol/Cmol] * qglc[Cmmol/gCDW-h] = q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "pyruvate production rate",
        "unit" => "mmol/gCDW h",
        "val" => TABLE_S2["Ypyr_glc"]["val"] .* qglc ./ 3,
        "err" => TABLE_S2["Ypyr_glc"]["err"] .* qglc ./ 3,
    )
    BUNDLE["uSUCC"] = Dict(
        # Y[Cmol/Cmol] * qglc[Cmmol/gCDW-h] = q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "succinate production rate",
        "unit" => "mmol/gCDW h",
        "val" => TABLE_S2["Ysucc_glc"]["val"] .* qglc ./ 4,
        "err" => TABLE_S2["Ysucc_glc"]["err"] .* qglc ./ 4,
    )
    BUNDLE["uLAC"] = Dict(
        # Y[Cmol/Cmol] * qglc[Cmmol/gCDW-h] = q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "lactate production rate",
        "unit" => "mmol/gCDW h",
        "val" => TABLE_S2["Ylac_glc"]["val"] .* qglc ./ 3,
        "err" => TABLE_S2["Ylac_glc"]["err"] .* qglc ./ 3,
    )
    BUNDLE["uFORM"] = Dict(
        # Y[Cmol/Cmol] * qglc[Cmmol/gCDW-h] = q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "formate production rate",
        "unit" => "mmol/gCDW h",
        "val" => TABLE_S2["Yform_glc"]["val"] .* qglc ./ 1,
        "err" => TABLE_S2["Yform_glc"]["err"] .* qglc ./ 1,
    )
    BUNDLE["uAC"] = Dict(
        # Y[Cmol/Cmol] * qglc[Cmmol/gCDW-h] = q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "acetate production rate",
        "unit" => "mmol/gCDW h",
        "val" => TABLE_S2["Yac_glc"]["val"] .* qglc ./ 2,
        "err" => TABLE_S2["Yac_glc"]["err"] .* qglc ./ 2,
    )
    BUNDLE["uCO2"] = Dict(
        # Y[Cmol/Cmol] * qglc[Cmmol/gCDW-h] = q[Cmmol/gCDW-h] / #C = q[mmol/gCDW-h]
        "name" => "acetate production rate",
        "unit" => "mmol/gCDW h",
        "val" => TABLE_S2["Yco2_glc"]["val"] .* qglc ./ 1,
        "err" => TABLE_S2["Yco2_glc"]["err"] .* qglc ./ 1,
    )

    # There negigeble glucose present in the medium at steady state
    # so we can compute X = c * D/qglc
    BUNDLE["X"] = Dict(
        # c[mmol/L] * D[1/h] / q[mmol/gCDW h] = X[gCDW/L]
        "name" => "cell concentration",
        "unit" => "gCDW/ L",
        "val" => abs.(MEDIUM["glc"] .* BUNDLE["D"]["val"] ./ BUNDLE["uGLC"]["val"]),
        "err" => zero(BUNDLE["uGLC"]["err"]) #TODO: fix this
    )

    # xi
    BUNDLE["xi"] = Dict(
        "name" => "dilution specific cell concentration",
        "unit" => "gCDW/ L h",
        "val" => BUNDLE["X"]["val"] ./ BUNDLE["D"]["val"],
        "err" => BUNDLE["X"]["err"] #TODO: fix this
    )

    # medium
    BUNDLE["cGLC"] = Dict(
        "name" => "glucose feed concentration",
        "unit" => "gCDW/ L h",
        "val" => fill(MEDIUM["glc"], 4),
        "err" => zeros(4)
    )

    BUNDLE
end

## ------------------------------------------------------------------
# API
const EXPS = 1:4
const msd_mets = ["GLC", "PYR", "SUCC", "LAC", "FORM", "AC"]

_get_val(id, dk) = BUNDLE[string(id)][dk]
_get_val(id, dk, exp::Int) = BUNDLE[string(id)][dk][exp]
function _get_val(id, dk, D::Float64)
    exp = findfirst(BUNDLE["D"][dk] .== D)
    isnothing(exp) && error("No experiment with D = $D")
    _get_val(id, dk, exp)
end
_get_val(id, dk, ref, dflt) =
    try; _get_val(id, dk, ref); catch err; dflt end

for fun in [:val, :err]
    dk = string(fun)
    @eval $fun(id) = _get_val(id, $dk)
    @eval $fun(id, ref) = _get_val(id, $dk, ref)
    @eval $fun(id, ref, dflt) = _get_val(id, $dk, ref, dflt)

    for p in [:u, :c]
        pstr = string(p)
        pfun = Symbol(p, fun)
        @eval $pfun(id, args...) = $fun(string($pstr, id), args...)
    end
end

name(id) = BUNDLE[string(id)]["name"]
unit(id) = BUNDLE[string(id)]["unit"]

ciD_X(id) = [cval(id, exp, 0.0) * val(:D, exp) / val(:X, exp) for exp in EXPS]
ciD_X(id, exp) = cval(id, exp, 0.0) * val(:D, exp) / val(:X, exp)


text:
@testset "mutate_node!" begin
    mtg = read_mtg("files/simple_plant.mtg", Dict)
    # Using a leaf node from the mtg:
    leaf_node = mtg.children[2].children[3].children[4].children[5]
    # Add a new attributes, x based on node field, y based on x and z using a function:
    @mutate_node!(mtg, x = length(node.name), y = node.x + 2, z = sum(node.y))
    @test mtg[:x] == 6
    @test mtg[:y] == mtg[:x] + 2
    @test mtg[:z] == sum(mtg[:y])

    # Same on a leaf:
    @mutate_node!(leaf_node, x = length(node.name), y = node.x + 3, z = sum(node.y))
    @test leaf_node[:x] == 6
    @test leaf_node[:y] == leaf_node[:x] + 3
    @test leaf_node[:z] == sum(leaf_node[:y])

    # Test by using node[:variable] format:
    @mutate_node!(leaf_node, node[:test] = node[:y])
    @test leaf_node[:test] == leaf_node[:y]

    # Test by using node fields instead of attributes:
    @mutate_node!(leaf_node, node.MTG.scale = 4)
    @test leaf_node.MTG.scale == 4

end

@testset "mutate_mtg!" begin
    mtg = read_mtg("files/simple_plant.mtg", Dict)
    # Using a leaf node from the mtg:
    leaf_node = mtg.children[2].children[3].children[4].children[5]

    # Add a new attribute based on the name:
    @mutate_mtg!(mtg, nametest = node.name * "_test")
    @test mtg[:nametest] == "node_1_test"
    @test leaf_node[:nametest] == "node_5_test"

    # Mutate with a filter:
    @mutate_mtg!(mtg, nametest = node.name * "_test2", scale = 3)
    @test mtg[:nametest] == "node_1_test"
    @test leaf_node[:nametest] == "node_5_test2"
end

@testset "append!" begin

    # On Dict attributes
    mtg = read_mtg("files/simple_plant.mtg", Dict)
    append!(mtg, (test = 3,))
    haskey(mtg.attributes, :test)

    @test haskey(mtg.attributes, :test)
    @test mtg[:test] == 3

    # On MutableNamedTuple attributes
    mtg = read_mtg("files/simple_plant.mtg", MutableNamedTuple)
    # Using a leaf node from the mtg:
    append!(mtg, (test = 3,))
    @test :test in keys(mtg.attributes)
    @test mtg[:test] == 3

    # On NamedTuple attributes
    mtg = read_mtg("files/simple_plant.mtg", NamedTuple)
    # Using a leaf node from the mtg:
    append!(mtg, (test = 3,))
    @test :test in keys(mtg.attributes)
    @test mtg[:test] == 3
end

@testset "pop!" begin
    # On Dict attributes
    mtg = read_mtg("files/simple_plant.mtg", Dict)
    pop!(mtg, :symbols)
    @test !haskey(mtg.attributes, :symbols)

    # On MutableNamedTuple attributes
    mtg = read_mtg("files/simple_plant.mtg", MutableNamedTuple)
    # Using a leaf node from the mtg:
    pop!(mtg, :symbols)
    @test !(:symbols in keys(mtg.attributes))

    # On NamedTuple attributes
    mtg = read_mtg("files/simple_plant.mtg", NamedTuple)
    # Using a leaf node from the mtg:
    pop!(mtg, :symbols)
    @test !(:symbols in keys(mtg.attributes))
end


text:
"""
    Fermi.PhysicalConstants

Module for storing physical constants and convertion factors used in computations.

# Functions:

    atomic_number   Given an element symbol, return the atomic number.
"""
module PhysicalConstants

export atomic_number

"""
    Fermi.PhysicalConstants.bohr_to_angstrom

Float64 object with the convertion factor from Bohr to Angstrom.
Source: [NIST](https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0)
"""
#bohr_to_angstrom = 0.529177210903
#psi4s
bohr_to_angstrom = 0.529177210670000

"""
    Fermi.PhysicalConstants.angstrom_to_bohr

Float64 object with the convertion factor from Angstrom to Bohr.
Computed as 1/Fermi.PhysicalConstants.bohr_to_angstrom
"""
angstrom_to_bohr = 1/bohr_to_angstrom

hartree_to_kcalmol = 627.509474
hartree_to_kjmol = 2625.49964
hartree_to_cm = 219474.6313708
hartree_to_eV = 27.21138505

"""
    Fermi.PhysicalConstants.atomic_number(atom::String)

Given a string with the element symbol, return its atomic number
"""
function atomic_number(atom::String)

    try
        return atom_num[atom]
    catch
        error("Invalid atom: $atom")
    end

end

atom_num = Dict(
    "H" => 1,
    "He" => 2,
    "Li" => 3,
    "Be" => 4,
    "B" => 5,
    "C" => 6,
    "N" => 7,
    "O" => 8,
    "F" => 9,
    "Ne" => 10,
    "Na" => 11,
    "Mg" => 12,
    "Al" => 13,
    "Si" => 14,
    "P" => 15,
    "S" => 16,
    "Cl" => 17,
    "Ar" => 18,
    "K" => 19,
    "Ca" => 20,
    "Sc" => 21,
    "Ti" => 22,
    "V" => 23,
    "Cr" => 24,
    "Mn" => 25,
    "Fe" => 26,
    "Co" => 27,
    "Ni" => 28,
    "Cu" => 29,
    "Zn" => 30,
    "Ga" => 31,
    "Ge" => 32,
    "As" => 33,
    "Se" => 34,
    "Br" => 35,
    "Kr" => 36,
    "Rb" => 37,
    "Sr" => 38,
    "Y" => 39,
    "Zr" => 40,
    "Nb" => 41,
    "Mo" => 42,
    "Tc" => 43,
    "Ru" => 44,
    "Rh" => 45,
    "Pd" => 46,
    "Ag" => 47,
    "Cd" => 48,
    "In" => 49,
    "Sn" => 50,
    "Sb" => 51,
    "Te" => 52,
    "I" => 53,
    "Xe" => 54,
    "Cs" => 55,
    "Ba" => 56,
    "La" => 57,
    "Ce" => 58,
    "Pr" => 59,
    "Nd" => 60,
    "Pm" => 61,
    "Sm" => 62,
    "Eu" => 63,
    "Gd" => 64,
    "Tb" => 65,
    "Dy" => 66,
    "Ho" => 67,
    "Er" => 68,
    "Tm" => 69,
    "Yb" => 70,
    "Lu" => 71,
    "Hf" => 72,
    "Ta" => 73,
    "W" => 74,
    "Re" => 75,
    "Os" => 76,
    "Ir" => 77,
    "Pt" => 78,
    "Au" => 79,
    "Hg" => 80,
    "Tl" => 81,
    "Pb" => 82,
    "Bi" => 83,
    "Po" => 84,
    "At" => 85,
    "Rn" => 86,
    "Fr" => 87,
    "Ra" => 88,
    "Ac" => 89,
    "Th" => 90,
    "Pa" => 91,
    "U" => 92,
    "Np" => 93,
    "Pu" => 94,
    "Am" => 95,
    "Cm" => 96,
    "Bk" => 97,
    "Cf" => 98,
    "Es" => 99,
    "Fm" => 100,
    "Md" => 101,
    "No" => 102,
    "Lr" => 103,
    "Rf" => 104,
    "Db" => 105,
    "Sg" => 106,
    "Bh" => 107,
    "Hs" => 108,
    "Mt" => 109,
    "Ds" => 110,
    "Rg" => 111,
    "Cn" => 112,
    "Nh" => 113,
    "Fl" => 114,
    "Mc" => 115,
    "Lv" => 116,
    "Ts" => 117,
    "Og" => 118
)

end #module


text:
using QuDiffEq
using OrdinaryDiffEq, Test
"""
Linearizing non linear ODE and solving it using QuLDE

du/dt = f
u0 -> inital condition for the ode

Jacobian J and vector b is calculated at every time step h.
Δu -> difference for fixed point

Equation input to QuLDE circuit : d(Δu)/dt = J * Δu + b
k -> order of Taylor Expansion in QuLDE circuit
Δu is added to previous value of u.
"""
function f(du,u,p,t)
    du[1] = -2*(u[2]^2)*u[1]
    du[2] = 3*(u[1]^(1.5)) - 0.1*u[2]
end

u0 = [0.2,0.1]
h = 0.1
k = 2
tspan = (0.0,0.8)
prob = ODEProblem(f,u0,tspan)

qsol = solve(prob,QuLDE(k),dt = h)
sol = solve(prob,Tsit5(),dt = h,adaptive = false)
r_out = transpose(hcat(sol.u...))

@test isapprox.(r_out,qsol, atol = 1e-3) |> all


text:
################################################################################
#
#    Module PotentialLearning.jl
#
################################################################################

module PotentialLearning

using GalacticOptim, Optim
using BlackBoxOptim
using Printf

export get_conf_params, generate_data, get_dft_data, learning_problem, learn,
       error_metrics, potential_energy, forces

abstract type PotentialLearningProblem end

include("./Potentials/Potentials.jl")
include("InputLoading.jl")
include("LearningData.jl")

"""
    learning_problem( dft_train_data::Vector{Float64},
                      ref_train_data::Vector{Float64},
                      params::Dict)

Creates a potential learning problem.
"""
function learning_problem(dft_train_data::Vector{Float64},
                          ref_train_data::Vector{Float64},
                          params::Dict)
    lp = Symbol(params["global"]["potential"])
    return @eval $lp($dft_train_data, $ref_train_data, $params)
end

"""
    learn(p::PotentialLearningProblem, params::Dict)

Fits the potentials, forces, and stresses against the DFT and reference data
using the configuration parameters.
"""
function learn(p::PotentialLearningProblem, params::Dict)
    if params["solver"]["name"] == "\\"
        p.β = p.A \ p.b
    elseif params["solver"]["name"] == "NelderMead"
        β0 = zeros(length(p.A[1,:]))
        prob = GalacticOptim.OptimizationProblem( (x, pars) -> error(x, p), β0, [])
        p.β = GalacticOptim.solve(prob, NelderMead(), maxiters=500)
    elseif params["solver"]["name"] == "BBO"
        β0 = zeros(p.cols)
        lb0 = -0.5ones(p.cols)
        ub0 = 0.5ones(p.cols)
        prob = GalacticOptim.OptimizationProblem( (x, pars) -> error(x, p), β0, [],
                                                  ub = ub0, lb = lb0)
        p.β = solve(prob, BBO(); maxiters=500)
    end
end

"""
    error_metrics(p::PotentialLearningProblem, val_data::Vector{Float64}, params::Dict)

Calculates error metrics: max_rel_error, mae, rmse
"""
function error_metrics(p::PotentialLearningProblem, val_data::Vector{Float64}, params::Dict)
    fit_forces = params["global"]["fit_forces"]
    rcutfac = params["global"]["rcutfac"]
    no_atoms_per_conf = params["global"]["no_atoms_per_conf"]
    no_train_atomic_conf = params["global"]["no_train_atomic_conf"]
    no_atomic_conf = params["global"]["no_atomic_conf"]
    no_val_atomic_conf = no_atomic_conf - no_train_atomic_conf

    global metrics = Dict()
    metrics["energy"] = Dict()

    # Energy metrics
    energies = val_data[1:no_val_atomic_conf]
    fitted_energies = [potential_energy(p, j) for j in no_train_atomic_conf+1:no_atomic_conf]
    metrics["energy"]["max_rel_error"] =
                    maximum(abs.(fitted_energies .- energies) ./ fitted_energies)
    metrics["energy"]["mae"] =
                    sum(abs.(fitted_energies .- energies)) / length(energies)
    metrics["energy"]["rmse"] =
                    sqrt(sum((fitted_energies .- energies).^2) / length(energies))

    # Force metrics
    if fit_forces
        forces_ = val_data[no_val_atomic_conf+1:end]
        fitted_forces = linearize([forces(p, j)
                                   for j in no_train_atomic_conf+1:no_atomic_conf])
        metrics["force"] = Dict()
        metrics["force"]["max_rel_error"] =
                    maximum(abs.((fitted_forces .- forces_) ./ forces_))
        metrics["force"]["mae"] =
                    sum(abs.(fitted_forces .- forces_)) / length(forces_)
        metrics["force"]["rmse"] =
                    sqrt(sum((fitted_forces .- forces_).^2) / length(forces_))
    end

    # TODO: Stress metrics

   return metrics
end

end


text:
load("Distributions")
using Distributions

load("MCMC")
using MCMC

function f(x)
  log(0.75 * pdf(Normal(0.0, 1.0), x) + 0.25 * pdf(Normal(10.0, 1.0), x))
end

function g(x)
  log(pdf(Normal(x,0,.5)))
end

function h(x)
  log(pdf(Normal(x,-1,.35)) + 2*pdf(Normal(x,1,.5)) + 1.5*pdf(Normal(x,2,.25)))
end

function mcmc(x::Float64, g::Function, sampler::Function, niter::Int64)
  xs = zeros(niter)
  gx = g(x)
  for iter = 1:niter
    x,gx = sampler(x,g,gx)
    xs[iter] = x
  end
  return xs
end

@elapsed samples = mcmc(0.0, f, slice_sampler, 50_000)
csvwrite("examples/results/slice.f.csv", samples)

@elapsed samples = mcmc(x0,g,slice_sampler,niter)
csvwrite("examples/results/slice.g.csv",samples)

@elapsed samples = mcmc(x0,h,slice_sampler,niter)
csvwrite("examples/results/slice.h.csv",samples)

@elapsed samples = mcmc(x0,g,mh_sampler,niter)
csvwrite("examples/results/mh.g.csv",samples)

@elapsed samples = mcmc(x0,h,mh_sampler,niter)
csvwrite("examples/results/mh.h.csv",samples)


text:
# ---
# title: 601. Human Traffic of Stadium
# id: problem601
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories:
# link: <https://leetcode.com/problems/human-traffic-of-stadium/description/>
# hidden: true
# ---
#
# Table: `Stadium`
#
#
#
#     +---------------+---------+
#     | Column Name   | Type    |
#     +---------------+---------+
#     | id            | int     |
#     | visit_date    | date    |
#     | people        | int     |
#     +---------------+---------+
#     visit_date is the primary key for this table.
#     Each row of this table contains the visit date and visit id to the stadium with the number of people during the visit.
#     No two rows will have the same visit_date, and as the id increases, the dates increase as well.
#
#
#
#
# Write an SQL query to display the records with three or more rows with
# **consecutive** `id`'s, and the number of people is greater than or equal to
# 100 for each.
#
# Return the result table ordered by `visit_date` in **ascending order**.
#
# The query result format is in the following example.
#
#
#
#
#
#     Stadium table:
#     +------+------------+-----------+
#     | id   | visit_date | people    |
#     +------+------------+-----------+
#     | 1    | 2017-01-01 | 10        |
#     | 2    | 2017-01-02 | 109       |
#     | 3    | 2017-01-03 | 150       |
#     | 4    | 2017-01-04 | 99        |
#     | 5    | 2017-01-05 | 145       |
#     | 6    | 2017-01-06 | 1455      |
#     | 7    | 2017-01-07 | 199       |
#     | 8    | 2017-01-09 | 188       |
#     +------+------------+-----------+
#
#     Result table:
#     +------+------------+-----------+
#     | id   | visit_date | people    |
#     +------+------------+-----------+
#     | 5    | 2017-01-05 | 145       |
#     | 6    | 2017-01-06 | 1455      |
#     | 7    | 2017-01-07 | 199       |
#     | 8    | 2017-01-09 | 188       |
#     +------+------------+-----------+
#     The four rows with ids 5, 6, 7, and 8 have consecutive ids and each of them has >= 100 people attended. Note that row 8 was included even though the visit_date was not the next day after row 7.
#     The rows with ids 2 and 3 are not included because we need at least three consecutive ids.
#
#
## @lc code=start
using LeetCode

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


text:
"""
    Neighbor

Represents a communications neighbor. Constructors are defined by the
backends and require a unique ID that represents the neighbor, the array
type to be used for staging (usually `Array` on the CPU and `CuArray`
for GPU code), the element type for the array, and the dimensions for
the send and receive buffers.
"""
struct Neighbor{B}
    pid::Int
    send_buf::B
    recv_buf::B

    function Neighbor{B}(pid, send_buf, recv_buf) where {B}
        new{B}(pid, send_buf, recv_buf)
    end
end
Neighbor(::Type{CC}, x...) where {CC <: AbstractCommsContext} =
    error("No `Neighbor` constructor defined for $CC")
pid(nbr::Neighbor) = nbr.pid
send_buffer(nbr::Neighbor) = nbr.send_buf
recv_buffer(nbr::Neighbor) = nbr.recv_buf

"""
    send_stage(nbr::Neighbor)

Return the staging area for data that is to be sent to the specified
neighbor.
"""
send_stage(nbr::Neighbor) = get_stage(nbr.send_buf)

"""
    recv_stage(nbr::Neighbor)

Return the staging area for data that has been received from the
specified neighbor.
"""
recv_stage(nbr::Neighbor) = get_stage(nbr.recv_buf)


text:
"""
```
imgb = binarize(MinimumIntermodes(), img)
binarize!(MinimumIntermodes(), img)
```

nder the assumption that the histogram is bimodal the images histogram is
iterativly smoothed untill two modes remain. The threshold is then set to the
minimum value between the two thresholds. Falls back on
`HistogramThresholding.UnimodalRosin()` if only one mode is found.

# Arguments

The function arguments are described in more detail below.

##  `img`

An `AbstractGray` image which is to be binarized into background (`Gray(0)`)
and foreground (`Gray(1)`) pixels.

## Reference

1. C. A. Glasbey, “An Analysis of Histogram-Based Thresholding Algorithms,” *CVGIP: Graphical Models and Image Processing*, vol. 55, no. 6, pp. 532–537, Nov. 1993. [doi:10.1006/cgip.1993.1040](https://doi.org/10.1006/cgip.1993.1040)
2. J. M. S. Prewitt and M. L. Mendelsohn, “THE ANALYSIS OF CELL IMAGES*,” *Annals of the New York Academy of Sciences*, vol. 128, no. 3, pp. 1035–1053, Dec. 2006. [doi:10.1111/j.1749-6632.1965.tb11715.x](https://doi.org/10.1111/j.1749-6632.1965.tb11715.x)
"""
function binarize(algorithm::MinimumIntermodes,  img::AbstractArray{T,2}) where T <: AbstractGray
  binarize!(algorithm, copy(img))
end

function binarize!(algorithm::MinimumIntermodes,  img::AbstractArray{T,2}) where T <: AbstractGray
  edges, counts = build_histogram(img,  256)
  t = find_threshold(HistogramThresholding.MinimumIntermodes(), counts[1:end], edges)
  map!(img,img) do val
    val < t ? Gray(0) : Gray(1)
  end
end


text:
export make_node, score

import StatsFuns: logπ
import SpecialFunctions: loggamma
import Distributions: Wishart

@node Wishart Stochastic [ out, (ν, aliases = [ df ]), (S, aliases = [ scale ]) ]

@average_energy Wishart (q_out::Any, q_ν::Any, q_S::Any) = begin
    d = dim(q_out)

    m_q_ν   = mean(q_ν)
    m_q_S   = mean(q_S)
    m_q_out = mean(q_out)

    return 0.5 * (
        m_q_ν * (logdet(m_q_S) + d * log(2)) -
        logdet(m_q_out) * (m_q_ν- d - 1.0) +
        tr(cholinv(m_q_S) * m_q_out) + 0.5 * d * (d - 1) * logπ
    ) + mapreduce(i -> loggamma(0.5 * (m_q_ν + 1.0 - i)), +, 1:d)
end


text:
import DataStructures.OrderedSet

flatten(token) = string(token)
flatten(token::AbstractString) = token
flatten(token::Array) = map(string, token)
flatten(token::Set) = map(string, collect(token))

# the following doesn't work in Julia v0.5
# flatten(token::Dict) = map(string, vcat(map(collect, token)...))
function flatten(token::Dict)
    r=AbstractString[]
    for (k,v) in token
        push!(r, string(k))
        push!(r, string(v))
    end
    r
end

function flatten(token::Tuple{T, U}...) where {T <: Number, U <: AbstractString}
    r=[]
    for item in token
        push!(r, item[1])
        push!(r, item[2])
    end
    r
end

flatten_command(command...) = vcat(map(flatten, command)...)

######## Type Conversions #########

convert_response(::Type{Float64}, response::T) where {T <: AbstractString} = parse(Float64, response)::Float64
convert_response(::Type{Float64}, response::T) where {T <: Real} = float(response)::Float64
convert_response(::Type{Bool}, response::AbstractString) = response == "OK" || response == "QUEUED" ? true : false
convert_response(::Type{Bool}, response::Integer) = response == 1 ? true : false
convert_response(::Type{Set{AbstractString}}, response) = Set{AbstractString}(response)
convert_response(::Type{OrderedSet{AbstractString}}, response) = OrderedSet{AbstractString}(response)

function convert_response(::Type{Dict{AbstractString, AbstractString}}, response)
    iseven(length(response)) || throw(ClientException("Response could not be converted to Dict"))
    retdict = Dict{AbstractString, AbstractString}()
    for i=1:2:length(response)
        retdict[response[i]] = response[i+1]
    end
    retdict
end

function convert_eval_response(::Any, response)
    return response
end

# import Base: ==
# ==(A::Union{T, Nothing}, B::Union{U, Nothing}) where {T<:AbstractString, U<:AbstractString} = A == B
# ==(A::Union{T, Nothing}, B::Union{U, Nothing}) where {T<:Number, U<:Number} = A == B

convert_response(::Type{AbstractString}, response) = string(response)
convert_response(::Type{Integer}, response) = response

function convert_response(::Type{Array{AbstractString, 1}}, response)
    r = Array{AbstractString, 1}()
    for item in response
        push!(r, item)
    end
    r
end

function convert_response(::Type{Union{T, Nothing}}, response) where {T<:Number}
    return response
end

function convert_response(::Type{Union{T, Nothing}}, response) where {T <: AbstractString}
    return response
end

# redundant
function convert_response(::Type{Array{Union{T, Nothing}, 1}}, response) where {T<:Number}
    if response == nothing
        Array{Union{T, Nothing}, 1}()
   else
        r = Array{Union{T, Nothing}, 1}()
        for item in response
            push!(r, tryparse(T, item))
        end
        r
    end
end

function convert_response(::Type{Array{Union{T, Nothing}, 1}}, response) where {T <: AbstractString}
    if response == nothing
        Array{Union{T, Nothing}, 1}()
   else
        r = Array{Union{T, Nothing}, 1}()
        for item in response
            push!(r, item)
        end
        r
    end
end

function open_transaction(conn::RedisConnection)
    t = TransactionConnection(conn)
    multi(t)
    t
end

function reset_transaction(conn::TransactionConnection)
    discard(conn)
    multi(conn)
end

function open_pipeline(conn::RedisConnection)
    PipelineConnection(conn)
end

function read_pipeline(conn::PipelineConnection)
    result = Any[]
    for i=1:conn.num_commands
        push!(result, read_reply(conn))
    end
    conn.num_commands = 0
    result
end

nullcb(err) = @debug err
function open_subscription(conn::RedisConnection, err_callback=nullcb)
    s = SubscriptionConnection(conn)
    @async subscription_loop(s, err_callback)
    s
end

function subscription_loop(conn::SubscriptionConnection, err_callback::Function)
    while is_connected(conn)
        try
            l = getline(conn.socket)
            reply = parseline(l, conn.socket)
            message = SubscriptionMessage(reply)
            if message.message_type == SubscriptionMessageType.Message
                conn.callbacks[message.channel](message.message)
            elseif message.message_type == SubscriptionMessageType.Pmessage
                conn.pcallbacks[message.channel](message.message)
            end
        catch err
            err_callback(err)
        end
    end
end

macro redisfunction(command::AbstractString, ret_type, args...)
    func_name = esc(Symbol(command))
    command = lstrip(command,'_')
    command = split(command, '_')

    if length(args) > 0
        return quote
            function $(func_name)(conn::RedisConnection, $(args...))
                response = execute_command(conn, flatten_command($(command...), $(args...)))
                convert_response($ret_type, response)
            end
            function $(func_name)(conn::TransactionConnection, $(args...))
                execute_command(conn, flatten_command($(command...), $(args...)))
            end
            function $(func_name)(conn::PipelineConnection, $(args...))
                execute_command_without_reply(conn, flatten_command($(command...), $(args...)))
                conn.num_commands += 1
            end
        end
    else
        return quote
            function $(func_name)(conn::RedisConnection)
                response = execute_command(conn, flatten_command($(command...)))
                convert_response($ret_type, response)
            end
            function $(func_name)(conn::TransactionConnection)
                execute_command(conn, flatten_command($(command...)))
            end
            function $(func_name)(conn::PipelineConnection)
                execute_command_without_reply(conn, flatten_command($(command...)))
                conn.num_commands += 1
            end
        end
    end
end

macro sentinelfunction(command, ret_type, args...)
    func_name = esc(Symbol(string("sentinel_", command)))
    return quote
        function $(func_name)(conn::SentinelConnection, $(args...))
            response = execute_command(conn, flatten_command("sentinel", $command, $(args...)))
            convert_response($ret_type, response)
        end
    end
end


text:
function calc_grad_Haux(pw::PWGrid, Potentials, Focc, evals, psi::Array{ComplexF64,2}, kT)
    Nstates = size(Focc)[1]
    g = zeros(ComplexF64,Nstates,Nstates)
    #
    denum_dmudeta = 0.0
    for ist = 1:Nstates
        denum_dmudeta = denum_dmudeta + Focc[ist]*(1 - Focc[ist])
    end
    #println("denum_dmudeta =", denum_dmudeta)
    #
    # Diagonal elements
    Hpsi = op_H(pw, Potentials, psi)
    dFdmu = 0.0 + im*0.0
    #
    for ist = 1:Nstates
        term1 = (dot( psi[:,ist], Hpsi[:,ist] ) - evals[ist])/kT*Focc[ist]*(1 - Focc[ist])
        #println("term1 = ", term1)
        g[ist,ist] = -term1
        dFdmu = dFdmu + term1
    end
    for ist = 1:Nstates
        g[ist,ist] = g[ist,ist] + Focc[ist]*(1 - Focc[ist])/denum_dmudeta * dFdmu
    end
    #
    # Off diagonals
    for ist = 1:Nstates
        for jst = (ist+1):Nstates
            pref = (Focc[jst] - Focc[ist])/(evals[jst] - evals[ist])
            #@printf("pref = %2d %2d %18.10f\n", ist, jst, pref)
            g[ist,jst] = pref*dot( psi[:,ist], Hpsi[:,jst] )
            g[jst,ist] = conj(g[ist,jst])
        end
    end
    #println("Pass 36 sum g Haux = ", sum(g))
    return g
end


text:
"""
    mutable struct Event{T <: SSDFloat}

Collection struct for individual events. This (mutable) struct is ment to be used to look at individual events,
not to process a hugh amount of events.
"""
mutable struct Event{T <: SSDFloat}
    locations::Union{Vector{<:AbstractCoordinatePoint{T}}, Missing}
    energies::Union{Vector{T}, Missing}
    drift_paths::Union{Vector{EHDriftPath{T}}, Missing}
    waveforms::Union{Vector{<:Any}, Missing}
    Event{T}() where {T <: SSDFloat} = new{T}(missing, missing, missing, missing)
end

function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity{T}} = ones(T, length(locations)))::Event{T} where {T <: SSDFloat}
    evt = Event{T}()
    evt.locations = locations
    evt.energies = to_internal_units(internal_energy_unit, energies)
    evt.waveforms = missing
    return evt
end

function Event(evt::NamedTuple{(:evtno, :detno, :thit, :edep, :pos),
         <:Tuple{
            Union{Integer, AbstractVector{<:Integer}},
            Union{Integer, AbstractVector{<:Integer}},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:AbstractVector{<:RealQuantity}}
        }}, T = missing)
    if ismissing(T) T = eltype(to_internal_units(internal_energy_unit, evt[:edep][:])) end

    event = Event(
        CartesianPoint{T}.(to_internal_units.(internal_length_unit, evt[:pos][:])),
        T.(to_internal_units.(internal_energy_unit, evt[:edep][:]))
    )

    return event
end

in(evt::Event, detector::SolidStateDetector) = all( pt -> pt in detector, evt.locations)
in(evt::Event, simulation::Simulation) = all( pt -> pt in simulation.detector, evt.locations)

function drift_charges!(event::Event{T}, sim::Simulation{T}; max_nsteps::Int = 1000, Δt::RealQuantity = 5u"ns", verbose::Bool = true)::Nothing where {T <: SSDFloat}
    event.drift_paths = drift_charges(sim, CartesianPoint.(event.locations), Δt = Δt, max_nsteps = max_nsteps, verbose = verbose)
    nothing
end
function get_signal!(event::Event{T}, sim::Simulation{T}, contact_id::Int; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}
    @assert !ismissing(event.drift_paths) "No drit path for this event. Use `drift_charges(event::Event, sim::Simulation` first."
    @assert !ismissing(sim.weighting_potentials[contact_id]) "No weighting potential for contact $(contact_id). Use `calculate_weighting_potential!(sim, contact_id)` first."
    if ismissing(event.waveforms)
        event.waveforms = Union{Missing, RadiationDetectorSignals.RDWaveform}[missing for i in eachindex(sim.detector.contacts)]
    end
    event.waveforms[contact_id] = get_signal(sim, event.drift_paths, event.energies, contact_id)
    nothing
end
function get_signals!(event::Event{T}, sim::Simulation{T}; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}
    @assert !ismissing(event.drift_paths) "No drit path for this event. Use `drift_charges(event::Event, sim::Simulation` first."
    if ismissing(event.waveforms)
        event.waveforms = Union{Missing, RadiationDetectorSignals.RDWaveform}[missing for i in eachindex(sim.detector.contacts)]
    end
    for contact in sim.detector.contacts
        @assert !ismissing(sim.weighting_potentials[contact.id]) "No weighting potential for contact $(contact.id). Use `calculate_weighting_potential!(sim, contact_id)` first."
        event.waveforms[contact.id] = get_signal(sim, event.drift_paths, event.energies, contact.id, Δt = Δt)
    end
    nothing
end

@recipe function f(wvs::Vector{Union{Missing, RadiationDetectorSignals.RDWaveform}})
    @series begin
        RadiationDetectorSignals.RDWaveform[wv for wv in skipmissing(wvs)]
    end
end

function simulate!(event::Event{T}, sim::Simulation{T}; max_nsteps::Int = 1000, Δt::RealQuantity = 5u"ns", verbose::Bool = true)::Nothing where {T <: SSDFloat}
    drift_charges!(event, sim, max_nsteps = max_nsteps, Δt = Δt, verbose = verbose)
    get_signals!(event, sim, Δt = Δt)
    nothing
end

function add_baseline_and_extend_tail(wv::RadiationDetectorSignals.RDWaveform{T,U,TV,UV}, n_baseline_samples::Int, total_waveform_length::Int) where {T,U,TV,UV}
    new_signal::Vector{eltype(UV)} = Vector{eltype(UV)}(undef, total_waveform_length)
    new_signal[1:n_baseline_samples] .= zero(eltype(UV))
    if length(wv.value) <= total_waveform_length - n_baseline_samples
        new_signal[n_baseline_samples+1:n_baseline_samples+length(wv.value)] = wv.value
        new_signal[n_baseline_samples+length(wv.value)+1:end] .= wv.value[end]
    else
        new_signal[n_baseline_samples+1:end] = wv.value[1:total_waveform_length - n_baseline_samples]
    end
    new_times = if TV <: AbstractRange
        range( zero(first(wv.time)), step = step(wv.time), length = total_waveform_length )
    else
        error("Not yet definted for timestamps of type `$(TV)`")
    end
    return RDWaveform( new_times, new_signal )
end


text:
using FactCheck

facts("MoralVector") do
    context("Create an Vector") do
        @fact MoralVector() --> not(nothing)
        @fact MoralVector(2) --> not(nothing)
        @fact MoralVector([1.,2,3]) --> not(nothing)
    end
end

facts("BasicSociety") do
    context("Create a BasicSociety") do
        @fact BasicSociety() --> not(nothing)
        # @fact BasicSociety(3) --> not(nothing)
        # @fact BasicSociety(rand(3,3)) --> not(nothing)
        # @fact BasicSociety(rand(3,2)) --> nothing
    end
end