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

using Markdown
using InteractiveUtils

# ╔═╡ 7697c424-7277-11eb-105d-150285a535c4
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add("Yao")
	Pkg.add("YaoPlots")
end

# ╔═╡ db7c91af-1621-4d79-a61c-2ff7b86db491
Pkg.add("Plots")

# ╔═╡ a49017b7-64b9-4519-9a91-457f2a303e1c
using Yao, YaoPlots

# ╔═╡ a120d2b5-e408-466e-9394-33e377f90de1
using Plots, LinearAlgebra

# ╔═╡ 427e33d9-d399-4e5d-9622-431e1fbd2b07
import Optimisers

# ╔═╡ c788534a-ee16-4972-91ce-b0b3ed2799fa
md"""
# Quick Start Guide

This guide is for people who:

- doesn't know much Julia but knows Python (or other programming languages)
- knows quantum computing basics
- want to try out programming quantum circuits quickly with Julia

the contents are organized as following

1. Julia Basics
2. Basic Concepts of Yao
3. How to write a Quantum Circuit Born Machine using Yao
"""

# ╔═╡ 01ebf798-c068-4115-9cb9-381f5d8e8b35
md"""
## Julia Basics

You can jump to next section if you already know how to use Julia. We will only be introduce some basic concepts of Julia and differences comparing to Python in this section, if you wish to learn this language more seriously please refer to the learning materials on the official website: [julialang.org/learning/](https://julialang.org/learning/). This short Julia tutorial is based on [MIT computational thinking course](https://computationalthinking.mit.edu/).

If you haven't installed Julia, please refer to the front page of Julia: [julialang.org](https://julialang.org/)

If you wish to run this notebook locally, you will need to install [Pluto](https://github.com/fonsp/Pluto.jl):

1. open your [Julia interactive session (known as REPL)](https://docs.julialang.org/en/v1/manual/getting-started/)
2. press `]` key in the REPL to use the package mode, then type the following command

```julia
pkg> add Pluto
```

you will now see Julia's package manager start downloading this package. After the download is finished, press `backspace` and run the following command.

```julia
julia> import Pluto

julia> Pluto.run()
```

now it should open a web page for you in the browser, choose the downloaded notebook file `quick-start.jl`, you are good to go!
"""

# ╔═╡ 41e660a2-0d51-4a05-aacd-a78728d77b1c
md"""
### Variables and Some Basic Types

In Julia, you can define a variable similar to how you define it in Python, e.g we can define a `x` using `=` (assignment)
"""

# ╔═╡ c4f6b254-07e4-4ed7-8870-452818cb6e67
x = 1

# ╔═╡ 3b1e55f8-1616-462f-8996-6c49b92b1f0a
md"""every variable has a type, you can check it using `typeof`"""

# ╔═╡ 8742c25c-fcea-4799-a91c-04094f5d0e0b
typeof(x)

# ╔═╡ 662685dd-96e5-4a18-ac48-9c3056926a5c
md"By default Julia displays the output of the last operation. (You can suppress the output by adding `;` (a semicolon) at the end.)"

# ╔═╡ 905f2ff2-29fe-4c65-b955-f1a6675ff965
md"""
### Functions

In Julia, you can also define short-form, one-line functions using `=` (assignment) similar to how you write things mathematically.
"""

# ╔═╡ 8569be06-1f4a-47b2-aa64-b65027a3d2b8
f(x) = 2x

# ╔═╡ ee81eb70-168b-4b98-8a89-bcb8d2a5eb38
md"Typing the function's name gives information about the function. To call it we must use parentheses:"

# ╔═╡ b6bdc0df-54a2-4c9a-a4e2-e3aee8fdcd05
f

# ╔═╡ a72f47f4-0c75-4433-8afd-468609d0bbaf
f(2)

# ╔═╡ 34f81e8d-7714-40ea-b546-2383d7813eed
md"For longer functions we use the following syntax with the `function` keyword and `end`:"

# ╔═╡ b5b567da-ec59-4a9e-a70e-a226a828eb50
function g(x, y)
	z = x + y
	return z^2
end

# ╔═╡ 04f8a8d5-1bbb-42df-a0d0-ed80aa0ce957
md"""
### Control Flows

In Julia, there are `for`, `if` and `while`, they look like the following
"""

# ╔═╡ 10cd61db-09b9-4876-a3c9-2114578da639
begin
	s = 0
	for i in 1:10
		s += 1
	end
end

# ╔═╡ a1aa126f-15f3-44a2-8acb-ebccd8f16cfd
md"here `begin...end` is used to wrap a code block as Pluto required, you don't have to write it outside Pluto, we can now check the value of `s` by typing it again"

# ╔═╡ 9c002ca1-28ac-451c-ada2-8e55b90d3e64
s

# ╔═╡ 6154997e-a121-430e-bafe-c794412d9478
md"Here, `1:10` is a **range** representing the numbers from 1 to 10:"

# ╔═╡ 2ae4588f-6b88-47b0-928a-36e6df4bd755
typeof(1:10)

# ╔═╡ 201c7fd9-5723-46ce-84fa-7cb400a02f7b
md"""
the if else statement looks like the following
"""

# ╔═╡ 766d92dd-3df7-43f5-853c-43c4194aa842
if s < 10
	# do something
elseif 10 < s < 13
	# do something
else
	# do something
end

# ╔═╡ c7a414ca-a016-4959-8e06-080907ff244c
md"""
### Matrix and Array

Julia carries its own `Array` type, if you use Python, it is similar to `numpy.array` in Python except:

1. index starts from 1
2. the multi-dimensional index is column-wise

You can also have list comprehension:
"""

# ╔═╡ 382f7095-a3fe-4e99-9a2f-1b745f41830f
[i for i in 1:10]

# ╔═╡ e19c3b51-8405-4c47-9e7d-df809d281b9b
md"it works for multi-dimensional case too:"

# ╔═╡ d18227df-56d5-4879-ac26-9f45d35ba144
[(i, j) for i in 1:10, j in 1:5]

# ╔═╡ 4537e253-4954-4d2c-bc33-18aee3a2bdf5
md"most functions follow the same convention as numpy or MATLAB, e.g you can create a random matrix using:"

# ╔═╡ 301a8694-0b01-4515-ab40-d0f8564cfc78
rand(5, 5)

# ╔═╡ 749bf936-eefc-4dd5-b495-939d82f496a0
md"""
if you have question about using a function, you can always type question mark `?` in your REPL following the function name

```julia
julia> ?rand
```
"""

# ╔═╡ 0571129d-bdbf-4ff4-b607-dacb6d6d3911
md"""
### Package Manager & Environments

Julia carries its own package manager, you can use it as a normal package:

```julia
julia> using Pkg
```

to install a pacakge, you can use

```julia
julia> Pkg.add("Yao")
```

to remove a pacakge, you can use

```julia
julia> Pkg.rm("Yao")
```

All Julia program **runs inside an environment**, it is the global environment by default. It is usually recommended to run your notebook in a local environment, so we won't hit **any version conflicts** between different packages. You can create and activate a new temperary environment and install Yao as following:
"""

# ╔═╡ 8cfe674b-a624-465f-9caa-f839ee83380e
md"""
# Yao Basics

First let's import the module by "using" it:
"""

# ╔═╡ 62c75c87-eaa5-4320-9b49-3c2f8a33397b
md"""
Yao represents quantum circuits/gates using **Yao Blocks**, they are a collection of Julia objects.

For example, you can implement a quantum Fourier transformation circuit as following:
"""

# ╔═╡ 763c233b-95cf-4a58-b3fb-ae3485a09f52
A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1))))

# ╔═╡ 22075649-4798-4a89-8a37-02d6be489f3f
B(n, k) = chain(n, j==k ? put(k=>H) : A(j, k) for j in k:n)

# ╔═╡ 250eac58-7c75-482d-8ce7-a423e281fee5
qft(n) = chain(B(n, k) for k in 1:n)

# ╔═╡ 35f0bb9a-6676-4afa-864e-14e15f56ff0d
YaoPlots.plot(qft(3))

# ╔═╡ 20d1d437-ceb3-4fb8-ab81-51a9aebe1bb3
md"""
here we use `plot` function to plot the generated quantum circuit, you can also use it to check what are the block `A` and block `B`.

The `chain` function is used to chain two blocks of same size together:
"""

# ╔═╡ cd837663-02bf-4540-9d75-1c0a456940c7
YaoPlots.plot(chain(X, Y, H))

# ╔═╡ 76d40c22-3307-4ff8-b0cd-19d8cb815c72
md"the `put` function is used to put a gate on a specific location, it thus creates a larger block"

# ╔═╡ 81a4bb41-5790-42de-bf8c-ee909cc816df
YaoPlots.plot(put(5, 2=>H))

# ╔═╡ 960629fa-1409-441d-a3c0-f5276d291c6c
md"""the control gates are defined using `control` block with another block as its input.

- the 1st argument is the number of qubits
- the 2nd argument is the controlled gate and its location
"""

# ╔═╡ a6aa360a-632a-46d7-b061-f48841042b13
YaoPlots.plot(control(5, 3, 2=>H))

# ╔═╡ 575c593a-1957-4d97-84b6-622b9a99fa42
md"the quantum blocks defined for a quantum circuit eventually form a tree-like structure, they are also printed in this way:"

# ╔═╡ da1f61c6-f17a-4944-8c23-05ce1af83ca2
qft(3)

# ╔═╡ 3a7e4835-555d-4938-bba7-a0c0339bda5e
md"now we can define registers, Yao provides a builtin register type `ArrayReg` for full amplitude simulation:"

# ╔═╡ 3123fc39-a783-4ac9-b8ca-9c60c72a8a73
@doc ArrayReg

# ╔═╡ 01afe151-a6aa-4121-8db7-66560d3eec94
md"you can also create some common states using convenient functions, such as `rand_state`, `zero_state`:"

# ╔═╡ b6cafdc0-e778-44e7-94af-01d60fef2f09
zero_state(5)

# ╔═╡ 2afa7e59-1974-42f8-8399-2fd8b7680536
r = rand_state(5)

# ╔═╡ ce9eea0e-7ae7-4e5d-a3a5-db658df5a797
md"you can measure the register using `measure`"

# ╔═╡ b1a009d1-b098-4b8a-9f77-c0a56fe9f309
measure(r; nshots=5)

# ╔═╡ eaa99d2c-c96d-4034-ba98-96585c46c15f
md"""
or a convenient function `expect` is provided to calcualte the expectation on a given operator
"""

# ╔═╡ 563acff1-3bf9-4808-b202-e0fd6b5f18c9
expect(sum(put(5, i=>X) for i in 1:5), r)

# ╔═╡ 62a8531d-2eb8-4a07-84ae-75c130b56b7b
md"here we use `kron` to construct a simple 5-site Hamiltonian ``\sum_{i=1}^5 X_i``:"

# ╔═╡ 5cffa9e8-fcdd-4500-9d48-62339a2857ed
hamiltonian = sum(put(5, i=>X) for i in 1:5)

# ╔═╡ 48e8742b-0c07-443e-845e-4b4883fc659f
md"now if we define a simple parameterized circuit using ``R_x(θ)``"

# ╔═╡ ceac3d5e-93ef-419d-80e7-af69b2e94864
circuit(θ) = chain(5, put(3=>Rx(θ)))

# ╔═╡ 663cdd78-6204-4307-98d2-7037723c6f24
md"we can use the following syntax to calculate the expectation:"

# ╔═╡ 987f381d-32b1-4b21-a80d-b633772b1a4b
expect(hamiltonian, r=>circuit(2.0))

# ╔═╡ fddc9185-6a54-49b1-888d-f4c5f227e065
md"and the gradient of θ:"

# ╔═╡ 64bbde8f-3d50-4bfe-8972-985ad083dba2
expect'(hamiltonian, r=>circuit(2.0))

# ╔═╡ 803052be-4c40-443d-b8fa-8c3398ae8094
md"""
## Implementing Quantum Circuit Born Machine

Yao is designed with variational quantum circuits in mind, and this section, we
will introduce how to use Yao for this kind of task by implementing a quantum
circuit born machine described in [Jin-Guo Liu, Lei Wang (2018)](https://arxiv.org/abs/1804.04168)

first we install one more package for plotting:
"""

# ╔═╡ 3b7ba2a8-8759-4ab1-adc1-1a77c389fb41
md"""
## Training Target

In this tutorial, we will ask the variational circuit to learn the most basic
distribution: a guassian distribution. It is defined as follows:

```math
f(x \left| \mu, \sigma^2\right) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
```

We implement it as `gaussian_pdf`:

"""

# ╔═╡ d3869fc6-4a69-4d11-8e03-0d3b6a224e75
function gaussian_pdf(x, μ::Real, σ::Real)
    pl = @. 1 / sqrt(2pi * σ^2) * exp(-(x - μ)^2 / (2 * σ^2))
    pl / sum(pl)
end

# ╔═╡ 6d5dc55e-4ada-4a90-8522-f0c270cd3f42
pg = gaussian_pdf(1:1<<6, 1<<5-0.5, 1<<4)

# ╔═╡ c122f2b6-72ef-4cc8-b14f-0d26160c7fea
md"We can plot the distribution, it looks like"

# ╔═╡ 7b42ad5d-0350-4b86-8a67-7fc041d9eff5
Plots.plot(pg)

# ╔═╡ ed4c3900-cd72-4012-9279-5e806d7325cb
md"""
## Create the Circuit
A quantum circuit born machine is composited by two different layers: **rotation layer** and **entangler layer**.

## Rotation Layer

Arbitrary rotation is built with **Rotation Gate** on **Z, X, Z** axis
with parameters.

```math
Rz(\theta) \cdot Rx(\theta) \cdot Rz(\theta)
```

Since our input will be a ``|0\dots 0\rangle`` state.
The first layer of arbitrary rotation can just
use ``Rx(\theta) \cdot Rz(\theta)`` and the last
layer of arbitrary rotation could just
use ``Rz(\theta)\cdot Rx(\theta)``

In **幺**, every Hilbert operator is a **block** type, this
includes all **quantum gates** and **quantum oracles**.
In general, operators appears in a quantum circuit
can be divided into **Composite Blocks** and **Primitive Blocks**.

We follow the low abstraction principle and
thus each block represents a certain approach
of calculation. The simplest **Composite Block**
is a **Chain Block**, which chains other blocks
(oracles) with the same number of qubits together.
It is just a simple mathematical composition of
operators with same size. e.g.

```math
\text{chain(X, Y, Z)} \iff X \cdot Y \cdot Z
```

We can construct an arbitrary rotation block by chain ``Rz``, ``Rx``, ``Rz`` together.
"""

# ╔═╡ 103d36e0-f3b8-4f1e-9a0d-aeabb83818d1
YaoPlots.plot(chain(Rz(0.0), Rx(0.0), Rz(0.0)))

# ╔═╡ 4d151c49-34a8-4eea-b6b2-edb9991d5fbc
md"""
`Rx`, `Rz` will construct new rotation gate,
which are just shorthands for `rot(X, 0.0)`, etc.

Then let's chain them up
"""

# ╔═╡ d78eeaf9-f110-41f1-82e3-b7f2ba84f098
layer(nbit::Int, x::Symbol) = layer(nbit, Val(x))

# ╔═╡ 45be3e84-db17-462d-92c5-ad0d7ae0294b
layer(nbit::Int, ::Val{:first}) = chain(nbit, put(i=>chain(Rx(0), Rz(0))) for i = 1:nbit)

# ╔═╡ 8e54478a-d207-40fb-a5d0-91136dd06f1e
md"""
We do not need to feed the first `n` parameter into `put` here.
All factory methods can be **lazy** evaluate **the first arguements**, which is the number of qubits.
It will return a lambda function that requires a single interger input.
The instance of desired block will only be constructed until all the information is filled.
When you filled all the information in somewhere of the declaration, 幺 will be able to infer the others.
We will now define the rest of rotation layers
"""

# ╔═╡ 2eb04052-86e8-4d09-9e70-1bfe77923085
layer(nbit::Int, ::Val{:last}) = chain(nbit, put(i=>chain(Rz(0), Rx(0))) for i = 1:nbit)

# ╔═╡ 86eda4a1-8f18-4412-bd61-6d64c150e82b
layer(nbit::Int, ::Val{:mid}) = chain(nbit, put(i=>chain(Rz(0), Rx(0), Rz(0))) for i = 1:nbit)

# ╔═╡ 68ad1aee-9517-4a08-8a64-6d981efd29c8
md"""
## Entangler

Another component of quantum circuit born machine are
several **CNOT** operators applied on different qubits.
"""

# ╔═╡ 714b624d-9d5c-46db-a4a9-71be7564bcdc
entangler(pairs) = chain(control(ctrl, target=>X) for (ctrl, target) in pairs)

# ╔═╡ 6c03277b-185f-4c6c-a88c-d4651b085a5a
YaoPlots.plot(entangler([1=>2, 3=>4])(5))

# ╔═╡ 777ef904-e44f-419e-96fd-b6f893dd8f8a
md"We can then define such a born machine"

# ╔═╡ e6a9ee42-719d-4181-be56-eecd9464057c
function build_circuit(n, nlayers, pairs)
    circuit = chain(n)
    push!(circuit, layer(n, :first))
    for i in 2:nlayers
        push!(circuit, cache(entangler(pairs)))
        push!(circuit, layer(n, :mid))
    end
    push!(circuit, cache(entangler(pairs)))
    push!(circuit, layer(n, :last))
    return circuit
end

# ╔═╡ 691d6352-9766-4a15-b442-04237f9e1ee0
YaoPlots.plot(build_circuit(5, 3, [1=>2, 2=>3, 3=>4, 4=>5, 5=>1]))

# ╔═╡ 25253e3b-c087-45d8-8583-3b399e43fdee
md"""
We use the method `cache` here to tag the entangler block that it
should be cached after its first run, because it is actually a
constant oracle. Let's see what will be constructed
"""

# ╔═╡ 24d3a673-c9bc-4e8e-adfd-27988b2a5b8a
md"""
## MMD Loss & Gradients

The MMD loss is describe below:

```math
\begin{aligned}
\mathcal{L} &= \left| \sum_{x} p \theta(x) \phi(x) - \sum_{x} \pi(x) \phi(x) \right|^2\\
            &= \langle K(x, y) \rangle_{x \sim p_{\theta}, y\sim p_{\theta}} - 2 \langle K(x, y) \rangle_{x\sim p_{\theta}, y\sim \pi} + \langle K(x, y) \rangle_{x\sim\pi, y\sim\pi}
\end{aligned}
```

We will use a squared exponential kernel here.
"""

# ╔═╡ cfa937f1-c861-4dc9-ad9a-7d16b60f0c04
begin
	struct RBFKernel
	    σ::Float64
	    m::Matrix{Float64}
	end

	function RBFKernel(σ::Float64, space)
	    dx2 = (space .- space').^2
	    return RBFKernel(σ, exp.(-1/2σ * dx2))
	end

	kexpect(κ::RBFKernel, x, y) = x' * κ.m * y
end

# ╔═╡ 43bc5e0f-083b-432f-ad20-0b0024fce80a
md"""
There are two different way to define the loss:

In simulation we can use the probability distribution of the state directly
"""

# ╔═╡ 47bfb83f-7b5c-4aee-9ecc-99c836143540
get_prob(qcbm) = probs(zero_state(nqubits(qcbm)) |> qcbm)

# ╔═╡ 52f516ef-07f5-4eb4-9a50-fc5c55a8378b
function loss(κ, c, target)
    p = get_prob(c) - target
    return kexpect(κ, p, p)
end

# ╔═╡ 8a684763-10e5-4206-9e64-83567ad909da
md"""
Or if you want to simulate the whole process with measurement (which is entirely
physical), you should define the loss with measurement results, for convenience
we directly use the simulated results as our loss

### Gradients

the gradient of MMD loss is

```math
\begin{aligned}
\frac{\partial \mathcal{L}}{\partial \theta^i_l} &= \langle K(x, y) \rangle_{x\sim p_{\theta^+}, y\sim p_{\theta}} - \langle K(x, y) \rangle_{x\sim p_{\theta}^-, y\sim p_{\theta}}\\
&- \langle K(x, y) \rangle _{x\sim p_{\theta^+}, y\sim\pi} + \langle K(x, y) \rangle_{x\sim p_{\theta^-}, y\sim\pi}
\end{aligned}
```

which can be implemented as
"""

# ╔═╡ 8e06a6e1-3fbf-42b5-ad48-123d89fa543b
function gradient(qcbm, κ, ptrain)
    n = nqubits(qcbm)
    prob = get_prob(qcbm)
    grad = zeros(Float64, nparameters(qcbm))

    count = 1
    for k in 1:2:length(qcbm), each_line in qcbm[k], gate in content(each_line)
        dispatch!(+, gate, π/2)
        prob_pos = probs(zero_state(n) |> qcbm)

        dispatch!(-, gate, π)
        prob_neg = probs(zero_state(n) |> qcbm)

        dispatch!(+, gate, π/2) # set back

        grad_pos = kexpect(κ, prob, prob_pos) - kexpect(κ, prob, prob_neg)
        grad_neg = kexpect(κ, ptrain, prob_pos) - kexpect(κ, ptrain, prob_neg)
        grad[count] = grad_pos - grad_neg
        count += 1
    end
    return grad
end

# ╔═╡ 43a17b19-16c3-4637-bcad-47bf0c562251
md"Now let's setup the training, we will use the ADAM optimizer from our [quantum algorithm zoo](https://github.com/QuantumBFS/QuAlgorithmZoo.jl)"

# ╔═╡ 0a80c152-1b02-44b2-b5a1-89b93e5239c6
qcbm = build_circuit(6, 10, [1=>2, 3=>4, 5=>6, 2=>3, 4=>5, 6=>1])

# ╔═╡ f7cfa8d9-7695-4ad8-930d-cf242059b400
YaoPlots.plot(dispatch!(qcbm, :random)) # initialize the parameters

# ╔═╡ 4dba452c-ef69-40ce-bd8f-c8a8b189a0a8
begin
	κ = RBFKernel(0.25, 0:2^6-1)
	opt = Optimisers.setup(Optimisers.ADAM(0.01), parameters(qcbm));
end

# ╔═╡ d01425cb-336b-4729-9ef0-7a83fea32f99
function train(qcbm, κ, opt, target)
    history = Float64[]
    for _ in 1:100
        push!(history, loss(κ, qcbm, target))
        ps = parameters(qcbm)
        Optimisers.update!(opt, ps, gradient(qcbm, κ, target))
        dispatch!(qcbm, ps)
    end
    return history
end

# ╔═╡ b03bfe4b-f71e-493e-8f7f-886f4ba8b3b6
history = train(qcbm, κ, opt, pg)

# ╔═╡ ee0f95f9-2800-4f45-b380-c19a1a32464e
md"""
we can see we "learn" this simple gausian distribution!
"""

# ╔═╡ a2d494f6-833b-44d9-80c1-0fc6a64c1e0b
trained_pg = probs(zero_state(nqubits(qcbm)) |> qcbm)

# ╔═╡ 883125f1-633d-43fc-8f07-e21b6c4b50cb
begin
	fig2 = Plots.plot(1:1<<6, trained_pg; label="trained")
	Plots.plot!(fig2, 1:1<<6, pg; label="target")
	title!("distribution")
	xlabel!("x"); ylabel!("p")
	fig2
end

# ╔═╡ d8faaef1-fbf0-432b-a9a1-33c91a963dcf
md"now let's plot the training history"

# ╔═╡ 5a5963c0-5270-469d-b3f5-b83fa1422bde
Plots.plot(history)

# ╔═╡ Cell order:
# ╟─c788534a-ee16-4972-91ce-b0b3ed2799fa
# ╟─01ebf798-c068-4115-9cb9-381f5d8e8b35
# ╟─41e660a2-0d51-4a05-aacd-a78728d77b1c
# ╠═c4f6b254-07e4-4ed7-8870-452818cb6e67
# ╟─3b1e55f8-1616-462f-8996-6c49b92b1f0a
# ╠═8742c25c-fcea-4799-a91c-04094f5d0e0b
# ╟─662685dd-96e5-4a18-ac48-9c3056926a5c
# ╟─905f2ff2-29fe-4c65-b955-f1a6675ff965
# ╠═8569be06-1f4a-47b2-aa64-b65027a3d2b8
# ╟─ee81eb70-168b-4b98-8a89-bcb8d2a5eb38
# ╠═b6bdc0df-54a2-4c9a-a4e2-e3aee8fdcd05
# ╠═a72f47f4-0c75-4433-8afd-468609d0bbaf
# ╟─34f81e8d-7714-40ea-b546-2383d7813eed
# ╠═b5b567da-ec59-4a9e-a70e-a226a828eb50
# ╟─04f8a8d5-1bbb-42df-a0d0-ed80aa0ce957
# ╠═10cd61db-09b9-4876-a3c9-2114578da639
# ╟─a1aa126f-15f3-44a2-8acb-ebccd8f16cfd
# ╠═9c002ca1-28ac-451c-ada2-8e55b90d3e64
# ╟─6154997e-a121-430e-bafe-c794412d9478
# ╠═2ae4588f-6b88-47b0-928a-36e6df4bd755
# ╟─201c7fd9-5723-46ce-84fa-7cb400a02f7b
# ╠═766d92dd-3df7-43f5-853c-43c4194aa842
# ╟─c7a414ca-a016-4959-8e06-080907ff244c
# ╠═382f7095-a3fe-4e99-9a2f-1b745f41830f
# ╟─e19c3b51-8405-4c47-9e7d-df809d281b9b
# ╠═d18227df-56d5-4879-ac26-9f45d35ba144
# ╟─4537e253-4954-4d2c-bc33-18aee3a2bdf5
# ╠═301a8694-0b01-4515-ab40-d0f8564cfc78
# ╟─749bf936-eefc-4dd5-b495-939d82f496a0
# ╟─0571129d-bdbf-4ff4-b607-dacb6d6d3911
# ╠═7697c424-7277-11eb-105d-150285a535c4
# ╟─8cfe674b-a624-465f-9caa-f839ee83380e
# ╠═a49017b7-64b9-4519-9a91-457f2a303e1c
# ╟─62c75c87-eaa5-4320-9b49-3c2f8a33397b
# ╠═763c233b-95cf-4a58-b3fb-ae3485a09f52
# ╠═22075649-4798-4a89-8a37-02d6be489f3f
# ╠═250eac58-7c75-482d-8ce7-a423e281fee5
# ╠═35f0bb9a-6676-4afa-864e-14e15f56ff0d
# ╟─20d1d437-ceb3-4fb8-ab81-51a9aebe1bb3
# ╠═cd837663-02bf-4540-9d75-1c0a456940c7
# ╟─76d40c22-3307-4ff8-b0cd-19d8cb815c72
# ╠═81a4bb41-5790-42de-bf8c-ee909cc816df
# ╟─960629fa-1409-441d-a3c0-f5276d291c6c
# ╠═a6aa360a-632a-46d7-b061-f48841042b13
# ╟─575c593a-1957-4d97-84b6-622b9a99fa42
# ╠═da1f61c6-f17a-4944-8c23-05ce1af83ca2
# ╟─3a7e4835-555d-4938-bba7-a0c0339bda5e
# ╠═3123fc39-a783-4ac9-b8ca-9c60c72a8a73
# ╟─01afe151-a6aa-4121-8db7-66560d3eec94
# ╠═b6cafdc0-e778-44e7-94af-01d60fef2f09
# ╠═2afa7e59-1974-42f8-8399-2fd8b7680536
# ╟─ce9eea0e-7ae7-4e5d-a3a5-db658df5a797
# ╠═b1a009d1-b098-4b8a-9f77-c0a56fe9f309
# ╟─eaa99d2c-c96d-4034-ba98-96585c46c15f
# ╠═563acff1-3bf9-4808-b202-e0fd6b5f18c9
# ╟─62a8531d-2eb8-4a07-84ae-75c130b56b7b
# ╠═5cffa9e8-fcdd-4500-9d48-62339a2857ed
# ╟─48e8742b-0c07-443e-845e-4b4883fc659f
# ╠═ceac3d5e-93ef-419d-80e7-af69b2e94864
# ╟─663cdd78-6204-4307-98d2-7037723c6f24
# ╠═987f381d-32b1-4b21-a80d-b633772b1a4b
# ╟─fddc9185-6a54-49b1-888d-f4c5f227e065
# ╠═64bbde8f-3d50-4bfe-8972-985ad083dba2
# ╟─803052be-4c40-443d-b8fa-8c3398ae8094
# ╠═db7c91af-1621-4d79-a61c-2ff7b86db491
# ╠═a120d2b5-e408-466e-9394-33e377f90de1
# ╟─3b7ba2a8-8759-4ab1-adc1-1a77c389fb41
# ╠═d3869fc6-4a69-4d11-8e03-0d3b6a224e75
# ╠═6d5dc55e-4ada-4a90-8522-f0c270cd3f42
# ╟─c122f2b6-72ef-4cc8-b14f-0d26160c7fea
# ╠═7b42ad5d-0350-4b86-8a67-7fc041d9eff5
# ╟─ed4c3900-cd72-4012-9279-5e806d7325cb
# ╠═103d36e0-f3b8-4f1e-9a0d-aeabb83818d1
# ╟─4d151c49-34a8-4eea-b6b2-edb9991d5fbc
# ╠═d78eeaf9-f110-41f1-82e3-b7f2ba84f098
# ╠═45be3e84-db17-462d-92c5-ad0d7ae0294b
# ╟─8e54478a-d207-40fb-a5d0-91136dd06f1e
# ╠═2eb04052-86e8-4d09-9e70-1bfe77923085
# ╠═86eda4a1-8f18-4412-bd61-6d64c150e82b
# ╟─68ad1aee-9517-4a08-8a64-6d981efd29c8
# ╠═714b624d-9d5c-46db-a4a9-71be7564bcdc
# ╠═6c03277b-185f-4c6c-a88c-d4651b085a5a
# ╟─777ef904-e44f-419e-96fd-b6f893dd8f8a
# ╠═e6a9ee42-719d-4181-be56-eecd9464057c
# ╠═691d6352-9766-4a15-b442-04237f9e1ee0
# ╟─25253e3b-c087-45d8-8583-3b399e43fdee
# ╟─24d3a673-c9bc-4e8e-adfd-27988b2a5b8a
# ╠═cfa937f1-c861-4dc9-ad9a-7d16b60f0c04
# ╟─43bc5e0f-083b-432f-ad20-0b0024fce80a
# ╠═47bfb83f-7b5c-4aee-9ecc-99c836143540
# ╠═52f516ef-07f5-4eb4-9a50-fc5c55a8378b
# ╟─8a684763-10e5-4206-9e64-83567ad909da
# ╠═8e06a6e1-3fbf-42b5-ad48-123d89fa543b
# ╟─43a17b19-16c3-4637-bcad-47bf0c562251
# ╠═427e33d9-d399-4e5d-9622-431e1fbd2b07
# ╠═0a80c152-1b02-44b2-b5a1-89b93e5239c6
# ╠═f7cfa8d9-7695-4ad8-930d-cf242059b400
# ╠═4dba452c-ef69-40ce-bd8f-c8a8b189a0a8
# ╠═d01425cb-336b-4729-9ef0-7a83fea32f99
# ╠═b03bfe4b-f71e-493e-8f7f-886f4ba8b3b6
# ╠═ee0f95f9-2800-4f45-b380-c19a1a32464e
# ╠═a2d494f6-833b-44d9-80c1-0fc6a64c1e0b
# ╠═883125f1-633d-43fc-8f07-e21b6c4b50cb
# ╠═d8faaef1-fbf0-432b-a9a1-33c91a963dcf
# ╠═5a5963c0-5270-469d-b3f5-b83fa1422bde


text:
# ### Function Contour plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

## ref https://jp.mathworks.com/matlabcentral/fileexchange/35283-matlab-plot-gallery-function-contour-plot?focused=7851904&tab=example
f(x,y) = sin(3x) * cos(x+y)

x = 0:0.01:5
y = 0:0.01:5
contour(x, y, f, xlabel="x", ylabel="y", fill=true)

ext: jl


text:
using LiteQTL # always your module
using Test # and of course the test package.

@testset "LiteQTL.jl" begin
    # Write your own tests here.
    # include("download_data.jl")
    # run(`Rscript --vanilla clean_data.R`)
    # run(`Rscript --vanilla rqtl2scan.R`)
    include("scan_test.jl")
end

# # @test "testing data io" begin include("data_io_test.jl") end
# @test "testing against R/qtl2" begin  end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libLAS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libLAS")
JLLWrappers.@generate_main_file("libLAS", UUID("84dbc9eb-108a-5969-a599-69afb4e6b6ac"))
end  # module libLAS_jll


text:
mutable struct Model
    snpsOld
    snpsfile
    phenotypesprefix
    phenotypesfilename
    phenotypes
    m
    maxbp
    nlayer1
    nlayer2
    losoInterval
    layerlimit
    inlayer1
    setsize
    notcomplete
    phenoOrder
    phenoOrder2
    genoOrder
    snps
    n
    k
    outputfilename
    printchain
    submodelcount
    futures
    weightsfilename
    weights
    weightOrder
    maxlayerlimit
    type
end

mutable struct SubModel
    m
    n
    k
    snpset
    phenotype
    phenotypesfilename
    K
    rng
    yty
    output_stats
    output_filename
    printchain
    id
    weights
    type
end

mutable struct BinaryMarkers{M<:Array{String,2}, I<:Array{String,1}, B<:BEDFile}
    markerInfo::M
    IID::I
    binarygenotypes::B
end

mutable struct Phenotypes{I<:Array{String,1}, T}
    IID::I
    trait::T # n x k
end

mutable struct Weights{I<:Array{String,1}, W<:Array{Float64,1}}
    IID::I
    rinverse::W
end

mutable struct BinMarkers{M<:Array{String,2}, I<:Array{String,1}, B<:Array{UInt8,2}}
    markerInfo::M
    IID::I
    binarygenotypes::B
end

mutable struct Markers{M<:Array{String,2}, I<:Array{String,1}, G<:Array{Float64,2}}
    markerInfo::M # m x 5 columns
    IID::I # n x 1
    genotypes::G # n x m
end

mutable struct Alphas
    markerNames
    alphas
    m
    k
end

mutable struct YCreator
    phenotypes
    alphas
    genselmarkers
    n
    m
    k
    ynew
    nlayer1
    nlayer2
    alphas_adj
    outputprefix
    interval
    chr
    chrlength
    region
    code
end

mutable struct GenselMarkers
    IID
    SNPID
    markers
end

mutable struct id_dataset
        iid
        header
        data
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
using Documenter, Weave

makedocs(modules = Weave)


text:
global const cwin = Ref{Int}(0)
global const windows = Dict{Int,af_window}()

# export figure, plot, image

function figure(n)
    if !haskey(windows, n)
        windows[n] = create_window(640, 480, "figure $n")
    else
        show(windows[n])
        if is_window_closed(windows[n])
            windows[n] = create_window(640, 480, "figure $n")
        end
    end
    show(windows[n])
    cwin[] = n
end

function plot(x, y)
    if cwin[] == 0
        figure(1)
    end
    draw_plot(windows[cwin[]], x, y, Ref(af_cell(0, 0, pointer(""), AF_COLORMAP_DEFAULT)))
    show(windows[cwin[]])
end

function image(img)
    if cwin[] == 0
        figure(1)
    end
    draw_image(windows[cwin[]], img, Ref(af_cell(0, 0, pointer(""), AF_COLORMAP_SPECTRUM)))
    show(windows[cwin[]])
end

function grid(wind::af_window,rows::Integer,cols::Integer)
    _error(ccall((:af_grid,af_lib),af_err,(af_window,Cint,Cint),wind,Cint(rows),Cint(cols)))
end


text:
module ReSeis

using Distributions
using DataFrames
using FFTW
using RecipesBase
using Distributed

import StatsBase: sample

# patterns
include("patterns.jl")
include("patternset.jl")

# simulators
include("simulators.jl")

# solvers
include("solvers.jl")

# problems
include("problems.jl")

# seismic utilities
include("seis_utils.jl")

# plot recipes
include("plotrecipes/empirical_pattern.jl")

export
    # pattern set
    patternset,
    sample,

    # empirical pattern
    reals,
    realcount,
    labels,
    stack,

    # seismic related
    Kennet,
    ricker,
    d2t
end


text:
function find_state_update_events(event_vector_dict, augmented_positive_times, time_ranges, start_value, end_value, global_var_dict, global_var_id, global_var_value)
  co_occurring_events = find_state_update_events_false_positives(event_vector_dict, augmented_positive_times, time_ranges, start_value, end_value, global_var_dict, global_var_id, global_var_value)
  matches = map(x -> (x[1], x[3]), filter(tuple -> tuple[2] == 0, co_occurring_events))
  sort(matches, by=x -> length(x))
end

function find_state_update_events_false_positives(orig_event_vector_dict, augmented_positive_times, time_ranges, start_value, end_value, global_var_dict, global_var_id, global_var_value, object_specific=false)
  @show global_var_dict

  if object_specific
    event_vector_dict = orig_event_vector_dict
  else
    event_vector_dict = Dict()

    for event in filter(e -> (orig_event_vector_dict[e] isa Array) && !occursin("globalVar", e), collect(keys(orig_event_vector_dict)))
      event_vector_dict[event] = orig_event_vector_dict[event]
    end

    for event in filter(e -> !(orig_event_vector_dict[e] isa Array) && !occursin("globalVar", e) && occursin("id) x", e), collect(keys(orig_event_vector_dict)))
      for i in keys(orig_event_vector_dict[event])
        event_vector_dict[replace(event, "id) x" => "id) $(i)")] = orig_event_vector_dict[event][i]
      end
    end
  end

  co_occurring_events = []
  @show collect(keys(event_vector_dict))
  for event in sort(collect(keys(event_vector_dict)), by=length)
    event_vector = event_vector_dict[event]
    event_times = findall(x -> x == 1, event_vector)
    if time_ranges == []
      if event_times == [] || augmented_positive_times == [] || (end_value == augmented_positive_times[1][2])
        push!(co_occurring_events, (event, 0, [], []))
        push!(co_occurring_events, ("(& $(event) (== globalVar$(global_var_id) $(start_value)))", 0, [], []))
      elseif start_value != augmented_positive_times[1][2]
        push!(co_occurring_events, ("(& $(event) (== globalVar$(global_var_id) $(start_value)))", 0, [], []))
      end
    else
      if !(0 in map(time_range -> length(findall(time -> ((time >= time_range[1]) && (time <= time_range[2])), event_times)) > 0, time_ranges))
        @show event
        # co-occurring event found
        # compute number of false positives
        co_occurring_times = [event_times[i] for i in vcat(map(time_range -> findall(time -> ((time >= time_range[1]) && (time <= time_range[2])), event_times), time_ranges)...)]
        false_positive_times = [time for time in event_times if !(time in co_occurring_times)]

        # find closest event in the future; if value matches end_value, not a false_positive event!
        desired_end_values = []
        desired_start_values = []
        for time in false_positive_times
          # handle end value

          if ((time == length(global_var_dict[global_var_id])) && (global_var_dict[global_var_id][time] < global_var_value)) || ((time < length(global_var_dict[global_var_id])) && (global_var_dict[global_var_id][time + 1] < global_var_value))
            desired_end_value = ((time + 1) > length(global_var_dict[global_var_id])) ? global_var_dict[global_var_id][time] : global_var_dict[global_var_id][time + 1]
            push!(desired_end_values, desired_end_value)
          elseif (time >= maximum(map(tuple -> tuple[1], augmented_positive_times)))
            desired_end_value = -1 # augmented_positive_times[end][2]
            push!(desired_end_values, desired_end_value)
          else
            future_augmented_positive_times = filter(tuple -> tuple[1] > time, augmented_positive_times)
            closest_future_augmented_time = future_augmented_positive_times[1]
            push!(desired_end_values, closest_future_augmented_time[2])
          end

          # handle start value
          if (time < minimum(map(tuple -> tuple[1], augmented_positive_times))) || global_var_dict[global_var_id][time] < global_var_value
            desired_start_value = global_var_dict[global_var_id][time]
            push!(desired_start_values, desired_start_value)
          elseif (time < minimum(map(tuple -> tuple[1], augmented_positive_times)))
            desired_start_value = augmented_positive_times[1][2] # -1
            push!(desired_start_values, desired_start_value)
          else
            earlier_augmented_positive_times = filter(tuple -> tuple[1] <= time, augmented_positive_times)
            closest_earlier_augmented_time = earlier_augmented_positive_times[end]
            push!(desired_start_values, closest_earlier_augmented_time[2])
          end

        end
        num_false_positives_with_effects = count(x -> x != end_value && x != -1, desired_end_values)
        false_positive_with_effects_times = [false_positive_times[i] for i in findall(v -> v != end_value && v != -1, desired_end_values)]
        no_effect_times = [false_positive_times[i] for i in findall(v -> v == end_value || v == -1, desired_end_values)]
        push!(co_occurring_events, (event,
                                    num_false_positives_with_effects,
                                    [co_occurring_times..., no_effect_times...],
                                    false_positive_with_effects_times))
        # (tuple structure: (event, # of false positives, list of correct event times, list of false positive event times ))

        zipped_values = [(desired_start_values[i], desired_end_values[i]) for i in 1:length(desired_end_values)]
        num_false_positives_with_effects_state = count(x -> !(x[2] in [end_value, -1]) && (x[1] == start_value), zipped_values)
        false_positives_with_effects_state_times = [false_positive_times[i] for i in findall(x -> !(x[2] in [end_value, -1]) && (x[1] == start_value), zipped_values)]
        push!(co_occurring_events, ("(& $(event) (== (prev globalVar$(global_var_id)) $(start_value)))",
                                    num_false_positives_with_effects_state,
                                    [co_occurring_times..., [false_positive_times[i] for i in findall(x -> (x[2] in [end_value, -1]) && (x[1] == start_value), zipped_values)]...],
                                    false_positives_with_effects_state_times))
      end
    end
  end
  @show co_occurring_events
  co_occurring_events = sort(co_occurring_events, by=x->x[2])

  # among events with minimum # of false positives, sort by length of event (i.e. so "left" appears before "& left (== globalVar1 1)")
  if co_occurring_events == []
    []
  else
    min_false_positives = co_occurring_events[1][2]
    min_false_positive_events = sort(filter(e -> e[2] == min_false_positives, co_occurring_events), by=x->length(x[1]))
    other_events = filter(e -> e[2] != min_false_positives, co_occurring_events)
    co_occurring_events = vcat(min_false_positive_events, other_events)
  end
end

function group_ranges(ranges)
  dict = Dict()

  for range in ranges
    start_time = range[1][1]
    start_value = range[1][2]

    end_time = range[2][1] - 1
    end_value = range[2][2]

    if !((start_value, end_value) in keys(dict))
      dict[(start_value, end_value)] = [range]
    else
      push!(dict[(start_value, end_value)], range)
    end
  end

  grouped_ranges = sort(map(k -> dict[k], collect(keys(dict))), by=group->group[1][1][1]) # sort by first start time in each group
  grouped_ranges
end

function group_ranges(ranges_dict::Dict)
  ranges = vcat(filter(x -> length(x) > 0, vcat(map(id -> map(range -> ((range[1][1], range[1][2], id), (range[2][1], range[2][2], id)), ranges_dict[id]), collect(keys(ranges_dict)))))...)
  group_ranges(ranges)
end

function find_state_update_events_object_specific(event_vector_dict, augmented_positive_times_dict, grouped_range, object_ids, object_mapping, max_state_value)
  # extract object-specific events
  events = filter(e -> !(event_vector_dict[e] isa AbstractArray) && sort(collect(keys(event_vector_dict[e]))) == sort(object_ids), collect(keys(event_vector_dict)))
  global_events = filter(e -> event_vector_dict[e] isa Array, collect(keys(event_vector_dict)))

  update_events_dict = Dict()
  for object_id in object_ids
    object_event_vector_dict = Dict()
    for event in events
      object_event_vector_dict[event] = event_vector_dict[event][object_id]
    end

    augmented_positive_times = augmented_positive_times_dict[object_id]

    time_ranges = map(range -> (range[1][1], range[2][1]), filter(r -> r[1][3] == object_id, grouped_range))

    start_value = grouped_range[1][1][2]
    end_value = grouped_range[1][2][2]

    if "field1" in map(x -> x[1], filter(obj -> !isnothing(obj), object_mapping[object_id])[1].type.custom_fields)
      custom_field_index = findall(field_tuple -> field_tuple[1] == "field1", filter(obj -> !isnothing(obj), object_mapping[object_ids[1]])[1].type.custom_fields)[1]
      object_var_dict = Dict(1 => map(obj -> isnothing(obj) ? -1 : obj.custom_field_values[custom_field_index], object_mapping[object_id])[1:end-1])
    else
      object_var_dict = Dict(1 => ones(Int, length(event_vector_dict[global_events[1]])))
    end
    object_var_value = max_state_value

    # @show object_id
    # @show object_event_vector_dict
    # @show augmented_positive_times
    # @show time_ranges
    # @show start_value
    # @show end_value
    update_events = find_state_update_events(object_event_vector_dict, augmented_positive_times, time_ranges, start_value, end_value, object_var_dict, 1, object_var_value)

    update_events_dict[object_id] = update_events
  end
  @show update_events_dict
  common_events = intersect(map(id -> map(x -> x[1], update_events_dict[id]), object_ids)...)
  if common_events == []
    # FAILURE CASE
    []
  else
    common_event = common_events[1]
    event_times = []
    for object_id in object_ids
      object_event_times = filter(event_tuple -> event_tuple[1] == common_event, update_events_dict[object_id])[1][2]
      augmented_times = map(time -> (time, object_id), object_event_times)
      event_times = vcat(event_times..., augmented_times...)
    end
    [(common_event, event_times)]
  end

end

function is_co_occurring(event, event_vector, update_function_times)
  event_times = findall(x -> x == 1, event_vector)
  if repr(sort(intersect(event_times, update_function_times))) == repr(sort(update_function_times))
    true
  else
    false
  end
end

function num_false_positives(event_vector, update_function_times)
  event_times = findall(x -> x == 1, event_vector)
  length([ time for time in event_times if !(time in update_function_times) ])
end


text:
module Bitcoin

using ECC, Base58
using SHA: sha1, sha256
using Ripemd: ripemd160
import Base.show
export Tx, TxIn, TxOut, Script
export address, wif, txparse, txserialize, txid, txfee, txsighash,
       scriptevaluate, txfetch, txverify, txsigninput

const SIGHASH_ALL = 1
const SIGHASH_NONE = 2
const SIGHASH_SINGLE = 3

include("helper.jl")
include("address.jl")
include("op.jl")
include("script.jl")
include("tx.jl")

end # module


text:
module ReST

end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # Julia is fast
#
# Very often, benchmarks are used to compare languages.  These benchmarks can
# lead to long discussions, first as to exactly what is being benchmarked and
# secondly what explains the differences.  These simple questions can sometimes
# get more complicated than you at first might imagine.
#
# The purpose of this notebook is for you to see a simple benchmark for
# yourself.
#
# (This material began life as a wonderful lecture by Steven Johnson at MIT:
# [Boxes and registers](https://github.com/stevengj/18S096/blob/master/lectures/lecture1/Boxes-and-registers.ipynb)).

#-

# # Outline of this notebook
#
# - Define the sum function
# - Implementations & benchmarking of sum in...
#     - Julia (built-in)
#     - Julia (hand-written)
#     - C (hand-written)
#     - python (built-in)
#     - python (numpy)
#     - python (hand-written)
# - Towards exploiting parallelism with Julia
#     - Allowing for floating point associativity
#     - Making use of four cores at once: built-in
#     - Making use of four cores at once: hand-written
# - Summary of benchmarks

#-

# # `sum`: An easy enough function to understand

#-

# Consider the  **sum** function `sum(a)`, which computes
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i,
# $$
# where $n$ is the length of `a`.

a = rand(10^7) # 1D vector of random numbers, uniform on [0,1)

#-

sum(a)

# The expected result is ~0.5 * 10^7, since the mean of each entry is 0.5

#-

# # Benchmarking a few ways in a few languages

@time sum(a)

#-

@time sum(a)

#-

@time sum(a)

# The `@time` macro can yield noisy results, so it's not our best choice for benchmarking!
#
# Luckily, Julia has a `BenchmarkTools.jl` package to make benchmarking easy and accurate:

#-

using BenchmarkTools

#-

@benchmark sum($a)

# # 1. Julia Built-in
#
# So that's the performance of Julia's built-in sum — but that could be doing any number of tricks to be fast, including not using Julia at all in the first place! Of course, it is indeed written in Julia, but would it perform if we write a naive implementation ourselves?

@which sum(a)

# Let's save these benchmark results to a dictionary so we can start keeping track of them and comparing them down the line.

j_bench = @benchmark sum($a)

#-

d = Dict()
d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# # 2. Julia (hand-written)

function mysum(A)
    s = 0.0
    for a in A
        s += a
    end
    return s
end

#-

j_bench_hand = @benchmark mysum($a)

#-

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# So that's about 2x slower than the builtin definition. We'll see why later on.
#
# But first: is this fast?  How would we know?  Let's compare it to some other languages...

#-

# #  3. The C language
#
# C is often considered the gold standard: difficult on the human, nice for the machine. Getting within a factor of 2 of C is often satisfying. Nonetheless, even within C, there are many kinds of optimizations possible that a naive C writer may or may not get the advantage of.
#
# The current author does not speak C, so he does not read the cell below, but is happy to know that you can put C code in a Julia session, compile it, and run it. Note that the `"""` wrap a multi-line string.

using Libdl
C_code = """
    #include <stddef.h>
    double c_sum(size_t n, double *X) {
        double s = 0.0;
        for (size_t i = 0; i < n; ++i) {
            s += X[i];
        }
        return s;
    }
"""

const Clib = tempname()   # make a temporary file

## compile to a shared library by piping C_code to gcc
## (works only if you have gcc installed):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

## define a Julia function that calls the C function:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

#-

c_sum(a)

#-

c_sum(a) ≈ sum(a) # type \approx and then <TAB> to get the ≈ symbolb

# We can now benchmark the C code directly from Julia:

c_bench = @benchmark c_sum($a)

#-

d["C"] = minimum(c_bench.times) / 1e6  # in milliseconds
d

# # 4. Python's built in `sum`

#-

# The `PyCall` package provides a Julia interface to Python:

using PyCall

#-

## get the Python built-in "sum" function:
pysum = pybuiltin("sum")

#-

pysum(a)

#-

pysum(a) ≈ sum(a)

#-

py_list_bench = @benchmark $pysum($a)

#-

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# # 5. Python: `numpy`
#
# `numpy` is an optimized C library, callable from Python.
# It may be installed within Julia as follows:

using Conda

#-

## Conda.add("numpy")

#-

numpy_sum = pyimport("numpy")["sum"]

py_numpy_bench = @benchmark $numpy_sum($a)

#-

numpy_sum(a)

#-

numpy_sum(a) ≈ sum(a)

#-

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# # 6. Python, hand-written

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

#-

py_hand = @benchmark $sum_py($a)

#-

sum_py(a)

#-

sum_py(a) ≈ sum(a)

#-

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# # Summary so far

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value; digits=1), 6, "."))
end

# We seem to have three different performance classes here: The numpy and Julia
# builtins are leading the pack, followed by the hand-written Julia and C
# definitions. Those seem to be about 2x slower.  And then we have the much much
# slower Python definitions over 100x slower.

#-

# # Exploiting parallelism with Julia

#-

# The fact that our hand-written Julia solution was almost an even multiple of
# 2x slower than the builtin solutions is a big clue: perhaps theres some sort
# of 2x parallelism going on here?
#
# (In fairness, there are ways to exploit parallelism in other languages, too,
# but for brevity we won't cover them)

#-

# # 7. Julia (allowing floating point associativity)

#-

# The `for` loop
#
# ```julia
# for a in A
#     s += a
# end
# ```
#
# defines a very strict _order_ to the summation: Julia follows exactly what you
# wrote and adds the elements of `A` to the result `s` in the order it iterates.
# Since floating point numbers aren't associative, a rearrangement here would
# change the answer — and Julia is loathe to give you different answer than
# the one you asked for.
#
# You can, however, tell Julia to relax that rule and allow for associativity
# with the `@fastmath` macro. This might allow Julia to rearrange the sum in an
# advantageous manner.

function mysum_fast(A)
    s = 0.0
    for a in A
        @fastmath s += a
    end
    s
end

#-

j_bench_hand_fast = @benchmark mysum_fast($a)

#-

mysum_fast(a)

#-

d["Julia hand-written fast"] = minimum(j_bench_hand_fast.times) / 1e6
d

# # 8. Distributed Julia (built-in)

#-

# We can take this one step further: nearly every modern computer these days has
# multiple cores. All the above solutions are working one core hard, but all the
# others are just sitting by idly. Let's put them to work!

using Distributed
using DistributedArrays
addprocs(4)
@sync @everywhere workers() include("/opt/julia-1.0/etc/julia/startup.jl") # Needed just for JuliaBox
@everywhere using DistributedArrays

#-

adist = distribute(a)
j_bench_dist = @benchmark sum($adist)

#-

d["Julia 4x built-in"] = minimum(j_bench_dist.times) / 1e6
d

# # 9. Distributed Julia (hand-written)
#
# Ok, that might be cheating, too — it's again just calling a library
# function. Is it possible to write distributed sum ourselves?

function mysum_dist(a::DArray)
    r = Array{Future}(undef, length(procs(a)))
    for (i, id) in enumerate(procs(a))
        r[i] = @spawnat id sum(localpart(a))
    end
    return sum(fetch.(r))
end

#-

j_bench_hand_dist = @benchmark mysum_dist($adist)

#-

d["Julia 4x hand-written"] = minimum(j_bench_hand_dist.times) / 1e6
d

# # Overall Summary

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value; digits=1), 6, "."))
end

# # Key take-aways
#
# * Julia allows for serial C-like performance, even with hand-written functions
# * Julia allows us to exploit many forms of parallelism to further improve performance. We demonstrated:
#     * Single-processor parallelism with SIMD
#     * Multi-process parallelism with DistributedArrays
# * But there are many other ways to express parallelism, too!

max_stars_repo_path: Courses/Parallel_Computing/015 Julia is fast.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
import Base: ==
using FerriteMeshParser
using Ferrite
using Test

# Define equality check between grids
==(a::Grid,b::Grid) = all([isequal(getfield.((a,b), (name,))...) for name in fieldnames(typeof(a))])

# Overload default_interpolation where not given in Ferrite
if !isdefined(Main, :SerendipityQuadrilateral)
    const SerendipityQuadrilateral = Cell{2,8,4}
end
Ferrite.default_interpolation(::Type{SerendipityQuadrilateral}) = Serendipity{2, RefCube, 2}()
Ferrite.vertices(c::SerendipityQuadrilateral) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
Ferrite.faces(c::SerendipityQuadrilateral) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))

@testset "CheckVolumes" begin
    unit_volume_files = ("2D_UnitArea_Linear", "2D_UnitArea_Quadratic",
                         "3D_UnitVolume_LinearHexahedron", "3D_UnitVolume_LinearTetrahedron",
                         "3D_UnitVolume_QuadraticHexahedron", "3D_UnitVolume_QuadraticTetrahedron")
    for base_name in unit_volume_files
        filename = joinpath(@__DIR__, "test_files", base_name*".inp")
        grid = get_ferrite_grid(filename)
        dh = MixedDofHandler(grid)
        unique_celltypes = unique(typeof.(grid.cells))

        fields = [Field(:u, Ferrite.default_interpolation(type), 1) for type in unique_celltypes]
        cellsets = [findall(x->isa(x,type), grid.cells) for type in unique_celltypes]
        fieldhandlers = [FieldHandler([field], Set(set)) for (field,set) in zip(fields, cellsets)]
        push!.((dh,), fieldhandlers)
        close!(dh)

        cv_vec = Any[]
        for type in unique_celltypes
            ip = Ferrite.default_interpolation(type)
            dim = Ferrite.getdim(ip)
            ref = Ferrite.getrefshape(ip)
            qr = QuadratureRule{dim, ref}(1)
            push!(cv_vec, CellScalarValues(qr, ip))
        end

        for (cellset, cv, type) in zip(cellsets, cv_vec, unique_celltypes)
            for cell in CellIterator(dh, cellset)
                reinit!(cv, cell)
                V = getdetJdV(cv, 1)
                volcheck = V ≈ 1.0
                !volcheck && println("Volume check failure for \"$base_name.inp\" with a cell of type \"$type\"")
                @test volcheck
            end
        end
    end
end

@testset "facesetgeneration" begin
    filename = joinpath(@__DIR__, "test_files", "compact_tension.inp")
    grid = get_ferrite_grid(filename)
    @test create_faceset(grid, getnodeset(grid, "Hole")) == create_faceset(grid, getnodeset(grid, "Hole"), getcellset(grid, "Hole"))    # Test that including cells doesn't change the created sets
end

@testset "exceptions" begin
    filename = joinpath(@__DIR__, "runtests.jl")
    @test_throws FerriteMeshParser.UndetectableMeshFormatError get_ferrite_grid(filename)

    filename = joinpath(@__DIR__, "test_files", "twoinstances.inp")
    @test_throws FerriteMeshParser.InvalidFileContent get_ferrite_grid(filename)

    filename = joinpath(@__DIR__, "test_files", "unsupported_element.inp")
    @test_throws FerriteMeshParser.UnsupportedElementType get_ferrite_grid(filename)

end

@testset "ordering" begin
    filename0 = joinpath(@__DIR__, "test_files", "2D_UnitArea_Linear.inp")
    grid0 = get_ferrite_grid(filename0)

    filename1 = joinpath(@__DIR__, "test_files", "2D_UnitArea_Linear_flipelementorder.inp")
    grid1 = get_ferrite_grid(filename1)

    filename2 = joinpath(@__DIR__, "test_files", "2D_UnitArea_Linear_perturbnodeorder.inp")
    grid2 = get_ferrite_grid(filename2)

    # comparing grid0==grid<i> doesn't work, so check each field:
    for key in fieldnames(typeof(grid0))
        @test getfield(grid0, key) == getfield(grid1, key)
        @test getfield(grid0, key) == getfield(grid2, key)
    end
end

@testset "caseinsensitive" begin
    filename = joinpath(@__DIR__, "test_files", "uppercase_test.inp")
    grid = get_ferrite_grid(filename)
    @test getncells(grid) == 2
    @test getnnodes(grid) == 6
end

@testset "sets" begin
    filename = joinpath(@__DIR__, "test_files", "generated_set.inp")
    grid = get_ferrite_grid(filename)
    @test getcellset(grid, "lower") == Set(1:8)
    @test getnodeset(grid, "lower") == Set((1,  2,  3,  4,  7,  8,  9, 10, 11, 12, 13, 14, 20, 21, 22))

    filename = joinpath(@__DIR__, "test_files", "2D_UnitArea_Quadratic.inp")
    grid = get_ferrite_grid(filename)
    @test getcellset(grid, "mysetname") == Set((1,))
end

@testset "custom_cells" begin
    filename = joinpath(@__DIR__, "test_files", "unsupported_element.inp")
    grid = get_ferrite_grid(filename; user_elements=Dict("CPS5R"=>Cell{2,4,4}))
    @test getcellset(grid, "CPS5R") == Set((1,))
    @test isa(getcells(grid)[1], Cell{2,4,4})
end

@testset "specify_file_type" begin
    filename = joinpath(@__DIR__, "test_files", "2D_UnitArea_Quadratic.inp")
    filename_cp = joinpath(@__DIR__, "tmp.txt")
    cp(filename, filename_cp; force=true)
    grid_cp = get_ferrite_grid(filename_cp; meshformat=FerriteMeshParser.AbaqusMeshFormat())
    grid = get_ferrite_grid(filename)
    @test grid_cp == grid
    rm(filename_cp)
end


text:
# specialize this function with required MatchEvent values
function analyzer(::Val{MatchEvent}, pst::PlayState) where {MatchEvent}
    nothing
end

function tablesize(stmt::SQLite.Stmt)::Int
    return Tables.rowtable(DBInterface.execute(stmt))[1][1]
end

function tabledata(stmt::SQLite.Stmt, args::Tuple{Int,Int,Int,Int})::Array{
        NamedTuple{(:id, :timestamp, :content),
        Tuple{String,Int64,Array{UInt8,1}}},1}
    return Tables.rowtable(DBInterface.execute(stmt, args))
end

function analyseDatabase(dbpath::String;
        readsize::Int = 1000, offset::Int = 0, limit::Int = 0,
        mindt::DateTime = DateTime(2000), maxdt::DateTime = now(UTC))

    # open database and compile sql statement
    if !isfile(dbpath) error("Database not found") end
    db = SQLite.DB(dbpath)
    selectsize = DBInterface.prepare(db, "SELECT COUNT(id) FROM records")
    selectdata = DBInterface.prepare(db, "SELECT * FROM records WHERE id IN
    (SELECT id FROM records LIMIT ? OFFSET ?) AND timestamp BETWEEN ? AND ?")
    mindate = Dates.datetime2epochms(mindt)
    maxdate = Dates.datetime2epochms(maxdt)

    # initialize main data structures
    if iszero(limit) limit = tablesize(selectsize) end
    playstates = [PlayState(undef) for i=1:Threads.nthreads()]

    # read and parse chunks of records
    while offset < limit

        # select required subtable
        table = tabledata(selectdata, (readsize, offset, mindate, maxdate))

        # parsing problems are fast when done in parallel
        Threads.@threads for i = 1:min(length(table), readsize)

            # get playstate and register log id
            st = playstates[Threads.threadid()]
            st.logid = table[i].id

            # decompress tenhou log contents into immutable string
            str = String(transcode(LZ4FrameDecompressor, table[i].content))

            status = 1
            strlen = sizeof(str)

            # read tags sequentially
            while status < strlen

                brk(c) = ('0' <= c <= '9' || c == ' ' || c == '/')

                tagbeg::Int = findnext('<', str, status)
                tagend::Int = findnext(brk, str, tagbeg)
                status::Int = findnext('>', str, tagend)

                tag = str[tagbeg+1:tagend-1]
                data = str[tagend:status-2]

                try
                    # parse using function dictionary
                    if isempty(tag)             ev = END(data, st)
                    elseif tag == "GO"          ev = GO(data, st)
                    elseif tag == "UN"          ev = UN(data, st)
                    elseif tag == "INIT"        ev = INIT(data, st)
                    elseif tag == "AGARI"       ev = AGA(data, st)
                    elseif tag == "RYUUKYOKU"   ev = RYU(data, st)
                    elseif tag == "T"           ev = DRAW(data, st, 1)
                    elseif tag == "U"           ev = DRAW(data, st, 2)
                    elseif tag == "V"           ev = DRAW(data, st, 3)
                    elseif tag == "W"           ev = DRAW(data, st, 4)
                    elseif tag == "D"           ev = DROP(data, st, 1)
                    elseif tag == "E"           ev = DROP(data, st, 2)
                    elseif tag == "F"           ev = DROP(data, st, 3)
                    elseif tag == "G"           ev = DROP(data, st, 4)
                    elseif tag == "N"           ev = MELD(data, st)
                    elseif tag == "REACH"       ev = RCH(data, st)
                    elseif tag == "DORA"        ev = DORA(data, st)
                    elseif tag == "BYE"         ev = BYE(data, st)
                    else                        ev = SKIP(data, st)
                    end

                    # and callback on any user event
                    if     ev == noevents       analyzer(Val(noevents), st)
                    elseif ev == matchset       analyzer(Val(matchset), st)
                    elseif ev == roundinit      analyzer(Val(roundinit), st)
                    elseif ev == roundwin       analyzer(Val(roundwin), st)
                    elseif ev == roundtie       analyzer(Val(roundtie), st)
                    elseif ev == matchend       analyzer(Val(matchend), st)
                    elseif ev == tiledraw       analyzer(Val(tiledraw), st)
                    elseif ev == tiledrop       analyzer(Val(tiledrop), st)
                    elseif ev == meldcall       analyzer(Val(meldcall), st)
                    elseif ev == riichicall     analyzer(Val(riichicall), st)
                    elseif ev == doraflip       analyzer(Val(doraflip), st)
                    elseif ev == playerdc       analyzer(Val(playerdc), st)
                    elseif ev == playerin       analyzer(Val(playerin), st)
                    else                        error("Unknown MatchEvent")
                    end

                catch ex
                    println("Error while parsing log: \t",
                            "http://tenhou.net/0/?log=", st.logid, "\t|\t",
                            " Round:", st.cycle, " Turn:", st.turn)
                    rethrow(ex)
                end

            end
        end

        # shift subtable offset
        offset = offset + readsize
    end

    return nothing
end

function analyseLog(dbpath::String, i::Int)

    # log indexing starts from 1 for consistency
    analyseDatabase(dbpath; readsize=1, offset=i-1, limit=i)

    return nothing
end


text:
using Random
using Distributions
using MultiAgentPathFinding
using MultiAgentAllocationTransit
using BenchmarkTools
using Statistics
using JSON
using JLD2
using Logging
global_logger(SimpleLogger(stderr, Logging.Warn))

rng = MersenneTwister(2345)

# Script arguments relating to transit files etc.
const city_params_file = "./data/sfmta/sf_params.toml"
const stop_coords_file = "./data/sfmta/stop_to_coords.json"
const trips_file = "./data/sfmta/trips.json"
const drone_params_file = "./data/drone_params.toml"
const bb_params_file = "./data/sfmta/sf_bb_params.toml"
const out_file = "./data/temp_mult_generic.json"
const city_travel_time_estimates = "./data/sfmta/sf_halton_tt_estimates.jld2"

# MAPF-TN params
const TRANSIT_CAP_RANGE = (2, 4)
const ECBS_WEIGHT = 1.1
const N_DEPOTS = 3

# Change this one
const N_AGENTS = 5

const N_SITES = 3*N_AGENTS

## Hard-code a bunch of depots and many more sites
const DEPOT1 = LatLonCoords((lat = 37.762892, lon = -122.472193))
const DEPOT2 = LatLonCoords((lat = 37.751751, lon = -122.410654))
const DEPOT3 = LatLonCoords((lat = 37.718779, lon = -122.462401))
depots = [DEPOT1, DEPOT2, DEPOT3]

city_params = parse_city_params(city_params_file)
bb_params = parse_city_params(bb_params_file)
drone_params = parse_drone_params(drone_params_file)

lat_dist = Uniform(bb_params.lat_start, bb_params.lat_end)
lon_dist = Uniform(bb_params.lon_start, bb_params.lon_end)

# Transit Graph Preprocessing
tg = load_transit_graph_latlong(stop_coords_file, trips_file, TRANSIT_CAP_RANGE, rng)
tg, stop_idx_to_trips =
                transit_graph_preprocessing(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean, drone_params)

# Load Halton stuff
@load city_travel_time_estimates halton_nn_tree city_halton_points travel_time_estimates

sites = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_SITES]
depot_sites = vcat(depots, sites)

# Load OTG stuff
otg = OffTransitGraph(depots = depots, sites = sites)
aug_trips_fws_dists = augmented_trip_meta_graph_fws_dists(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                                                      length(depots), length(sites),
                                                      vcat(depots, sites),
                                                      drone_params)
state_graph, depot_sites_to_vtx, trip_to_vtx_range = setup_state_graph(tg, otg)

# Set the cost function using the wrapper
env = MAPFTransitEnv(off_transit_graph = otg, transit_graph = tg, state_graph = state_graph,
                     agent_states = AgentState[], depot_sites_to_vtx = depot_sites_to_vtx, trip_to_vtx_range = trip_to_vtx_range,
                     stop_idx_to_trips = stop_idx_to_trips, aug_trips_fws_dists = aug_trips_fws_dists,
                     drone_params = drone_params, dist_fn = MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                     curr_site_points = [], threshold_global_conflicts = 10)

cost_fn(i, j) = allocation_cost_wrapper_estimate(env, ECBS_WEIGHT, N_DEPOTS, N_SITES,
                                                                     halton_nn_tree, travel_time_estimates, i, j)

agent_tours = task_allocation(N_DEPOTS, N_SITES, N_AGENTS,
                              depot_sites, cost_fn)
@show agent_tours

agent_tasks = get_agent_task_set(agent_tours, N_DEPOTS, N_SITES)

true_n_agents = length(agent_tasks)
@show true_n_agents

# Now reset agent_states
env.agent_states = [AgentState(task=agt_task) for agt_task in agent_tasks]
env.curr_site_points = zeros(Int64, true_n_agents)
env.plan_ref_times = zeros(true_n_agents)

initial_states = Vector{MAPFTransitVertexState}(undef, true_n_agents)
for i = 1:true_n_agents
 initial_states[i] = env.state_graph.vertices[env.depot_sites_to_vtx[string("d-", agent_tasks[i].origin)]]
end

solver = ECBSSolver{MAPFTransitVertexState,MAPFTransitAction,Float64,Makespan,MAPFTransitConflict,MAPFTransitConstraints,MAPFTransitEnv}(env = env, weight = ECBS_WEIGHT)

solution = search!(solver, initial_states)

env_copy = deepcopy(env)
solution_copy = deepcopy(solution)

did_replan_indiv, el_time_indiv = replan_individual!(env, solution, N_DEPOTS, N_SITES, agent_tours, ECBS_WEIGHT)

did_replan_collec, el_time_collec, new_soln_collec = replan_collective!(env_copy, solution_copy, N_DEPOTS, N_SITES, agent_tours, ECBS_WEIGHT)

@show el_time_collec


text:
# src/BoundedRationals.jl
# Rational Numbers with Bounded Denominators

__precompile__(true)

"""```
module BoundedRationals
```
Rational Numbers with Bounded Denominators
See https://github.com/bhgomes/BoundedRationals.jl for more details.
"""
module BoundedRationals

include("core.jl")
include("impl.jl")

end  # module BoundedRationals


text:
struct Biped{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    g
    μ

    # torso
    l_torso
    d_torso
    m_torso
    J_torso

    # leg 1
        # thigh
    l_thigh1
    d_thigh1
    m_thigh1
    J_thigh1

        # calf
    l_calf1
    d_calf1
    m_calf1
    J_calf1

    # leg 2
        # thigh
    l_thigh2
    d_thigh2
    m_thigh2
    J_thigh2

        # calf
    l_calf2
    d_calf2
    m_calf2
    J_calf2

    # joint limits
    qL
    qU

    # torque limits
    uL
    uU

    nq
    nu
    nc
    nf
    nb
    ns

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# Dimensions
nq = 2 + 5                # configuration dimension
nu = 4                    # control dimension
nc = 2                    # number of contact points
nf = 2                    # number of parameters for friction cone
nb = nc * nf
ns = 1

# World parameters
μ = 1.0      # coefficient of friction
g = 9.81     # gravity

# Model parameters
m_torso = 0.5 + 0.48 * 2.0
m_thigh = 0.8112
m_calf = 0.3037

J_torso = 0.0029
J_thigh = 0.00709
J_calf = 0.00398

l_torso = 0.15 + 0.15
l_thigh = 0.2755
l_calf = 0.308

d_torso = 0.0342
d_thigh = 0.2176
d_calf = 0.1445

n = 2 * nq
m = nu + nc + nb + nc + nb + ns
d = nq

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

qL = -Inf * ones(nq)
qU = Inf * ones(nq)

uL = -100.0 * ones(nu)
uU = 100.0 * ones(nu)

function kinematics_1(model::Biped, q; body = :torso, mode = :ee)
	x = q[1]
	z = q[2]

	if body == :torso
		l = model.l_torso
		d = model.d_torso
		θ = q[3]
		if mode == :ee
			return [x - l * sin(θ); z + l * cos(θ)]
		elseif mode == :com
			return [x - d * sin(θ); z + d * cos(θ)]
		end
	elseif body == :thigh_1
		l = model.l_thigh1
		d = model.d_thigh1
		θ = q[4]
	elseif body == :thigh_2
		l = model.l_thigh2
		d = model.d_thigh2
		θ = q[6]
	else
		@error "incorrect body specification"
	end

	if mode == :ee
		return [x + l * sin(θ); z - l * cos(θ)]
	elseif mode == :com
		return [x + d * sin(θ); z - d * cos(θ)]
	else
		@error "incorrect mode specification"
	end
end

function jacobian_1(model::Biped, q; body = :torso, mode = :ee)
	jac = zeros(eltype(q), 2, model.nq)
	jac[1, 1] = 1.0
	jac[2, 2] = 1.0
	if body == :torso
		r = mode == :ee ? model.l_torso : model.d_torso
		θ = q[3]
		jac[1, 3] = -r * cos(θ)
		jac[2, 3] = -r * sin(θ)
	elseif body == :thigh_1
		r = mode == :ee ? model.l_thigh1 : model.d_thigh1
		θ = q[4]
		jac[1, 4] = r * cos(θ)
		jac[2, 4] = r * sin(θ)
	elseif body == :thigh_2
		r = mode == :ee ? model.l_thigh2 : model.d_thigh2
		θ = q[6]
		jac[1, 6] = r * cos(θ)
		jac[2, 6] = r * sin(θ)
	else
		@error "incorrect body specification"
	end

	return jac
end

function kinematics_2(model::Biped, q; body = :calf_1, mode = :ee)

	if body == :calf_1
		p = kinematics_1(model, q, body = :thigh_1, mode = :ee)

		θb = q[5]

		lb = model.l_calf1
		db = model.d_calf1
	elseif body == :calf_2
		p = kinematics_1(model, q, body = :thigh_2, mode = :ee)

		θb = q[7]

		lb = model.l_calf2
		db = model.d_calf2
	else
		@error "incorrect body specification"
	end

	if mode == :ee
		return p + [lb * sin(θb); -1.0 * lb * cos(θb)]
	elseif mode == :com
		return p + [db * sin(θb); -1.0 * db * cos(θb)]
	else
		@error "incorrect mode specification"
	end
end

function jacobian_2(model::Biped, q; body = :calf_1, mode = :ee)

	if body == :calf_1
		jac = jacobian_1(model, q, body = :thigh_1, mode = :ee)

		θb = q[5]

		r = mode == :ee ? model.l_calf1 : model.d_calf1

		jac[1, 5] += r * cos(θb)
		jac[2, 5] += r * sin(θb)
	elseif body == :calf_2
		jac = jacobian_1(model, q, body = :thigh_2, mode = :ee)

		θb = q[7]

		r = mode == :ee ? model.l_calf2 : model.d_calf2

		jac[1, 7] += r * cos(θb)
		jac[2, 7] += r * sin(θb)
	else
		@error "incorrect body specification"
	end

	return jac
end

# Lagrangian

function lagrangian(model::Biped, q, q̇)
	L = 0.0

	# torso
	p_torso = kinematics_1(model, q, body = :torso, mode = :com)
	J_torso = jacobian_1(model, q, body = :torso, mode = :com)
	v_torso = J_torso * q̇

	L += 0.5 * model.m_torso * transpose(v_torso) * v_torso
	L += 0.5 * model.J_torso * q̇[3]^2.0
	L -= model.m_torso * model.g * p_torso[2]

	# thigh 1
	p_thigh_1 = kinematics_1(model, q, body = :thigh_1, mode = :com)
	J_thigh_1 = jacobian_1(model, q, body = :thigh_1, mode = :com)
	v_thigh_1 = J_thigh_1 * q̇

	L += 0.5 * model.m_thigh1 * transpose(v_thigh_1) * v_thigh_1
	L += 0.5 * model.J_thigh1 * q̇[4]^2.0
	L -= model.m_thigh1 * model.g * p_thigh_1[2]

	# leg 1
	p_calf_1 = kinematics_2(model, q, body = :calf_1, mode = :com)
	J_calf_1 = jacobian_2(model, q, body = :calf_1, mode = :com)
	v_calf_1 = J_calf_1 * q̇

	L += 0.5 * model.m_calf1 * transpose(v_calf_1) * v_calf_1
	L += 0.5 * model.J_calf1 * q̇[5]^2.0
	L -= model.m_calf1 * model.g * p_calf_1[2]

	# thigh 2
	p_thigh_2 = kinematics_1(model, q, body = :thigh_2, mode = :com)
	J_thigh_2 = jacobian_1(model, q, body = :thigh_2, mode = :com)
	v_thigh_2 = J_thigh_2 * q̇

	L += 0.5 * model.m_thigh2 * transpose(v_thigh_2) * v_thigh_2
	L += 0.5 * model.J_thigh2 * q̇[6]^2.0
	L -= model.m_thigh2 * model.g * p_thigh_2[2]

	# leg 2
	p_calf_2 = kinematics_2(model, q, body = :calf_2, mode = :com)
	J_calf_2 = jacobian_2(model, q, body = :calf_2, mode = :com)
	v_calf_2 = J_calf_2 * q̇

	L += 0.5 * model.m_calf2 * transpose(v_calf_2) * v_calf_2
	L += 0.5 * model.J_calf2 * q̇[7]^2.0
	L -= model.m_calf2 * model.g * p_calf_2[2]

	return L
end

function _dLdq(model::Biped, q, q̇)
	Lq(x) = lagrangian(model, x, q̇)
	ForwardDiff.gradient(Lq, q)
end

function _dLdq̇(model::Biped, q, q̇)
	Lq̇(x) = lagrangian(model, q, x)
	ForwardDiff.gradient(Lq̇, q̇)
end

# Methods
function M_func(model::Biped, q)
	M = Diagonal([0.0, 0.0, model.J_torso, model.J_thigh1, model.J_calf1, model.J_thigh2, model.J_calf2])

	# torso
	J_torso = jacobian_1(model, q, body = :torso, mode = :com)
	M += model.m_torso * J_torso' * J_torso

	# thigh 1
	J_thigh_1 = jacobian_1(model, q, body = :thigh_1, mode = :com)
	M += model.m_thigh1 * J_thigh_1' * J_thigh_1

	# leg 1
	J_calf_1 = jacobian_2(model, q, body = :calf_1, mode = :com)
	M += model.m_calf1 * J_calf_1' * J_calf_1

	# thigh 2
	J_thigh_2 = jacobian_1(model, q, body = :thigh_2, mode = :com)
	M += model.m_thigh2 * J_thigh_2' * J_thigh_2

	# leg 2
	J_calf_2 = jacobian_2(model, q, body = :calf_2, mode = :com)
	M += model.m_calf2 * J_calf_2' * J_calf_2

	return M
end

function _C_func(model::Biped, q, q̇)
	tmp_q(z) = _dLdq̇(model, z, q̇)
	tmp_q̇(z) = _dLdq̇(model, q, z)

	ForwardDiff.jacobian(tmp_q, q) * q̇ - _dLdq(model, q, q̇)
end

function ϕ_func(model::Biped, q)
	p_calf_1 = kinematics_2(model, q, body = :calf_1, mode = :ee)
	p_calf_2 = kinematics_2(model, q, body = :calf_2, mode = :ee)

	@SVector [p_calf_1[2], p_calf_2[2]]
end

function B_func(model::Biped, q)
	@SMatrix [0.0 0.0 0.0 1.0 0.0 0.0 0.0;
			  0.0 0.0 0.0 0.0 1.0 0.0 0.0;
			  0.0 0.0 0.0 0.0 0.0 1.0 0.0;
			  0.0 0.0 0.0 0.0 0.0 0.0 1.0]
end

function N_func(model::Biped, q)
	J_calf_1 = jacobian_2(model, q, body = :calf_1, mode = :ee)
	J_calf_2 = jacobian_2(model, q, body = :calf_2, mode = :ee)

	return [view(J_calf_1, 2:2, :);
			view(J_calf_2, 2:2, :)]
end

function _P_func(model::Biped, q)
	J_calf_1 = jacobian_2(model, q, body = :calf_1, mode = :ee)
	J_calf_2 = jacobian_2(model, q, body = :calf_2, mode = :ee)

	return [view(J_calf_1, 1:1, :);
			view(J_calf_2, 1:1, :)]
end

function P_func(model::Biped, q)
	J_calf_1 = jacobian_2(model, q, body = :calf_1, mode = :ee)
	J_calf_2 = jacobian_2(model, q, body = :calf_2, mode = :ee)
	map = [1.0; -1.0]

	return [map * view(J_calf_1, 1:1, :);
			map * view(J_calf_2, 1:1, :)]
end

function friction_cone(model::Biped, u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	return @SVector [model.μ * λ[1] - sum(view(b, 1:2)),
					 model.μ * λ[2] - sum(view(b, 3:4))]
end

function maximum_dissipation(model::Biped{Discrete, FixedTime}, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	ψ = view(u, model.idx_ψ)
	ψ_stack = [ψ[1] * ones(2); ψ[2] * ones(2)]
	η = view(u, model.idx_η)
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function maximum_dissipation(model::Biped{Discrete, FreeTime}, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	ψ = view(u, model.idx_ψ)
	ψ_stack = [ψ[1] * ones(2); ψ[2] * ones(2)]
	η = view(u, model.idx_η)
	h = u[end]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function no_slip(model::Biped{Discrete, FixedTime}, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	λ = view(u, model.idx_λ)
	s = view(u, model.idx_s)
	return s[1] - (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function no_slip(model::Biped{Discrete, FreeTime}, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	λ = view(u, model.idx_λ)
	s = view(u, model.idx_s)
	h = u[end]
	return s[1] - (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function fd(model::Biped{Discrete, FixedTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

    SVector{14}([q2⁺ - q2⁻;
    ((1.0 / h) * (M_func(model, q1) * (SVector{7}(q2⁺) - SVector{7}(q1))
    - M_func(model, q2⁺) * (SVector{7}(q3) - SVector{7}(q2⁺)))
    + h * (transpose(B_func(model, q3)) * SVector{4}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{2}(λ)
    + transpose(P_func(model, q3)) * SVector{4}(b))
    - h * C_func(model, q3, (q3 - q2⁺) / h)
	+ h * w)])
end

function fd(model::Biped{Discrete, FreeTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	h = u[end]

    SVector{14}([q2⁺ - q2⁻;
    ((1.0 / h) * (M_func(model, q1) * (SVector{7}(q2⁺) - SVector{7}(q1))
    - M_func(model, q2⁺) * (SVector{7}(q3) - SVector{7}(q2⁺)))
    + h * (transpose(B_func(model, q3)) * SVector{4}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{2}(λ)
    + transpose(P_func(model, q3)) * SVector{4}(b))
    - h * C_func(model, q3, (q3 - q2⁺) / h)
	+ h * w)])
end

model = Biped{Discrete, FixedTime}(n, m, d,
			  g, μ,
			  l_torso, d_torso, m_torso, J_torso,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_calf, d_calf, m_calf, J_calf,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_calf, d_calf, m_calf, J_calf,
			  qL, qU,
			  uL, uU,
			  nq,
			  nu,
			  nc,
			  nf,
			  nb,
			  ns,
			  idx_u,
			  idx_λ,
			  idx_b,
			  idx_ψ,
			  idx_η,
			  idx_s)

#NOTE: if a new model is instantiated, re-run the lines below
@variables z_sym[1:model.n]
l(z) = lagrangian(model, view(z, 1:model.nq), view(z, model.nq .+ (1:model.nq)))
_l = simplify.(l(z_sym))
_dL = ModelingToolkit.gradient(_l, z_sym)
_dLq = view(_dL, 1:model.nq)
_dLq̇ = view(_dL, model.nq .+ (1:model.nq))
_ddL = ModelingToolkit.sparsehessian(_l, z_sym)
_ddLq̇q = view(_ddL, model.nq .+ (1:model.nq), 1:model.nq)

dL = eval(ModelingToolkit.build_function(_dL, z_sym)[1])
dLq = eval(ModelingToolkit.build_function(_dLq, z_sym)[1])
dLq̇ = eval(ModelingToolkit.build_function(_dLq̇, z_sym)[1])
ddLq̇q = eval(ModelingToolkit.build_function(_ddLq̇q, z_sym)[1])
ddL = eval(ModelingToolkit.build_function(_ddL, z_sym)[1])

function C_func(model::Biped, q, q̇)
	ddLq̇q([q; q̇]) * q̇ - dLq([q; q̇])
end

qq = rand(model.nq)
vv = rand(model.nq)
norm(C_func(model, qq, vv) - _C_func(model, qq, vv))
dL(qq)

# visualization
function visualize!(vis, model::Biped, q;
      r = 0.035, Δt = 0.1)

	default_background!(vis)

	torso = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l_torso),
		convert(Float32, 0.025))
	setobject!(vis["torso"], torso,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	thigh_1 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_thigh1),
		convert(Float32, 0.025))
	setobject!(vis["thigh1"], thigh_1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	calf_1 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_calf1),
		convert(Float32, 0.025))
	setobject!(vis["calf1"], calf_1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	thigh_2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_thigh2),
		convert(Float32, 0.025))
	setobject!(vis["thigh2"], thigh_2,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	calf_2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_calf2),
		convert(Float32, 0.025))
	setobject!(vis["calf2"], calf_2,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	setobject!(vis["foot1"], Sphere(Point3f0(0.0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
	setobject!(vis["foot2"], Sphere(Point3f0(0.0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
	setobject!(vis["knee1"], Sphere(Point3f0(0.0),
		convert(Float32, 0.025)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
	setobject!(vis["knee2"], Sphere(Point3f0(0.0),
		convert(Float32, 0.025)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
	setobject!(vis["hip"], Sphere(Point3f0(0.0),
		convert(Float32, 0.025)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
	setobject!(vis["torso_top"], Sphere(Point3f0(0.0),
		convert(Float32, 0.025)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q)
	p_shift = [0.0; 0.0; r]
	for t = 1:T
		MeshCat.atframe(anim, t) do
			p = [q[t][1]; 0.0; q[t][2]] + p_shift

			k_torso = kinematics_1(model, q[t], body = :torso, mode = :ee)
			p_torso = [k_torso[1], 0.0, k_torso[2]] + p_shift

			k_thigh_1 = kinematics_1(model, q[t], body = :thigh_1, mode = :ee)
			p_thigh_1 = [k_thigh_1[1], 0.0, k_thigh_1[2]] + p_shift

			k_calf_1 = kinematics_2(model, q[t], body = :calf_1, mode = :ee)
			p_calf_1 = [k_calf_1[1], 0.0, k_calf_1[2]] + p_shift

			k_thigh_2 = kinematics_1(model, q[t], body = :thigh_2, mode = :ee)
			p_thigh_2 = [k_thigh_2[1], 0.0, k_thigh_2[2]] + p_shift

			k_calf_2 = kinematics_2(model, q[t], body = :calf_2, mode = :ee)
			p_calf_2 = [k_calf_2[1], 0.0, k_calf_2[2]] + p_shift

			settransform!(vis["thigh1"], cable_transform(p, p_thigh_1))
			settransform!(vis["calf1"], cable_transform(p_thigh_1, p_calf_1))
			settransform!(vis["thigh2"], cable_transform(p, p_thigh_2))
			settransform!(vis["calf2"], cable_transform(p_thigh_2, p_calf_2))
			settransform!(vis["torso"], cable_transform(p_torso,p))
			settransform!(vis["foot1"], Translation(p_calf_1))
			settransform!(vis["foot2"], Translation(p_calf_2))
			settransform!(vis["knee1"], Translation(p_thigh_1))
			settransform!(vis["knee2"], Translation(p_thigh_2))
			settransform!(vis["hip"], Translation(p))
			settransform!(vis["torso_top"], Translation(p_torso))

		end
	end

	settransform!(vis["/Cameras/default"],
	   compose(Translation(0.0 , 1.0 , -1.0), LinearMap(RotZ(-pi / 2.0))))

	MeshCat.setanimation!(vis, anim)
end

function initial_configuration(model, θ_torso, θ_thigh_1, θ_leg_1, θ_thigh_2)
    q1 = zeros(model.nq)
    q1[3] = θ_torso
    q1[4] = θ_thigh_1
    q1[5] = θ_leg_1
    z1 = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])
    q1[6] = θ_thigh_2
    q1[7] = -1.0 * acos((z1 - model.l_thigh2 * cos(q1[6])) / model.l_calf2)
    q1[2] = z1

    p1 = kinematics_2(model, q1, body = :calf_1, mode = :ee)
    p2 = kinematics_2(model, q1, body = :calf_2, mode = :ee)
    @show stride = abs(p1[1] - p2[1])

    q1[1] = -1.0 * p1[1]

    qM = copy(q1)
    qM[4] = q1[6]
    qM[5] = q1[7]
    qM[6] = q1[4]
    qM[7] = q1[5]
    qM[1] = abs(p2[1])

    pM_1 = kinematics_2(model, qM, body = :calf_1, mode = :ee)
    pM_2 = kinematics_2(model, qM, body = :calf_2, mode = :ee)

    qT = copy(q1)
    qT[1] = 2 * stride

    pT_1 = kinematics_2(model, qT, body = :calf_1, mode = :ee)
    pT_2 = kinematics_2(model, qT, body = :calf_2, mode = :ee)

    return q1, qM, qT
end


text:
# ============================================================
# Governing Equations
# ============================================================

include("eq_scalar.jl")
include("eq_advection.jl")
include("eq_euler.jl")
include("eq_bgk.jl")


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_semiinfinite/algorithms/sip_hybrid.jl
# FUTURE FEATURE... NOT CURRENTLY FUNCTIONAL
# Defines the SIP-hybrid algorithm which implements Algorithm #2 of Djelassi,
# Hatim, and Alexander Mitsos. "A hybrid discretization algorithm with guaranteed
# feasibility for the global solution of semi-infinite programs."
# Journal of Global Optimization 68.2 (2017): 227-253.
#############################################################################

"""
    SIPHybrid

Specifies that the SIPHybrid algorithm which implements Algorithm #2 of Djelassi,
Hatim, and Alexander Mitsos. "A hybrid discretization algorithm with guaranteed
feasibility for the global solution of semi-infinite programs." Journal of Global
Optimization 68.2 (2017): 227-253 should be used.
"""
struct SIPHybrid <: AbstractSIPAlgo end

function set_tolerance!(t::DefaultExt, alg::SIPHybrid, s::LowerLevel1, m::JuMP.Model, sr::SIPSubResult, i::Int)
    set_tolerance_inner!(t, alg, s, m, sr.eps_l[i])
end
function set_tolerance!(t::DefaultExt, alg::SIPHybrid, s::LowerLevel2, m::JuMP.Model, sr::SIPSubResult, i::Int)
    set_tolerance_inner!(t, alg, s, m, sr.eps_u[i])
end
function set_tolerance!(t::DefaultExt, alg::SIPHybrid, s::LowerProblem, m::JuMP.Model, sr::SIPSubResult, i::Int)
    set_tolerance_inner!(t, alg, s, m, sr.lbd.tol)
end
function set_tolerance!(t::DefaultExt, alg::SIPHybrid, s::UpperProblem, m::JuMP.Model, sr::SIPSubResult, i::Int)
    set_tolerance_inner!(t, alg, s, m, sr.ubd.tol)
end

function get_disc_set(t::ExtensionType, alg::SIPHybrid, s::S, sr::SIPSubResult,
                      i::Int) where S <: Union{LowerProblem,UpperProblem,ResProblem}
    sr.disc_l[i]
end

function sip_solve!(t::ExtensionType, alg::SIPHybrid, buffer::SIPSubResult, prob::SIPProblem,
                    result::SIPResult, cb::SIPCallback)

    verb = prob.verbosity

    # begin main solution loop
    @label main_iteration

    # solve lower bounding problem and check feasibility
    sip_bnd!(t, alg, LowerProblem(), buffer, result, prob, cb)
    result.lower_bound = buffer.lbd.obj_val
    if !buffer.lbd.feas
        result.feasibility = false
        println("Terminated: lower bounding problem infeasible.")
        @goto main_end
    end
    print_summary!(LowerProblem(), verb, buffer)

    # solve inner program and update lower discretization set
    is_llp1_nonpositive = true
    for i = 1:prob.nSIP
        sip_llp!(t, alg, LowerLevel1(), result, buffer, prob, cb, i)
        buffer.disc_l_buffer .= buffer.llp1.sol
        print_summary!(LowerLevel1(), verb, buffer, i)
        if buffer.llp1.obj_bnd <= 0.0
            continue
        elseif buffer.llp1.obj_val > 0.0
            push!(buffer.disc_l[i], deepcopy(buffer.disc_l_buffer))
            is_llp1_nonpositive = false
        else
            buffer.eps_l[i] = (buffer.llp1.obj_bnd - buffer.llp1.obj_val)/buffer.r_l
            is_llp1_nonpositive = false
        end
    end

    # if the lower problem is feasible then it's solution is the optimal value
    if is_llp1_nonpositive
        result.upper_bound = buffer.lbd.obj_val
        result.xsol .= buffer.lbd.sol
        result.feasibility = true
        @goto main_end
    end

    # solve upper bounding problem, if feasible solve lower level problem,
    # and potentially update upper discretization set
    @label upper_problem
    sip_bnd!(t, alg, UpperProblem(), buffer, result, prob, cb)
    print_summary!(UpperProblem(), verb, buffer)
    if buffer.ubd.feas
        is_llp2_nonpositive = true
        for i = 1:prob.nSIP
            sip_llp!(t, alg, LowerLevel2(), result, buffer, prob, cb, i)
            buffer.disc_l_buffer .= buffer.llp2.sol
            print_summary!(LowerLevel2(), verb, buffer, i)
            if buffer.llp2.obj_bnd <= 0.0
                buffer.eps_g[i] /= buffer.r_g
                continue
            else
                push!(buffer.disc_l[i], deepcopy(buffer.disc_l_buffer))
                is_llp2_nonpositive = false
            end
        end
        if is_llp2_nonpositive
            if buffer.ubd.obj_val <= result.upper_bound
                result.upper_bound = buffer.ubd.obj_val
                result.xsol .= buffer.ubd.sol
            end
        end
    else
        buffer.eps_g ./= buffer.r_g
    end
    check_convergence(result, prob.abs_tolerance, verb) && @goto main_end

    # solve restriction problem updating lower and upper bound as appropriate
    @label res_problem
    sip_res!(t, alg, buffer, result, prob, cb)
    if buffer.res.obj_bnd < 0
        result.lower_bound = buffer.res.obj_val
        buffer.res.obj_val = 0.5*(result.upper_bound + result.lower_bound)
        result.res_iteration_number += 1
        @goto res_problem
    elseif buffer.res.obj_val > 0
        is_llp3_nonpositive = false
        for i = 1:prob.nSIP
            sip_llp!(t, alg, LowerLevel3(), result, buffer, prob, cb, i)
            buffer.disc_l_buffer .= buffer.llp3.sol
            print_summary!(LowerLevel3(), verb, buffer, i)
            if buffer.llp3.obj_bnd <= 0.0
                if buffer.res.obj_bnd/buffer.r_g < buffer.eps_g[i]
                    buffer.eps_g[i] = buffer.res.obj_bnd/buffer.r_g
                end
            elseif (buffer.llp3.obj_bnd > 0.0) &&
                   (result.res_iteration_number < prob.res_iteration_limit)
                is_llp3_nonpositive = false
                push!(buffer.disc_l[i], deepcopy(buffer.disc_l_buffer))
                result.res_iteration_number += 1
                @goto res_problem
            else
                result.res_iteration_number = 1
                @goto main_iteration
            end
        end
        if is_llp3_nonpositive
            if buffer.res.obj_val <= result.upper_bound
                result.upper_bound = buffer.res.obj_val
                result.xsol = buffer.res.sol
                result.res_iteration_number = 1
                @goto upper_problem
            end
        end
    else
        result.res_iteration_number = 1
        @goto main_iteration
    end

    # print iteration information and advance
    print_int!(verb, prob, result, buffer.r_g)
    result.iteration_number += 1
    result.iteration_number < prob.iteration_limit && @goto main_iteration

    @label main_end
    return nothing
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("classical.jl")
using .Classical
using Plots

r = Classical.ring(64;xf=4.0)
fig = Plots.plot()

for b in [0.0,1.0,2.0]#,-0.5,1.0,2.0,3.0]
    r.β = b
    for j in 1:1000000
        Classical.iterate!(r)
    end
    Plots.plot!(r.x_[1:24],r.σ_[1:24];legend=false)
    display(fig)
end

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
type KronSumMat{T} <: AbstractKronMat{T}
    terms::Vector{AbstractMatrix{T}}
end

function KronSumMat(terms...)
    T = promote_type([eltype(term) for term in terms]...)
    KronSumMat(AbstractMatrix{T}[terms...])
end

#ok, evidently there is going to be too much dry violation here.. should correct before proceeding too far.


text:
@testset "Macro1" begin
    @bind begin
        x = 5.6
        y = 10.1
    end
    expr = x^2 + 10*x*y + 4*y^2
    res = seval(expr, (x,y))
    @test res == 10.0
end

@testset "ops1" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    @bind env begin
        x = 0.56
        y = 0.101
    end
    expr = exp(x^2 + 10*x*y + 4*y^2)
    res = seval(expr, (x,y), env, cache)
    @test res == seval(exp(x^2 + 10*x*y + 4*y^2) * (10*x + 8*y) * (2*x + 10*y) + exp(x^2 + 10*x*y + 4*y^2) * 10, env)
end

@testset "Promotion1" begin
    a = symbolic(:a)
    b = symbolic(:b, Int)
    x = [a, b]
    @test typeof(b) <: SymbolicVariable{Int}
    @test typeof(x[1]) == SymbolicVariable{Float64}
    @test typeof(x[2]) == SymbolicVariable{Float64}
end

@testset "Promotion2" begin
    a = symbolic(:a)
    x = [a, 1.0]
    @test typeof(x[1]) == SymbolicVariable{Float64}
    @test typeof(x[2]) == SymbolicValue{Float64}
end

@testset "Promotion3" begin
    a = symbolic(:a)
    x = [a, 1]
    @test typeof(x[1]) == SymbolicVariable{Float64}
    @test typeof(x[2]) == SymbolicValue{Float64}
end

@testset "Promotion4" begin
    a = symbolic(:a)
    x = a + 1
    println(x)
    @test typeof(x) == SymbolicExpression{Float64}
end

@testset "Var1" begin
    @bind begin
        x = 0.8
        y = 0.9
    end
    z = x + y
    result = seval(z)
    @test result == 0.8 + 0.9
end

@testset "Var3" begin
    @bind begin
        x = 0.8
        y = 0.9
    end
    z = x + y
    result = seval(z, x)
    @test result == 1.0
end


text:
module EMLHeatExchangerDetailed
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...heat_exchangers.EMLHEX_Engine
	import EMLtypes.length
	include("HeatExchangerDetailed/ShellandTubesBasic.jl")
	include("HeatExchangerDetailed/ShellandTubes_NTU.jl")
	include("HeatExchangerDetailed/ShellandTubes_LMTD.jl")
end


text:
module CODATA2018

using PhysicalConstants, Measurements, Unitful, Roots

import Unitful: Ω, A, C, F, Hz, J, kg, K, m, mol, N, Pa, s, T, W
import PhysicalConstants: @constant, @derived_constant

## Helper functions for the Wien displacement law constants
_λf(x) = (x - 5) * exp(x) + 5
_Dλf(x) = (x - 4) * exp(x)
const _λx0 = 4.965114231744276
_λx(T) = find_zero((_λf, _Dλf), T(_λx0), Roots.Newton())

_νf(x) = (x - 3) * exp(x) + 3
_Dνf(x) = (x - 2) * exp(x)
const _νx0 = 2.8214393721220787
_νx(T) = find_zero((_νf, _Dνf), T(_νx0), Roots.Newton())
## end

@constant(FineStructureConstant, α, "Fine-structure constant", 7.297_352_5693e-3,
          BigFloat(72_973_525_693)/BigFloat(10_000_000_000_000), Unitful.NoUnits,
          1.1e-12, BigFloat(11)/BigFloat(10_000_000_000_000), "CODATA 2018")
@constant(BohrRadius, a_0, "Bohr radius", 5.291_772_109_03e-11,
          BigFloat(529_177_210_903)/BigFloat(10_000_000_000_000_000_000_000), m,
          8e-21, BigFloat(80)/BigFloat(10_000_000_000_000_000_000_000), "CODATA 2018")
@constant(StandardAtmosphere, atm, "Standard atmosphere", 101_325.0, BigFloat(101_325), Pa,
          0.0, BigFloat(0), "CODATA 2018")
@constant(SpeedOfLightInVacuum, c_0, "Speed of light in vacuum", 299_792_458.0,
          BigFloat(299_792_458.0), m / s, 0.0, BigFloat(0), "CODATA 2018")
@constant(VacuumMagneticPermeability, µ_0, "Vacuum magnetic permeability", 1.256_637_062_12e-6,
          BigFloat(125_663_706_212)/BigFloat(100_000_000_000_000_000), N * A^-2, 1.9e-16,
          BigFloat(19)/BigFloat(100_000_000_000_000_000), "CODATA 2018")
@constant(VacuumElectricPermittivity, ε_0, "Vacuum electric permittivity", 8.854_187_8128e-12,
          BigFloat(88_541_878_128)/BigFloat(10_000_000_000_000_000_000_000), F * m^-1,
          1.3e-21, BigFloat(13)/BigFloat(10_000_000_000_000_000_000_000), "CODATA 2018")
@constant(ElementaryCharge, e, "Elementary charge", 1.602_176_634e-19,
          big(1_602_176_634)/big(10_000_000_000_000_000_000_000_000_000),
          C, 0.0, BigFloat(0.0), "CODATA 2018")
@constant(NewtonianConstantOfGravitation, G, "Newtonian constant of gravitation",
          6.674_30e-11, big(667_430)/big(10_000_000_000_000_000), m^3 * kg^-1 * s^-2,
          1.5e-15, big(15)/big(10_000_000_000_000_000), "CODATA 2018")
@constant(StandardAccelerationOfGravitation, g_n, "Standard acceleration of gravitation",
          9.806_65, big(980_665)/big(100_000), m * s^-2, 0, 0, "CODATA 2018")
@constant(PlanckConstant, h, "Planck constant", 6.626_070_15e-34,
          6_626_070_15/1_000_000_000_000_000_000_000_000_000_000_000_000_000_000,
          J * s, 0.0, BigFloat(0.0), "CODATA 2018")
@derived_constant(ReducedPlanckConstant, ħ, "Reduced Planck constant",
                  1.0545718176461565e-34, ustrip(big(h))/(2 * big(pi)), J * s,
                  measurement(h)/2pi, measurement(BigFloat, h)/(2 * big(pi)), "CODATA 2018")
@constant(BoltzmannConstant, k_B, "Boltzmann constant", 1.380_649e-23,
          BigFloat(1_380_649)/BigFloat(100_000_000_000_000_000_000_000_000_000), J * K^-1,
          0.0, BigFloat(0.0), "CODATA 2018")
@constant(BohrMagneton, µ_B, "Bohr magneton", 9.274_010_0783e-24,
          BigFloat(92_740_100_783)/BigFloat(10_000_000_000_000_000_000_000_000_000_000_000),
          J * T^-1, 2.8e-33,
          BigFloat(28)/BigFloat(100_00_000_000_000_000_000_000_000_000_000_000), "CODATA 2018")
@constant(ElectronMass, m_e, "Electron mass", 9.109_383_7015e-31,
          BigFloat(91_093_837_015)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000_000),
          kg, 2.8e-40,
          BigFloat(28)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(NeutronMass, m_n, "Neutron mass", 1.674_927_498_04e-27,
          BigFloat(167_492_749_804)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000),
          kg, 9.5e-37,
          BigFloat(95)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(ProtonMass, m_p, "Proton mass", 1.672_621_923_69e-27,
          BigFloat(167_262_192_369)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000),
          kg, 5.1e-37,
          BigFloat(51)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(AtomicMassConstant, m_u, "Atomic mass constant", 1.660_539_066_60e-27,
          BigFloat(166_053_906_660)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000),
          kg, 5.0e-37,
          BigFloat(50)/BigFloat(100_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(AvogadroConstant, N_A, "Avogadro constant", 6.022_140_76e23,
          BigFloat(602_214_076_000_000_000_000_000), mol^-1,
          0.0, BigFloat(0.0), "CODATA 2018")
@derived_constant(MolarGasConstant, R, "Molar gas constant", 8.314_462_618_153_24,
                  ustrip(big(N_A) * big(k_B)), J * mol^-1 * K^-1,
                  measurement(N_A) * measurement(k_B),
                  measurement(BigFloat, N_A) * measurement(BigFloat, k_B), "CODATA 2018")
@constant(RydbergConstant, R_∞, "Rydberg constant", 10_973_731.568_160,
          BigFloat(10_973_731_568_160)/BigFloat(1_000_000), m^-1,
          2.1e-5, BigFloat(21)/BigFloat(1_000_000), "CODATA 2018")
@derived_constant(StefanBoltzmannConstant, σ, "Stefan-Boltzmann constant", 5.670_374_419_184_4294e-8,
                  ustrip(2 * big(pi) ^ 5 * big(k_B) ^ 4 / (15 * big(h) ^ 3 * big(c_0) ^ 2)), W * m^-2 * K^-4,
                  (2 * pi ^ 5 * measurement(k_B) ^ 4) / (15 * measurement(h) ^ 3 * measurement(c_0) ^ 2),
                  (2 * big(pi) ^ 5 * measurement(BigFloat, k_B) ^ 4) / (15 * measurement(BigFloat, h) ^ 3 * measurement(BigFloat, c_0) ^ 2),
                  "CODATA 2018")
@constant(ThomsonCrossSection, σ_e, "Thomson cross section", 6.652_458_7321e-29,
          BigFloat(66_524_587_321)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          m^2, 6.0e-38,
          BigFloat(60)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@derived_constant(WienWavelengthDisplacementLawConstant, b, "Wien wavelength displacement law constant",
                  2.897_771_955_185_1727e-3,
                  ustrip(big(h) * big(c_0) / (_λx(BigFloat) * big(k_B))), m * K,
                  measurement(h) * measurement(c_0) / (_λx0 * measurement(k_B)),
                  measurement(BigFloat, h) * measurement(BigFloat, c_0) / (_λx(BigFloat) * measurement(BigFloat, k_B)),
                  "CODATA 2018")
@derived_constant(WienFrequencyDisplacementLawConstant, b′, "Wien frequency displacement law constant",
                  5.878_925_757_646_825e10,
                  ustrip(_νx(BigFloat) * big(k_B) / big(h)), Hz / K,
                  _νx0 * measurement(k_B) / measurement(h),
                  _νx(BigFloat) * measurement(BigFloat, k_B) / measurement(BigFloat, h), "CODATA 2018")
@constant(CharacteristicImpedanceOfVacuum, Z_0, "Characteristic impedance of vacuum",
          376.730_313_668, BigFloat(376_730_313_668)/BigFloat(1000_000_000), Ω, 5.7e-8,
          BigFloat(57) / BigFloat(1000_000_000), "CODATA 2018")

end # module CODATA2018

max_issues_repo_name: JuliaPhysics/Constants

max_forks_repo_name: JuliaPhysics/Constants


text:
using Test, Algorithms, Documenter
DocMeta.setdocmeta!(Algorithms, :DocTestSetup, :(using Algorithms); recursive=true)
doctest(Algorithms)
#@test Maths.abs_value(5)==10