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

using Markdown
using InteractiveUtils

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

# ╔═╡ 8d7c0eb0-84ca-11eb-0018-bdb2e389ee52
begin
	using PlutoUI
	using Distributions
	using Statistics
	using Plots
	using Random
	using JLD
	using Images
	using LsqFit
	md"_hidden code cell_"
end

# ╔═╡ e8d31f7c-88aa-11eb-05b7-31ff25d50afb
html"""
Term Paper 3

This work is a submission for the course Modelling Complex Systems(IDC-631) for the year 2020-21 Spring (Even) Semester by Dhruva Sambrani, MS18163
"""

# ╔═╡ 0d6813c6-88ae-11eb-3f39-b59ffd71ffd5
md"# Abstract

This paper mainly attempts to recreate Kuramoto's model of coupled oscillators as defined in [^strogatz].
This paper consists of four sections-

- §1. Describe and implement the model as in the paper
- §2. Check that the model evolves as we expect it to and make pretty animations
- §3. Show that simulated results agree with certain analytical properties
- §4. Study the robustness of the model to different distributions

"

# ╔═╡ 831e6f04-8a55-11eb-35da-2f77e5dd68a5
md"## How to read the paper

This paper is made using the Pluto[^pluto] Package which provides an interface which can show both code and theory side by side, which will hopefully enhance the understanding of the work and also provide a reproducable implementation without having to go into the details of it. Those who wish to read the paper for the theory and results can easily skip over the code and still understand everything. Those who are attempting to recreate the work can use the code as a guide. However, in the interest of keeping the code aligned to the topic at hand, certain auxillary functions have been hidden. The content of these hidden cells can be found in the executable version of this notebook [here](https://github.com/DhruvaSambrani/idc621/blob/main/TermPaper3.jl)
"

# ╔═╡ bd9326c2-8954-11eb-003e-b52f7107f199
TableOfContents()

# ╔═╡ 3effc012-8a56-11eb-1cd2-475e917c969b
md"# §0 Packages used

The Plots.jl[^plots] package is used for plotting and the Statistics[^stats] Package for some statistical analysis. Distributions.jl[^distributions] is used to sample the . PlutoUI.jl[^plutoui] is used for some UI elements such as the ToC. JLD[^jld] is used for storing data between runs. Images.jl[^images] is used for Gaussian smoothing. LsqFit.jl[^lsqfit] is used to fit the $\kappa_c$ vs $\sigma$ function."

# ╔═╡ 27e115ae-88b0-11eb-1d7d-57ae925c4c9f
md"# §1 - The Kuramoto Model

This is a model which attempts to simulate coupled oscillators. Consider a set of $n$ oscillators $A_{i\ |\ 1\le i\le n}$. Each $A_i$ has an angular position and internal angular velocity $\theta_i$ and $\omega_i$. Then, the most general coupled oscillator equation would be

$$\dot{\theta}_j = \omega_j + f(\theta_i\dots)_{1\le i\le n}$$

If the coupling is weak, and symmetric, this can be reduced to

$$\dot{\theta}_j = \omega_j + \kappa \sum_{i=1}^{n} \sin(\theta_j-\theta_i)$$

Which can be recast into the following equation

$$\dot{\theta}_j = \omega_j + \kappa r \sin(\Psi-\theta_i)$$
$$\text{where}$$
$$re^{\iota\Psi} = \text{mean}(e^{\iota\theta_i}\dots)_{1\le i \le n}$$
"

# ╔═╡ bc07ce24-88b5-11eb-2052-53543e8ef9fc
md"## Agent

Let us define a `struct`ure called Agent which holds two variables R and V, which correspond to $\theta$ and $\omega$ respectively.
"

# ╔═╡ a34a7c40-84ca-11eb-3433-41615dfc1c08
mutable struct Agent
	R::Float64
	V::Float64
end

# ╔═╡ 56ce942c-88b8-11eb-08fc-0970f8de6359
md"Now let's define how an agent updates, given κ, r, Ψ and dt, as defined in the equations before"

# ╔═╡ c4b72810-84ca-11eb-1855-23703a18b3e4
function update!(agent::Agent, κr::Float64, Phi::Float64, dt::Float64)
	agent.R += (agent.V + κr * sin(Phi - agent.R)) * dt
	# θ = θ + dθ/dt * Δt = θ + (ω + κr⋅sin(Ψ-θ))⋅dt
end

# ╔═╡ 9a1dec4e-88e7-11eb-0736-6b6e6b24627f
md"## Helper Functions"

# ╔═╡ c12bde70-84ca-11eb-1fe2-634e16d69909
begin
	Statistics.mean(arr::Array{Agent})::Complex = mean(map(arr) do ag; cis(ag.R); end)
	init(n::Int, σ::Float64) = Agent.(1., rand(Normal(0., σ), n))
end

# ╔═╡ 9fa5bf26-89b8-11eb-0636-ed01bb9b689b
init(n::Int, d::Distribution) = Agent.(1., rand(d, n))

# ╔═╡ a7d91872-88e7-11eb-2799-3502fc01e720
md"## Evolution and Simulation

As always, we define a single step evolution along with a multistep simulate function.

The evolve function takes a list of agents and updates them according to the dynamics explained above.
The simulate function sets the parameters and records certain variables for plotting and further analysis
"

# ╔═╡ cb3e73a0-84ca-11eb-3c68-8fc3fc5cf239
function evolve!(agents::Array{Agent}, κ::Float64, dt::Float64)
	meanvec = mean(agents)
	update!.(agents, κ*abs(meanvec), angle(meanvec), dt)
	meanvec
end

# ╔═╡ d00b64fe-84ca-11eb-1e68-f73bbca9e373
begin
	function simulate(steps::Int, n::Int, σ::Float64, κ::Float64, time::Int64)
		dt = time/steps
		agents = init(n, σ)
		means = Array{Complex,1}(undef, Int(steps))
		for step in 1:steps
			means[step] = evolve!(agents, κ, dt)
		end
		means
	end
	function simulate(steps::Int, agents::Array{Agent,1}, κ::Float64, time::Int64)
		dt = time/steps
		agents = deepcopy(agents)
		means = Array{Complex,1}(undef, Int(steps))
		for step in 1:steps
			means[step] = evolve!(agents, κ, dt)
		end
		means
	end
end

# ╔═╡ 393ed800-88e9-11eb-2160-edef4c90ee74
md"# §2 Visualizing the model

Let us first make sure that the model is running correctly. We can either visualize the positions of the oscillators with time, or we can plot $r$ vs time and make sure that the results match what is given in the paper. We will do both.
"

# ╔═╡ 86408954-88f9-11eb-2ee7-85532daf659d
md"## $r_\infty$ depends on κ"

# ╔═╡ 8c9b5c8e-84d3-11eb-0f9d-47c844a46a68
begin
	no_corr_means = simulate(100_000, 1000, 3., 1., 40)
	corr_means = simulate(100_000, 1000, 3., 7., 40)
end;

# ╔═╡ c618ba40-84d2-11eb-3990-93a786c8f8a9
begin
	plot(abs.(corr_means), ylims=(0,1), label="κ=7")
	plot!(abs.(no_corr_means), ylims=(0,1), label="κ=1")
	plot!(legend=:right, xlabel="Timestep", ylabel="r", title="r vs t")
end

# ╔═╡ cbb28294-88f9-11eb-1b77-fdacf61b8c9b
md"## Animated visualization of the dynamics

For this, we define a `simulate2` function which is of the same form as the previous function, but it also records the agents for each time step.
Then we plot the locations of the oscillators and the mean vector for each time step and make an animation.
Due to computational limitations, only 100 particles for 10000 timesteps are taken, and animation is made for every 20th frame. Also, for reproducability, the random seed is fixed. The seed was selected so that the results were clearer to show. However all observations are observed for any number of unseeded runs. To test, simply comment out the seed line.
"

# ╔═╡ f564ad60-88ea-11eb-31c5-a9f36ba1c93b
function simulate2(steps::Int, n::Int, σ::Float64, κ::Float64, time::Int64)
	dt = time/steps
	agents = init(n, σ)
	means = Array{Complex,1}(undef, Int(steps))
	recd = Array{Array{Agent, 1}, 1}(undef, Int(steps))
	for step in 1:steps
		means[step] = evolve!(agents, κ, dt)
		recd[step] = deepcopy(agents)
	end
	means, recd
end

# ╔═╡ bef4689a-88ec-11eb-2297-3bfe711d504b
begin
	Random.seed!(1901)
	ms, ags = simulate2(10000, 100, 3., 7., 40)
end;

# ╔═╡ 2a5fcef8-88ed-11eb-36e9-a743d38881f6
begin
	cmpl = ags .|> a->map(a) do ag
		cis(ag.R)
	end
	anim = @animate for i in 1:20:10000
		plot(real.(cmpl[i]), imag.(cmpl[i]), st=:scatter, xlims=(-1,1), ylims=(-1,1), aspect_ratio=1)
		plot!([0,real(ms[i])], [0, imag(ms[i])], width=2, label="Mean Vector")
	end
	md"_hidden code cell_"
end

# ╔═╡ 586e0782-88f3-11eb-2f30-03211c7a70a2
gif(anim, fps=3)

# ╔═╡ f2a5c968-88fb-11eb-3a39-2b95d2a8d978
begin
	plot(abs.(ms), ylim=(0,1), label="r", title="r vs Timestep", xlabel="Timestep", ylabel="r")
end

# ╔═╡ 222d4b00-88f4-11eb-2f5b-bd320a59b2ab
begin
	rs = ags .|> a->map(a) do ag; ag.R; end
	vs = (rs[2:end] .- rs[1:end-1]) .* 250
	md"_hidden code cell_"
end

# ╔═╡ 9537f2f8-88f4-11eb-091d-b5827649cb40
plot(1:10:10000, hcat(vs...)'[1:10:10000, :], legend=false, xaxis="Timestep", yaxis="Velocity", title="Effective Velocity vs time")

# ╔═╡ def3ccb4-8949-11eb-1ab1-a3403dc946ad
md"# §3 Simulations agree with analysis

In section 3.3 of [^strogatz], the author mentions the following observations of the simulations."

# ╔═╡ d2d048a8-8956-11eb-0536-972068f6862f
md"## r∞ vs κ

- For all $K$ less than a certain threshold $K_c$, the oscillators act as if they were uncoupled: the phases become uniformly distributed around the circle, starting from any initial condition. Then $r(t)$ decays to a tiny jitter of size $\mathcal O(N^{−1/2})$
- When $K > K_c$, this incoherent state becomes unstable and $r(t)$ grows exponentially, reflecting the nucleation of a small cluster of oscillators that are mutually synchronized, thereby generating a collective oscillation. Eventually $r(t)$ saturates at some level $r_\infty < 1$, though still with $\mathcal O(N^{−1/2})$ fluctuations.

Both of these claims are shown to be true in the first plot. For a \"high\" $K(=7)$, $r$ settles with jitter at 1, and for a \"low\" $K(=1)$, $r$ settles with jitter at 0.
"

# ╔═╡ 963c13ba-8957-11eb-3d45-f1d14f7b1c87
md"
## Individual Oscillators

- The population splits into two groups: the oscillators near the center of the frequency distribution lock together at the mean frequency and co-rotate with the average phase ψ(t), those in the tails run near their natural frequencies and drift relative to the synchronized cluster.
- This mixed state is often called partially synchronized. With further increases in K, more and more oscillators are recruited into the synchronized cluster, and $r_\infty$ tends to 1

Both of these claims are also shown to be true in the animation and in the plots of the effective velocities. In the animation, we clearly see that most of the oscillators are in a clump, but there are a few that drift away from the mean and stay out of the cluster. In the velcity plot, this is almost clearly evident. Also, if looked at carefully, it shows that the particles which are not in sync are the ones which where far away from the mean in the initial distribution"

# ╔═╡ afbd4280-8957-11eb-12a7-6762ab68086d
md"## κc vs g(ω)

The paper also analyses $\kappa_c$ as a function of $g(\omega)$. While Kuramoto and Nishikawa had previously assumed that the drifters would not have any effect on the coherence of the oscillators, they later recanted their stance and assumed that the drifters **did** in fact have an effect on the evolution of the coherence. Probably if they had the luxury of fast computers and easy computer languages, they would have seen that the drifters did in fact have an effect on $r(t)$. In the animation above, there are small but fast changes in the angle and the amplitude of the mean vector, and these are directly related to the movement of the drifters. Obviously, a simulation of only 100 oscillators is obviously not in the thermodynamic limit and the effect on $r$ is more apparent.

None the less, we can be reasonably sure that Kuramoto and Nishikawa's second attempt (§6.2 in [^strogatz]) was more accurate, and by extension so was Strogatz's own analysis. From their analysis,

$$\kappa_c = \frac{2}{\pi g(0)}$$

In the case of the gaussian distribution,

$$g(0) = \frac{1}{\sigma \sqrt{2\pi}}$$

Hence,

$$\kappa_c = \frac{2\sigma\sqrt{2\pi}}{\pi} = \sqrt\frac{8}{\pi}\sigma$$
"

# ╔═╡ ff877638-8967-11eb-3c88-7b20d5b80377
md"### Computational implementation

We define a function `simulate_for_sigma` which searches by bisection for $\kappa$s which are useful for our plot. We limit the difference of $\kappa$ and difference in $r_{\infty}$ for the sake of computational time. $r_{\infty}$ is evaluated by taking a mean over the last $10\%$ timesteps.

In the ideal case, we'd average $r_\infty$ over many initial conditions to make sure that there are infact no artifacts of the initial conditions. This is especially true, near the bifurcation, where the transition times can be very long. However, we take the easier path by simply box-smoothing over the data (a tool often employed by statisticians when they get noisy data).
"

# ╔═╡ 441a57e4-889a-11eb-27b8-4bdcc4bc2826
get_rinf(arr::Array{Complex}) = mean(abs.(last(arr, 10_000)))

# ╔═╡ 1d3dcba2-8894-11eb-1429-3d90997bbc68
function simulate_for_sigma(σ::Float64, κmn::Float64, κmx::Float64, κdiff::Float64, rdiff::Float64)
	function inner_function!(d::Dict{Float64, Float64}, κmn::Float64, κmx::Float64)
		if ((κmx - κmn) < κdiff) || ((d[κmx] - d[κmn]) < rdiff)
			return
		end
		κ = (κmx + κmn)/2
		d[κ] = get_rinf(simulate(100_000, 1000, σ, κ, 100))
		inner_function!(d, κ, κmx)
		inner_function!(d, κmn, κ)
	end
	d = Dict{Float64, Float64}()
	d[κmn] = get_rinf(simulate(100_000, 1000, σ, κmn, 100))
	d[κmx] = get_rinf(simulate(100_000, 1000, σ, κmx, 100))
	inner_function!(d, κmn, κmx)
	return d
end

# ╔═╡ 8937d90e-88a0-11eb-051a-1564ca42a34a
boxsmooth(arr, b) = [mean(arr[i-b:i+b]) for i in b+1:length(arr)-b]

# ╔═╡ 42f81e66-8a78-11eb-3ee9-cbfa2e4a0cce
md"Run new simulation? $(@bind runsim CheckBox(false))"

# ╔═╡ 0778d0b8-8962-11eb-3d45-ef790f6f1701
begin
	ks = 0
	rinfs = 0
	if runsim
		d = simulate_for_sigma(3, 0., 10., 0.01, 0.01)
		ks = sort(keys(d))
		rinfs = map(ks) do k; d[k]; end
	else
		d = load("save_sigma_3_nice.jld")
		ks, rinfs = d["ks"], d["rinfs"]
	end
	plot(boxsmooth(ks, 3), boxsmooth(rinfs, 3), label="Data")
	vline!([√(8/π)*3], label="Theoretical prediction")
	vline!([1.56*3], label="\\kappa_{c} from fitting below")
	plot!(xlabel="\\kappa", ylabel="\\ r_∞", title="\\kappa vs r_∞", legend=:right)
end

# ╔═╡ 206c24ea-896b-11eb-1391-89083f16f768
md"### Stability of κc over σ

Using the powerful computers at our disposal, we can run the simulation on multiple $\kappa$s and $\sigma$s and find the $r_\infty$ for each of these. Then, we can analyse the results to find $\kappa_c$ for each $\sigma$ and fit it to a straight line and see if the slope is equal to $\sqrt{\frac{8}{\pi}}$. Like before, we smooth our data, but this time we convolve the data with a 2dimensional Gaussian blur filter.

We measure the points of phase transition by taking a derivative along $\kappa$, and picking points where the derivative is above a certain threshold(=0.09).
"

# ╔═╡ b6ce7d86-8a77-11eb-3b8e-b36f29d2839b
md"""
#### Code to generate data
```julia
function graph(κs::Array{Float64}, σs::Array{Float64})
	arr = zeros(Float64, length(κs), length(σs))
	for (i, κ) in enumerate(κs), (j, σ) in enumerate(σs)
		arr[i,j] = mean(
				abs,
				@view simulate(Int(1e5), 1000, σ, κ, 1e-4)[end - 1000:end]
		)
	end
	arr
end
κs, σs = collect(0.0:0.1:2.5), collect(0.0:0.1:2.5)
zs = graph(κs, σs)
save("save.jld", "zs", zs, "ks", κs, "ss", σs)
```
"""

# ╔═╡ a028f8a0-8972-11eb-05dc-3103777aa93a
begin
	kcs_data = load("save.jld")
	kcszs, kcsκs, kcsσs = kcs_data["zs"], kcs_data["ks"], kcs_data["ss"]
	kcsnzs = imfilter(kcszs, Kernel.gaussian(0.7))
	deriv = kcsnzs[2:end, :] - kcsnzs[1:end-1, :]
	j = findall(deriv .> 0.09)

	ys = map(j) do i; kcsκs[i[1]]; end
	xs = map(j) do i; kcsσs[i[2]]; end

	@. model(x, m) = m*x
	fit = curve_fit(model, xs, ys, [0.])
end;

# ╔═╡ d61d1470-8989-11eb-1ddb-cff31dd6d747
md"_Interactive plot below_"

# ╔═╡ 4ad38dd0-896c-11eb-2da5-335c39437209
let
	plotlyjs()
	p = plot(
		kcsκs, kcsσs, kcsnzs',
		st=:surface, zlims=(0, 1),
		xlabel="κ" , ylabel="σ", zlabel="r∞",
		title="\$R_∞(\\sigma, \\kappa)\$", extra_plot_kwargs = KW(:include_mathjax => "cdn")
	)
	gr()
	p
end

# ╔═╡ e29ac97e-8989-11eb-372b-3be6756949d8
md"_Interactive plot below_"

# ╔═╡ 1f26e4dc-8973-11eb-3428-ff96adccd9d6
let
	plotlyjs()
	p = plot(
		kcsκs, kcsσs, deriv',
		zlims=(0,0.5), st=:surface,
		xlabel="κ" , ylabel="σ", zlabel="r∞",
		title="\$\\dfrac{dr_\\infty}{d\\kappa}(\\sigma, \\kappa)\$", extra_plot_kwargs = KW(:include_mathjax => "cdn")
	)
	gr()
	p
end

# ╔═╡ 97c0c8a4-8973-11eb-31c1-a3e993ec94b6
begin
	plot(0:0.5:2.5, model(0:0.5:2.5, fit.param), label="fit",
		yaxis="κ", xaxis="σ",
		xlims=(0,2.5), ylims=(0, 2.5))
	plot!(xs, ys, st=:scatter, label = "data",
		legend=:bottomright, title="\\kappa_{c} vs \\sigma")
end

# ╔═╡ b9a7c122-8980-11eb-1dba-a790aec65322
Markdown.parse("\$\$\\kappa_c = $(round(fit.param[1], sigdigits=3))\\cdot\\sigma\$\$")

# ╔═╡ 6b896608-898f-11eb-2c33-ff41bc3b1c87
md"### Observations

The single run and the extended analysis both show the same result that the $\kappa_c$ follows the analysis in the paper closely, however, there is a small disparity. Any finite simulation will inevitably not have the sharp nature of the bifurcation. In §4, we run the same simulations for other distributions, and see the effect of the kind of distribution on this.
"

# ╔═╡ e0a1ebd2-89a9-11eb-15b0-f5e35ecd8ffc
md"# §4 Robustness of results to distributions

As mentioned in the previous subsection, let's run the same simulation with some other distributions"

# ╔═╡ a64c11a4-89b1-11eb-2050-f5267b9c57c1
kappas = 0:0.5:10

# ╔═╡ a1d5c076-89ad-11eb-1530-cbd02ee0a46f
md"## Gaussian Distribution

We've already run this before, but let's do it again.

$$g(\omega) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{\omega^2}{2\sigma^2}\right];\  g(0) = \frac{1}{\sqrt{2\pi\sigma^2}};\ \kappa_c(3) = \sqrt\frac{8}{\pi}*3 = 1.596$$
"

# ╔═╡ e88b4f7a-89af-11eb-1b9f-c9cfc4510d86
gaussian_rs = map(kappas) do kappa
	mean(abs.(last(simulate(100_000, 1000, 3., kappa, 40), 1000)))
end;

# ╔═╡ 78e8340a-89b7-11eb-219c-135ab851ccc2
begin
	plot(kappas, gaussian_rs, xlabel="\\kappa", ylabel="\\ r_∞")
	vline!([sqrt(8/pi)*3], title="Gaussian")
end

# ╔═╡ 12d46c62-89b7-11eb-1d5e-1717a0643ad1
md"## Lorenz-Cauchy distributions

$$g(\omega) = \frac{1}{\pi \sigma \left(1 + \left(\frac{\omega}{\sigma} \right)^2 \right)};\ g(0) = \frac{1}{\pi\sigma};\ \kappa_c(3) = 6$$
"

# ╔═╡ f3e77840-89b8-11eb-0d5a-b3501386bb19
begin
	cauchy_kappas = 3:0.5:12
	cauchy_rs_1 = map(cauchy_kappas) do kappa
		mean(abs.(last(simulate(50_000, init(1000, Cauchy(0,3)), kappa, 40), 1000)))
	end;
	cauchy_rs_2 = map(cauchy_kappas) do kappa
		mean(abs.(last(simulate(50_000, init(1000, Cauchy(0,3)), kappa, 40), 1000)))
	end;
	cauchy_rs_3 = map(cauchy_kappas) do kappa
		mean(abs.(last(simulate(50_000, init(1000, Cauchy(0,3)), kappa, 40), 1000)))
	end;
	cauchy_rs = (cauchy_rs_1 + cauchy_rs_2 + cauchy_rs_3)/3
end

# ╔═╡ ef8ad520-89b9-11eb-126a-db54fc2e65e9
begin
	plot(cauchy_kappas, cauchy_rs, xlabel="\\kappa", ylabel="\\ r_∞")
	vline!([6], title="Lorenz-Cauchy", ylim=(0,1))
end

# ╔═╡ 7df01b22-89c4-11eb-1422-e5aa9bd7eac6
md"## Reflected exponential decay

$$g(\omega) = \frac{1}{2a}\exp\left(-\frac{|x|}{a}\right);\ g(0) = \frac{1}{2a};\ \kappa_c(3) = \frac{4\pi}{3} = 4.189$$
"

# ╔═╡ cbdf714a-8a44-11eb-0f4d-d30ab28b7ab7
begin
	function rejectionSample(f, N)
		i, x_list = 1, zeros(N)
		while(i < N + 1)
			x=20*rand() - 10
			y=rand()
			if(f(x)>y)
				x_list[i] = x
				i += 1
			end
		end
		return x_list
	end
	red(a)= x -> exp(-abs(x)/a)/2a
end

# ╔═╡ 123f09de-8a45-11eb-13e5-fdac2199b036
red_ags = Agent.(0.0, rejectionSample(red(3), 1000));

# ╔═╡ cd5c8b66-8a49-11eb-3cf5-75946228434e
begin
	red_ags_2 = deepcopy(red_ags)
	red_rs = map(kappas) do kappa
		mean(abs.(last(simulate(100_000, red_ags_2, kappa, 40), 1000)))
	end
end;

# ╔═╡ 23872ebe-8a4b-11eb-0e7d-2d7500b874fb
begin
	plot(kappas, red_rs, xlabel="\\kappa", ylabel="\\ r_∞")
	vline!([4pi/3], title="Reflected Exponential Decay")
end

# ╔═╡ ca49a058-8a4e-11eb-1445-417745db689d
md"## Observations

Like the Gaussian, the other distributions also have a slight discrepancy from the expected value. However by a visual analysis, we can say that the theoretical results are close to what we would get in the infinite limit. For a more accurate analysis, one could fit the function as mentioned in the paper to the data, and then look at the value of $\kappa_c$.

Another observation that we can see is that distributions that are tail heavy are much more prone to variations in the partially coherent part of the graph. In order to reduce this, we average results over 3 runs with random initial conditions for the Lorenz Cauchy distribution
"

# ╔═╡ 30c5c5d4-8a57-11eb-248c-a7961f0d4590
md"# References
[^strogatz]: Steven H. Strogatz, From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators, Physica D: Nonlinear Phenomena, Volume 143, Issues 1–4, [https://doi.org/10.1016/S0167-2789(00)00094-4](https://doi.org/10.1016/S0167-2789(00)00094-4.)
[^julia]: Julia: A Fresh Approach to Numerical Computing. Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah. (2017) SIAM Review, 59: 65–98. doi: [10.1137/141000671](https://dx.doi.org/10.1137/141000671)
[^plots]: Plots - powerful convenience for visualization in Julia. Thomas Breloff. [docs.juliaplots.org](docs.juliaplots.org)
[^pluto]: Pluto - Reactive Notebooks for Julia. Fons van der Plas and Mikołaj Bochenski. [https://github.com/fonsp/Pluto.jl](https://github.com/fonsp/Pluto.jl)
[^plutoui]: PlutoUI. Fons van der Plas. [https://github.com/fonsp/PlutoUI.jl](https://github.com/fonsp/PlutoUI.jl)
[^stats]: Statistics - Julia stdlib for Statistics. [https://docs.julialang.org/en/v1/stdlib/Statistics/](https://docs.julialang.org/en/v1/stdlib/Statistics/)
[^images]: JuliaImages: image processing and machine vision for Julia. [https://juliaimages.org/stable/](https://juliaimages.org/stable/)
[^lsqfit]: LsqFit.jl: Basic least-squares fitting in pure Julia under an MIT license. [https://github.com/JuliaNLSolvers/LsqFit.jl](https://github.com/JuliaNLSolvers/LsqFit.jl)
[^distributions]: 1.Dahua Lin. JuliaStats/Distributions.jl: v0.24.15. (2021). [doi:10.5281/zenodo.4577319](https://doi.org/10.5281/zenodo.2647458)
[^jld]: JLD.jl. JuliaIO. [https://github.com/JuliaIO/JLD.jl](https://github.com/JuliaIO/JLD.jl)
"

# ╔═╡ Cell order:
# ╟─e8d31f7c-88aa-11eb-05b7-31ff25d50afb
# ╟─0d6813c6-88ae-11eb-3f39-b59ffd71ffd5
# ╟─831e6f04-8a55-11eb-35da-2f77e5dd68a5
# ╟─bd9326c2-8954-11eb-003e-b52f7107f199
# ╟─3effc012-8a56-11eb-1cd2-475e917c969b
# ╟─8d7c0eb0-84ca-11eb-0018-bdb2e389ee52
# ╟─27e115ae-88b0-11eb-1d7d-57ae925c4c9f
# ╟─bc07ce24-88b5-11eb-2052-53543e8ef9fc
# ╠═a34a7c40-84ca-11eb-3433-41615dfc1c08
# ╟─56ce942c-88b8-11eb-08fc-0970f8de6359
# ╠═c4b72810-84ca-11eb-1855-23703a18b3e4
# ╟─9a1dec4e-88e7-11eb-0736-6b6e6b24627f
# ╠═c12bde70-84ca-11eb-1fe2-634e16d69909
# ╠═9fa5bf26-89b8-11eb-0636-ed01bb9b689b
# ╟─a7d91872-88e7-11eb-2799-3502fc01e720
# ╠═cb3e73a0-84ca-11eb-3c68-8fc3fc5cf239
# ╠═d00b64fe-84ca-11eb-1e68-f73bbca9e373
# ╟─393ed800-88e9-11eb-2160-edef4c90ee74
# ╟─86408954-88f9-11eb-2ee7-85532daf659d
# ╠═8c9b5c8e-84d3-11eb-0f9d-47c844a46a68
# ╟─c618ba40-84d2-11eb-3990-93a786c8f8a9
# ╟─cbb28294-88f9-11eb-1b77-fdacf61b8c9b
# ╟─f564ad60-88ea-11eb-31c5-a9f36ba1c93b
# ╠═bef4689a-88ec-11eb-2297-3bfe711d504b
# ╟─2a5fcef8-88ed-11eb-36e9-a743d38881f6
# ╟─586e0782-88f3-11eb-2f30-03211c7a70a2
# ╟─f2a5c968-88fb-11eb-3a39-2b95d2a8d978
# ╟─222d4b00-88f4-11eb-2f5b-bd320a59b2ab
# ╟─9537f2f8-88f4-11eb-091d-b5827649cb40
# ╟─def3ccb4-8949-11eb-1ab1-a3403dc946ad
# ╟─d2d048a8-8956-11eb-0536-972068f6862f
# ╟─963c13ba-8957-11eb-3d45-f1d14f7b1c87
# ╟─afbd4280-8957-11eb-12a7-6762ab68086d
# ╟─ff877638-8967-11eb-3c88-7b20d5b80377
# ╠═441a57e4-889a-11eb-27b8-4bdcc4bc2826
# ╠═1d3dcba2-8894-11eb-1429-3d90997bbc68
# ╠═8937d90e-88a0-11eb-051a-1564ca42a34a
# ╟─42f81e66-8a78-11eb-3ee9-cbfa2e4a0cce
# ╠═0778d0b8-8962-11eb-3d45-ef790f6f1701
# ╟─206c24ea-896b-11eb-1391-89083f16f768
# ╟─b6ce7d86-8a77-11eb-3b8e-b36f29d2839b
# ╠═a028f8a0-8972-11eb-05dc-3103777aa93a
# ╟─d61d1470-8989-11eb-1ddb-cff31dd6d747
# ╟─4ad38dd0-896c-11eb-2da5-335c39437209
# ╟─e29ac97e-8989-11eb-372b-3be6756949d8
# ╟─1f26e4dc-8973-11eb-3428-ff96adccd9d6
# ╟─97c0c8a4-8973-11eb-31c1-a3e993ec94b6
# ╟─b9a7c122-8980-11eb-1dba-a790aec65322
# ╟─6b896608-898f-11eb-2c33-ff41bc3b1c87
# ╟─e0a1ebd2-89a9-11eb-15b0-f5e35ecd8ffc
# ╠═a64c11a4-89b1-11eb-2050-f5267b9c57c1
# ╟─a1d5c076-89ad-11eb-1530-cbd02ee0a46f
# ╠═e88b4f7a-89af-11eb-1b9f-c9cfc4510d86
# ╟─78e8340a-89b7-11eb-219c-135ab851ccc2
# ╟─12d46c62-89b7-11eb-1d5e-1717a0643ad1
# ╠═f3e77840-89b8-11eb-0d5a-b3501386bb19
# ╟─ef8ad520-89b9-11eb-126a-db54fc2e65e9
# ╟─7df01b22-89c4-11eb-1422-e5aa9bd7eac6
# ╟─cbdf714a-8a44-11eb-0f4d-d30ab28b7ab7
# ╠═123f09de-8a45-11eb-13e5-fdac2199b036
# ╠═cd5c8b66-8a49-11eb-3cf5-75946228434e
# ╟─23872ebe-8a4b-11eb-0e7d-2d7500b874fb
# ╟─ca49a058-8a4e-11eb-1445-417745db689d
# ╟─30c5c5d4-8a57-11eb-248c-a7961f0d4590


text:
# this file was generated by this code
# using SnoopCompile
# inf_timing = @snoopi tmin=0.01 include("test/runtests.jl")
# pc = SnoopCompile.parcel(inf_timing)
# SnoopCompile.write("/tmp/precompile", pc)
# I have no idea what's going on here...

# Use
#    @warnpcfail precompile(args...)
# if you want to be warned when a precompile directive fails
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(circlepath)),NamedTuple{(:reversepath,), Tuple{Bool}},typeof(circlepath),Point,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(crescent)),NamedTuple{(:vertices, :reversepath), Tuple{Bool, Bool}},typeof(crescent),Point,Int64,Point,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(crescent)),NamedTuple{(:vertices, :reversepath, :steps), Tuple{Bool, Bool, Int64}},typeof(crescent),Point,Int64,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(crescent)),NamedTuple{(:vertices,), Tuple{Bool}},typeof(crescent),Point,Float64,Float64})
    Base.precompile(Tuple{Core.kwftype(typeof(ellipse)),NamedTuple{(:stepvalue, :vertices, :reversepath), Tuple{Float64, Bool, Bool}},typeof(ellipse),Point,Point,Point,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(epitrochoid)),NamedTuple{(:action, :close, :stepby, :vertices), Tuple{Symbol, Bool, Float64, Bool}},typeof(epitrochoid),Int64,Int64,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(epitrochoid)),NamedTuple{(:close, :stepby, :vertices), Tuple{Bool, Float64, Bool}},typeof(epitrochoid),Int64,Int64,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(hypotrochoid)),NamedTuple{(:close, :stepby, :vertices), Tuple{Bool, Float64, Bool}},typeof(hypotrochoid),Int64,Int64,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(ngon)),NamedTuple{(:action,), Tuple{Symbol}},typeof(ngon),Int64,Int64,Int64,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(ngon)),NamedTuple{(:vertices, :reversepath), Tuple{Bool, Bool}},typeof(ngon),Point,Int64,Int64,Float64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(ngon)),NamedTuple{(:vertices,), Tuple{Bool}},typeof(ngon),Point,Int64,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(polycross)),NamedTuple{(:vertices, :reversepath), Tuple{Bool, Bool}},typeof(polycross),Point,Int64,Int64,Float64,Float64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(prettypoly)),NamedTuple{(:action, :close), Tuple{Symbol, Bool}},typeof(prettypoly),Vector{Point}})
    Base.precompile(Tuple{Core.kwftype(typeof(prettypoly)),NamedTuple{(:action,), Tuple{Symbol}},typeof(prettypoly),BoundingBox})
    Base.precompile(Tuple{Core.kwftype(typeof(prettypoly)),NamedTuple{(:close,), Tuple{Bool}},typeof(prettypoly),Vector{Point},Symbol,Function})
    Base.precompile(Tuple{Core.kwftype(typeof(prettypoly)),NamedTuple{(:close,), Tuple{Bool}},typeof(prettypoly),Vector{Point},Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(spiral)),NamedTuple{(:stepby, :period, :vertices, :log), Tuple{Float64, Float64, Bool, Bool}},typeof(spiral),Int64,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(squircle)),NamedTuple{(:rt, :vertices, :stepby, :reversepath), Tuple{Float64, Bool, Float64, Bool}},typeof(squircle),Point,Int64,Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(star)),NamedTuple{(:action, :vertices, :reversepath), Tuple{Symbol, Bool, Bool}},typeof(star),Int64,Int64,Int64,Int64,Float64,Float64})
    Base.precompile(Tuple{Core.kwftype(typeof(star)),NamedTuple{(:vertices, :reversepath), Tuple{Bool, Bool}},typeof(star),Point,Int64,Int64,Float64,Irrational{:π},Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(textoutlines)),NamedTuple{(:action, :halign, :valign, :startnewpath), Tuple{Symbol, Symbol, Symbol, Bool}},typeof(textoutlines),String,Point})
    Base.precompile(Tuple{Type{BoundingBox}})
    Base.precompile(Tuple{Type{Drawing}})
    Base.precompile(Tuple{typeof(arc2r),Point,Point,Point,Symbol})
    Base.precompile(Tuple{typeof(arc2sagitta),Point,Point,Int64,Symbol})
    Base.precompile(Tuple{typeof(box),Point,Int64,Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(center3pts),Point,Point,Point})
    Base.precompile(Tuple{typeof(ngon),Int64,Int64,Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(ngon),Point,Float64,Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(ngon),Point,Int64,Int64})
    Base.precompile(Tuple{typeof(sector),Float64,Float64,Float64,Float64,Int64,Symbol})
    Base.precompile(Tuple{typeof(sector),Point,Int64,Int64,Int64,Float64,Int64,Symbol})
    Base.precompile(Tuple{typeof(snapshot)})
    let fbody = try __lookup_kwbody__(which(arrow, (Point,Point,Vector{Int64},Symbol,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Int64,Int64,Float64,Bool,Bool,Bool,Float64,Function,Nothing,typeof(arrow),Point,Point,Vector{Int64},Symbol,))
        end
    end
end

ext: jl


text:
#winreg.jl: Register AudioStick as a Windows Explorer contextual menu.

const STRING_SYNCWITHAUDIOSTICK = "Synchronize with AudioStick"

#==
===============================================================================#
function generate_inf_string()
	infqt(s) = string("\"\"", s, "\"\"") #.inf escaped quote (2 are necessary)
	module_root = dirname(dirname(abspath(@__FILE__)))
	juliacmd = infqt(joinpath(Sys.BINDIR, Base.julia_exename()))
	path_audiostick = infqt(joinpath(module_root, "run.jl"))
	#Compensate for shell bug by adding " " before quote
	#(Interprets last "\" of "Z:\" as an escape char!)
	path_dflttgt = infqt(DFLT_TGTROOT * " ")

	inf_file = """
[version]
signature="\$CHICAGO\$"

[DefaultInstall]
AddReg = Explore.AddReg

[Strings]
ACTION_TEXT = "$STRING_SYNCWITHAUDIOSTICK"
ACTION_CMD = "$juliacmd --color=yes $path_audiostick ""%1"" $path_dflttgt"

[Explore.AddReg]
HKCR,SystemFileAssociations\\.m3u\\shell\\AudioStick,,,%ACTION_TEXT%
HKCR,SystemFileAssociations\\.m3u\\shell\\AudioStick\\command,,,%ACTION_CMD%
HKCR,SystemFileAssociations\\.m3u8\\shell\\AudioStick,,,%ACTION_TEXT%
HKCR,SystemFileAssociations\\.m3u8\\shell\\AudioStick\\command,,,%ACTION_CMD%
"""
end

#Create .inf file to register AudioStick as a Windows Explorer context menu
function writeinf(filepath::String)
	filepath = abspath(filepath)
	filename = basename(filepath)
	@info("Writing $filepath...")
	open(filepath, "w") do f
		write(f, generate_inf_string())
	end

	println("""\n
You can now install the AudioStick context menu from the auto-generated
.inf file:
   1) In Windows Explorer, right-click on:
         $filename
   2) Click "Install".

You will then be able to synchronize playlists from the windows Explorer:
   1) In Windows Explorer, right-click any .m3u/.m3u8 file.
	2) Click "$STRING_SYNCWITHAUDIOSTICK".
""")
end
writeinf() = writeinf("AudioStick_AddWinExContextMenu.inf")

#Last line

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick


text:
module CellMapApplyTests

using Test
using Gridap
using Gridap.CachedArrays
using TensorValues

using ..CellValuesMocks
using ..MapsMocks

l = 10

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)
r = evaluate(m,p)

cm = TestIterCellValue(m,l)
cp = TestIterCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]

test_iter_cell_map(cm,cp,rm)

cm2 = apply(-,cm,broadcast=true)
rm = [ CachedArray(-r) for i in 1:l]

test_iter_cell_map(cm2,cp,rm)

cm = TestIndexCellValue(m,l)
cp = TestIndexCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]

test_index_cell_map_with_index_arg(cm,cp,rm)

ca2 = evaluate(cm,cp)

test_index_cell_array(ca2,rm)

cm = TestIterCellValue(m,l)
cp = TestIterCellValue(p,l)
rm = [ CachedArray(r.-p) for i in 1:l]
cm2 = apply(-,cm,cp,broadcast=true)
test_iter_cell_map(cm2,cp,rm)

a = 2
cm = TestIterCellValue(m,l)
ca = TestIterCellValue(a,l)
rm = [ CachedArray(r.-a) for i in 1:l]
cm2 = apply(-,cm,ca,broadcast=true)
test_iter_cell_map(cm2,cp,rm)

cm = TestIndexCellValue(m,l)
cp = TestIndexCellValue(p,l)
rm = [ CachedArray(r.-p) for i in 1:l]
cm2 = apply(-,cm,cp,broadcast=true)
test_index_cell_map_with_index_arg(cm2,cp,rm)

a = 2
cm = TestIndexCellValue(m,l)
ca = TestIndexCellValue(a,l)
rm = [ CachedArray(r.-a) for i in 1:l]
cm2 = apply(-,cm,ca,broadcast=true)
test_index_cell_map_with_index_arg(cm2,cp,rm)

a = VectorValue(10.1,10.3)
m = MockMap(a)
r = evaluate(m,p)

cm = TestIterCellValue(m,l)
cp = TestIterCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]
test_iter_cell_map(cm,cp,rm)

cm2 = apply(-,cm,broadcast=true)
rm = [ CachedArray(-r) for i in 1:l]
test_iter_cell_map(cm2,cp,rm)

cm = TestIndexCellValue(m,l)
cp = TestIndexCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]
test_index_cell_map(cm,cp,rm)

cm2 = apply(-,cm,broadcast=true)
rm = [ CachedArray(-r) for i in 1:l]
test_index_cell_map_with_index_arg(cm2,cp,rm)

end # module


text:
"""
This submodule contains the HTTP interface to the AlphaVantage API.

The `AlphaVantageClient` struct owns the data needed for the HTTP requests.

A `Base.get` method takes as arguments the above client, a data type parameter that inherits from
AVFunction (see av_types.jl), and a Dict of query parameters. It then executes the HTTP request
as well as dispatches to appropriate pre/postprocessor functions (see av_handlers.jl).
"""

"""
HTTP client for connecting to AlphaVantage API.

```
AlphaVantageClient(;scheme, host, key)
```

Appropriate defaults are set for `scheme` and `host` fields. You must provide an API key via a named parameter
or by creating an ALPHA_VANTAGE_API_KEY environment variable prior to launching your Julia session. See
https://www.alphavantage.co/support/#api-key to obtain an API key. If a key is not provided, the default is
to use "demo" which may have restricted support.
"""
struct AlphaVantageClient
  scheme::String
  host::String
  key::String
end

"Retrieve API key from environment, user argument, or use the default demo key."
_api_key(k=get(ENV, ENV_VARIABLE_API_KEY, "")) = if k == "" DEFAULT_API_KEY else k end

AlphaVantageClient(;scheme=HTTP_SCHEME, host=HTTP_HOST, key=_api_key()) = AlphaVantageClient(scheme, host, key)
AlphaVantageClient(k::String) = AlphaVantageClient(;key=k)

"Dynamically creates and evaluates an HTTP.URI call for a given API request."
_create_uri(c::AlphaVantageClient, q::Dict{String}) = string(HTTP.URI(scheme=c.scheme, host=c.host, query=q))

# `Base.get` method dispatching
function Base.get(c::AlphaVantageClient, params)
  params["apikey"] = c.key
  req = HTTP.request("GET", _create_uri(c, params))

  if req.status != 200
      throw("HTTP request returned status code $(req.status)")
  end

  deepcopy(req.body)
end

function Base.get(c::AlphaVantageClient, f::A where A <: AVFunction, params::Dict{String})
  params["function"] = string(f)
  params = preprocess(f, params)
  resp = get(c, params)
  resp = postprocess(f, resp)
  resp
end

function Base.get(c::AlphaVantageClient, f::DataType, params::Dict{String})
  get(c, f(), params)
end


text:
using EE263
using Test

@testset "EE263.jl" begin
    # Write your tests here.
end


text:
web = "http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/b37.vcf/"
vcf = "1kg.phase3.v5a.vcf.gz"
tbi = "1kg.phase3.v5a.vcf.gz.tbi"

for chr in 1:22
    println("downloading chromosome $chr data")
    vcffile = web * "chr$chr." * vcf
    tbifile = web * "chr$chr." * tbi
    isfile(vcffile) || download(vcffile, "chr$chr." * vcf)
    isfile(tbifile) || download(tbifile, "chr$chr." * tbi)
end

vcffile = web * "chrX." * vcf
tbifile = web * "chrX." * tbi
isfile(vcffile) || download(vcffile, "chrX." * vcf)
isfile(tbifile) || download(tbifile, "chrX." * tbi)

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
function sr_step(peps, h)
    vars = linsolve(x->apply_smatrix(peps, x), -im .* fvec(peps, h), variables(peps), GMRES())
    load_variables!(copy(peps), vars)
end

# TODO: this needs to be fixed
function apply_smatrix(peps::PEPS, x)
    g = TensorAD.gradient(x->inner_product(conj(peps), x), peps)
    vg = variables(g)
    res = zero(vg)
    # replace target tensor by x_i
    for i=1:nsite(peps)
        qeqs = copy(peps)
        qeqs.tensors[i] = x.tensors[i]
        g_i = TensorAD.gradient(x->inner_product(conj(qeqs), x), peps)
        res .+= variables(g_i)
    end
    return res .- vg .* (vg' * x)
end

# im*L₂
function iloss2(h, peps, variables)
    pl = load_variables(peps, variables)
    pl = pl * inv.(norm(pl))
    pr = peps * inv.(norm(peps))
    return real(expect(h, conj(pl), pr))
end
function fvec(peps::PEPS, h)
    variables = TensorNetworkEvolve.variables(peps)
    return -im*TensorAD.gradient(x->iloss2(h, wrap_difftensor(peps), x), DiffTensor(variables))[1]
end
function wrap_difftensor(peps::PEPS)
    return replace_tensors(peps, DiffTensor.(alltensors(peps)))
end

# @non_differentiable OMEinsum.optimize_greedy(code, size_dict)
# @non_differentiable replace(vec, pairs...)

# # zygote patch
# using Zygote
# function Zygote.accum(x::Vector, y::Tuple)
#     Zygote.accum.(x, y)
# end


text:
"""
    CompatHelper(;
        file="$(contractuser(default_file("github", "workflows", "CompatHelper.yml")))",
        destination="CompatHelper.yml",
    )

Integrates your packages with [CompatHelper](https://github.com/bcbi/CompatHelper.jl) via GitHub Actions.

## Keyword Arguments
- `file::AbstractString`: Template file for the workflow file.
- `destination::AbstractString`: Destination of the workflow file,
  relative to `.github/workflows`.
"""
@with_kw_noshow struct CompatHelper <: BasicPlugin
    file::String = default_file("github", "workflows", "CompatHelper.yml")
    destination::String = "CompatHelper.yml"
end

source(p::CompatHelper) = p.file
destination(p::CompatHelper) = joinpath(".github", "workflows", p.destination)
tags(::CompatHelper) = "<<", ">>"

view(p::CompatHelper, t::Template, ::AbstractString) = Dict(
    "VERSION" => format_version(max(v"1.2", t.julia)),
)


text:
using EventSimulation
using Distributions
using Printf
using Random

# Objectives of the example:
# * show how SimQueue object can be used
# * handling interruptions of actions

mutable struct Customer
    name::String
    arrival_time::Float64
    renege_a::Action{Float64} # this is initially left #undef

    function Customer(name, arrival_time)
        x = new()
        x.name = name
        x.arrival_time = arrival_time
        return x
    end
end

mutable struct BankState <: AbstractState
    all_cust::Int # how many customers we want to process
    n::Int # how many customers were processed
    reneged_count::Int # number of customers that reneged
    arrival::Exponential{Float64}
    patience::Uniform{Float64}
    cashing::Exponential{Float64}
    counter::SimQueue{Customer}
    report::Bool
end

function source(s::Scheduler)
    s.state.n += 1
    c = Customer(string("Customer", s.state.n), s.now)
    s.state.report && @printf("%6.4f %s: Here I am\n", s.now, c.name)
    provide!(s, s.state.counter, c)
    # we remember the action to remove it from event_queue if needed
    c.renege_a = register!(s, x -> renege(x, c), rand(s.state.patience))
    if s.state.n < s.state.all_cust
        register!(s, source, rand(s.state.arrival))
    end
end

function renege(s, c)
    s.state.report && @printf("%6.4f %s: RENEGED after %6.4f\n",
                              s.now, c.name, s.now-c.arrival_time)
    s.state.reneged_count += 1
    # customer reneged and leaves the queue
    withdraw!(s.state.counter, c) || error("should be succesffull")
end

function serve(s, c)
    s.state.report && @printf("%6.4f %s: waited %6.4f\n",
                              s.now, c.name, s.now-c.arrival_time)
    # customer got to service so we interrupt renege action
    interrupt!(s, c.renege_a) || error("should be succesfull")
    register!(s, x -> finish(x, c), rand(s.state.cashing))
end

function finish(s, c)
    s.state.report && @printf("%6.4f %s: Finished\n", s.now, c.name)
    request!(s, s.state.counter, serve)
end

function run(ccount, report)
    bs = BankState(ccount, 0, 0, Exponential(10.0),
                   Uniform(1.0, 3.0), Exponential(12.0),
                   SimQueue{Customer}(), report)
    report && println("Bank renege")
    s = Scheduler(bs)
    request!(s, bs.counter, serve)
    source(s)
    go!(s, Inf)
    println("Renege probability: ", bs.reneged_count / bs.all_cust)
end

Random.seed!(1)
@time run(10, true)
@time run(1_000_000, false)


text:
# Various functions regarding graph connectivity

export components, num_components, is_connected, spanning_forest, random_spanning_forest
export find_path, dist, diam, is_cut_edge, is_acyclic, wiener_index
export eccentricity, radius, graph_center

"""
`components(G)` returns the vertex sets of the connected components of
`G` (as a `Partition`).
"""
function components(G::SimpleGraph{T}) where {T}
    if cache_check(G, :components)
        return cache_recall(G, :components)
    end

    P = Partition(G.V)
    for e in G.E
        (u, v) = e
        merge_parts!(P, u, v)
    end
    cache_save(G, :components, P)
    return P
end

"""
`num_components(G)` returns the number of connected components in `G`.
"""
function num_components(G::SimpleGraph{T})::Int where {T}
    if cache_check(G, :num_components)
        return cache_recall(G, :num_components)
    end
    result = num_parts(components(G))
    cache_save(G, :num_components, result)
    return result
end

# determine if the graph is connected
"""
`is_connected(G)` determines if `G` is connected.
"""
function is_connected(G::SimpleGraph{T}) where {T}
    return num_components(G) <= 1
end

# create a spanning forest of G, i.e., a maximal acyclic spanning
# subgraph. If G is connected, this is a tree.

"""
`spanning_forest(G)` creates a maximal acyclic subgraph of `G`.
"""
function spanning_forest(G::SimpleGraph{T}) where {T}
    if cache_check(G, :spanning_forest)
        return cache_recall(G, :spanning_forest)
    end
    H = SimpleGraph{T}()
    if NV(G) == 0
        return H
    end

    for v in G.V
        add!(H, v)
    end

    P = Partition(G.V)

    for e in G.E
        (u, v) = e
        if in_same_part(P, u, v)
            continue
        end
        add!(H, u, v)
        merge_parts!(P, u, v)
        if num_parts(P) == 1
            break
        end
    end
    cache_save(G, :spanning_forest, H)
    return H
end

"""
`random_spanning_forest(G)` returns a random spanning forest of the graph `G`.
Each component of the result spans a component of `G`.

This differs from `spanning_forest` in that repeated invocations of this function
can return different results.
"""
function random_spanning_forest(G::SimpleGraph)
    VT = eltype(G)
    ET = Tuple{VT,VT}

    T = SimpleGraph{VT}()  # output

    VV = vlist(G)
    EE = elist(G)

    for v in VV
        add!(T, v)
    end

    wt = Dict{ET,Float64}()
    for e in EE
        wt[e] = rand()
    end

    wt_list = [wt[e] for e in EE]
    p = sortperm(wt_list)
    E_list = [EE[j] for j in p]

    P = Partition(G.V)

    for e in E_list
        a, b = e
        if !in_same_part(P, a, b)
            add!(T, a, b)
            merge_parts!(P, a, b)
        end
    end
    return T
end

"""
`find_path(G,s,t)` finds a shortest path from `s` to `t`. If no
such path exists, an empty list is returned.
"""
function find_path(G::AbstractSimpleGraph, s, t)
    T = eltype(G)
    if ~has(G, s) || ~has(G, t)
        error("Source and/or target vertex is not in this graph")
    end
    if s == t
        #   result = Array{T}(undef,1)
        #   result[1] = s
        return [s]
    end

    # set up a queue for vertex exploration
    Q = Queue{T}()
    enqueue!(Q, s)

    # set up trace-back dictionary
    tracer = Dict{T,T}()
    tracer[s] = s

    while length(Q) > 0
        v = dequeue!(Q)
        Nv = G.N[v]
        for w in Nv
            if haskey(tracer, w)
                continue
            end
            tracer[w] = v
            enqueue!(Q, w)

            if w == t  # success!
                path = Array{T}(undef, 1)
                path[1] = t
                while path[1] != s
                    v = tracer[path[1]]
                    pushfirst!(path, v)
                end
                return path

            end
        end
    end
    return T[]   # return empty array if no path found
end

# functions for finding distances between vertices in a graph. the
# distance between two vertices is the number of edges in a shortest
# path between them. if there is no such path, the distance is
# undefined (or infinite) but since we want the return value to be an
# Int, we use -1 to signal this.

# find the distance between specified vertices

"""
`dist(G,u,v)` finds the length of a shortest path from `u` to `v` in
`G`. Returns `-1` if no such path exists.

`dist(G,u)` finds the distance from `u` to all other vertices in
`G`. Result is returned as a `Dict`.

`dist(G)` finds all pairs of distances in `G`. Result is a `Dict`
whose `[u,v]` entry is the distance from `u` to `v`.
"""
function dist(G::AbstractSimpleGraph, u, v)
    if !has(G, u) || !has(G, v)
        error("One or both of $u and $v are not vertices of this graph")
    end
    if typeof(G) <: SimpleGraph && cache_check(G, :dist)
        d = cache_recall(G, :dist)
        return d[u, v]
    end
    if u == v
        return 0
    end
    return length(find_path(G, u, v)) - 1
end

# find all distances from a given vertex
function dist(G::AbstractSimpleGraph, v)
    T = eltype(G)
    d = Dict{T,Int}()
    if !has(G, v)
        error("Given vertex is not in this graph")
    end

    d[v] = 0
    Q = Queue{T}()
    enqueue!(Q, v)

    while length(Q) > 0
        w = dequeue!(Q)  # get 1st vertex in the queue
        Nw = G.N[w]
        for x in Nw
            if !haskey(d, x)
                d[x] = d[w] + 1
                enqueue!(Q, x)
            end
        end
    end

    # record -1 for any vertices we missed
    for x in G.V
        if !haskey(d, x)
            d[x] = -1
        end
    end

    return d
end

# find all distances between all vertices
function dist(G::AbstractSimpleGraph)
    if cache_check(G, :dist)
        return cache_recall(G, :dist)
    end

    T = eltype(G)
    dd = Dict{Tuple{T,T},Int}()
    vtcs = vlist(G)

    for v in vtcs
        d = dist(G, v)
        for w in vtcs
            dd[(v, w)] = d[w]
        end
    end
    cache_save(G, :dist, dd)
    return dd
end

"""
`eccentricity(G,v)` returns the eccentricty of vertex `v`
in the graph `G`. This is the maximum distance from `v`
to another vertex (or -1 if the graph is not connected).
"""
function eccentricity(G::SimpleGraph, v)
    if !has(G, v)
        error("$v is not a vertex of this graph")
    end
    d = collect(values(dist(G, v)))
    if minimum(d) < 0
        return -1
    end
    return maximum(d)
end

"""
`graph_center(G)` returns the set of vertices of a `SimpleGraph` with minimum
eccentricities.
"""
function graph_center(G::SimpleGraph)::Set
    if cache_check(G, :graph_center)
        return cache_recall(G, :graph_center)
    end
    if G.cache_flag
        dist(G) # force all pairs distance computation
    end

    xtable = Dict((v, eccentricity(G, v)) for v in G.V)
    min_r = minimum(values(xtable))

    if min_r < 0
        return deepcopy(G.V)
    end

    A = Set(v for v in keys(xtable) if xtable[v] == min_r)
    cache_save(G, :graph_center, A)
    return A
end

"""
`radius(G)` returns the radius of the graph `G`. This is the
minimum `eccentricity` of a vertex of `G` (or -1 if the graph
is not connected).
"""
function radius(G::SimpleGraph)
    if cache_check(G, :radius)
        return cache_recall(G, :radius)
    end
    D = dist_matrix(G)
    if minimum(D) < 0
        return -1
    end
    r = minimum(maximum(D, dims = 1))
    cache_save(G, :radius, r)
    return r
end

"""
`wiener_index(G)` is the sum of the distances between vertices in `G`.
Returns -1 if `G` is not connected.
"""
function wiener_index(G::SimpleGraph)::Int
    if is_connected(G)
        return div(sum(values(dist(G))), 2)
    end
    return -1
end

# Calculate the diameter of a graph, but return -1 if not connected.

"""
`diam(G)` returns the diameter of `G` or `-1` if `G` is not connected.
"""
function diam(G::SimpleGraph)
    if is_connected(G)
        return maximum(values(dist(G)))
    end
    return -1
end

# Determine if a given edge in a graph is a cut edge. If there is no
# such edge in the graph, an error is raised.

"""
`is_cut_edge(G,u,v)` [or `is_cut_edge(G,e)`] determins if `(u,v)` [or
`e`] is a cut edge of `G`.
"""
function is_cut_edge(G::SimpleGraph, u, v)
    if !has(G, u, v)
        error("No such edge in this graph")
    end

    SimpleGraphs.delete!(G, u, v)
    P = find_path(G, u, v)
    result = false
    if length(P) == 0
        result = true
    end
    add!(G, u, v)
    return result
end

# When called as is_cut_edge(G,e), we assume e is a tuple or list
# whose first two entries are the end points of the edge
function is_cut_edge(G::SimpleGraph, e)
    return is_cut_edge(G, e[1], e[2])
end

"""
`is_acyclic(G)` returns `true` if `G` has no cycles and `false`
otherwise.
"""
function is_acyclic(G::SimpleGraph)
    n = NV(G)
    m = NE(G)
    c = num_components(G)
    return m == n - c
end


text:
######################################################################
# Input Abstraction Layer
# -----
# Licensed under MIT License
using BitFlags
using GLFW

export InputAction, Press, Release, Repeat
@enum InputAction Press = Integer(GLFW.PRESS) Release = Integer(GLFW.RELEASE) Repeat = Integer(GLFW.REPEAT)
InputAction(action::GLFW.Action) = InputAction(Integer(action))

export ModifierKey, NoMod, ShiftMod, ControlMod, AltMod, SuperMod, CapsLockMod, NumLockMod
@bitflag ModifierKey::UInt16 NoMod = 0 ShiftMod ControlMod AltMod SuperMod CapsLockMod NumLockMod

export pollevents, waitevents, post_empty_event
pollevents() = GLFW.PollEvents()
waitevents() = GLFW.WaitEvents()
waitevents(timeout::AbstractFloat) = GLFW.WaitEvents(timeout)
post_empty_event() = GLFW.PostEmptyEvent()

include("./input.keyboard.jl")
include("./input.mouse.jl")
# It appears GLFW.jl does not support the latest version of GLFW? Gamepad support is incomplete.
# include("./input.gamepad.jl")

# ===== Events =====

"""Called when a key is pressed down."""
function on_key_press(wnd::Window, key::Key, scancode::Integer, modifiers::ModifierKey) end
"""Called when a key is released."""
function on_key_release(wnd::Window, key::Key, scancode::Integer, modifiers::ModifierKey) end
"""Called when a key is being held, triggering continuous textual input."""
function on_key_repeat(wnd::Window, key::Key, scancode::Integer, modifiers::ModifierKey) end
"""Called when a key results in a Unicode character input."""
function on_receive_char(wnd::Window, char::Char) end
"""Called when the mouse cursor is moved while within the window's boundaries."""
function on_mouse_move(wnd::Window, xpos::AbstractFloat, ypos::AbstractFloat) end
"""Called when the mouse cursor enters the window's boundaries."""
function on_mouse_enter(wnd::Window) end
"""Called when the mouse cursor leaves the window's boundaries."""
function on_mouse_leave(wnd::Window) end
"""Called when a mouse button is pressed down."""
function on_mouse_press(wnd::Window, button::MouseButton, modifiers::ModifierKey) end
"""Called when a mouse button is released."""
function on_mouse_release(wnd::Window, button::MouseButton, modifiers::ModifierKey) end
"""Called when a new gamepad is connected."""

function register_default_input_callbacks(wnd::Window)
    register_key_callback(wnd) do _::Window, key::Key, scancode::Integer, action::InputAction, modifiers::ModifierKey
        if action === Press
            on_key_press(wnd, key, scancode, modifiers)
        elseif action === Release
            on_key_release(wnd, key, scancode, modifiers)
        elseif action === Repeat
            on_key_repeat(wnd, key, scancode, modifiers)
        else
            error("unknown key action $action")
        end
    end
    register_text_callback(wnd) do _::Window, char::Char; on_receive_char(wnd, char) end
    register_cursor_pos_callback(wnd) do _::Window, xpos::AbstractFloat, ypos::AbstractFloat; on_mouse_move(wnd, xpos, ypos) end
    register_cursor_enter_callback(wnd) do _::Window, entered::Bool
        if entered
            on_mouse_enter(wnd)
        else
            on_mouse_leave(wnd)
        end
    end
    register_mouse_button_callback(wnd) do _::Window, button::MouseButton, action::InputAction, modifiers::ModifierKey
        @assert action !== Repeat
        if action === Press
            on_mouse_press(wnd, button, modifiers)
        else
            on_mouse_release(wnd, button, modifiers)
        end
    end
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVCProjStd

using Test, RLinearAlgebra, LinearAlgebra, Random

@testset "LSVC Projection Standard -- Procedural" begin
    # Supertype and aliases
    @test supertype(LinSysVecColProjStd) == LinSysVecColProjection
    @test CoordinateDescent == LinSysVecColProjStd
    @test GaussSeidel == LinSysVecColProjStd

    # Verify that the residual projection is zero
    Random.seed!(1010)

    A = rand(5,10)
    x = rand(10)
    b = A * x

    rsub = GaussSeidel()

    for i = 1:10
        @test let
            # Initialization of iteration
            z = rand(10)

            # Search Direction
            e = zeros(10)
            e[i] = 1.0

            # Update
            RLinearAlgebra.rsubsolve!(rsub, z, (e, A, dot(A * e, A * z - b)), i)

            abs(dot(A[:,i], A * z - b)) < 1e-15
        end
    end
end
end # End module


text:
"""
Delay coordinates embedding, `Dataset` structure and neighhborhoods.
Basic package used in the ecosystem of DynamicalSystems.jl.
"""
module DelayEmbeddings

include("dataset.jl")
include("subdataset.jl")
include("embeddings.jl")
include("utils.jl")
include("neighborhoods.jl")

include("traditional_de/estimate_delay.jl")
include("traditional_de/estimate_dimension.jl")
include("traditional_de/automated.jl")

include("unified_de/pecora.jl")
include("unified_de/uzal_cost.jl")
include("unified_de/MDOP.jl")
include("unified_de/garcia_almeida.jl")
include("unified_de/pecuzal.jl")

include("deprecate.jl")

end


text:
using JSON
using Random
using CUDA
using BenchmarkTools
using Statistics
using Logging
using Dates
using DataStructures
include("optimizers/optimizer.jl")
include("environments/collect_points_env.jl")
include("brains/continuous_time_rnn.jl")

configuration = JSON.parsefile("configurations/CMA_ES_Deap_CTRNN_Dense.json",dicttype=DataStructures.OrderedDict)

number_generations = configuration["number_generations"]
number_validation_runs = configuration["number_validation_runs"]
number_rounds = configuration["number_rounds"]
maximum_env_seed = configuration["maximum_env_seed"]
environment = configuration["environment"]
brain = configuration["brain"]
optimizer = configuration["optimizer"]
brain_cfg = ContinuousTimeRNN(brain["delta_t"],brain["number_neurons"],separated,brain["clipping_range_min"],brain["clipping_range_max"],brain["alpha"])
number_inputs = get_number_inputs()
number_outputs = get_number_outputs()
environment_cfg = CollectPoints(environment["maze_columns"],environment["maze_rows"],environment["maze_cell_size"],environment["agent_radius"],environment["point_radius"],environment["agent_movement_range"],environment["reward_per_collected_positive_point"],environment["reward_per_collected_negative_point"],environment["number_time_steps"],number_inputs,number_outputs)
number_individuals = optimizer["population_size"]
brain_state = generate_brain_state(number_inputs,number_outputs,brain)
free_parameters = get_individual_size(brain_state)

optimizer = inititalize_optimizer(free_parameters,optimizer)

for generation in 1:number_generations
    individuals = fill(0.0f0,number_individuals,free_parameters)
    genomes = ask(optimizer)
    for i in 1:number_individuals
        for j in 1:free_parameters
            individuals[i,j] = (genomes[i])[j]
        end
    end
    rewards_training = 100 .* rand(Float32,number_individuals)
    tell(optimizer,rewards_training)
end


text:
module BilayerNetworks

using AlgebraicPetri
using Catlab
using Catlab.CategoricalAlgebra
using Catlab.Graphics
import Catlab.CategoricalAlgebra: migrate!

export ThBilayerNetwork, AbstractBilayerNetwork, BilayerNetwork,
    ThLabelledBilayerNetwork, AbstractLabelledBilayerNetwork, LabelledBilayerNetwork

@present ThBilayerNetwork(FreeSchema) begin
    (Qin, Qout, Win, Wn, Wa, Box)::Ob
    arg::Hom(Win, Qin)
    call::Hom(Win, Box)
    influx::Hom(Wa, Box)
    infusion::Hom(Wa, Qout)
    efflux::Hom(Wn, Box)
    effusion::Hom(Wn, Qout)
end

@present ThLabelledBilayerNetwork <: ThBilayerNetwork begin
    Name::Data

    parameter::Attr(Box, Name)
    variable::Attr(Qin, Name)
    tanvar::Attr(Qout, Name)
end

const AbstractBilayerNetwork = AbstractACSetType(ThBilayerNetwork)
const BilayerNetwork = ACSetType(ThBilayerNetwork)

const AbstractLabelledBilayerNetwork = AbstractACSetType(ThLabelledBilayerNetwork)
const LabelledBilayerNetwork = ACSetType(ThLabelledBilayerNetwork){Symbol}

function migrate!(bn::AbstractBilayerNetwork, pn::AbstractPetriNet)
    migrate!(bn, pn,
             Dict(:Qin=>:S, :Qout=>:S, :Box=>:T, :Win=>:I, :Wn=>:I, :Wa=>:O),
             Dict(:arg=>:is,
                  :call=>:it,
                  :efflux=>:it,
                  :effusion=>:is,
                  :influx=>:ot,
                  :infusion=>:os))
end

function migrate!(bn::AbstractLabelledBilayerNetwork, pn::AbstractLabelledPetriNet)
    migrate!(bn, pn,
             Dict(:Qin=>:S, :Qout=>:S, :Box=>:T, :Win=>:I, :Wn=>:I, :Wa=>:O, :Name=>:Name),
             Dict(:arg=>:is,
                  :call=>:it,
                  :efflux=>:it,
                  :effusion=>:is,
                  :influx=>:ot,
                  :infusion=>:os,
                  :parameter=>:tname,
                  :variable=>:sname,
                  :tanvar=>:sname))
end

function migrate!(pn::AbstractPetriNet, bn::AbstractBilayerNetwork)
    migrate!(pn,bn,
         Dict(:S=>:Qin, :T=>:Box, :I=>:Win, :O=>:Wa),
         Dict(:is=>:arg,
              :it=>:call,
              :ot=>:influx,
              :os=>:infusion))
end

function migrate!(pn::AbstractLabelledPetriNet, bn::AbstractLabelledBilayerNetwork)
    migrate!(pn,bn,
         Dict(:S=>:Qin, :T=>:Box, :I=>:Win, :O=>:Wa, :Name=>:Name),
         Dict(:is=>:arg,
              :it=>:call,
              :ot=>:influx,
              :os=>:infusion,
              :tname=>:parameter,
              :sname=>:variable))
end

function propertygraph(bn::AbstractBilayerNetwork;
                       prog::AbstractString="dot", graph_attrs::AbstractDict=Dict(),
                       node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
                       node_labels::Bool=true, edge_labels::Bool=false)

    Qinshift = nparts(bn, :Box)
    Qoutshift = nparts(bn, :Box) + nparts(bn, :Qin)

    colors = ["#E28F41", "#6C9AC3", "#80BE63", "#A8DCD9", "#E0E5CD"]

    # vertex partition helpers V = coproduct(Box, Qin, Qout)
    isbox(v) = v <=Qinshift
    isqin(v) = Qinshift < v <= Qoutshift
    isqout(v) = Qoutshift < v

    # edge partition helpers E = coproduct(Win, W, Wn)
    calledge(e) = e <= nparts(bn, :Win)
    influxedge(e) = nparts(bn, :Win) < e <= nparts(bn, :Win) + nparts(bn, :Wa)
    effluxedge(e) = nparts(bn, :Win) + nparts(bn, :Wa) < e

    node_labeler(v) = if isbox(v)
        Dict(:label=>"p[$v]", :style=>"rounded, filled", :width=>"0.75", :height=>"0.5", :fixedsize=>"true", :shape=>"rectangle", :color=>colors[1])
    elseif isqin(v)
        return Dict(:label=>"u[$(v-Qinshift)]", :style=>"filled", :color=>colors[2])
    elseif isqout(v)
        return Dict(:label=>"u[$(v-Qoutshift)]'", :style=>"filled", :color=>colors[3])
    else
        return Dict(:label=>"$v", :style=>"filled", :color=>colors[4])
    end

    edge_labeler(e) = if calledge(e)
        return Dict(:label => edge_labels ? string(e) : "", :color=>"black", :style=>"dashed")
    elseif effluxedge(e)
        return Dict(:label => edge_labels ? string(e) : "", :color=>"red")
    else
        return Dict(:label => edge_labels ? string(e) : "", :color=>"black")
    end

    g′ = Catlab.Graphs.Graph()
    add_parts!(g′, :E, nparts(bn, :Win) + nparts(bn, :Wa) + nparts(bn, :Wn))
    add_parts!(g′, :V, nparts(bn, :Box) + nparts(bn, :Qin) + nparts(bn, :Qout))

    srcs = vcat(bn[:arg].+Qinshift, bn[:influx], bn[:efflux])
    tgts = vcat(bn[:call], bn[:infusion].+Qoutshift, bn[:effusion].+Qoutshift)

    g′[:, :src] = srcs
    g′[:, :tgt] = tgts
    Catlab.Graphs.PropertyGraph{Any}(g′, node_labeler, edge_labeler;
                                     prog = prog,
                                     graph = merge!(Dict(:rankdir => "TB"), graph_attrs),
                                     node = merge!(Graphics.GraphvizGraphs.default_node_attrs(node_labels), node_attrs),
                                     edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
                                     )
end

function GraphvizGraphs.to_graphviz(g::AbstractBilayerNetwork; kwargs...)
    to_graphviz(propertygraph(g; kwargs...))
end
function propertygraph(bn::AbstractLabelledBilayerNetwork;
                       prog::AbstractString="dot", graph_attrs::AbstractDict=Dict(),
                       node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
                       node_labels::Bool=true, edge_labels::Bool=false)

    Qinshift = nparts(bn, :Box)
    Qoutshift = nparts(bn, :Box) + nparts(bn, :Qin)

    colors = ["#E28F41", "#6C9AC3", "#80BE63", "#A8DCD9", "#E0E5CD"]

    # vertex partition helpers V = coproduct(Box, Qin, Qout)
    isbox(v) = v <=Qinshift
    isqin(v) = Qinshift < v <= Qoutshift
    isqout(v) = Qoutshift < v

    # edge partition helpers E = coproduct(Win, W, Wn)
    calledge(e) = e <= nparts(bn, :Win)
    influxedge(e) = nparts(bn, :Win) < e <= nparts(bn, :Win) + nparts(bn, :Wa)
    effluxedge(e) = nparts(bn, :Win) + nparts(bn, :Wa) < e

    node_labeler(v) = if isbox(v)
        Dict(:label=>"$(bn[v, :parameter])", :style=>"rounded, filled", :width=>"0.75", :height=>"0.5", :fixedsize=>"true", :shape=>"rectangle", :color=>colors[1])
    elseif isqin(v)
        return Dict(:label=>"$(bn[v-Qinshift, :variable])", :style=>"filled", :color=>colors[2])
    elseif isqout(v)
        return Dict(:label=>"$(bn[v-Qoutshift, :tanvar])'", :style=>"filled", :color=>colors[3])
    else
        return Dict(:label=>"$v", :style=>"filled", :color=>colors[4])
    end

    edge_labeler(e) = if calledge(e)
        return Dict(:label => edge_labels ? string(e) : "", :color=>"black", :style=>"dashed")
    elseif effluxedge(e)
        return Dict(:label => edge_labels ? string(e) : "", :color=>"red")
    else
        return Dict(:label => edge_labels ? string(e) : "", :color=>"black")
    end

    g′ = Catlab.Graphs.Graph()
    add_parts!(g′, :E, nparts(bn, :Win) + nparts(bn, :Wa) + nparts(bn, :Wn))
    add_parts!(g′, :V, nparts(bn, :Box) + nparts(bn, :Qin) + nparts(bn, :Qout))

    srcs = vcat(bn[:arg].+Qinshift, bn[:influx], bn[:efflux])
    tgts = vcat(bn[:call], bn[:infusion].+Qoutshift, bn[:effusion].+Qoutshift)

    g′[:, :src] = srcs
    g′[:, :tgt] = tgts
    Catlab.Graphs.PropertyGraph{Any}(g′, node_labeler, edge_labeler;
                                     prog = prog,
                                     graph = merge!(Dict(:rankdir => "TB"), graph_attrs),
                                     node = merge!(Graphics.GraphvizGraphs.default_node_attrs(node_labels), node_attrs),
                                     edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
                                     )
end

function GraphvizGraphs.to_graphviz(g::AbstractLabelledBilayerNetwork; kwargs...)
    to_graphviz(propertygraph(g; kwargs...))
end

function evaluate!(du, ϕ, bn::AbstractLabelledBilayerNetwork, state; params...)
    du.= 0.0
    ϕ .= 1.0
    for i in parts(bn, :Win)
        ϕ[bn[i, :call]] *= state[bn[i, :arg]]
    end
    for i in parts(bn, :Box)
        ϕ[i] *= params[bn[i, :parameter]]
    end

    for i in parts(bn, :Wn)
        du[bn[i, :effusion]] -= ϕ[bn[i, :efflux]]
    end
    for i in parts(bn, :Wa)
        du[bn[i, :infusion]] += ϕ[bn[i, :influx]]
    end
    return du
end

evaluate(bn::AbstractLabelledBilayerNetwork, state; params...) = evaluate!(zeros(length(state)), ones(nparts(bn, :Box)), bn, state; params...)

function paramexps(bn::AbstractLabelledBilayerNetwork, params::Symbol)
    map(parts(bn, :Box)) do i
        p = bn[i, :parameter]
        :(ϕ[$i] *= params[$(Meta.quot(p))])
    end
end

function paramexps(bn::AbstractLabelledBilayerNetwork, params)
    map(parts(bn, :Box)) do i
        p = bn[i, :parameter]
        β = params[p]
        :(ϕ[$i] *= $β)
    end
end

funcwrap(du, ϕ, state, params::Symbol, body::Expr) = :(f!($du, $ϕ, $state, $params, t) = $body)
funcwrap(du, ϕ, state, params, body::Expr) = :(f!($du, $ϕ, $state, t) = $body)

function compile(bn::Union{AbstractLabelledBilayerNetwork, AbstractBilayerNetwork}, du::Symbol, ϕ::Symbol, state::Symbol, params)
    body = quote
        $du.= 0.0
        $ϕ .= 1.0
    end

    ϕs = map(parts(bn, :Win)) do i
        j = bn[i, :arg]
        k = bn[i, :call]
        :(ϕ[$k] *= $state[$j])
    end
    append!(body.args, ϕs)
    ps = paramexps(bn, params)
    append!(body.args, ps)

    effs = map(parts(bn, :Wn)) do i
        j = bn[i, :efflux]
        k = bn[i, :effusion]
        :(du[$k] -= ϕ[$j])
    end
    append!(body.args, effs)

    infs = map(parts(bn, :Wa)) do i
        j = bn[i,:influx]
        k = bn[i,:infusion]
        :(du[$k] += ϕ[$j])
    end
    append!(body.args, infs)
    push!(body.args, :(return $du))
    return funcwrap(du, ϕ, state, params, body)
end

compile(bn, du, ϕ, state; params...) = compile(bn, du, ϕ, state, params)

end


text:
mutable struct Config

    # phenotype config
    input_nodes::Int64
    output_nodes::Int64
    initial_hidden_nodes::Int64
    max_weight::Float64
    min_weight::Float64
    nn_activation::Symbol
    #weight_stdev::Float64

    # GA config
    pop_size::Int64
    max_fitness_threshold::Float64
    prob_addnode::Float64
    prob_removenode::Float64
    prob_addlayer::Float64
    prob_mutatebias::Float64
    bias_mutation_power::Float64
    prob_mutate_weight::Float64 # dynamic mutation rate (future release)
    weight_mutation_power::Float64
    prob_structural_mutation::Float64
    prob_enable_mutation::Float64
    elitism::Bool

    # genotype compatibility
    compatibility_threshold::Float64
    compatibility_change::Float64
    excess_coeficient::Float64
    disjoint_coeficient::Float64
    weight_coeficient::Float64

    # species
    species_size::Int64
    survival_threshold::Float64 # only the best 20% for each species is allowed to mate
    old_threshold::Int64
    youth_threshold::Int64
    old_penalty::Float64    # always in (0,1)
    youth_boost::Float64    # always in (1,2)
    max_stagnation::Int64
    best_max_stagnation::Int64
    limit_species::Int64

    # evaluation
    dataset_id::Int64
    fitness_function_id::Int64
    normalise_input::Bool
    balance_dataset::Bool

    # checkpoint
    checkpoint_filename::String

    function Config(params::Dict{String,String})

        new(
            # phenotype
            parse(Int64, params["input_nodes"]),
            parse(Int64, params["output_nodes"]),
            parse(Int64, params["initial_hidden_nodes"]),
            parse(Float64, params["max_weight"]),
            parse(Float64, params["min_weight"]),
            Meta.parse(params["nn_activation"]),
            #float(params["weight_stdev"]),

            # GA
            parse(Int64, params["pop_size"]),
            parse(Float64, params["max_fitness_threshold"]),
            parse(Float64, params["prob_addnode"]),
            parse(Float64, params["prob_removenode"]),
            parse(Float64, params["prob_addlayer"]),
            parse(Float64, params["prob_mutatebias"]),
            parse(Float64, params["bias_mutation_power"]),
            parse(Float64, params["prob_mutate_weight"]),
            parse(Float64, params["weight_mutation_power"]),
            parse(Float64, params["prob_structural_mutation"]),
            parse(Float64, params["prob_enable_mutation"]),
            parse(Bool, params["elitism"]),

            # genotype compatibility
            parse(Float64, params["compatibility_threshold"]),
            parse(Float64, params["compatibility_change"]),
            parse(Float64, params["excess_coeficient"]),
            parse(Float64, params["disjoint_coeficient"]),
            parse(Float64, params["weight_coeficient"]),

            # species
            parse(Int64, params["species_size"]),
            parse(Float64, params["survival_threshold"]),
            parse(Int64, params["old_threshold"]),
            parse(Int64, params["youth_threshold"]),
            parse(Float64, params["old_penalty"]),
            parse(Float64, params["youth_boost"]),
            parse(Int64, params["max_stagnation"]),
            parse(Int64, params["best_max_stagnation"]),
            parse(Int64, params["limit_species"]),

            # evaluation
            parse(Int64, params["dataset_id"]),
            parse(Int64, params["fitness_function_id"]),
            parse(Bool, params["normalise_input"]),
            parse(Bool, params["balance_dataset"]),

            # checkpoint
            params["checkpoint_filename"]
         )
    end
end

function loadConfig(file::String)
    str = open(f->read(f, String), file)

    str = replace(str, r"\r(\n)?" => '\n')
    ls = split(str, "\n")

    ls = filter(l->length(l)>0 && l[1] != '#', ls)
    lsMap = map(x->split(x,'='),ls)
    params = Dict{String,String}()
    for i = 1:length(lsMap)
#         println((i,lsMap[i][1],lsMap[i][2]))
        params[rstrip(lsMap[i][1])] = lstrip(lsMap[i][2])
    end

    return params
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
function node_f(net::FastNet,c,rp)
    @inbounds nodepos=net.ncstart[c]+rp
    @inbounds net.nid[nodepos]
end

function node_f(net::FastNet,rp)
    @inbounds net.nid[rp]
end

function nodestate_f(net::FastNet,nid)
    @inbounds net.nstate[nid]
end

function nodestate_f!(net::FastNet,nid,cls)
    @inbounds begin
        if net.nstate[nid]===cls
            return nothing
        end
        while net.nstate[nid]<cls
            upstate_node!(net,nid)
        end
        while  net.nstate[nid]>cls
            destate_node!(net,nid)
        end
        cur=net.nout[nid]
        while cur!=0
            dst=net.kdst[cur]
            dc=net.nstate[dst]
            nlc=net.ttable[cls,dc]
            linkstate_f!(net,cur,nlc)
            cur=net.knexts[cur]
        end
        cur=net.nin[nid]
        while cur!=0
            src=net.ksrc[cur]
            sc=net.nstate[src]
            nlc=net.ttable[sc,cls]
            linkstate_f!(net,cur,nlc)
            cur=net.knextd[cur]
        end
    end
    nothing
end

function makenode_f!(net::FastNet,cls)
    nid=node_f(net,net.C,1)
    nodestate_f!(net,nid,cls)
    nid
end

function destroynode_f!(net::FastNet,nid)
    @inbounds begin
        cur=net.nout[nid]
        while (cur!=0)
            destroylink_f!(net,cur)
            cur=net.nout[nid]
        end
        cur=net.nin[nid]
        while (cur!=0)
            destroylink_f!(net,cur)
            cur=net.nin[nid]
        end
        nodestate_f!(net,nid,net.C)
    end
    nothing
end

function countnodes_f(net::FastNet)
    ret=0
    nc=net.C-1
    for i=1:nc
        @inbounds ret+=net.nclen[i]
    end
    ret
end

function countnodes_f(net::FastNet,cls::Integer)
    @inbounds net.nclen[cls]
end

function nodecounts_f(net::FastNet)
    n=net.C-1
    @inbounds net.nclen[1:n]
end

function countnodes_f(net::FastNet,cls::Union{Vector,Tuple})
    ret=0
    for c in cls
        @inbounds ret+=net.nclen[c]
    end
    ret
end

function firstlinkin_f(net::FastNet,nid)
    net.nin[nid]
end

function firstlinkout_f(net::FastNet,nid)
    net.nout[nid]
end

function degree_f(net::FastNet,nid)
    indegree_f(net,nid)+outdegree_f(net,nid)
end

function indegree_f(net::FastNet,nid)
    ret=0
    @inbounds begin
        cur=net.nin[nid]
        while cur!=0
             ret+=1
             cur=net.knextd[cur]
        end
    end
    ret
end

function outdegree_f(net::FastNet,nid)
    ret=0
    @inbounds begin
        cur=net.nout[nid]
        while cur!=0
             ret+=1
             cur=net.knexts[cur]
        end
    end
    ret
end

function nodeexists_f(net::FastNet,nid)
    @inbounds net.nstate[nid]!=net.C
end

function makenodes_f!(net::FastNet,N,cls)
    for n=1:N
        makenode!(net,cls)
    end
    nothing
end

function randomnode_f(net::FastNet)
    node_f(net,rand(net.rng,1:countnodes_f(net)))
end

function randomnode_f(net::FastNet,cls::Integer)
    node_f(net,cls,rand(net.rng,1:countnodes_f(net,cls)))
end

function randomnode_f(net::FastNet,cls::Union{Array,Tuple})
    tot=countnodes_f(net,cls)
    r=rand(net.rng,1:tot)
    old=r
    i=1
    cr=old-countnodes_f(net,cls[i])
    while cr>0
        i+=1
        old=cr
        cr=old-countnodes_f(net,cls[i])
    end
    node_f(net,cls[i],old)
end

function adjacent_f(net::FastNet,nida,nidb)
    @inbounds begin
        link=net.firstlinkout[nida]
        while (link!==0)
            if net.kdst[link]===nidb
                return link
            end
            link=net.knexts[link]
        end
    end
    link=net.firstlinkin[nida]
    while (link!==0)
        if net.ksrc[link]===nidb
            return link
        end
        link=net.knextd[link]
    end
    0
end


text:
using Luxor

p1 = [5, 3]
p2 = [5, -3]

r = p1 + t * (p2 - p1)

max_stars_repo_name: TheCedarPrince/DataExplorers

max_issues_repo_name: TheCedarPrince/DataExplorers

max_forks_repo_name: TheCedarPrince/DataExplorers


text:
"""
Test script forOptStoic/MinFlux procedure.
Model construction and package imports can be found in `testmodel_construction.jl`.
Output is expected to be as follows:
Dict{Any, Any} with 10 entries:
  "r5"      => -1.0
  "r2"      => 1.0
  "EX_m6"   => 2.0
  "biomass" => -2.0
  "EX_co1"  => -2.0
  "EX_m1"   => -1.0
  "EX_h"    => 2.0
  "EX_co2"  => 2.0
  "r3"      => 1.0
  "r4"      => 1.0
"""
include("testmodel_construction.jl")

minflux_bounds = MinFlux_preprocessing(m, reaction_dG)
substrates = ["m1"]
targets = ["m6"]
cofactors = ["co1", "co2", "h"] #allow cofactors

bounds = Dict("substrates" => -1, "reactants" => (-15, 15))

os = build_OptStoic_model(m, substrates, targets, cofactors, Gurobi, energy_dc = formation_dG, variable_bounds = bounds)
os = optimize_OptStoic(os)
os_solution = solution_OptStoic(os)

if !isnothing(os_solution)
    m, mf = build_MinFlux_model(m, os_solution, minflux_bounds, optimizer = Gurobi)
    mf = MinFlux(mf)
    MinFlux_solution(m, mf)
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
module TestMLBaseWrapper

using Test
using TSML
using DataFrames

@testset "MLBase transformers" begin

  @testset "StandardScaler transforms features" begin
    features = [
      5 10;
      -5 0;
      0 5;
    ] |> DataFrame
    labels = [
      "x";
      "y";
      "z";
    ]
    expected_transformed = [
      1.0 1.0;
      -1.0 -1.0;
      0.0 0.0;
    ] |> DataFrame
    standard_scaler = StandardScaler()
    fit!(standard_scaler, features, labels)
    transformed = transform!(standard_scaler, features)
    @test (transformed .== expected_transformed) |> Matrix |> sum == 6
  end

end

end # module


text:
using Base.Test
using BSBP

@testset "Simulation" begin
    score,gap,rtime,ncount = Simulation()
    @test mean(score) >= 7
    @test mean(gap) == 0
    @test mean(rtime) <= 2
    @test mean(ncount) <= 2
end

@testset "TransportationMode" begin
    score,gap,rtime,ncount = TransportationMode()
    @test score == 756
    @test gap <= 40
    @test rtime <= 100
    @test ncount >= 2000
end

@testset "TransportationModeCV" begin
    score,gap,rtime,ncount = TransportationModeCV()
    @test mean(score) >=605
    @test mean(gap) <= 40
    @test mean(rtime) <= 100
    @test mean(ncount) >= 2000
end


text:
module Utilities

include("plushalf.jl")

end # module


text:
struct DeviceEnergyBalanceConstraintSpec
    constraint_name::Symbol
    energy_variable::Symbol
    initial_condition::Type{<:InitialConditionType}
    pin_variable_names::Vector{Symbol}
    pout_variable_names::Vector{Symbol}
    parameter_name::Union{Nothing, String}
    forecast_label::Union{Nothing, String}
    multiplier_func::Union{Nothing, Function}
    constraint_func::Function
end

function DeviceEnergyBalanceConstraintSpec(;
    constraint_name::Symbol,
    energy_variable::Symbol,
    initial_condition::Type{<:InitialConditionType},
    constraint_func::Function,
    pin_variable_names::Vector{Symbol} = Vector{Symbol}(),
    pout_variable_names::Vector{Symbol} = Vector{Symbol}(),
    parameter_name::Union{Nothing, String} = nothing,
    forecast_label::Union{Nothing, String} = nothing,
    multiplier_func::Union{Nothing, Function} = nothing,
)
    return DeviceEnergyBalanceConstraintSpec(
        constraint_name,
        energy_variable,
        initial_condition,
        pin_variable_names,
        pout_variable_names,
        parameter_name,
        forecast_label,
        multiplier_func,
        constraint_func,
    )
end

function add_constraints!(
    optimization_container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {
    T <: EnergyBalanceConstraint,
    U <: VariableType,
    V <: PSY.Device,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    use_parameters = model_has_parameters(optimization_container)
    use_forecasts = model_uses_forecasts(optimization_container)
    # @assert !(use_parameters && !use_forecasts)
    spec = DeviceEnergyBalanceConstraintSpec(
        T,
        U,
        V,
        W,
        X,
        feedforward,
        use_parameters,
        use_forecasts,
    )
    energy_balance_constraints!(optimization_container, devices, model, feedforward, spec)
end

function energy_balance_constraints!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    spec::DeviceEnergyBalanceConstraintSpec,
) where {T <: PSY.Device, U <: AbstractDeviceFormulation}
    initial_conditions =
        get_initial_conditions(optimization_container, spec.initial_condition, T)
    _apply_energy_balance_constraint_spec!(
        optimization_container,
        spec,
        devices,
        initial_conditions,
        model,
        feedforward,
    )
end

struct DeviceEnergyBalanceConstraintSpecInternal
    constraint_infos::Vector{<:EnergyBalanceConstraintInfo}
    constraint_name::Symbol
    energy_variable::Symbol
    pin_variable_names::Vector{Symbol}
    pout_variable_names::Vector{Symbol}
    param_reference::Union{Nothing, UpdateRef}
end

function _apply_energy_balance_constraint_spec!(
    optimization_container,
    spec,
    devices::IS.FlattenIteratorWrapper{T},
    initial_conditions::Vector{InitialCondition},
    model,
    ff_affected_variables,
) where {T <: PSY.Device}
    constraint_infos = Vector{EnergyBalanceConstraintInfo}(undef, length(devices))
    for (ix, ic) in enumerate(initial_conditions)
        device = ic.device
        dev_name = PSY.get_name(device)
        if !isnothing(spec.forecast_label)
            ts_vector = get_time_series(optimization_container, device, spec.forecast_label)
            multiplier = spec.multiplier_func(device)
            constraint_info = EnergyBalanceConstraintInfo(
                dev_name,
                get_efficiency(device, spec.initial_condition),
                ic,
                multiplier,
                ts_vector,
            )
        else
            constraint_info = EnergyBalanceConstraintInfo(;
                component_name = dev_name,
                efficiency_data = get_efficiency(device, spec.initial_condition),
                ic_energy = ic,
            )
        end
        constraint_infos[ix] = constraint_info
    end

    spec.constraint_func(
        optimization_container,
        DeviceEnergyBalanceConstraintSpecInternal(
            constraint_infos,
            spec.constraint_name,
            spec.energy_variable,
            spec.pin_variable_names,
            spec.pout_variable_names,
            spec.parameter_name === nothing ? nothing :
            UpdateRef{T}(spec.parameter_name, spec.forecast_label),
        ),
    )
    return
end

@doc raw"""
Constructs multi-timestep constraint from initial condition, efficiency data, and variable tuple
# Constraints
If t = 1:
``` varenergy[name, 1] == initial_conditions[ix].value + varin[name, 1]*eff_in*fraction_of_hour - varout[name, 1]*fraction_of_hour/eff_out ```
If t > 1:
``` varenergy[name, t] == varenergy[name, t-1] + varin[name, t]*eff_in*fraction_of_hour - varout[name, t]*fraction_of_hour/eff_out ```
# LaTeX
`` x^{energy}_1 == x^{energy}_{init} + frhr \eta^{in} x^{in}_1 - \frac{frhr}{\eta^{out}} x^{out}_1, \text{ for } t = 1 ``
`` x^{energy}_t == x^{energy}_{t-1} + frhr \eta^{in} x^{in}_t - \frac{frhr}{\eta^{out}} x^{out}_t, \forall t \geq 2 ``
# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* inputs::Vector{DeviceEnergyBalanceConstraintSpecInternal} : stores constraint information
"""
function energy_balance!(
    optimization_container::OptimizationContainer,
    inputs::DeviceEnergyBalanceConstraintSpecInternal,
)
    time_steps = model_time_steps(optimization_container)
    resolution = model_resolution(optimization_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR
    names = [get_component_name(x) for x in inputs.constraint_infos]

    varenergy = get_variable(optimization_container, inputs.energy_variable)
    varin = [get_variable(optimization_container, var) for var in inputs.pin_variable_names]
    varout =
        [get_variable(optimization_container, var) for var in inputs.pout_variable_names]

    constraint = add_cons_container!(
        optimization_container,
        inputs.constraint_name,
        names,
        time_steps,
    )

    for info in inputs.constraint_infos
        eff_in = info.efficiency_data.in
        eff_out = info.efficiency_data.out
        name = get_component_name(info)
        # Create the PGAE outside of the constraint definition
        !isnothing(info.timeseries) ? ts_value = info.timeseries[1] * info.multiplier :
        ts_value = 0.0
        expr = JuMP.AffExpr(0.0)
        for var in varin
            JuMP.add_to_expression!(expr, var[name, 1], eff_in * fraction_of_hour)
        end
        for var in varout
            JuMP.add_to_expression!(expr, var[name, 1], -1.0 * fraction_of_hour / eff_out)
        end
        constraint[name, 1] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            varenergy[name, 1] == info.ic_energy.value + expr + ts_value
        )
    end

    for t in time_steps[2:end], info in inputs.constraint_infos
        eff_in = info.efficiency_data.in
        eff_out = info.efficiency_data.out
        name = get_component_name(info)
        !isnothing(info.timeseries) ? ts_value = info.timeseries[t] * info.multiplier :
        ts_value = 0.0
        expr = JuMP.AffExpr(0.0, varenergy[name, t - 1] => 1.0)
        for var in varin
            JuMP.add_to_expression!(expr, var[name, t], eff_in * fraction_of_hour)
        end
        for var in varout
            JuMP.add_to_expression!(expr, var[name, t], -1.0 * fraction_of_hour / eff_out)
        end
        constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            varenergy[name, t] == expr + ts_value
        )
    end

    return
end

@doc raw"""
Constructs multi-timestep constraint from initial condition, efficiency data, and variable tuple
# Constraints
If t = 1:
``` varenergy[name, 1] == initial_conditions[ix].value + varin[name, 1]*eff_in*fraction_of_hour - varout[name, 1]*fraction_of_hour/eff_out ```
If t > 1:
``` varenergy[name, t] == varenergy[name, t-1] + varin[name, t]*eff_in*fraction_of_hour - varout[name, t]*fraction_of_hour/eff_out ```
# LaTeX
`` x^{energy}_1 == x^{energy}_{init} + frhr \eta^{in} x^{in}_1 - \frac{frhr}{\eta^{out}} x^{out}_1, \text{ for } t = 1 ``
`` x^{energy}_t == x^{energy}_{t-1} + frhr \eta^{in} x^{in}_t - \frac{frhr}{\eta^{out}} x^{out}_t, \forall t \geq 2 ``
# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* inputs::Vector{DeviceEnergyBalanceConstraintSpecInternal} : stores constraint information
"""
function energy_balance_param!(
    optimization_container::OptimizationContainer,
    inputs::DeviceEnergyBalanceConstraintSpecInternal,
)
    time_steps = model_time_steps(optimization_container)
    resolution = model_resolution(optimization_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR
    names = [get_component_name(x) for x in inputs.constraint_infos]
    has_parameter_data = !isnothing(inputs.param_reference)

    varenergy = get_variable(optimization_container, inputs.energy_variable)
    varin = [get_variable(optimization_container, var) for var in inputs.pin_variable_names]
    varout =
        [get_variable(optimization_container, var) for var in inputs.pout_variable_names]

    if has_parameter_data
        container = add_param_container!(
            optimization_container,
            inputs.param_reference,
            names,
            time_steps,
        )
        multiplier = get_multiplier_array(container)
        param = get_parameter_array(container)
    end

    constraint = add_cons_container!(
        optimization_container,
        inputs.constraint_name,
        names,
        time_steps,
    )

    for info in inputs.constraint_infos
        eff_in = info.efficiency_data.in
        eff_out = info.efficiency_data.out
        name = get_component_name(info)
        # Create the PGAE outside of the constraint definition
        expr = zero(PGAE)
        for var in varin
            JuMP.add_to_expression!(expr, var[name, 1], eff_in * fraction_of_hour)
        end
        for var in varout
            JuMP.add_to_expression!(expr, var[name, 1], -1 * fraction_of_hour / eff_out)
        end
        if has_parameter_data
            multiplier[name, 1] = info.multiplier
            param[name, 1] =
                PJ.add_parameter(optimization_container.JuMPmodel, info.timeseries[1])
            JuMP.add_to_expression!(expr, param[name, 1], multiplier[name, 1])
        end
        constraint[name, 1] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            varenergy[name, 1] == info.ic_energy.value + expr
        )
    end

    for t in time_steps[2:end], info in inputs.constraint_infos
        eff_in = info.efficiency_data.in
        eff_out = info.efficiency_data.out
        name = get_component_name(info)
        expr = expr = zero(PGAE)
        JuMP.add_to_expression!(expr, varenergy[name, t - 1])
        for var in varin
            JuMP.add_to_expression!(expr, var[name, t], eff_in * fraction_of_hour)
        end
        for var in varout
            JuMP.add_to_expression!(expr, var[name, t], -1 * fraction_of_hour / eff_out)
        end
        if has_parameter_data
            multiplier[name, t] = info.multiplier
            param[name, t] =
                PJ.add_parameter(optimization_container.JuMPmodel, info.timeseries[t])
            JuMP.add_to_expression!(expr, param[name, t], multiplier[name, t])
        end
        constraint[name, t] =
            JuMP.@constraint(optimization_container.JuMPmodel, varenergy[name, t] == expr)
    end

    return
end