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

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

# ╔═╡ 1ef174fa-16f0-11eb-328a-afc201effd2f
using Pkg, Printf

# ╔═╡ 55cfdab8-d792-11ea-271f-e7383e19997c
using PlutoUI;

# ╔═╡ 9e509f80-d485-11ea-0044-c5b7e750aacb
using NiLang

# ╔═╡ 37ed073a-d492-11ea-156f-1fb155128d0f
using Zygote, BenchmarkTools

# ╔═╡ 4d75f302-d492-11ea-31b9-bbbdb43f344e
using NiLang.AD

# ╔═╡ 627ea2fb-6530-4ea0-98ee-66be3db54411
html"""

  <svg class="octicon octicon-mark-github v-align-middle" height="32" viewBox="0 0 16 16" version="1.1" width="32" aria-hidden="true"><path fill-rule="evenodd" d="M8 0C3.58 0 0 3.58 0 8c0 3.54 2.29 6.53 5.47 7.59.4.07.55-.17.55-.38 0-.19-.01-.82-.01-1.49-2.01.37-2.53-.49-2.69-.94-.09-.23-.48-.94-.82-1.13-.28-.15-.68-.52-.01-.53.63-.01 1.08.58 1.23.82.72 1.21 1.87.87 2.33.66.07-.52.28-.87.51-1.07-1.78-.2-3.64-.89-3.64-3.95 0-.87.31-1.59.82-2.15-.08-.2-.36-1.02.08-2.12 0 0 .67-.21 2.2.82.64-.18 1.32-.27 2-.27.68 0 1.36.09 2 .27 1.53-1.04 2.2-.82 2.2-.82.44 1.1.16 1.92.08 2.12.51.56.82 1.27.82 2.15 0 3.07-1.87 3.75-3.65 3.95.29.25.54.73.54 1.48 0 1.07-.01 1.93-.01 2.2 0 .21.15.46.55.38A8.013 8.013 0 0016 8c0-4.42-3.58-8-8-8z"></path></svg>

Download this notebook

"""

# ╔═╡ 94b2b962-e02a-11ea-09a5-81b3226891ed
md"""# 连猩猩都能懂的可逆编程
### (Reversible programming made simple)
[https://github.com/JuliaReverse/NiLangTutorial/](https://github.com/JuliaReverse/NiLangTutorial/)

$(html"
")

**Jinguo Liu** (github: [GiggleLiu](https://github.com/GiggleLiu/))

*Postdoc, Institute of physics, Chinese academy of sciences* (when doing this project)

*Consultant, QuEra Computing* (current)

*Postdoc, Havard* (soon)
"""

# ╔═╡ a5ee60c8-e02a-11ea-3512-7f481e499f23
md"""
# Table of Contents
1. Reversible programming basics
2. Differentiate everything with a reversible programming language
4. Real world applications and benchmarks
"""

# ╔═╡ a11c4b60-d77d-11ea-1afe-1f2ab9621f42
md"""
## In this talk,
We use the reversible eDSL [NiLang](https://github.com/GiggleLiu/NiLang.jl) is a [Julia](https://julialang.org/) as our reversible programming tool.

A package that can differentiate everything.

![NiLang](https://raw.githubusercontent.com/GiggleLiu/NiLang.jl/master/docs/src/asset/logo3.png)

Authors:
[GiggleLiu](https://github.com/GiggleLiu), [Taine Zhao](https://github.com/thautwarm)
"""

# ╔═╡ e54a1be6-d485-11ea-0262-034c56e0fda8
md"""
## Sec I. Reversible programming basic

### Reversible function definition

A reversible function `f` is defined as
```julia
(~f)(f(x, y, z...)...) == (x, y, z...)
```
"""

# ╔═╡ d1628f08-ddfb-11ea-241a-c7e6c1a22212
md"""
##  Example 1: reversible adder
```math
\begin{align}
f &: x, y → x+y, y\\
{\small \mathrel{\sim}}f &: x, y → x-y, y
\end{align}
```
"""

# ╔═╡ 278ac6b6-e02c-11ea-1354-cd7ecd1099be
md"The reversible macro `@i` defines two functions, the function itself and its inverse."

# ╔═╡ a28d38be-d486-11ea-2c40-a377b74a05c1
@i function reversible_plus(x, y)
	x += y
end

# ╔═╡ e93f0bf6-d487-11ea-1baa-21d51ddb4a20
reversible_plus(2.0, 3.0)

# ╔═╡ fc932606-d487-11ea-303e-75ca8b7a02f6
(~reversible_plus)(5.0, 3.0)

# ╔═╡ e3d2b23a-ddfb-11ea-0f5e-e72ed299bb45
md"## The difference to a regular programming language"

# ╔═╡ a961e048-ddf2-11ea-0262-6d19eb82b36b
md"**Comment 1**: The return statement is not allowed, a reversible function returns input arguments directly."

# ╔═╡ 2d22f504-ddf1-11ea-28ec-5de6f4ee79bb
md"**Comment 2**: Every operation is reversible. `+=` is considered as reversible for integers and floating point numbers in NiLang, although for floating point numbers, there are *rounding errors*."

# ╔═╡ 7d08ac24-e143-11ea-2085-539fd9e35889
md"### A case where `+=` is not reversible"

# ╔═╡ 9fcdd77c-e0df-11ea-09e6-49a2861137e5
let
	x, y = 1e-20, 1e20
	x += y
	x -= y
	(x, y)
end

# ╔═╡ 0a1a8594-ddfc-11ea-119a-1997c86cd91b
md"""
## Use this function
"""

# ╔═╡ 0b4edb1a-ddf0-11ea-220c-91f2df7452e7
@i function reversible_plus2(x, y)
	reversible_plus(x, y)  # equivalent to `reversible_plus(x, y)`
	reversible_plus(x, y)
end

# ╔═╡ f875ecd6-ddef-11ea-22a1-619809d15b37
md"**Comment**: Inside a reversible function definition, a statement changes a variable *inplace*"

# ╔═╡ e7557bee-e0cc-11ea-1788-411e759b4766
reversible_plus2(2.0, 3.0)

# ╔═╡ cd7b2a2e-ddf5-11ea-04c4-f7583bbb5a53
md"A statement can be **uncalled** with `~`"

# ╔═╡ bc98a824-ddf5-11ea-1a6a-1f795452d3d0
@i function do_nothing(x, y)
	reversible_plus(x, y)
	~reversible_plus(x, y)  # uncall the expression
end

# ╔═╡ 05f8b91c-e0cd-11ea-09e3-f3c5c0e07e63
do_nothing(2.0, 3.0)

# ╔═╡ ac302844-e07b-11ea-35dd-e3e06054401b
md"## Example 2: Compute $x^5$"

# ╔═╡ b722e098-e07b-11ea-3483-01360fb6954e
@i function naive_power5(y, x::T) where T
	y = one(T)   # error 1: `=` is not reversible
	for i=1:5
		y *= x   # error 2: `*=` is not reversible
	end
end

# ╔═╡ bf8b722c-dfa4-11ea-196a-719802bc23c5
md"""
## Compute $x^5$ reversibly
"""

# ╔═╡ 330edc28-dfac-11ea-35a5-3144c4afbfcf
md"note: `*=` is not reversible for usual number systems"

# ╔═╡ 0a679e04-dfa7-11ea-0288-a1fa490c4387
@i function power5(x5, x4, x3, x2, x1, x)
	x1 += x
	x2 += x1 * x
	x3 += x2 * x
	x4 += x3 * x
	x5 += x4 * x
end

# ╔═╡ cc32cae8-dfab-11ea-0d0b-c70ea8de720a
power5(0.0, 0.0, 0.0, 0.0, 0.0, 2.0)

# ╔═╡ b4240c16-dfac-11ea-3a40-33c54436e3a3
md"# Don't make me so many input arguments!"

# ╔═╡ ade52358-dfac-11ea-2dd3-d3a691e7a8a2
@i function power5_twoinputs(x5, x::T) where T
	x1 ← zero(T)
	x2 ← zero(T)
	x3 ← zero(T)
	x4 ← zero(T)
	x1 += x
	x2 += x1 * x
	x3 += x2 * x
	x4 += x3 * x

	x5 += x4 * x

	x4 -= x3 * x
	x3 -= x2 * x
	x2 -= x1 * x
	x1 -= x
	x4 → zero(T)
	x3 → zero(T)
	x2 → zero(T)
	x1 → zero(T)
end

# ╔═╡ d86e2e5e-dfab-11ea-0053-6d52f1164bc5
power5_twoinputs(0.0, 2.0)

# ╔═╡ 7951b9ec-e030-11ea-32ee-b1de49378186
md"""
**Comment**:
`n ← zero(T)` is the variable allocation operation. It means
```
if n is defined
	error
else
	n = zero(T)
end
```
Its inverse is `n → zero(T)`. It means
```
@assert n == zero(T)
deallocate(n)
```
"""

# ╔═╡ 6bc97f5e-dfad-11ea-0c43-e30b6620e6e8
md"# Shorter: compute-copy-uncompute"

# ╔═╡ 80d24e9e-dfad-11ea-1dae-49568d534f10
@i function power5_twoinputs_shorter(x5, x::T) where T
	@routine begin  # compute
		@zeros T x1 x2 x3 x4
		x1 += x
		x2 += x1 * x
		x3 += x2 * x
		x4 += x3 * x
	end

	x5 += x4 * x   # copy

	~@routine    # uncompute
end

# ╔═╡ a8092b18-dfad-11ea-0989-474f37d05f73
power5_twoinputs_shorter(0.0, 2.0)

# ╔═╡ 43f0c2fc-e030-11ea-25d9-b323e6496a35
md"""**Comment**:
```
@routine statement
~@routine
```

is equivalent to
```
statement
~(statement)
```
This is the famous `compute-copy-uncompute` design pattern in reversible computing. Check this [reference](https://epubs.siam.org/doi/10.1137/0219046).
"""

# ╔═╡ b4ad5830-dfad-11ea-0057-055dda8cc9be
md"# How to compute x^1000?"

# ╔═╡ cf576d38-dfad-11ea-2682-7bd540db44a5
@i function power1000(x1000, x::T) where T
	@routine begin
		xs ← zeros(T, 1000)
		xs[1] += 1
		for i=2:1000
			xs[i] += xs[i-1] * x
		end
	end

	x1000 += xs[1000] * x

	~@routine
end

# ╔═╡ 35fff53c-dfae-11ea-3602-918a17d5a5fa
power1000(0.0, 1.001)

# ╔═╡ 9b9b5328-e030-11ea-1d00-f3341572734a
html"""

For loop

	Forward

	for i=start:step:stop
		# do something
	end

	Reverse

	for i=stop:-step:start
		# undo something
	end

"""

# ╔═╡ f3b87892-e080-11ea-353d-8d81c52cf9ac
md"### Sometimes, a for loop can break down"

# ╔═╡ b27a3974-e030-11ea-0bcd-7f7035d55165
@i function power1000_bad(x1000, x::T) where T
	@routine begin
		xs ← zeros(T, 1000)
		xs[1] += 1
		for i=2:length(xs)
			xs[i] += xs[i-1] * x
			PUSH!(xs, @const zero(T))
		end
	end

	x1000 += xs[1000] * x

	~@routine
end

# ╔═╡ e5d47096-e030-11ea-1e87-5b9b1dbecfe0
power1000_bad(0.0, 1.001)

# ╔═╡ 9c62289a-dfae-11ea-0fe0-b1cb80a87704
md"#  Don't allocate for me!"

# ╔═╡ 88838bce-dfaf-11ea-1a72-7d15629cfcb0
md"""
Multipling two unsigned logarithmic numbers `x = exp(lx)` and `y = exp(ly)`
```
x * y = exp(lx) * exp(ly) = exp(lx + ly)
```
"""

# ╔═╡ a593f970-dfae-11ea-2d79-876030850dee
@i function power1000_noalloc(x1000, x::T) where T
	if x!= 0
		@routine begin
			absx ← zero(T)
			lx ← one(ULogarithmic{T})
			lx1000 ← one(ULogarithmic{T})
			absx += abs(x)
			lx *= convert(absx)
			for i=1:1000
				lx1000 *= lx
			end
		end
		x1000 += convert(lx1000)
		~@routine
	end
end

# ╔═╡ f448548e-dfaf-11ea-05c0-d5d177683445
power1000_noalloc(0.0, 1.001)

# ╔═╡ 65cd13ca-e031-11ea-3fc6-977792eb5f8c
html"""

If statement

	Forward

	if (precondition, postcondition)
		# do A
	else
		# do B
	end

	Reverse

	if (postcondition, precondition)
		# undo A
	else
		# undo B
	end

"""

# ╔═╡ 53c02100-e08f-11ea-1f5d-8b2311b095d2
md"![](https://user-images.githubusercontent.com/6257240/116341762-78a31080-a7af-11eb-8376-d2ba0bf2b454.png)"

# ╔═╡ 75751b24-e0b8-11ea-2b37-9d138121345c
md"### You should not do"

# ╔═╡ 76b84de4-e031-11ea-0bcf-39b86a6b4552
@i function break_if(x)
	if x%2 == 1
		x += 1
	else
		x -= 1
	end
end

# ╔═╡ b1984d24-e031-11ea-3b13-3bd0119a2bcb
break_if(3)

# ╔═╡ 7f163d82-e0b8-11ea-2fe7-332bb4dee586
md"### You should do"

# ╔═╡ ddc6329e-e031-11ea-0e6e-e7332fa26e22
@i function happy_if(x)
	if (x%2 == 1, x%2 == 0)
		x += 1
	else
		x -= 1
	end
end

# ╔═╡ f3d5e1b0-e031-11ea-1a90-7bed88e28bad
happy_if(3)

# ╔═╡ ab67419a-dfae-11ea-27ba-09321303ad62
md"""# Wrap up

1. reversible arithmetic instructions `+=` and `-=`, besides, we have `SWAP`, `NEG`, `INC` and `ROT` et. al.
2. inverse statement `~`
3. there is no "`=`" operation in reversible computing, use "`←`" to allocate a new variable, and use "`→`" to deallocate an pre-emptied variable.
4. compute-uncompute macro `@routine` and `~@routine`
5. reversible control flow: `for` loop and `if` statement, the `while` statement is also available.
6. logarithmic number is reversible under `*=` and `/=`
"""

# ╔═╡ d5c2efbc-d779-11ea-11ad-1f5873b95628
md"""
![yeah](https://media1.tenor.com/images/40147f2eac14c0a7f18c34ecba73fa34/tenor.gif?itemid=7805520)
"""
# ![yeah](https://pic.chinesefontdesign.com/uploads/2017/03/chinesefontdesign.com_2017-03-07_08-19-24.gif)

# ╔═╡ 30af9642-e084-11ea-1f92-b52abfddcf06
md"# Sec II. Automatic differentiation in NiLang

### References
* Nextjournal [https://nextjournal.com/giggle/reverse-checkpointing](https://nextjournal.com/giggle/reverse-checkpointing)

* arXiv: 2003.04617
"

# ╔═╡ db1fab1c-e084-11ea-0bf0-b1fbe9e74b3f
html"""
Automatic differentiation?
"""

# ╔═╡ e1370f80-e0bc-11ea-2a90-d50cc762cbcb
md"When we start learning AD, we start by learning the backward rules of the matrix multiplication"

# ╔═╡ 3098411c-e0bc-11ea-2754-eb0afbd663de
function mymul!(out::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
	@assert size(A, 2) == size(B, 1) && size(out) == (size(A, 1), size(B, 2))
	for k=1:size(B, 2)
		for j=1:size(B, 1)
			for i=size(A, 1)
				@inbounds out[i, k] += A[i, j] * B[j, k]
			end
		end
	end
	return out
end

# ╔═╡ 3d0150ee-e0bd-11ea-0a5a-339465b496dc
md"Then, we learning how to use chain rule to chain different utilities."

# ╔═╡ 8016ff94-e0bc-11ea-3b9e-4f0676587edf
md"##### But wait! Why don't we start from the backward rules of `+` and `*`, then use the chain rule to derive the backward rule for matrix multiplication?"

# ╔═╡ 99108ace-e0bc-11ea-2744-d1b18db50ae1
md"# They are different"

# ╔═╡ b2337f26-e0bb-11ea-3da0-9507c35101ae
md"""
### Domain-specific autodiff (DS-AD)
* **Tensor**Flow
* **PyTorch**
* **Jax**
* **Flux (Zygote backended)**

### General Purposed autodiff (GP-AD)
* **Tapenade**
* **NiLang**
"""

# ╔═╡ 48db515c-e084-11ea-2eec-018b8545fa34
md"## Traditional AD uses checkpointing
Checkpoint every 100 steps. Blue and yellow objects are computing and re-computing. Here states 1 and state 101 are cached. Blue objects are computing, and yellow ones are re-computing. The state 100 is the desired state.
"

# ╔═╡ f531f556-e083-11ea-2f7e-77e110d6c53a
md"![](https://nextjournal.com/data/Qmes4v3ic2VrYQt6W9mWu4p6W53Gd1DmbDcYCuafbwTe7Y?filename=checkpointing.png&content-type=image/png)"

# ╔═╡ 62643fbc-e084-11ea-1b1f-39b87ff32b9e
md"## Reverse Computing
Reversible computing approach to free up memories (a) when no operations are reversible. (b) when all operations are reversible. Blue and yellow diamonds are reversible operations executed in forward and backward directions, red cubics are garbage variables.
"

# ╔═╡ 0bf54b08-e084-11ea-3d11-7be65f3ec022
md"![](https://nextjournal.com/data/QmPsgm4Z4mqVw2h2eC3RkGf96xTQp13KE9rdzmPeUe5KWN?filename=reversecomputing.png&content-type=image/png)"

# ╔═╡ 15f7c60a-e08e-11ea-31ea-a5cd055644db
md"## Difference Explained"

# ╔═╡ 55a3a260-d48e-11ea-06e2-1b7bd7bba6f5
md"""
![adprog](https://github.com/GiggleLiu/NiLang.jl/raw/master/docs/src/asset/adprog.png)
"""

# ╔═╡ 38014ad0-e08e-11ea-1905-198038ab7e5f
md"# Obtaining the gradient of norm in Zygote"

# ╔═╡ 2e6fe4da-d79d-11ea-1e90-f5215190395c
md"**Obtaining the gradient of the norm function**"

# ╔═╡ 6560c28c-e08e-11ea-1094-d333b88071ce
function regular_norm(x::AbstractArray{T}) where T
	res = zero(T)
	for i=1:length(x)
		@inbounds res += x[i]^2
	end
	return sqrt(res)
end

# ╔═╡ 744dd3c6-d492-11ea-0ed5-0fe02f99db1f
@benchmark Zygote.gradient($regular_norm, $(randn(1000))) seconds=1

# ╔═╡ f72246f8-e08e-11ea-3aa0-53f47a64f3e9
md"## The reversible counterpart"

# ╔═╡ f025e454-e08e-11ea-20d6-d139b9a6b301
@i function reversible_norm(res, y, x::AbstractArray{T}) where {T}
	for i=1:length(x)
		@inbounds y += x[i]^2
	end
	res += sqrt(y)
end

# ╔═╡ 8fedd65a-e08e-11ea-27f4-03bf9ed65875
let x = randn(1000)
	@assert Zygote.gradient(regular_norm, x)[1] ≈ NiLang.AD.gradient(reversible_norm, (0.0, 0.0, x), iloss=1)[3]
end

# ╔═╡ 8ad60dc0-d492-11ea-2cb3-1750b39ddf86
@benchmark NiLang.AD.gradient($reversible_norm, (0.0, 0.0, $(randn(1000))), iloss=1)

# ╔═╡ 7bab4614-d77e-11ea-037c-8d1f432fc3b8
md"""
![yeah](https://media1.tenor.com/images/40147f2eac14c0a7f18c34ecba73fa34/tenor.gif?itemid=7805520)
"""
# ![yeah](https://pic.chinesefontdesign.com/uploads/2017/03/chinesefontdesign.com_2017-03-07_08-19-24.gif)

# ╔═╡ fcca27ba-d4a4-11ea-213a-c3e2305869f1
#**1. The bundle adjustment jacobian benchmark**
#$(LocalResource("ba-origin.png"))
#![ba](https://github.com/JuliaReverse/NiBundleAdjustment.jl/raw/master/benchmarks/adbench.png)

#**2. The Gaussian mixture model benchmark**
#$(LocalResource("gmm-origin.png"))
#![gmm](https://github.com/JuliaReverse/NiGaussianMixture.jl/raw/master/benchmarks/adbench.png)

md"""
# Sec III. Applications in real world and benchmarks
"""

# ╔═╡ 519dc834-e092-11ea-2151-57ef23810b84
md"""
## 1. Bundle Adjustment (Jacobian)
![bundle adjustment](https://encrypted-tbn0.gstatic.com/images?q=tbn%3AANd9GcRgpGSCWRjHDDaIYQX5ejhMvyKY_GFhynVoQg&usqp=CAU)

*Srajer, Filip, Zuzana Kukelova, and Andrew Fitzgibbon. "A benchmark of selected algorithmic differentiation tools on some problems in computer vision and machine learning." Optimization Methods and Software 33.4-6 (2018): 889-906.*

### Benchmarks
**Devices**
* CPU: Intel(R) Xeon(R) Gold 6230 CPU @ 2.10GHz
* GPU: Nvidia Titan V.

**Github Repos**
* [https://github.com/microsoft/ADBench](https://github.com/microsoft/ADBench)
* [https://github.com/JuliaReverse/NiBundleAdjustment.jl](https://github.com/JuliaReverse/NiBundleAdjustment.jl)
"""

# ╔═╡ c89108f0-e092-11ea-0fe2-efad85008b28
html"""

"""

# ╔═╡ 2ec4c700-e093-11ea-06ff-47d2c21a068f
md"""##### NiLang.AD and Tapenade
![](https://user-images.githubusercontent.com/6257240/116341804-907a9480-a7af-11eb-934f-7eb94803f5f2.png)"""

# ╔═╡ 474aa228-e092-11ea-042b-bdfaeb99f16f
md"""
## 2. Gaussian Mixture Model (Gradient)
![gmm](https://prateekvjoshi.files.wordpress.com/2013/06/multimodal.jpg)
"""

# ╔═╡ 2baaff10-d56c-11ea-2a23-bfa3a7ae2e4b
md"""
### Benchmarks
*Srajer, Filip, Zuzana Kukelova, and Andrew Fitzgibbon. "A benchmark of selected algorithmic differentiation tools on some problems in computer vision and machine learning." Optimization Methods and Software 33.4-6 (2018): 889-906.*

**Devices**
* CPU: Intel(R) Xeon(R) Gold 6230 CPU @ 2.10GHz

**Github Repos**
* [https://github.com/microsoft/ADBench](https://github.com/microsoft/ADBench)
* [https://github.com/JuliaReverse/NiGaussianMixture.jl](https://github.com/JuliaReverse/NiGaussianMixture.jl)
"""

# ╔═╡ 102fbf2e-d56b-11ea-189d-c78d56c0a924
html"""

Results from the original benchmark

"""

# ╔═╡ cc0d5622-d788-11ea-19cd-3bf6864d9263
md"""##### Including NiLang.AD
![](https://github.com/JuliaReverse/NiLangTutorial/blob/master/notebooks/asset/benchmarks_gmm.png?raw=true)"""

# ╔═╡ a1646ef0-e091-11ea-00f1-e7c246e191ff
md"## 3. Solve the memory wall problem in machine learning"

# ╔═╡ b18b3ae8-e091-11ea-24a1-e968b70b217c
html"""
Learning a ring distribution with NICE network, before and after training

References

arXiv: 1410.8516

NiLang's documentation

"""

# ╔═╡ bf3774de-e091-11ea-3372-ef56452158e6
md"""
## 4. Solve the spinglass ground state configuration
Obtaining the optimal configuration of a spinglass problem on a $28 \times 28$ square lattice.

![](https://user-images.githubusercontent.com/6257240/116342088-067efb80-a7b0-11eb-935e-0b5e29010a22.png)

##### References
Jin-Guo Liu, Lei Wang, Pan Zhang, **arXiv 2008.06888**
"""

# ╔═╡ c8e4f7a6-e091-11ea-24a3-4399635a41a5
md"""
## 5. Optimizing problems in finance
Gradient based optimization of Sharpe rate.

600x acceleration comparing with using pure Zygote.

##### References
* Han Li's Github repo: [https://github.com/HanLi123/NiLang](https://github.com/HanLi123/NiLang) and his Zhihu blog [猴子掷骰子](https://zhuanlan.zhihu.com/c_1092471228488634368).
"""

# ╔═╡ bc872296-e09f-11ea-143b-9bfd5e52b14f
md"""## 6. Accelerate the performance critical part of variational mean field

[https://github.com/quantumlang/NiLangTest/pull/1](https://github.com/quantumlang/NiLangTest/pull/1)

600x acceleration comparing with using pure Zygote.
"""

# ╔═╡ e7b21fce-e091-11ea-180c-7b42e00598a9
md"""# Thank you!
Special thanks to my collaborator **Taine Zhao** and (ex-)advisor **Lei Wang**.

QuEra computing (a quantum computing company located in Boston) is hiring people.
"""

# ╔═╡ 7c79975c-d789-11ea-30b1-67ff05418cdb
md"""
![yeah](https://media1.tenor.com/images/40147f2eac14c0a7f18c34ecba73fa34/tenor.gif?itemid=7805520)
"""
# ![yeah](https://pic.chinesefontdesign.com/uploads/2017/03/chinesefontdesign.com_2017-03-07_08-19-24.gif)

# ╔═╡ 5f1c3f6c-d48b-11ea-3eb0-357fd3ece4fc
md"""
## Sec IV. More about number systems

* Integers are reversible under (`+=`, `-=`).
* Floating point number system is **irreversible** under (`+=`, `-=`) and (`*=`, `/=`).
* [Fixedpoint number system](https://github.com/JuliaMath/FixedPointNumbers.jl) are reversible under (`+=`, `-=`)
* [Logarithmic number system](https://github.com/cjdoris/LogarithmicNumbers.jl) is reversible under (`*=`, `/=`)
"""

# ╔═╡ 11ddebfe-d488-11ea-223a-e9403f6ec8de
md"""
##### Example 1: Affine transformation with rounding error

```julia
y = A * x + b
```
"""

# ╔═╡ 030e592e-d488-11ea-060d-97a3bb6353b7
@i function reversible_affine!(y!::AbstractVector{T}, W::AbstractMatrix{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T
    @safe @assert size(W) == (length(y!), length(x)) && length(b) == length(y!)
    for j=1:size(W, 2)
        for i=1:size(W, 1)
            @inbounds y![i] += W[i,j]*x[j]
        end
    end
    for i=1:size(W, 1)
        @inbounds y![i] += b[i]
    end
end

# ╔═╡ c8d26856-d48a-11ea-3cd3-1124cd172f3a
begin
	W = randn(10, 10)
	b = randn(10)
	x = randn(10)
end;

# ╔═╡ 37c4394e-d489-11ea-174c-b13bdddbe741
yout, Wout, bout, xout = reversible_affine!(zeros(10), W, b, x)

# ╔═╡ fef54688-d48a-11ea-340b-295b88d21382
# should be restored to 0, but not!
yin, Win, bin, xin = (~reversible_affine!)(yout, Wout, bout, xout)

# ╔═╡ 259a2852-d48c-11ea-0f01-b9634850e09d
md"""
### Reversible arithmetic functions

Computing basic functions like `power`, `exp` and `besselj` is not trivial for reversible programming.
There is no efficient constant memory algorithm using pure fixed point numbers only.
"""

# ╔═╡ f06fb004-d79f-11ea-0d60-8151019bf8c7
md"""
##### Example 2: Computing power function
To compute `x ^ n` reversiblly with fixed point numbers,
we need to either allocate a vector of size $O(n)$ or suffer from polynomial time overhead. It does not show the advantage to checkpointing.
"""

# ╔═╡ 26a8a42c-d7a1-11ea-24a3-45bc6e0674ea
@i function i_power_cache(y!::T, x::T, n::Int) where T
    @routine @invcheckoff begin
        cache ← zeros(T, n)  # allocate a buffer of size n
		cache[1] += x
        for i=2:n
            cache[i] += cache[i-1] * x
        end
    end

    y! += cache[n]

    ~@routine  # uncompute cache
end

# ╔═╡ 399552c4-d7a1-11ea-36bb-ad5ca42043cb
# To check the function
i_power_cache(Fixed43(0.0), Fixed43(0.99), 100)

# ╔═╡ 4bb19760-d7bf-11ea-12ed-4d9e4efb3482
md"""
##### Example 3: reversible thinker, the logarithmic number approach

With **logarithmic numbers**, we can still utilize reversibility. Fixed point numbers and logarithmic numbers can be converted via "a fast binary logarithm algorithm".

##### References
* [1] C. S. Turner, "A Fast Binary Logarithm Algorithm", IEEE Signal Processing Mag., pp. 124,140, Sep. 2010.
"""

# ╔═╡ 5a8ba8f4-d493-11ea-1839-8ba81f86799d
@i function i_power_lognumber(y::T, x::T, n::Int) where T
    @routine @invcheckoff begin
        lx ← one(ULogarithmic{T})
        ly ← one(ULogarithmic{T})
        ## convert `x` to a logarithmic number
        ## Here, `*=` is reversible for log numbers
        lx *= convert(x)
        for i=1:n
            ly *= lx
        end
    end

    ## convert back to fixed point numbers
    y += convert(ly)

    ~@routine
end

# ╔═╡ a625a922-d493-11ea-1fe9-bdd4a694cde0
# To check the function
i_power_lognumber(Fixed43(0.0), Fixed43(0.99), 100)

# ╔═╡ 4fd20ed2-d7a2-11ea-206e-13799234913f
md"**Less allocation, better speed**"

# ╔═╡ 692dfb44-d7a1-11ea-00da-af6550bc0622
@benchmark i_power_cache(Fixed43(0.0), Fixed43(0.99), 100)

# ╔═╡ 7e4ee09c-d7a1-11ea-0e56-c1921012bc30
@benchmark i_power_lognumber(Fixed43(0.0), Fixed43(0.99), 100)

# ╔═╡ 4c209bbe-d7b1-11ea-0628-33eb8d664f5b
md"""##### Example 4: The first kind Bessel function computed with Taylor expansion
```math
J_\nu(z) = \sum\limits_{n=0}^{\infty} \frac{(z/2)^\nu}{\Gamma(k+1)\Gamma(k+\nu+1)} (-z^2/4)^{n}
```

"""

# ╔═╡ fd44a3d4-d7a4-11ea-24ea-09456ff2c53d
@i function ibesselj(y!::T, ν, z::T; atol=1e-8) where T
	if z == 0
		if ν == 0
			out! += 1
		end
	else
		@routine @invcheckoff begin
			k ← 0
			@ones ULogarithmic{T} lz halfz halfz_power_2 s
			@zeros T out_anc
			lz *= convert(z)
			halfz *= lz / 2
			halfz_power_2 *= halfz ^ 2
			# s *= (z/2)^ν/ factorial(ν)
			s *= halfz ^ ν
			for i=1:ν
				s /= i
			end
			out_anc += convert(s)
			while (s.log > -25, k!=0) # upto precision e^-25
				k += 1
				# s *= 1 / k / (k+ν) * (z/2)^2
				@routine begin
					@zeros Int kkv kv
					kv += k+ ν
					kkv += kv*k
				end
				s *= halfz_power_2 / kkv
				if k%2 == 0
					out_anc += convert(s)
				else
					out_anc -= convert(s)
				end
				~@routine
			end
		end
		y! += out_anc
		~@routine
	end
end

# ╔═╡ 84272664-d7b7-11ea-2e37-dffd2023d8d6
md"z = $(@bind z Slider(0:0.01:10; default=1.0))"

# ╔═╡ 900e2ea4-d7b8-11ea-3511-6f12d95e638a
begin
	y = ibesselj(Fixed43(0.0), 2, Fixed43(z))[1]
	gz = NiLang.AD.gradient(Val(1), ibesselj, (Fixed43(0.0), 2, Fixed43(z)))[3]
end;

# ╔═╡ d76be888-d7b4-11ea-2989-2174682ead76
let
	md"""
| ``z`` | ``y`` | ``\partial y/\partial z`` |
| ----  | ----- | -------- |
| $(@sprintf "%.5f" z) | $(@sprintf "%.5f" y) | $(@sprintf "%.5f" gz) |
"""
end

# ╔═╡ 85c9edcc-d789-11ea-14c8-71697cd6a047
md"""
![yeah](https://media1.tenor.com/images/40147f2eac14c0a7f18c34ecba73fa34/tenor.gif?itemid=7805520)
"""
# ![yeah](https://pic.chinesefontdesign.com/uploads/2017/03/chinesefontdesign.com_2017-03-07_08-19-24.gif)

# ╔═╡ Cell order:
# ╟─1ef174fa-16f0-11eb-328a-afc201effd2f
# ╟─627ea2fb-6530-4ea0-98ee-66be3db54411
# ╟─94b2b962-e02a-11ea-09a5-81b3226891ed
# ╟─a5ee60c8-e02a-11ea-3512-7f481e499f23
# ╟─a11c4b60-d77d-11ea-1afe-1f2ab9621f42
# ╟─e54a1be6-d485-11ea-0262-034c56e0fda8
# ╟─55cfdab8-d792-11ea-271f-e7383e19997c
# ╟─d1628f08-ddfb-11ea-241a-c7e6c1a22212
# ╠═9e509f80-d485-11ea-0044-c5b7e750aacb
# ╟─278ac6b6-e02c-11ea-1354-cd7ecd1099be
# ╠═a28d38be-d486-11ea-2c40-a377b74a05c1
# ╠═e93f0bf6-d487-11ea-1baa-21d51ddb4a20
# ╠═fc932606-d487-11ea-303e-75ca8b7a02f6
# ╟─e3d2b23a-ddfb-11ea-0f5e-e72ed299bb45
# ╟─a961e048-ddf2-11ea-0262-6d19eb82b36b
# ╟─2d22f504-ddf1-11ea-28ec-5de6f4ee79bb
# ╟─7d08ac24-e143-11ea-2085-539fd9e35889
# ╠═9fcdd77c-e0df-11ea-09e6-49a2861137e5
# ╟─0a1a8594-ddfc-11ea-119a-1997c86cd91b
# ╠═0b4edb1a-ddf0-11ea-220c-91f2df7452e7
# ╟─f875ecd6-ddef-11ea-22a1-619809d15b37
# ╠═e7557bee-e0cc-11ea-1788-411e759b4766
# ╟─cd7b2a2e-ddf5-11ea-04c4-f7583bbb5a53
# ╠═bc98a824-ddf5-11ea-1a6a-1f795452d3d0
# ╠═05f8b91c-e0cd-11ea-09e3-f3c5c0e07e63
# ╟─ac302844-e07b-11ea-35dd-e3e06054401b
# ╠═b722e098-e07b-11ea-3483-01360fb6954e
# ╟─bf8b722c-dfa4-11ea-196a-719802bc23c5
# ╟─330edc28-dfac-11ea-35a5-3144c4afbfcf
# ╠═0a679e04-dfa7-11ea-0288-a1fa490c4387
# ╠═cc32cae8-dfab-11ea-0d0b-c70ea8de720a
# ╟─b4240c16-dfac-11ea-3a40-33c54436e3a3
# ╠═ade52358-dfac-11ea-2dd3-d3a691e7a8a2
# ╠═d86e2e5e-dfab-11ea-0053-6d52f1164bc5
# ╟─7951b9ec-e030-11ea-32ee-b1de49378186
# ╟─6bc97f5e-dfad-11ea-0c43-e30b6620e6e8
# ╠═80d24e9e-dfad-11ea-1dae-49568d534f10
# ╠═a8092b18-dfad-11ea-0989-474f37d05f73
# ╟─43f0c2fc-e030-11ea-25d9-b323e6496a35
# ╟─b4ad5830-dfad-11ea-0057-055dda8cc9be
# ╠═cf576d38-dfad-11ea-2682-7bd540db44a5
# ╠═35fff53c-dfae-11ea-3602-918a17d5a5fa
# ╟─9b9b5328-e030-11ea-1d00-f3341572734a
# ╟─f3b87892-e080-11ea-353d-8d81c52cf9ac
# ╠═b27a3974-e030-11ea-0bcd-7f7035d55165
# ╠═e5d47096-e030-11ea-1e87-5b9b1dbecfe0
# ╟─9c62289a-dfae-11ea-0fe0-b1cb80a87704
# ╟─88838bce-dfaf-11ea-1a72-7d15629cfcb0
# ╠═a593f970-dfae-11ea-2d79-876030850dee
# ╠═f448548e-dfaf-11ea-05c0-d5d177683445
# ╟─65cd13ca-e031-11ea-3fc6-977792eb5f8c
# ╟─53c02100-e08f-11ea-1f5d-8b2311b095d2
# ╟─75751b24-e0b8-11ea-2b37-9d138121345c
# ╠═76b84de4-e031-11ea-0bcf-39b86a6b4552
# ╠═b1984d24-e031-11ea-3b13-3bd0119a2bcb
# ╟─7f163d82-e0b8-11ea-2fe7-332bb4dee586
# ╠═ddc6329e-e031-11ea-0e6e-e7332fa26e22
# ╠═f3d5e1b0-e031-11ea-1a90-7bed88e28bad
# ╟─ab67419a-dfae-11ea-27ba-09321303ad62
# ╟─d5c2efbc-d779-11ea-11ad-1f5873b95628
# ╟─30af9642-e084-11ea-1f92-b52abfddcf06
# ╟─db1fab1c-e084-11ea-0bf0-b1fbe9e74b3f
# ╟─e1370f80-e0bc-11ea-2a90-d50cc762cbcb
# ╠═3098411c-e0bc-11ea-2754-eb0afbd663de
# ╟─3d0150ee-e0bd-11ea-0a5a-339465b496dc
# ╟─8016ff94-e0bc-11ea-3b9e-4f0676587edf
# ╟─99108ace-e0bc-11ea-2744-d1b18db50ae1
# ╟─b2337f26-e0bb-11ea-3da0-9507c35101ae
# ╟─48db515c-e084-11ea-2eec-018b8545fa34
# ╟─f531f556-e083-11ea-2f7e-77e110d6c53a
# ╟─62643fbc-e084-11ea-1b1f-39b87ff32b9e
# ╟─0bf54b08-e084-11ea-3d11-7be65f3ec022
# ╟─15f7c60a-e08e-11ea-31ea-a5cd055644db
# ╟─55a3a260-d48e-11ea-06e2-1b7bd7bba6f5
# ╟─38014ad0-e08e-11ea-1905-198038ab7e5f
# ╟─2e6fe4da-d79d-11ea-1e90-f5215190395c
# ╠═6560c28c-e08e-11ea-1094-d333b88071ce
# ╠═37ed073a-d492-11ea-156f-1fb155128d0f
# ╠═744dd3c6-d492-11ea-0ed5-0fe02f99db1f
# ╟─f72246f8-e08e-11ea-3aa0-53f47a64f3e9
# ╠═f025e454-e08e-11ea-20d6-d139b9a6b301
# ╠═4d75f302-d492-11ea-31b9-bbbdb43f344e
# ╠═8fedd65a-e08e-11ea-27f4-03bf9ed65875
# ╠═8ad60dc0-d492-11ea-2cb3-1750b39ddf86
# ╟─7bab4614-d77e-11ea-037c-8d1f432fc3b8
# ╟─fcca27ba-d4a4-11ea-213a-c3e2305869f1
# ╟─519dc834-e092-11ea-2151-57ef23810b84
# ╟─c89108f0-e092-11ea-0fe2-efad85008b28
# ╟─2ec4c700-e093-11ea-06ff-47d2c21a068f
# ╟─474aa228-e092-11ea-042b-bdfaeb99f16f
# ╟─2baaff10-d56c-11ea-2a23-bfa3a7ae2e4b
# ╟─102fbf2e-d56b-11ea-189d-c78d56c0a924
# ╟─cc0d5622-d788-11ea-19cd-3bf6864d9263
# ╟─a1646ef0-e091-11ea-00f1-e7c246e191ff
# ╟─b18b3ae8-e091-11ea-24a1-e968b70b217c
# ╟─bf3774de-e091-11ea-3372-ef56452158e6
# ╟─c8e4f7a6-e091-11ea-24a3-4399635a41a5
# ╟─bc872296-e09f-11ea-143b-9bfd5e52b14f
# ╟─e7b21fce-e091-11ea-180c-7b42e00598a9
# ╟─7c79975c-d789-11ea-30b1-67ff05418cdb
# ╟─5f1c3f6c-d48b-11ea-3eb0-357fd3ece4fc
# ╟─11ddebfe-d488-11ea-223a-e9403f6ec8de
# ╠═030e592e-d488-11ea-060d-97a3bb6353b7
# ╠═c8d26856-d48a-11ea-3cd3-1124cd172f3a
# ╠═37c4394e-d489-11ea-174c-b13bdddbe741
# ╠═fef54688-d48a-11ea-340b-295b88d21382
# ╟─259a2852-d48c-11ea-0f01-b9634850e09d
# ╟─f06fb004-d79f-11ea-0d60-8151019bf8c7
# ╠═26a8a42c-d7a1-11ea-24a3-45bc6e0674ea
# ╠═399552c4-d7a1-11ea-36bb-ad5ca42043cb
# ╟─4bb19760-d7bf-11ea-12ed-4d9e4efb3482
# ╠═5a8ba8f4-d493-11ea-1839-8ba81f86799d
# ╠═a625a922-d493-11ea-1fe9-bdd4a694cde0
# ╟─4fd20ed2-d7a2-11ea-206e-13799234913f
# ╠═692dfb44-d7a1-11ea-00da-af6550bc0622
# ╠═7e4ee09c-d7a1-11ea-0e56-c1921012bc30
# ╟─4c209bbe-d7b1-11ea-0628-33eb8d664f5b
# ╠═fd44a3d4-d7a4-11ea-24ea-09456ff2c53d
# ╟─84272664-d7b7-11ea-2e37-dffd2023d8d6
# ╠═900e2ea4-d7b8-11ea-3511-6f12d95e638a
# ╟─d76be888-d7b4-11ea-2989-2174682ead76
# ╟─85c9edcc-d789-11ea-14c8-71697cd6a047


text:
function UnivariateProcessNoise(arg0::RealMatrix, arg1::LOFType, arg2::PositionAngle, arg3::Vector{UnivariateFunction}, arg4::Vector{UnivariateFunction})
    return UnivariateProcessNoise((RealMatrix, LOFType, PositionAngle, Vector{UnivariateFunction}, Vector{UnivariateFunction}), arg0, arg1, arg2, arg3, arg4)
end

function get_lof_cartesian_orbital_parameters_evolution(obj::UnivariateProcessNoise)
    return jcall(obj, "getLofCartesianOrbitalParametersEvolution", Vector{UnivariateFunction}, ())
end

function get_lof_type(obj::UnivariateProcessNoise)
    return jcall(obj, "getLofType", LOFType, ())
end

function get_position_angle(obj::UnivariateProcessNoise)
    return jcall(obj, "getPositionAngle", PositionAngle, ())
end

function get_process_noise_matrix(obj::UnivariateProcessNoise, arg0::SpacecraftState, arg1::SpacecraftState)
    return jcall(obj, "getProcessNoiseMatrix", RealMatrix, (SpacecraftState, SpacecraftState), arg0, arg1)
end

function get_propagation_parameters_evolution(obj::UnivariateProcessNoise)
    return jcall(obj, "getPropagationParametersEvolution", Vector{UnivariateFunction}, ())
end

ext: jl


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using Legendre, GLMakie
using GeometryBasics, LinearAlgebra, StatsBase
using Makie: get_dim, surface_normals
GLMakie.activate!() # HIDE
let
    function Y(θ, ϕ, l, m)
        if m < 0
            return (-1)^m * √2 * Nlm(l, abs(m)) * Plm(l, abs(m), cos(θ)) * sin(abs(m)*ϕ)
        elseif m == 0
            return sqrt((2*l+1)/4π)*Plm(l, m, cos(θ))
        else
            return (-1)^m * √2 * Nlm(l, m) * Plm(l, m, cos(θ)) * cos(m*ϕ)
        end
    end
    function getMesh(x,y,z)
        positions = vec(map(CartesianIndices(z)) do i
        GeometryBasics.Point{3, Float32}(
            get_dim(x, i, 1, size(z)),
            get_dim(y, i, 2, size(z)),
            z[i])
        end)
        faces = decompose(GLTriangleFace, Rect2D(0f0, 0f0, 1f0, 1f0), size(z))
        normals = surface_normals(x, y, z)
        vertices = GeometryBasics.meta(positions; normals=normals)
        meshObj = GeometryBasics.Mesh(vertices, faces)
        meshObj
    end
    # Grids of polar and azimuthal angles
    θ = LinRange(0, π, 200)
    ϕ = LinRange(0, 2π, 200)
    x = [sin(θ)*sin(ϕ) for θ in θ, ϕ in ϕ]
    y = [sin(θ)*cos(ϕ) for θ in θ, ϕ in ϕ]
    z = [cos(θ)        for θ in θ, ϕ in ϕ]
    l = 4
    m = 1
    ambient =  Vec3f0(0.75, 0.75, 0.75)
    cmap = (:dodgerblue, :white)
    fig = Figure(resolution = (1200, 800), backgroundcolor = :black)

    Ygrid = [Y(θ, ϕ, l, m) for θ in θ, ϕ in ϕ]
    Ylm = abs.(Ygrid)
    Ygrid2 = vec(Ygrid)

    ax1 = Axis3(fig, aspect = :data)
    ax2 = Axis3(fig, aspect = :data)
    pltobj1 = mesh!(ax1, getMesh(x, y, z), color = Ygrid2, colormap = cmap,
        ambient = ambient)
    pltobj2 = mesh!(ax2, getMesh(Ylm .* x, Ylm .* y, Ylm .* z), color = Ygrid2,
        colormap = cmap,  ambient = ambient)
    hidespines!(ax1)
    hidedecorations!(ax1)
    hidespines!(ax2)
    hidedecorations!(ax2)
    cbar = Colorbar(fig, pltobj1, label = "Yₗₘ(θ,ϕ)", tickwidth = 2,
        labelcolor = :white, ticklabelcolor = :white, tickcolor = :grey,
        width = 25,ticklabelsize = 30,labelsize = 30,ticksize=25,tickalign = 1)
    fig[1,1] = ax1
    fig[1,2] = ax2
    fig[1,3] = cbar
    fig[0,1:2] = Label(fig, "Tesseral Spherical Harmonics l = $(l), m = $(m)",
        textsize = 30, color = (:white, 0.85))
    fig
    colgap!(fig.layout, 2)
    fig
    save(joinpath(@__DIR__, "output", "tesseralSphericalH.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie", "Legendre","GeometryBasics", "LinearAlgebra", "StatsBase", "Makie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
function useLockModel(conf, counts, workloadName)
	tps,J,f,K,regions,D,S0,S,g = prepareLockModel(1, conf, counts[1,:], workloadName);
	nRows = size(counts, 1);
	all_R =zeros(nRows, convert(Int,J));
	all_T_total =zeros(nRows, 1);
	all_M_total =zeros(nRows, 1);
	all_Vp =zeros(nRows, 9);
	all_V =zeros(nRows, 9);
	all_W =zeros(nRows, 1);
	all_Pcon =zeros(nRows, convert(Int,J));
	all_totalWaits =zeros(nRows, convert(Int,J));
	all_TimeSpentWaiting =zeros(nRows, convert(Int,J));
	all_LocksBeingHeld =zeros(nRows, convert(Int,J));

	for row=1:nRows
		tps,J,f,K,regions,D,S0,S,g = prepareLockModel(1, conf, counts[row,:], workloadName);
		filter = f ./ f;
		for j=1:length(filter)
			if isnan(filter[j])
				filter[j] = 0;
			end
		end

		R,T_total,M_total,Vp,V,W,U,Pcon,L = basicThomasian(tps, J, f, K, regions, D, S0, S, g);
		#% R_j=latency of tran type J
		#% T_total = givenT; % the requested TPS
		#% M_total=total number of transactions in the system
		#% Vp_i= mean waiting time w.r.t. the active transactions in DBR i
		#% V_i=blocking time when transactions encounter a lock conflict in DBR i
		#% W=mean waiting time per lock conflict
		#% U_jn=mean delay incurred by trans type j when they encounter a lock conflict at step n
		#% Pcon_j=prob of lock conflict PER LOCK request for transactions of type j
		#% L_ji= of locks held by trans of type j in DBR i

		all_R[row,:]=R .* filter;
		if(length(T_total)==1)
			T_total=T_total[1]
		end
		all_T_total[row,:]= T_total;
		if(length(M_total)==1)
			M_total=M_total[1]
		end
		all_M_total[row,:] = M_total;
		all_Vp[row,:] = Vp;
		all_V[row,:] = V;
		if(length(W)==1)
			W=W[1]
		end
		all_W[row,:] = W;
		all_Pcon[row,:] = Pcon .* filter;
		all_totalWaits[row,:] = Pcon .* K .* f * T_total[1] .* filter;
		all_TimeSpentWaiting[row,:] = Pcon .* K .* f * T_total[1] * W .* filter;

		all_LocksBeingHeld[row, :] = sum((f * T_total[1]) * L, 2)[1];
	end

	all_totalWaits2 = zeros(size(all_totalWaits,1),1)
	all_TimeSpentWaiting2 = zeros(size(all_TimeSpentWaiting,1),1)
	all_LocksBeingHeld2 = zeros(size(all_LocksBeingHeld,1),1)
	for i=1:size(all_totalWaits,1)
		all_totalWaits2[i] = dot(vec(all_totalWaits[i,:]),vec(counts[i,:]))
	end
	for i=1:size(all_TimeSpentWaiting,1)
		all_TimeSpentWaiting2[i] = dot(vec(all_TimeSpentWaiting[i,:]),vec(counts[i,:]))
	end
	for i=1:size(all_LocksBeingHeld,1)
		all_LocksBeingHeld2[i] = dot(vec(all_LocksBeingHeld[i,:]),vec(counts[i,:]))
	end

	return typeOfPredictions(all_R,all_T_total,all_M_total,all_Vp, all_V,all_W, all_Pcon,  all_totalWaits2,all_TimeSpentWaiting2,all_LocksBeingHeld2)

end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
module PAPI

using NUMA_jll

try
    include(joinpath(dirname(@__DIR__), "deps","deps.jl"))
catch e
    error("PAPI.jl not properly configured, please run `Pkg.build(\"PAPI\")`.")
end

const REFCOUNT = Ref(zero(UInt))
function deref_shutdown()
    REFCOUNT[] -= 1
    if REFCOUNT[] == 0
        # no objects to be finalized
        ccall((:PAPI_shutdown, libpapi), Cvoid, ())
    end
end

include("constants.jl")
include("error.jl")
include("events.jl")
include("eventsets.jl")
include("counters.jl")
include("components.jl")
include("sampling.jl")
include("prettyprint.jl")
include("serialization.jl")
include("numa.jl")

function __init__()
    papi_current_version = (papi_version.major << 24) | (papi_version.minor << 16)

    # init the library and make sure that some counters are available
    rv = ccall((:PAPI_library_init, libpapi), Cint, (Cint,), papi_current_version)
    if rv != papi_current_version
      if rv > 0
        error("PAPI library version mismatch!")
      else
        throw(PAPIError(rv))
      end
    end

    if num_counters() == 0
        error("PAPI init error: No counters are available on the current system")
    end

    REFCOUNT[] += 1
    atexit(deref_shutdown)
end

"""
Get the number of hardware counters available on the system

`PAPI.num_counters()` initializes the PAPI library if necessary.

`PAPI_num_counters()`` returns the number of hardware counters available on the system.
"""
function num_counters()
    errcode = ccall((:PAPI_num_hwctrs, libpapi), Cint, ())
    if errcode < 0
        throw(PAPIError(errcode))
    else
        Int(errcode)
    end
end

export PAPIError, num_counters, start_counters, read_counters!, accum_counters!, stop_counters, stop_counters!
export @profile, @sample, sample, profile, @numaprofile, numaprofile
export name_to_event, event_to_name, @event_str
export find_component, exists, available_presets, available_native, Event, Counts, EventSet
export load, save
end # module


text:
function train(rng, m::AbstractDPM, train, predict; iter = 2000, warmup = iter ÷ 2, thin = 1)
    gammachain = [zeros(Bool, m.G) for _ = 1:(iter - warmup) ÷ thin]
    fchain = [zeros(length(predict.y)) for _ = 1:(iter - warmup) ÷ thin]
    for t in 1:iter
        update!(rng, m, train)
        if t > warmup && ((t - warmup) % thin == 0)
            t0 = (t - warmup) ÷ thin
            gammachain[t0] = m.gamma[:]
            for i = 1:length(predict.y)
                fchain[t0][i] = predlik(m, train, predict, i)
            end
        end
    end
    return DPM_MCMC_Chain(gammachain, fchain)
end

function predlik(m::AbstractDPM, train, predict, i::Int)
    k̄ = first(passive_clusters(m))
    A = active_clusters(m)
    n = cluster_sizes(m)
    N = length(train.y)
    α = dp_mass(m)
    ans = 0.0
    for k in A
        ans += exp(logpredlik(m, train, predict, i, k)) * n[k] / (N + α)
    end
    ans += exp(logpredlik(m, train, predict, i, k̄)) * α / (N + α)
    return(ans)
end


text:
module StaticUnivariatePolynomials

export
    AbstractPolynomial,
    Polynomial,
    BernsteinPolynomial

using Base: tail

include("util.jl")
include("abstract_polynomial.jl")
include("monomial_basis.jl")
include("bernstein_basis.jl")

end # module


text:
#### Update forcing

####
#### Types
####

abstract type ForcingType end

struct NoForcing <: ForcingType end

struct StandardForcing{FT} <: ForcingType
    apply_subsidence::Bool
    apply_coriolis::Bool
    coriolis_param::FT
    StandardForcing(;
        apply_subsidence,
        apply_coriolis,
        coriolis_param::FT,
    ) where {FT} = new{FT}(apply_subsidence, apply_coriolis, coriolis_param)
end

####
#### Helper functions
####

function coriolis_force!(aux, grid, q, coriolis_param)
    gm, en, ud, sd, al = allcombinations(q)
    for k in over_elems_real(grid)
        aux[:forcing_u, k, gm] -= coriolis_param * (aux[:vg, k] - q[:v, k, gm])
        aux[:forcing_v, k, gm] += coriolis_param * (aux[:ug, k] - q[:u, k, gm])
    end
end

####
#### Update forcing methods
####

"""
    update_forcing!

Define forcing fields in the `aux` state.
"""
function update_forcing! end

function update_forcing!(
    aux::StateVec,
    q::StateVec,
    grid::Grid,
    params,
    ::NoForcing,
) end

function update_forcing!(
    aux::StateVec,
    q::StateVec,
    grid::Grid,
    params,
    forcing::StandardForcing,
)
    gm, en, ud, sd, al = allcombinations(q)
    @unpack param_set = params

    for k in over_elems_real(grid)
        # Apply large-scale horizontal advection tendencies
        aux[:forcing_θ_liq, k, gm] +=
            aux[:dTdt, k] / TD.exner_given_pressure(param_set, aux[:p_0, k])
        aux[:forcing_q_tot, k, gm] += aux[:dqtdt, k]
    end
    if forcing.apply_subsidence
        for k in over_elems_real(grid)
            # Apply large-scale subsidence tendencies
            aux[:forcing_θ_liq, k, gm] -=
                grad(q[:θ_liq, Dual(k), gm], grid) * aux[:subsidence, k]
            aux[:forcing_q_tot, k, gm] -=
                grad(q[:q_tot, Dual(k), gm], grid) * aux[:subsidence, k]
        end
    end

    if forcing.apply_coriolis
        coriolis_force!(aux, grid, q, forcing.coriolis_param)
    end
end


text:
struct StanResult{M,R,C}
  model::M
  return_code::R
  chains::C
end

function Base.show(io::IO, mime::MIME"text/plain", res::StanResult)
  show(io, mime, res.chains)
end

struct StanODEData
end

function generate_priors(n,priors)
  priors_string = ""
  if priors===nothing
    for i in 1:n
      priors_string = string(priors_string,"theta_$i ~ normal(0, 1)", " ; ")
    end
  else
    for i in 1:n
      priors_string = string(priors_string,"theta_$i ~ ",stan_string(priors[i]),";")
    end
  end
  priors_string
end

function generate_theta(n,priors)
  theta = ""
  for i in 1:n
    upper_bound = ""
    lower_bound = ""
    if !isnothing(priors) && maximum(priors[i]) != Inf
      upper_bound = string("upper=",maximum(priors[i]))
    end
    if !isnothing(priors) && minimum(priors[i]) != -Inf
      lower_bound = string("lower=",minimum(priors[i]))
    end
    if lower_bound != "" && upper_bound != ""
      theta = string(theta,"real","<$lower_bound",",","$upper_bound>"," theta_$i",";")
    elseif lower_bound != ""
      theta = string(theta,"real","<$lower_bound",">"," theta_$i",";")
    elseif upper_bound != ""
      theta = string(theta,"real","<","$upper_bound>"," theta_$i",";")
    else
      theta = string(theta,"real"," theta_$i",";")
    end
  end
  return theta
end

function stan_inference(prob::DiffEqBase.DEProblem,t,data,priors = nothing,
                            stanmodel = nothing;alg=:rk45,
                            num_samples=1000, num_warmup=1000, reltol=1e-3,
                            abstol=1e-6, maxiter=Int(1e5),likelihood=Normal,
                            vars=(StanODEData(),InverseGamma(3,3)),nchains=[1],
                            sample_u0 = false, save_idxs = nothing, diffeq_string = nothing,
                            printsummary = true, output_format = :mcmcchains)

  save_idxs !== nothing && length(save_idxs) == 1 ? save_idxs = save_idxs[1] : save_idxs = save_idxs
  length_of_y = length(prob.u0)
  save_idxs = something(save_idxs, 1:length_of_y)
  length_of_params = length(vars)
  if isnothing(diffeq_string)
    sys = ModelingToolkit.modelingtoolkitize(prob)
    length_of_parameter = length(ModelingToolkit.parameters(sys)) + sample_u0 * length(save_idxs)
  else
    length_of_parameter = length(prob.p) + sample_u0 * length(save_idxs)
  end

  if stanmodel === nothing
    if alg ==:adams
      algorithm = "ode_adams_tol"
    elseif alg ==:rk45
      algorithm = "ode_rk45_tol"
    elseif alg == :bdf
      algorithm = "ode_bdf_tol"
    else
      error("The choices for alg are :adams, :rk45, or :bdf")
    end
    hyper_params = ""
    tuple_hyper_params = ""
    setup_params = ""
    thetas = ""
    theta_names = ""
    theta_string = generate_theta(length_of_parameter,priors)
    for i in 1:length_of_parameter
      thetas = string(thetas,"real theta_$i",";")
      theta_names = string(theta_names,"theta_$i",",")
    end
    for i in 1:length_of_params
      if isa(vars[i],StanODEData)
        tuple_hyper_params = string(tuple_hyper_params,"u_hat[t,$save_idxs]",",")
      else
        dist = stan_string(vars[i])
        hyper_params = string(hyper_params,"sigma$(i-1) ~ $dist;")
        tuple_hyper_params = string(tuple_hyper_params,"sigma$(i-1)",",")
        setup_params = string(setup_params,"row_vector<lower=0>[$(length(save_idxs))] sigma$(i-1);")
      end
    end
    tuple_hyper_params = tuple_hyper_params[1:length(tuple_hyper_params)-1]
    priors_string = string(generate_priors(length_of_parameter,priors))
    stan_likelihood = stan_string(likelihood)
    if sample_u0
      nu = length(save_idxs)
      dv_names_ind = findfirst("$nu", theta_names)[1]
      if nu < length(prob.u0)
        u0 = ""
        for u_ in prob.u0[nu+1:length(prob.u0)]
          u0 = u0*string(u_)
        end
        integral_string = "u_hat = $algorithm(sho, [$(theta_names[1:dv_names_ind]),$u0]', t0, ts, $reltol, $abstol, $maxiter, $(rstrip(theta_names[dv_names_ind+2:end],',')));"
      else
        integral_string = "u_hat = $algorithm(sho, [$(theta_names[1:dv_names_ind])]', t0, ts, $reltol, $abstol, $maxiter, $(rstrip(theta_names[dv_names_ind+2:end],',')));"
      end
    else
      integral_string = "u_hat = $algorithm(sho, u0, t0, ts, $reltol, $abstol, $maxiter, $(rstrip(theta_names,',')));"
    end
    binsearch_string = """
      int bin_search(real x, int min_val, int max_val){
        int range = (max_val - min_val + 1) / 2;
        int mid_pt = min_val + range;
        int out;
        while (range > 0) {
            if (x == mid_pt) {
                out = mid_pt;
                range = 0;
            } else {
                range = (range + 1) / 2;
                mid_pt = x > mid_pt ? mid_pt + range: mid_pt - range;
            }
        }
        return out;
    }
    """
    if isnothing(diffeq_string)
      diffeq_string = ModelingToolkit.build_function(
        ModelingToolkit.equations(sys),ModelingToolkit.states(sys),
        ModelingToolkit.parameters(sys),ModelingToolkit.independent_variable(sys),
          fname = :sho,
          target = ModelingToolkit.StanTarget()
      )
    end
    parameter_estimation_model = "
    functions {
      $binsearch_string
      $diffeq_string
    }
    data {
      vector[$length_of_y] u0;
      int<lower=1> T;
      real internal_var___u[T,$(length(save_idxs))];
      real t0;
      real ts[T];
    }
    parameters {
      $setup_params
      $theta_string
    }
    model{
      vector[$length_of_y] u_hat[T];
      $hyper_params
      $priors_string
      $integral_string
      for (t in 1:T){
        internal_var___u[t,:] ~ $stan_likelihood($tuple_hyper_params);
        }
    }
    "
    stanmodel = StanSample.SampleModel("parameter_estimation_model", parameter_estimation_model, nchains; printsummary = printsummary, method = StanSample.Sample(;num_samples = num_samples, num_warmup = num_warmup))
  end
  parameter_estimation_data = Dict("u0"=>prob.u0, "T" => length(t), "internal_var___u" => view(data, :, 1:length(t))', "t0" => prob.tspan[1], "ts" => t)
  rc = stan_sample(stanmodel; data = parameter_estimation_data)
  if success(rc)
    return StanResult(stanmodel, rc, read_samples(stanmodel, output_format))
  else
     rc.err
  end
end


text:
# Author: Luaksz Bratos

include("diff.jl")
using Test

# dane testowe do zadania 1
x1 = [3.0, 1.0, 5.0, 6.0]
f1 = [1.0, -3.0, 2.0, 4.0]
ans1 = [1.0, 2.0, -3.0 / 8.0, 7.0 / 40.0]

x2 = [-1.0, 0.0, 1.0, 2.0, 3.0]
f2 = [2.0, 1.0, 2.0, -7.0, 10.0]
ans2 = [2.0, -1.0, 1.0, -2.0, 2.0]

# dane testowe do zadania 2
fx2 = Diff.ilorazyRoznicowe(x2, f2)

# dane testowe do zadania 3
ans3 = [1.0, 6.0, -1.0, -6.0, 2.0]

@testset "Ilorazy roznicowe" begin
    @test isapprox(Diff.ilorazyRoznicowe(x1, f1), ans1)
    @test isapprox(Diff.ilorazyRoznicowe(x2, f2), ans2)
end

@testset "Wartosc wielomianu Newtona" begin
    @test isapprox(Diff.warNewton(x2, fx2, 2.0), -7.0)
    @test isapprox(Diff.warNewton(x2, fx2, 10.0), 13961.0)
    @test isapprox(Diff.warNewton(x2, fx2, 5.0), 506.0)
    @test isapprox(Diff.warNewton(x2, fx2, -3.0), 298.0)
end

@testset "Wspolczynniki w postaci normalnej z postaci Newtona" begin
    @test isapprox(Diff.naturalna(x2, fx2), ans3)
end


text:
using Documenter, ProjectManagement

include("create_demos.jl")
create_demofiles()

makedocs(;
    modules=[ProjectManagement],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Demos" => [
            replace(file[1:end-2], "_"=>" ") => joinpath("demos", file)
            for file in readdir(joinpath(@__DIR__, "src", "demos"))
        ],
        "API" => "API.md",
    ],
    repo="https://github.com/oxinabox/ProjectManagement.jl/blob/{commit}{path}#L{line}",
    sitename="ProjectManagement.jl",
    authors="Lyndon White",
)

deploydocs(;
    repo="github.com/oxinabox/ProjectManagement.jl",
    push_preview=true,
)


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

using Images: save, load, Colors, clamp01nan, Gray, imresize
using ImageDraw: draw, Polygon, Point

#=
    notify_user(message::AbstractString) -> AbstractString

A function to pretty print a message to the user

# Arguments

- `message::AbstractString`: Some message to print

# Returns
- `A::AbstractString`: A message to print to the user
=#
function notify_user(message::AbstractString)
    return println("\033[1;34m===>\033[0;38m\033[1;38m\t$message\033[0;38m")
end

#=
    filtered_ls(path::AbstractString) -> Array{String, 1}

A function to filter the output of readdir

# Arguments

- `path::AbstractString`: Some path to folder

# Returns

- `Array{String, 1}`: An array of filtered files in the path
=#
function filtered_ls(path::AbstractString)
    return filter!(f -> ! occursin(r".*\.DS_Store", f), readdir(path, join=true, sort=false))
end

"""
    load_image(image_path::AbstractString) -> AbstractArray
Loads an image as gray_scale

# Arguments
- `image_path::AbstractString`: Path to an image

# Returns

 - `AbstractArray`: An array of floating point values representing the image
"""
function load_image(
    image_path::AbstractString;
    scale::Bool=false,
    scale_to::Tuple=(200,200)
    )

    img = load(image_path)
    img = convert(Array{Float64}, Gray.(img))
    img = scale ? imresize(img, scale_to) : img

    return to_integral_image(img)
end

"""
    determine_feature_size(
        pos_training_path::AbstractString,
        neg_training_path::AbstractString
    ) -> Tuple{Integer, Integer, Integer, Integer, Tuple{Integer, Integer}}

Takes images and finds the best feature size for the image size.

# Arguments

- `pos_training_path::AbstractString`: the path to the positive training images
- `neg_training_path::AbstractString`: the path to the negative training images

# Returns

- `max_feature_width::Integer`: the maximum width of the feature
- `max_feature_height::Integer`: the maximum height of the feature
- `min_feature_height::Integer`: the minimum height of the feature
- `min_feature_width::Integer`: the minimum width of the feature
- `min_size_img::Tuple{Integer, Integer}`: the minimum-sized image in the image directories
"""
function determine_feature_size(
    pos_training_path::AbstractString,
    neg_training_path::AbstractString;
    scale::Bool=false,
    scale_to::Tuple=(200, 200)
)

    min_feature_height = 0
    min_feature_width = 0
    max_feature_height = 0
    max_feature_width = 0

    min_size_img = (0, 0)
    sizes = []

    for picture_dir in[pos_training_path, neg_training_path]
        for picture in filtered_ls(picture_dir)
            img = load_image(picture, scale=scale, scale_to=scale_to)
            new_size = size(img)
            sizes = push!(sizes, new_size)
        end
    end

    min_size_img = minimum(sizes)

    max_feature_height = Int(round(min_size_img[2]*(10/19)))
    max_feature_width = Int(round(min_size_img[1]*(10/19)))
    min_feature_height = Int(round(max_feature_height - max_feature_height*(2/max_feature_height)))
    min_feature_width = Int(round(max_feature_width - max_feature_width*(2/max_feature_width)))

    return max_feature_width, max_feature_height, min_feature_height, min_feature_width, min_size_img

end

@doc raw"""
    ensemble_vote(int_img::AbstractArray, classifiers::AbstractArray) -> Integer

Classifies given integral image (Abstract Array) using given classifiers.  I.e., if the sum of all classifier votes is greater 0, the image is classified positively (1); else it is classified negatively (0). The threshold is 0, because votes can be +1 or -1.

That is, the final strong classifier is

```math
h(x) = \begin{cases}
1&\text{if }\sum_{t=1}^{T}\alpha_{th_{t(x)}}\geq\frac{1}{2}\sum_{t=1}^{T}\alpha_t\\
0&\text{otherwise}
\end{cases}
\text{ where }\alpha_t = \log{\left(\frac{1}{\beta_t}\right)}
```

# Arguments

- `int_img::AbstractArray`: Integral image to be classified
- `classifiers::Array{HaarLikeObject, 1}`: List of classifiers

# Returns

- `vote::Integer`
    1       ⟺ sum of classifier votes > 0
    0       otherwise
"""
function ensemble_vote(int_img::Matrix, classifiers::Array{HaarLikeObject, 1})
    return sum(c -> get_vote(c, int_img), classifiers) ≥ zero(Int8) ? one(Int8) : zero(Int8)
end

"""
    ensemble_vote_all(int_imgs::AbstractArray, classifiers::AbstractArray) -> AbstractArray
Classifies given integral image (Abstract Array) using given classifiers.  I.e., if the sum of all classifier votes is greater 0, the image is classified positively (1); else it is classified negatively (0). The threshold is 0, because votes can be +1 or -1.

# Arguments
- `int_img::AbstractArray`: Integral image to be classified
- `classifiers::Array{HaarLikeObject, 1}`: List of classifiers

# Returns

`votes::AbstractArray`: A list of assigned votes (see ensemble_vote).
"""
function ensemble_vote_all(
    image_path::AbstractString,
    classifiers::Array{HaarLikeObject, 1};
    scale::Bool=false,
    scale_to::Tuple=(200, 200)
    )

    return [ensemble_vote(load_image(i, scale=scale, scale_to=scale_to), classifiers) for i in filtered_ls(image_path)]
end

"""
    get_faceness(feature, int_img::AbstractArray) -> Number

Get facelikeness for a given feature.

# Arguments

- `feature::HaarLikeObject`: given Haar-like feature (parameterised replacement of Python's `self`)
- `int_img::AbstractArray`: Integral image array

# Returns

- `score::Number`: Score for given feature
"""
function get_faceness(feature::HaarLikeObject, int_img::Matrix{T}) where T
    score, faceness = get_score(feature, int_img)

    return (feature.weight * score) < (feature.polarity * feature.threshold) ? faceness : zero(T)
end

#=
    reconstruct(classifiers::AbstractArray, img_size::Tuple) -> AbstractArray

Creates an image by putting all given classifiers on top of each other producing an archetype of the learned class of object.

# Arguments

- `classifiers::Array{HaarLikeObject, 1}`: List of classifiers
- `img_size::Tuple{Integer, Integer}`: Tuple of width and height

# Returns

- `result::AbstractArray`: Reconstructed image
=#
function reconstruct(classifiers::AbstractArray, img_size::Tuple)
    image = zeros(img_size)

    for c in classifiers
        # map polarity: -1 -> 0, 1 -> 1
        polarity = ((1 + c.polarity)^2)/4
        if c.feature_type == feature_types["two_vertical"]
            for x in 1:c.width
                sign = polarity
                for y in 1:c.height
                    if y >= c.height/2
                        sign = mod((sign + 1), 2)
                    end
                    image[c.top_left[2] + y, c.top_left[1] + x] += 1 * sign * c.weight
                end
            end
        elseif c.feature_type == feature_types["two_horizontal"]
            sign = polarity
            for x in 1:c.width
                if x >= c.width/2
                    sign = mod((sign + 1), 2)
                end
                for y in 1:c.height
                    image[c.top_left[1] + x, c.top_left[2] + y] += 1 * sign * c.weight
                end
            end
        elseif c.feature_type == feature_types["three_horizontal"]
            sign = polarity
            for x in 1:c.width
                if iszero(mod(x, c.width/3))
                    sign = mod((sign + 1), 2)
                end
                for y in 1:c.height
                    image[c.top_left[1] + x, c.top_left[2] + y] += 1 * sign * c.weight
                end
            end
        elseif c.feature_type == feature_types["three_vertical"]
            for x in 1:c.width
                sign = polarity
                for y in 1:c.height
                    if iszero(mod(x, c.height/3))
                        sign = mod((sign + 1), 2)
                    end
                    image[c.top_left[1] + x, c.top_left[2] + y] += 1 * sign * c.weight
                end
            end
        elseif c.feature_type == feature_types["four"]
            sign = polarity
            for x in 1:c.width
                if iszero(mod(x, c.width/2))
                    sign = mod((sign + 1), 2)
                end
                for y in 1:c.height
                    if iszero(mod(x, c.height/2))
                        sign = mod((sign + 1), 2)
                    end
                    image[c.top_left[1] + x, c.top_left[2] + y] += 1 * sign * c.weight
                end
            end
        end
    end # end for c in classifiers
    # image .-= minimum(image) # equivalent to `min(image...)`
    # image ./= maximum(image)
    # image .*= 255
    #
    # image = replace!(image, NaN=>0.0) # change NaN to white (not that there should be any NaN values)
    #
    return image
end

#=
    get_random_image(
        face_path::AbstractString,
        non_face_path::AbstractString="",
        non_faces::Bool=false
    ) -> AbstractString

Chooses a random image from a given two directories.

# Arguments

- `face_path::AbstractString`: The path to the faces directory
- `non_face_path::AbstractString`: The path to the non-faces directory

# Returns

- `file_name::AbstractString`: The path to the file randomly chosen
=#
function get_random_image(
    face_path::AbstractString;
    non_face_path::AbstractString=string(),
    non_faces::Bool=false
)
    file_name = string()

    if non_faces
        face = rand(Bool)
        file_name = rand(filter!(f -> ! occursin(r".*\.DS_Store", f), readdir(face ? face_path : non_face_path, join=true)))
    else
        file_name = rand(filter!(f -> ! occursin(r".*\.DS_Store", f), readdir(face_path, join=true)))
    end

    return file_name
end

#=
    scale_box(
        top_left::Tuple{Integer, Integer},
        bottom_right::Tuple{Integer, Integer},
        genisis_size::Tuple{Integer, Integer},
        img_size::Tuple{Integer, Integer}
    ) -> NTuple{::Tuple{Integer, Integer}, 4}

Scales the bounding box around classifiers if the image we are pasting it on is a different size to the original image.

# Arguments

- `top_left::Tuple{Integer, Integer}`: the top left of the Haar-like feature
- `bottom_right::Tuple{Integer, Integer}`: the bottom right of the Haar-like feature
- `genisis_size::Tuple{Integer, Integer}`: the size of the test images
- `img_size::Tuple{Integer, Integer}`: the size of the image which we are pasting the bounding box on top of

# Returns

- `top_left::Tuple{Integer, Integer},`: new top left of box after scaling
- `bottom_left::Tuple{Integer, Integer},`: new bottom left of box after scaling
- `bottom_right::Tuple{Integer, Integer},`: new bottom right of box after scaling
- `top_right::Tuple{Integer, Integer},`: new top right of box after scaling
=#
function scale_box(
    top_left::Tuple{Integer, Integer},
    bottom_right::Tuple{Integer, Integer},
    genisis_size::Tuple{Integer, Integer},
    img_size::Tuple{Integer, Integer}
)
    T = typeof(first(top_left))
    image_ratio = (img_size[1]/genisis_size[1], img_size[2]/genisis_size[2])

    bottom_left = (top_left[1], bottom_right[2])
    top_right = (bottom_right[1], top_left[2])

    top_left = convert.(T, round.(top_left .* image_ratio))
    bottom_right = convert.(T, round.(bottom_right .* image_ratio))
    bottom_left = convert.(T, round.(bottom_left .* image_ratio))
    top_right = convert.(T, round.(top_right .* image_ratio))

    return top_left, bottom_left, bottom_right, top_right
end

#=
    generate_validation_image(image_path::AbstractString, classifiers::AbstractArray) -> AbstractArray

Generates a bounding box around the face of a random image.

# Arguments

- `image_path::AbstractString`: The path to images
- `classifiers::Array{HaarLikeObject, 1}`: List of classifiers/haar like features

# Returns

- `validation_image::AbstractArray`: The new image with a bounding box
=#
function generate_validation_image(image_path::AbstractString, classifiers::Array{HaarLikeObject, 1})

    # === THIS FUNCTION IS A WORK IN PROGRESS ===

    img = load_image(image_path)
    img_size = size(img)

    top_lefts = [c.top_left for c in classifiers]
    bottom_rights = [c.bottom_right for c in classifiers]
    x_coords = vcat([x[1] for x in top_lefts], [x[1] for x in bottom_rights])
    y_coords = vcat([y[2] for y in top_lefts], [y[2] for y in bottom_rights])
    min_x, max_x = extrema(x_coords)
    min_y, max_y = extrema(y_coords)
    top_left = min_x, min_y
    bottom_right = max_x, max_y

    box_dimensions = scale_box(top_left, bottom_right, (19, 19), img_size)

    return draw(load(image_path), Polygon([Point(box_dimensions[1]), Point(box_dimensions[2]), Point(box_dimensions[3]), Point(box_dimensions[4])]))
end


text:
# EMalgorithm() currently does not use Jacobians. If profiling shows that most of the time is spent in the M-step rather than teh E-step,
# this might something useful to do.

"""
Generic EM algorithm.

	em(model,data,thetai,L=1000)

Maximization of the likelihood by EM algorithm. The E-step is handled by `estep()`'s
forward-backward algorithm and the M-step is a generic numerical optimization step.
`thetai` is a starting value and `L` a maximum numer of iterations.

"""
function em(model::DynamicDiscreteModel,data,thetai=rand(dim(model)),L=1000,tol=1e-10)
	dx,dy=size(model.mu)
	w=zeros(dx,dy,dx,dy)
	llks=zeros(L)
	theta=copy(thetai)
	coef!(model,theta)

	nonzeroindicies=find(model.m.!=0.0)

	llks[1]=loglikelihood(model,data)
	l=1
	go=true
	while go
		l+=1
		w[:]=0
		estep(model,data,w)
		function ff(xtheta)
			coef!(model,xtheta)
			#w[] has zero at lest wherever model.m has zeros, and maybe more
			#we still must be careful of not calling log(0)
			#this is what nonzeroindicies is for
			res=0.0
			for i in nonzeroindicies
				res+=w[i]*log(model.m[i])
			end
			-res
		end
		ret = Optim.optimize(ff, theta,method=:cg,iterations=L)
		theta[:]=ret.minimum
		coef!(model,theta)
		llks[l]=loglikelihood(model,data)
		go=(abs(llks[l]-llks[l-1])>tol && l<L)
	end
	println()
	println(" estimating dynamic discrete model via EM algorithm...")
	println("  $l iterations, final log-likelihood: $(round(llks[l],4))")
	# plot(llks[1:l])
	theta
end


text:
using PartedArrays, Test, ForwardDiff
using BenchmarkTools
using DocStringExtensions

"Sense of a constraint (inequality / equality / null)"
abstract type ConstraintType end
"Inequality constraints"
abstract type Equality <: ConstraintType end
"Equality constraints"
abstract type Inequality <: ConstraintType end
"An empty constraint"
abstract type Null <: ConstraintType end

abstract type GeneralConstraint end
abstract type AbstractConstraint{S<:ConstraintType} <: GeneralConstraint end

import Base.+
function +(C::Vector{<:GeneralConstraint}, con::GeneralConstraint)
    C_ = append!(GeneralConstraint[], C)
    push!(C_, con)
end

function +(C::Vector{<:GeneralConstraint}, C2::Vector{<:GeneralConstraint})
    C_ = append!(GeneralConstraint[], C)
    append!(C_, C2)
end
+(con1::GeneralConstraint, con2::GeneralConstraint) = [con1] + con2

"$(SIGNATURES) Return the type of the constraint (Inequality or Equality)"
type(::AbstractConstraint{S}) where S = S

function con_methods(f::Function)
    term = hasmethod(f, (AbstractVector, AbstractVector))
    stage = hasmethod(f, (Vector, Vector, Vector))
    if term && stage
        :all
    elseif stage
        :stage
    elseif term
        :terminal
    else
        ArgumentError("The function doesn't match the expected signatures")
    end
end

"""$(TYPEDEF) General nonlinear vector-valued constraint with `p` entries at a single timestep for `n` states and `m` controls.
# Constructors
```julia
Constraint{S}(confun, ∇confun, n, m, p, label; inds, term, inputs)
Constraint{S}(confun, ∇confun, n,    p, label; inds, term, inputs)  # Terminal constraint
Constraint{S}(confun,          n, m, p, label; inds, term, inputs)
Constraint{S}(confun,          n,    p, label; inds, term, inputs)  # Terminal constraint
```

# Arguments
* `confun`: Inplace constraint function, either `confun!(v, x, u)` or `confun!(v, xN)`.
* `∇confun`: Inplace constraint jacobian function, either `∇confun!(V, x, u)` or `∇confun!(V, xN)`
* `n,m,p`: Number of state, controls, and constraints
* `label`: Symbol that unique identifies the constraint
* `inds`: Specifies the indices of `x` and `u` that are passed into `confun`. Default highly recommended.
* `term`: Where the constraint is applicable. One of `(:all, :stage, :terminal)`. Detected automatically based on defined methods of `confun` but can be manually given.
* `inputs`: Which inputs are actually used. One of `(:xu, :x, :u)`. Should be specified.

"""
struct Constraint{S} <: AbstractConstraint{S}
    c::Function
    ∇c::Function
    p::Int
    label::Symbol
    inds::Vector{Vector{Int}}
    type::Symbol
    inputs::Symbol

    function Constraint{S}(c::Function, ∇c::Function, p::Int, label::Symbol,
            inds::Vector{Vector{Int}}, type::Symbol, inputs::Symbol) where S <: ConstraintType
        if type ∉ [:all,:stage,:terminal]
            ArgumentError(string(type) * " is not a supported constraint type")
        end
        new{S}(c,∇c,p,label,inds,type,inputs)
    end
end

"$(TYPEDEF) Create a stage-wise constraint, using ForwardDiff to generate the Jacobian"
function Constraint{S}(c::Function, n::Int, m::Int, p::Int, label::Symbol;
        inds=[collect(1:n), collect(1:m)], term=con_methods(c), inputs=:xu) where S<:ConstraintType
    ∇c,c_aug = generate_jacobian(c,n,m,p)
    Constraint{S}(c, ∇c, p, label, inds, term, inputs)
end

"Create a terminal constraint using ForwardDiff"
function Constraint{S}(c::Function, n::Int, p::Int, label::Symbol;
        inds=[collect(1:n), collect(1:m)], term=:terminal, inputs=:x) where S<:ConstraintType
    m = 0
    ∇c,c_aug = generate_jacobian(c,n,m,p)
    Constraint{S}(c, ∇c, p, label, inds, term, inputs)
end

"$(TYPEDEF) Create a constraint, providing an analytical Jacobian"
function Constraint{S}(c::Function, ∇c::Function, n::Int, m::Int, p::Int, label::Symbol;
        inds=[collect(1:n), collect(1:m)], term=con_methods(c), inputs=:xu) where S<:ConstraintType
    Constraint{S}(c, ∇c, p, label, inds, term, inputs)
end

"Create a terminal constraint with analytical Jacobian"
function Constraint{S}(c::Function, ∇c::Function, n::Int, p::Int, label::Symbol;
        inds=[collect(1:n), collect(1:0)], term=:terminal, inputs=:x) where S<:ConstraintType
    Constraint{S}(c, ∇c, p, label, inds, term, inputs)
end

function Constraint()
    Constraint{Null}(x->nothing, x->nothing, 0, 0, 0, :null, term, :null)
end

evaluate!(v::AbstractVector, con::Constraint, x::AbstractVector, u::AbstractVector) = is_stage(con) ? con.c(v, x, u) : nothing
evaluate!(v::AbstractVector, con::Constraint, x::AbstractVector) = is_terminal(con) ? con.c(v,x) : nothing
jacobian!(V::AbstractMatrix, con::Constraint, x::AbstractVector, u::AbstractVector) = is_stage(con) ? con.∇c(V, x, u) : nothing
jacobian!(V::AbstractMatrix, con::Constraint, x::AbstractVector) = is_terminal(con) ? con.∇c(V,x) : nothing

is_terminal(con::Constraint) = con.type != :stage
is_stage(con::Constraint) = con.type != :terminal
function Base.length(con::Constraint, type=:stage)
    if type == :stage && is_stage(con)
        return con.p
    elseif type == :terminal && is_terminal(con)
        return con.p
    else
        return 0
    end
end

"""$(TYPEDEF) Linear bound constraint on states and controls
# Constructors
```julia
BoundConstraint(n, m; x_min, x_max, u_min, u_max)
```
Any of the bounds can be ±∞. The bound can also be specifed as a single scalar, which applies the bound to all state/controls.
"""
struct BoundConstraint{T} <: AbstractConstraint{Inequality}
    x_max::Vector{T}
    x_min::Vector{T}
    u_max::Vector{T}
    u_min::Vector{T}
    jac::SparseMatrixCSC{T,Int}
    label::Symbol
    active::NamedTuple{(:x_max, :u_max, :x_min, :u_min, :all, :x_all),NTuple{6,BitArray{1}}}
    inds::Vector{Vector{Int}}
    part::NamedTuple{(:x_max, :u_max, :x_min, :u_min),NTuple{4,UnitRange{Int64}}}
end

"""$(SIGNATURES) Create a stage bound constraint
Will default to bounds at infinity. "trim" will remove any bounds at infinity from the constraint function.
"""
function BoundConstraint(n::Int,m::Int; x_min=ones(n)*-Inf, x_max=ones(n)*Inf,
                                        u_min=ones(m)*-Inf, u_max=ones(m)*Inf, trim::Bool=true)
     # Validate bounds
     u_max, u_min = _validate_bounds(u_max,u_min,m)
     x_max, x_min = _validate_bounds(x_max,x_min,n)

     # Specify stacking
     part = create_partition((n,m,n,m),(:x_max,:u_max,:x_min,:u_min))

     # Pre-allocate jacobian
     jac_bnd = [Diagonal(1.0I,n+m); -Diagonal(1.0I,n+m)]
     jac_A = view(jac_bnd,:,part.x_max)
     jac_B = view(jac_bnd,:,part.u_max)

     inds2 = [collect(1:n), collect(1:m)]

     # Specify which controls and states get used  TODO: Allow this as an input

     if trim
         active = (x_max=isfinite.(x_max),u_max=isfinite.(u_max),
                   x_min=isfinite.(x_min),u_min=isfinite.(u_min))
         lengths = [count(val) for val in values(active)]
         part = create_partition(Tuple(lengths),keys(active))

         active_all = vcat(values(active)...)
         active = merge(active, (all=active_all, x_all=[active.x_max; falses(m); active.x_min; falses(m)]))
         ∇c_trim(C,x,u) = copyto!(C,jac_bnd[active_all,:])
     else
         active = (x_max=trues(n),u_max=trues(m),
                   x_min=trues(n),u_min=trues(m),
                   all=trues(2(n+m)), x_all=[trues(n); falses(m); trues(m); falses(m)])
     end
     return BoundConstraint(x_max, x_min, u_max, u_min, jac_bnd, :bound, active, inds2, part)
end

"""```julia
combine(::BoundConstraint, ::BoundConstraint)
```
Stack two bound constraints. Useful when the state is augmented.
"""
function combine(bnd1::BoundConstraint, bnd2::BoundConstraint)
    n1,m1 = length(bnd1.x_max), length(bnd1.u_max)
    n2,m2 = length(bnd2.x_max), length(bnd2.u_max)
    bnd = BoundConstraint(n1+n2, m1+m2,
                          x_max=[bnd1.x_max; bnd2.x_max],
                          x_min=[bnd1.x_min; bnd2.x_min],
                          u_max=[bnd1.u_max; bnd2.u_max],
                          u_min=[bnd1.u_min; bnd2.u_min])
end

"""```julia
vcat(bnd1::BoundConstraint, bnd2::BoundConstraint)
```
Stack two bound constraints. Useful when the state is augmented. Equivalent to combine(bnd1, bnd2).
"""
Base.vcat(bnd1::BoundConstraint, bnd2::BoundConstraint) = combine(bnd1, bnd2)

function evaluate!(v::AbstractVector, bnd::BoundConstraint, x::AbstractVector, u::AbstractVector)
    inds = bnd.part
    active = bnd.active
    v[inds.x_max] = (x - bnd.x_max)[active.x_max]
    v[inds.u_max] = (u - bnd.u_max)[active.u_max]
    v[inds.x_min] = (bnd.x_min - x)[active.x_min]
    v[inds.u_min] = (bnd.u_min - u)[active.u_min]
end

function evaluate!(v::AbstractVector, bnd::BoundConstraint, x::AbstractVector)
    inds = bnd.part
    active = bnd.active
    m = count(active.u_max)
    v[inds.x_max     ] = (x - bnd.x_max)[active.x_max]
    v[inds.x_min .- m] = (bnd.x_min - x)[active.x_min]
end

function jacobian!(V::AbstractMatrix, bnd::BoundConstraint, x::AbstractVector, u::AbstractVector)
    n,m = length(bnd.x_max), length(bnd.u_max)
    # TODO: avoid concatenation
    copyto!(V, bnd.jac[bnd.active.all,:])
end

function jacobian!(V::AbstractMatrix, bnd::BoundConstraint, x::AbstractVector)
    copyto!(V, bnd.jac[bnd.active.x_all,1:length(x)])
end

function Base.length(bnd::BoundConstraint, type=:stage)
    n,m = length(bnd.x_max), length(bnd.u_max)
    if type == :stage
        count(bnd.active.all)
    elseif type == :terminal
        count(bnd.active.x_max) + count(bnd.active.x_min)
    end
end

is_terminal(::BoundConstraint) = true
is_stage(::BoundConstraint) = true

"""
$(SIGNATURES)
Check max/min bounds for state and control.

Converts scalar bounds to vectors of appropriate size and checks that lengths
are equal and bounds do not result in an empty set (i.e. max > min).

# Arguments
* n: number of elements in the vector (n for states and m for controls)
"""
function _validate_bounds(max,min,n::Int)

    if min isa Real
        min = ones(n)*min
    end
    if max isa Real
        max = ones(n)*max
    end
    if length(max) != length(min)
        throw(DimensionMismatch("u_max and u_min must have equal length"))
    end
    if ~all(max .>= min)
        throw(ArgumentError("u_max must be greater than u_min"))
    end
    if length(max) != n
        throw(DimensionMismatch("limit of length $(length(max)) doesn't match expected length of $n"))
    end
    return max, min
end

"""$(SIGNATURES)
A constraint where x,y positions of the state must remain a distance r from a circle centered at x_obs
Assumes x,y are the first two dimensions of the state vector
"""
function planar_obstacle_constraint(n, m, x_obs, r_obs, label=:obstacle)
    c(v,x,u) = v[1] = circle_constraint(x, x_obs, r_obs)
    c(v,x) = circle_constraint(x, x_obs, r_obs)
    Constraint{Inequality}(c, n, m, 1, label)
end

"""```julia
goal_constraint(xf)
```
Creates a terminal equality constraint specifying the goal. All states must be specified.
"""
function goal_constraint(xf::Vector{T}) where T
    n = length(xf)
    terminal_constraint(v,xN) = copyto!(v,xN-xf)
    terminal_jacobian(C,xN) = copyto!(C,Diagonal(I,n))
    Constraint{Equality}(terminal_constraint, terminal_jacobian, n, :goal, [collect(1:n),collect(1:0)], :terminal, :x)
end

function infeasible_constraints(n::Int, m::Int)
    idx_inf = (n+m) .+ (1:n)
    u_inf = m .+ (1:n)
    ∇inf = zeros(n,2n+m)
    ∇inf[:,idx_inf] = Diagonal(1.0I,n)
    inf_con(v,x,u) = copyto!(v, u)
    inf_jac(C,x,u) = copyto!(C, ∇inf)
    Constraint{Equality}(inf_con, inf_jac, n, :infeasible, [collect(1:n), collect(u_inf)], :stage, :u)
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


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

# TODO: LegendHDF5File, LegendHDF5Input, LegendHDF5Output
# with Base.read/write, open/close, etc., and atomic file names.

const default_compression = ("shuffle", (), "deflate", 3)

const datatype_regexp = r"""^(([A-Za-z_]*)(<([0-9,]*)>)?)(\{(.*)\})?$"""
const arraydims_regexp = r"""^<([0-9,]*)>$"""

function _eldatatype_from_string(s::Union{Nothing,AbstractString}) where L_dset
    if s == nothing || s == ""
        RealQuantity
    else
        datatype_from_string(s)
    end
end

Base.@propagate_inbounds _tuple_droplast(x::NTuple{N,Any}, ::Val{M}) where {N,M} =
    Base.ntuple(i -> x[i], Val{N-M}())

_namedtuple_type(members::AbstractVector{<:AbstractString}) = NamedTuple{(Symbol.(members)...,)}

function datatype_from_string(s::AbstractString) where L_dset
    if s == "real"
        RealQuantity
    elseif s == "bool"
        Bool
    elseif s == "string"
        String
    elseif s == "symbol"
        Symbol
    elseif haskey(_datatype_dict, s)
        _datatype_dict[s]
    else
        m = match(datatype_regexp, s)
        m isa Nothing && throw(ErrorException("Invalid datatype string \"$s\""))
        tp = m[2]
        content = m[6]
        if tp == "struct"
            _namedtuple_type(split(content, ","))
        elseif tp == "table"
            TypedTables.Table{<:_namedtuple_type(split(content, ","))}
        elseif tp == "ntuple"
            T = _eldatatype_from_string(content)
            (NTuple{N,<:T} where N)
        else
            dims = parse.(Int, split(m[4], ","))
            eltp = content
            T = _eldatatype_from_string(eltp)
            if tp == "array_of_equalsized_arrays"
                length(dims) == 2 || throw(ErrorException("Invalid dims $dims for datatype \"$tp\""))
                N = dims[1]; M = dims[2]
                # T isa Array || throw(ErrorException("Datatype \"$tp\" with \"array\" nested in \"array_of_equalsized_arrays\" currently not supported"))
                AbstractArrayOfSimilarArrays{<:T,M,N}
            elseif tp == "array_of_encoded_arrays"
                length(dims) == 2 || throw(ErrorException("Invalid dims $dims for datatype \"$tp\""))
                N = dims[1]; M = dims[2]
                N != 2 || throw(ErrorException("Only one-dimensional arrays of encoded arrays are supported"))
                VectorOfEncodedArrays{<:T,dims[1]}
            elseif tp == "fixedsize_array"
                length(dims) == 1 || throw(ErrorException("Invalid dims $dims for datatype \"$tp\""))
                N = dims[1]
                N == 1  || throw(ErrorException("Datatype fixedsize_array with $dims dims currently not supported\"$tp\""))
                # T <: RealQuantity || throw(ErrorException("Element type \"$eltp\" in datatype \"$tp\" currently not supported"))
                StaticArray{Tuple{L},<:T,1} where {L}
            elseif tp == "array"
                length(dims) == 1 || throw(ErrorException("Invalid dims $dims for datatype \"$tp\""))
                N = dims[1]
                AbstractArray{<:T,N}
            else
                throw(ErrorException("Unknown datatype \"$tp\""))
            end
        end
    end
end

function _inner_datatype_to_string(::Type{T}) where T
    s = datatype_to_string(T)
    isempty(s) ? "" : "{$s}"
end

datatype_to_string(::Type{<:RealQuantity}) = "real"

datatype_to_string(::Type{Bool}) = "bool"

datatype_to_string(T::Type{<:Enum{U}}) where {U} = "enum{"*join(broadcast(x -> "$(string(x))=$(U(x))", instances.(T)), ",")*"}"

datatype_to_string(::Type{<:AbstractString}) = "string"

datatype_to_string(::Type{<:Symbol}) = "symbol"

datatype_to_string(::Type{NTuple{N,T}}) where {N,T} = "ntuple$(_inner_datatype_to_string(T))"

datatype_to_string(::Type{<:AbstractArray{T,N}}) where {T,N} =
    "array<$N>$(_inner_datatype_to_string(T))"

datatype_to_string(::Type{<:StaticArray{TPL,T,N}}) where {TPL,T<:RealQuantity,N} =
    "fixedsize_array<$N>$(_inner_datatype_to_string(T))"

datatype_to_string(::Type{<:ArrayOfSimilarArrays{T,M,N}}) where {T,M,N} =
    "array_of_equalsized_arrays<$N,$M>$(_inner_datatype_to_string(T))"

datatype_to_string(::Type{<:VectorOfEncodedArrays{T,N}}) where {T,N} =
    "array_of_encoded_arrays<1,$N>$(_inner_datatype_to_string(T))"

datatype_to_string(::Type{<:NamedTuple{K}}) where K = "struct{$(join(K,","))}"

datatype_to_string(::Type{<:TypedTables.Table{<:NamedTuple{K}}}) where K = "table{$(join(K,","))}"

function _eltype(dset::HDF5.HDF5Dataset)
    dtype = HDF5.datatype(dset)
    try
        HDF5.hdf5_to_julia_eltype(dtype)
    finally
        close(dtype)
    end
end

function _cumulative_length(A::VectorOfArrays)
    elem_ptr = ArraysOfArrays.internal_element_ptr(A)
    elem_ptr[(firstindex(elem_ptr) + 1):end] .- first(elem_ptr)
end

function _cumulative_length(A::AbstractVector{<:AbstractArray})
    # ToDo: improve implementation
    cumsum(length.(eachindex.(A)))
end

function _element_ptrs(clen::Vector{<:Integer})
    vcat([1], Int.(clen) .+ 1)
end

# _without_units(A::Array, unit::Unitful.Unitlike) = ustrip.(uconvert.(unit, A))
# _with_units(A::Array, unit::Unitful.Unitlike) = A * unit

function hasattribute(
    obj::Union{HDF5.HDF5Dataset, HDF5.HDF5.DataFile}, key::Symbol
) where {T<:Union{AbstractString,Real}}
    key_str = String(key)
    attributes = HDF5.attrs(obj)
    HDF5.exists(attributes, key_str)
end

function getattribute(
    obj::Union{HDF5.HDF5Dataset, HDF5.HDF5.DataFile}, key::Symbol, ::Type{T}
) where {T<:Union{AbstractString,Real}}
    key_str = String(key)
    attributes = HDF5.attrs(obj)
    x = read(attributes[key_str])
    x isa T ? x : convert(T, x)
end

function getattribute(
    obj::Union{HDF5.HDF5Dataset, HDF5.HDF5.DataFile}, key::Symbol, default_value::T
) where {T<:Union{AbstractString,Real}}
    if hasattribute(obj, key)
        getattribute(obj, key, T)
    else
        default_value
    end
end

function setattribute!(
    obj::Union{HDF5.HDF5Dataset, HDF5.HDF5Group}, key::Symbol,
    value::Real
)
    HDF5.attrs(obj)[String(key)] = value
    nothing
end

function setattribute!(
    obj::Union{HDF5.HDF5Dataset, HDF5.HDF5Group}, key::Symbol,
    value::AbstractString
)
    # Write variable-length string for h5py compatibility (see https://github.com/h5py/h5py/issues/585).
    s_arr = Array{String,0}(undef)
    s_arr[] = convert(String, value)
    HDF5.attrs(obj)[String(key)] = s_arr
    nothing
end

LegendDataTypes.getunits(dset::HDF5.HDF5Dataset) = units_from_string(getattribute(dset, :units, ""))

function LegendDataTypes.setunits!(dset::HDF5.HDF5Dataset, units::Unitful.Unitlike)
    ustr = units_to_string(units)
    # @debug "setunits!($(_infostring(dset)), \"$ustr\")"
    setattribute!(dset, :units, ustr)
end

default_datatype(dset::HDF5.HDF5Dataset) = AbstractArray{<:RealQuantity,length(size(dset))}
default_datatype(df::HDF5.DataFile) = NamedTuple{(Symbol.(names(df))...,)}

_size(dset::HDF5.HDF5Dataset) = size(dset)
_size(df::HDF5.DataFile) = ()

# HDF5.DataFile
# HDF5.HDF5Dataset
_infostring(x::HDF5.HDF5Group) = "group \"$(HDF5.name(x))\""
_infostring(x::HDF5.HDF5Dataset) = "dataset \"$(HDF5.name(x))\" with size $(size(x)) of $(eltype(x))"

function getdatatype(input::Union{HDF5.HDF5Dataset, HDF5.DataFile})
    dtstr = getattribute(input, :datatype, "")
    dt = isempty(dtstr) ? default_datatype(input) : datatype_from_string(dtstr)
    # @debug "getdatatype($(_infostring(input))) = $dt"
    dt
end

function setdatatype!(output::Union{HDF5.HDF5Dataset, HDF5.DataFile}, datatype::Type)
    dtstr = datatype_to_string(datatype)
    # @debug "setdatatype!($(_infostring(output)), \"$dtstr\")"
    setattribute!(output, :datatype, dtstr)
end

function getinfo(dset::HDF5.HDF5Dataset)
    (
        datatype = getattribute(dset, :datatype, ""),
        units = getunits(dset),
        cumsum_length_ds = getattribute(dset, :cumsum_length_ds, ""),
    )
end

function _getcontent_impl(dset::HDF5.HDF5Dataset, idxs::NTuple{N,Colon}, axs::NTuple{N}) where {N}
    read(dset)
end

function _getcontent_impl(dset::HDF5.HDF5Dataset, idxs::Tuple{Integer}, axs::NTuple{0}) where {N}
    idxs == (1,) || Base.throw_boundserror(dset, idxs)
    read(dset)
end

function _getcontent_impl(dset::HDF5.HDF5Dataset, idxs::NTuple{N,Any}, axs::NTuple{N}) where {N}
    canonical_idxs = Base.to_indices(dset, axs, idxs)
    isinbounds = Base.checkbounds_indices(Bool, axs, idxs)
    isinbounds || Base.throw_boundserror(dset, idxs)
    dset[canonical_idxs...]
end

function getcontent(dset::HDF5.HDF5Dataset, idxs::Tuple = axes(dset))
    _getcontent_impl(dset, idxs, axes(dset))
end

function LegendDataTypes.readdata(input::Union{HDF5.HDF5Dataset, HDF5.DataFile}, name::AbstractString)
    datatype = getdatatype(input[name])
    readdata(input, name, datatype)
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::Union{T,AbstractArray{T}},
    fulldatatype::DataType = typeof(x)
) where {T<:RealQuantity}
    units = unit(eltype(x))
    # @debug "name" name
    if units == NoUnits
        # @debug "without units"
        #output[name, "shuffle", (), "deflate", 3] = x
        output[name] = x
    else
        # @debug "with units"
        output[name] = ustrip.(x)
        setunits!(output[name], units)
    end
    if fulldatatype != Nothing
        setdatatype!(output[name], fulldatatype)
    end
    nothing
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    ::Type{<:Union{RealQuantity,AbstractArray}}
)
    dset = input[name]
    data = getcontent(dset)#::AT
    units = getunits(dset)
    if units == NoUnits
        data
    else
        data * units
    end
end

# Write Bool arrays as arrays of UInt8 for h5py compatibility (see
# https://github.com/h5py/h5py/issues/641). Issue may be fixed in recent
# versions of h5py (see https://github.com/h5py/h5py/pull/821).
function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::Union{Bool,AbstractArray{Bool}},
    fulldatatype::DataType = typeof(x)
)
    data = UInt8.(x)
    output[name] = data
    if fulldatatype != Nothing
        setdatatype!(output[name], fulldatatype)
    end
    nothing
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    ::Type{<:Union{Bool,AbstractArray{<:Bool}}}
)
    dset = input[name]
    units = getunits(dset)
    units == NoUnits || throw(ErrorExceptions("Can't interpret dataset with units as Bool values"))
    data = getcontent(dset)
    # Broadcast will return BitArray, map would return Array{UInt8}:
    (x -> x > 0).(data)
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractArray{<:Enum{T}},
    fulldatatype::DataType = typeof(x)
) where {T}
    writedata(output, name, reinterpret(T, x), fulldatatype)
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:AbstractArray{<:Enum}}
)
    ET = AT.body.parameters[1].ub
    data = readdata(input, name, AbstractArray{RealQuantity})
    ET.(data)
end

function _flatview_of_array_of_ntuple(A::AbstractArray{TPL,N}) where {L,T,N,TPL<:NTuple{L,T}}
    reshape(reinterpret(T, A), L, size(A)...)
end

function _flatview_to_array_of_ntuple(A::AbstractArray{T,N}, TPL::Type{NTuple{L,T}}) where {T,N,L}
    size_A = size(A)
    sz_out = Base.tail(size_A)
    N_out = length(sz_out)
    size_A[1] == L || throw(DimensionMismatch("Length $L of NTuple type does not match first dimension of array of size $size_A"))
    tmp = reshape(reinterpret(TPL, A), sz_out...)
    convert(Array{NTuple{L,T},N_out}, tmp)
end

_array_of_ntuple_innersize(A::AbstractArray{NTuple{L,T}}) where {L,T} = L

function _array_of_ntuple_innerconv(::Type{T}, A::AbstractArray{NTuple{L,U},N}) where {T,N,L,U}
    convert(Array{NTuple{L,T},N}, A)
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractArray{T,N},
    fulldatatype::DataType = typeof(x)
) where {L,T<:NTuple{L,RealQuantity},N}
    writedata(output, name, _flatview_of_array_of_ntuple(x), fulldatatype)
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:AbstractArray{<:NTuple}}
)
    N = length(size(input[name])) - 1
    data = readdata(input, name, AbstractArray{RealQuantity,N})
    SV = AT.var.ub
    L = size(data, 1)
    if SV isa DataType
        L_expected = SV.parameters[1].parameters[1][1]
        L_expected == L || throw(ErrorException("Trying to read array of NTuples of length $L_expected, but inner dimension of data has length $L"))
    end
    _flatview_to_array_of_ntuple(data, NTuple{L,eltype(data)})
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractArray{T,N},
    fulldatatype::DataType = typeof(x)
) where {TPL,T<:StaticArray{TPL,<:RealQuantity},N}
    writedata(output, name, flatview(x), fulldatatype)
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:AbstractArray{<:StaticVector}}
)
    N = length(size(input[name])) - 1
    data = readdata(input, name, AbstractArray{RealQuantity,N})
    SV = AT.var.ub
    L = size(data, 1)
    if SV isa DataType
        L_expected = SV.parameters[1].parameters[1][1]
        L_expected == L || throw(ErrorException("Trying to read array of static vectors of length $L_expected, but inner dimension of data has length $L"))
    end
    nestedview(data, SVector{L})
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractString,
    fulldatatype::DataType = typeof(x)
)
    output[name] = x
    if fulldatatype != Nothing
        setdatatype!(output[name], fulldatatype)
    end
    nothing
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    ::Type{<:String}
)
    dset = input[name]
    read(dset)#::AT
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::Symbol,
    fulldatatype::DataType = typeof(x)
)
    writedata(output, name, String(x), fulldatatype)
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    ::Type{<:Symbol}
)
    Symbol(readdata(input, name, String))
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractArray{<:AbstractArray{T,M},N},
    fulldatatype::DataType = typeof(x)
) where {T,M,N}
    N == 1 || throw(ArgumentError("Output of multi-dimensional arrays of arrays to HDF5 is not supported"))
    # ToDo: Support vectors of multi-dimensional arrays
    M ==1 || throw(ArgumentError("Output of vectors of multi-dimensional arrays to HDF5 is not supported"))
    writedata(output, "$(name)/flattened_data", flatview(x))
    writedata(output, "$(name)/cumulative_length", _cumulative_length(x))
    setdatatype!(output[name], fulldatatype)
    nothing
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:AbstractArray{<:AbstractArray}}
)
    data = readdata(input, "$name/flattened_data")
    clen = readdata(input, "$name/cumulative_length")
    data_vec = VectorOfVectors(data, _element_ptrs(clen))
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::VectorOfEncodedArrays{T,1,C},
    fulldatatype::DataType = typeof(x)
) where {T,N,C}
    codec = first(x.codec)
    @inbounds for c in x.codec
        c != codec && throw("Can't write VectorOfEncodedArrays that has non-uniform codec parameters")
    end

    writedata(output, "$name/encoded_data", x.encoded)
    writedata(output, "$name/decoded_size", x.size)

    dset = output[name]

    codec_name = LegendDataTypes.array_codecs[C]
    setattribute!(dset, :codec, String(codec_name))
    write_to_properties!(setattribute!, dset, codec)

    setdatatype!(dset, fulldatatype)
    nothing
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:VectorOfEncodedArrays}
)
    data_vec = readdata(input, "$name/encoded_data")
    size_vec_in = readdata(input, "$name/decoded_size")
    N = _array_of_ntuple_innersize(size_vec_in)
    size_vec = _array_of_ntuple_innerconv(Int, size_vec_in)

    SV = AT.var.ub
    if SV isa DataType
        N_expected = SV.parameters[1].parameters[1][1]
        N_expected == N || throw(ErrorException("Trying to read a vector of encoded arrays with $N_expected dimensions, but data indicates $N dimensions"))
    end

    dset = input[name]

    codec_name = Symbol(getattribute(dset, :codec, String))
    C = LegendDataTypes.array_codecs[codec_name]
    codec = read_from_properties(getattribute, dset, C)
    n = length(data_vec)
    codec_vec = StructArray(fill(codec, n)) # ToDo: Improve

    StructArray{EncodedArray{Int32,N,C,Vector{UInt8}}}(
        (
            codec_vec,
            size_vec,
            data_vec
        )
    )
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::AbstractArrayOfSimilarArrays{T,M,N},
    fulldatatype::DataType = typeof(x)
) where {T,M,N}
    writedata(output, name, flatview(x), fulldatatype)
    nothing
end

# Hack:
function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:AbstractVectorOfSimilarVectors}
)
    nestedview(readdata(input, name, AbstractArray{<:RealQuantity}))
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::NamedTuple,
    fulldatatype::DataType = typeof(x)
)
    for k in keys(x)
        y = x[k]
        writedata(output, "$name/$k", y, typeof(y))
    end
    setdatatype!(output[name], fulldatatype)
    nothing
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:NamedTuple}
)
    syms = AT.body.parameters[1]
    NamedTuple{syms}(map(k -> readdata(input, "$name/$k"), syms))
end

function LegendDataTypes.writedata(
    output::HDF5.DataFile, name::AbstractString,
    x::Any,
    fulldatatype::DataType = typeof(x)
)
    # @debug("writedata(\"$(HDF5.name(output))\", \"$name\", $fulldatatype")
    Tables.istable(x) || throw(ArgumentError("Value to write, of type $(typeof(x)),is not a table"))
    cols = Tables.columns(x)
    writedata(output, name, cols, fulldatatype)
end

function LegendDataTypes.readdata(
    input::HDF5.DataFile, name::AbstractString,
    AT::Type{<:TypedTables.Table}
)
    # Hack:
    TypedTables.Table(readdata(input, name, NamedTuple{AT.var.ub.body.parameters[1]}))
end


text:
module TestMOIWrapper

using Test
using KNITRO

const MOI = KNITRO.MOI

function runtests()
    for name in names(@__MODULE__; all=true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

const TEST_CONFIG = MOI.Test.Config(
    atol=1e-4,
    rtol=1e-4,
    optimal_status=MOI.LOCALLY_SOLVED,
    exclude=Any[
        MOI.ConstraintBasisStatus,
        MOI.DualObjectiveValue,
        MOI.ObjectiveBound,
        MOI.ListOfConstraintTypesPresent,
        MOI.ConstraintFunction,
        MOI.ObjectiveFunction,
    ],
)

const MOI_BASE_EXCLUDED = String[
    # KNITRO does not support problem's modification
    "test_modification",
    # KNITRO does not support delete
    "_delete_",
    # KNITRO returns LOCALLY_INFEASIBLE, not INFEASIBLE
    "INFEASIBLE",
    # MODEL
    "test_model_copy_to_", # TODO: No Exception thrown when we copy BadConstraintModel
    "test_model_ScalarFunctionConstantNotZero",  # RequirementUnmet: _supports(config, MOI.ConstraintFunction)
    # VARIABLE
    "test_variable_get_VariableIndex", # Knitro does not support get(::, ::MathOptInterface.VariableIndex, ::String)
    # SOLVE
    "test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_",
    # CONSTRAINTS
    "test_constraint_ZeroOne", # Wrong solution returned
    "test_constraint_VectorAffineFunction_duplicate", # Knitro does not support get(::, ::MathOptInterface.VariableIndex, ::String)
    "test_constraint_PrimalStart_DualStart_SecondOrderCone", #TODO : bug in conic interface
    # LINEAR
    "test_linear_transform", # Knitro does not support model transform
    "test_linear_Semicontinuous_integration", # Wrong return status
    "test_linear_Semiinteger_integration", # Wrong return status
    "test_linear_integration_Interval", # TODO ? Knitro 13.0 converges to 9.99999* instead of 10.0
    # QUADRATIC
    "test_quadratic_Integer_SecondOrderCone", # MOI.get(model, MOI.TerminationStatus()) == MOI.OTHER_LIMIT
    # CONIC
    "test_conic", # TODO: solve issues with conic interface
]

function test_MOI_Test_cached()
    bridged = MOI.Bridges.full_bridge_optimizer(KNITRO.Optimizer(), Float64)
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        bridged,
    )
    MOI.set(model, MOI.Silent(), true)
    excluded = copy(MOI_BASE_EXCLUDED)
    append!(excluded, [
        "test_quadratic_nonhomogeneous", # Knitro diverges on second problem solved
    ])
    return MOI.Test.runtests(model, TEST_CONFIG; exclude=excluded)
end

function test_MOI_Test_bridged()
    model = MOI.Bridges.full_bridge_optimizer(KNITRO.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)
    excluded = copy(MOI_BASE_EXCLUDED)
    append!(
        excluded,
        [
            "test_add_constrained_variables_vector", # Knitro does not support getting MOI.ConstraintSet
            "test_basic", # TODO: Need better support for names
            "test_model", # TODO: Need better support for names
            "test_objective_set_via_modify", # KNITRO does not support getting MOI.ListOfModelAttributesSet
            "test_objective_get_ObjectiveFunction_ScalarAffineFunction", # KNITRO does not support getting MOI.ObjectiveFunction
            "test_objective_ObjectiveFunction_VariableIndex", # KNITRO does not support getting MOI.ObjectiveFunctionType
            "test_quadratic_duplicate_terms", # Knitro does not support getting MOI.ObjectiveFunction / MOI.ConstraintFunction
            "test_quadratic_integration", # Knitro does not support getting ObjectiveFunction / MOI.ConstraintFunction
            "test_constraint_get_ConstraintIndex", # Knitro does not support get(::, ::MathOptInterface.VariableIndex, ::String)
        ],
    )
    MOI.Test.runtests(model, TEST_CONFIG; exclude=excluded)
    return
end

end

TestMOIWrapper.runtests()