text:
using JuMP
using Cbc

struct DenseSync
    N::Int64
    p::Int64
    μ::Rational
    x::Array
end

Base.show(io::IO, ds::DenseSync) = print(io, "DenseSync{\n system size: ", ds.N, ", p-twisted state: ", ds.p, ", \n μ = ", ds.μ, " ≃ ", float(ds.μ), "\n x = ", ds.x, "\n}")

function max_connectivity(N::Int64, p::Int64, ε = 10^(-5), limit_sec = 60)
    function L_kl(k, l)
        return cos(2π*p*l/N)*(-1+cos(2*π*k*l/N))
    end
    L = [L_kl(k,l) for k in 1:N-1, l in 1:N-1]
    # model
    model = Model(Cbc.Optimizer)
    set_time_limit_sec(model, limit_sec)
    set_silent(model)
    # variable
    @variable(model, x[i=1:N-1], Bin)
    # maximizing function
    ex = @expression(model, sum(x))
    @objective(model, Max, ex)
    # constraint
    @constraint(model, L*x .≤ -ε) # negative eigenvalues
    @constraint(model, [i=1:(N-1)÷2], x[i]==x[N-i]) # undirected constraint
    # solve!!
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
        μ = objective_value(model)
        μ = Int(μ)//(N-1)
        ans = [Int(round(a)) for a in value.(x)]
        time = solve_time(model)
        @info "maximum connectivity obtained by direct optimization" N p μ=string(μ, " ≃ ",round(μ, digits=7)) time=string(round(time, digits=3), "[s]")
        return DenseSync(N, p, μ, ans)
    else
        μ = -1
        @info "maximum connectivity obtained by direct optimization" N p μ
        return DenseSync(N, p, μ, zeros(Int, N-1))
    end
end

function exact_sol(N::Int64, p::Int64)
    @assert gcd(N,p)==p
    m = p
    n = N ÷ m
    @assert n ≥ 5
    b(x) = -cos(x)*(1-cos(x))
    # compute critical index kc and number of indecies additonally set to one
    kc, add = begin
        xs = [2π*l/n for l in 1:n-1]
        ys = b.(xs)
        sk = cumsum(ys)
        critical_index = findfirst(x -> x ≥ 0, sk)
        additioal = 2*Int(ceil(-m*sk[critical_index-1]/(sk[critical_index]-sk[critical_index-1])-1))
        critical_index, additioal
    end
    # calcilate maximum connectivity
    connects = 2*m*(kc-1) + m-1 + additioal
    μ = connects//(N-1)
    # next, we obtain optimal solution
    x = zeros(Int, N-1)
    # step 1
    for i in 1:m
        for j in 1:kc-1
            x[(i-1)*n + j] = 1
            x[i*n - kc + j] = 1
        end
    end
    # step 2
    for i in 1:m-1
        x[i*n] = 1
    end
    # step 3
    candidates = Int[]
    for i in 1:m
        append!(candidates, (i-1)*n+kc)
        append!(candidates, i*n-kc)
    end
    for i in 1:add÷2
        x[candidates[i]] = 1
        x[N-candidates[i]] = 1
    end
    @info "maximum connectivity obtained by Theorem 1" N p μ=string(μ, " ≃ ",round(μ, digits=7))
    return DenseSync(N, p, μ, x)
end

max_stars_repo_name: yonesuke/DenseSync

max_issues_repo_name: yonesuke/DenseSync

max_forks_repo_name: yonesuke/DenseSync


text:
function n_depth_increments(depths)
  difference = depths[2:1:end] - depths[1:1:end-1]

  n_increments = 0
  for diff in difference
    if diff > 0
      n_increments += 1
    end
  end
  return n_increments
end

function three_measurement_window(depths)
  window = depths[1:1:end-2] + depths[2:1:end-1] + depths[3:1:end]
  difference = window[2:1:end] - window[1:1:end-1]

  n_increments = 0
  for diff in difference
    if diff > 0
      n_increments += 1
    end
  end
  return n_increments
end

# Test case
# Each entry is the sea floor depth in increasing distance from sub.
test_input = [199, 200, 208, 210, 200, 207, 240, 269, 260, 263]
n_test = n_depth_increments(test_input)
println("Number of increments in the test case: $n_test")

n_test_three_window = three_measurement_window(test_input)
println("Number of increments in the three window test case: $n_test_three_window")

# Application to actual input.
open("input.txt", "r") do f
  data = [parse(Int, s) for s in readlines(f)]
  n_increments = n_depth_increments(data)
  println("Number of increments: $n_increments")

  n_three_window = three_measurement_window(data)
  println("Number of increments in the three window case: $n_three_window")
end

ext: jl


text:
"""
This is a generalised DEB model. It was developed for plant modelling, but can potentially
model any organisms and symbioses.

This model can also be run in microclimates provided by the NicheMapr R package, and
can use wide a range of photosynthesis and stomatal conductance formulations from
[Photosynthesis.jl](https://github.com/rafaqz/Photosynthesis.jl).

It is also an in-progress attempt at using Julia's multiple-dispatch methods to abstract and generalise DEB theory and maintain a short, maintainable codebase
for multiplt models - potentially any organism.  Code is adapted from the original [DEBtool](https://github.com/add-my-pet/DEBtool_M) plant model by Bas Kooijman.  """
module DynamicEnergyBudgets

using Unitful,
      OrdinaryDiffEq,
      DocStringExtensions,
      MacroTools,
      Distributions,
      LabelledArrays,
      SimpleRoots,
      Mixers,
      FieldMetadata,
      FieldDefaults,
      Photosynthesis,
      Microclimate,
      Apply,
      Flatten

using Unitful: °C, K, Pa, kPa, MPa, J, kJ, W, L, g, kg, cm, m, s, hr, d, mol, mmol, μmol, σ, R
using Base: tail

import FieldDefaults: get_default
import FieldMetadata: @prior, @default, @description, @units, @limits, @logscaled, @flattenable, @plottable, @selectable,
                      prior, default, description, units, limits, logscaled, flattenable, plottable, selectable
import Photosynthesis: potential_dependence

export tempcorr,
       rate_bracket,
       rate_formula,
       catabolic_fluxes,
       half_saturation,
       stoich_merge,
       synthesizing_unit

export integrate,
       get_state1_names,
       init_state,
       describe,
       paramrange,
       runmodel!,
       debmodel!,
       apply_environment!

export AbstractAssimilation,
       AbstractNAssim, NH4_NO3Assim, NAssim, ConstantNAssim,
       AbstractCAssim, C3Photosynthesis,
       KooijmanPhotosynthesis, KooijmanSLAPhotosynthesis, KooijmanWaterPotentialPhotosynthesis,
       KooijmanNH4_NO3Assim, FvCBPhotosynthesis, ConstantCAssim

export AbstractSynthesizingUnit, ParallelComplementarySU, MinimumRuleSU, KfamilySU

export AbstractRate, SimpleRate, FZeroRate

export AbstractStateFeedback, LosslessAutophagy, StructuralLossAutophagy, DissipativeAutophagy

export AbstractTemperatureCorrection, TempCorr, TempCorrLower, TempCorrLowerUpper, ParentTardieu

export AbstractScaling, KooijmanArea

export AbstractMaturity, Maturity

export AbstractMaintenance, Maintenance

export AbstractDEBCore, DEBCore

export AbstractGermination, ThresholdGermination

export AbstractProduction, Production

export AbstractShape, Isomorph, V0morph, V1morph, V1V0morph, Plantmorph

export AbstractRejection, LosslessRejection, DissipativeRejection

export AbstractCatabolism, Catabolism, CatabolismCN, CatabolismCNshared, CatabolismCNE

export AbstractTranslocation, LosslessMultipleTranslocation, DissipativeMultipleTranslocation,
       LosslessTranslocation, DissipativeTranslocation

export AbstractAllometry, Allometry, SqrtAllometry, FixedAllometry

export AbstractParams, Params, SharedParams,
       Vars, CarbonVars, NitrogenVars, Records,
       AbstractOrgan, Organ,
       AbstractOrganism, Plant

# Auto docstrings
@template TYPES =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

# Field metadata columns
@chain columns @description @logscaled @limits @prior @units @udefault_kw

const STATELEN = 6

const STATE = (:P, :V, :M, :C, :N, :E)
const STATE1 = (:EE, :CN, :C, :N, :E)
const TRANS = (:ass, :gro, :mai, :mat, :rej, :tra, :fbk)
const TRANS1 = (:ctb, :rej, :los)
const BI_XTOL = 1e-10
const BI_MAXITER = 100

include("traits.jl")
include("components/synthesizing_units.jl")
include("components/temperature_correction.jl")
include("components/rate.jl")
include("components/shape.jl")
include("components/allometry.jl")
include("components/autophagy.jl")
include("components/assimilation.jl")
include("components/germination.jl")
include("components/maturity.jl")
include("components/production.jl")
include("components/translocation.jl")
include("components/catabolism.jl")
include("components/maintenance.jl")
include("types.jl")
include("environment.jl")
include("getters.jl")
include("functions.jl")
include("setup.jl")
include("model.jl")
include("aliases.jl")

end # module


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

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

# ╔═╡ 45f1331e-59b1-11eb-2bc0-c138b1421b55
begin
	using Pkg
	Pkg.activate(".")
	Pkg.instantiate()
end

# ╔═╡ fc1a9e20-56aa-11eb-29b1-3be8043c6ede
begin
	using Random
	Random.seed!(0)
	using Zygote
	using LinearAlgebra
	using Plots
	using PlutoUI
	using Images
	using LaTeXStrings
end

# ╔═╡ 60a23536-59d0-11eb-368a-5588b71a6093
md"# Curve fitting using AD"

# ╔═╡ 832e2c5a-59ce-11eb-32d7-85b05440f919
md"## Problem Formulation"

# ╔═╡ e35400cc-56ad-11eb-2794-592020c6ac7e
md"### The underlying mathematics of the problem"

# ╔═╡ 761fc420-57af-11eb-24b3-9728b06c482e
md"We are trying to fit a function $f_{org}(x)$ to a polynomial function $f_{ML}(x)$"

# ╔═╡ 021f6506-57b1-11eb-2179-7b6f153fce2f
md"$$f_{ML}(x_i) = \displaystyle\sum\limits_{j=0}^nw_jx_i^j = w_0 + x_iw_1 + x_i^2w_2+\ldots+x_i^nw_n$$"

# ╔═╡ df0e42cc-59be-11eb-367e-cf25ab5f105f
md"so that $f_{org}(x_i) \approx f_{ML}(x_i)$"

# ╔═╡ 2730398c-59bf-11eb-060e-276c0138395e
md"We define $X, \textbf{w}, \textbf{y}$ for a more concise and all-round understanding of the problem,"

# ╔═╡ e50c3cc4-59b7-11eb-1479-8bba33d251fd
md"$\begin{gather}
X=
\begin{bmatrix}
1 & x_1 & \ldots & x_1^n\\
\vdots & & \ddots &\vdots \\
1 & x_m & & x_m^n
\end{bmatrix}
\quad\quad\quad
\textbf{w} = \begin{bmatrix}
w_o \\
w_1 \\
\vdots\\
w_n
\end{bmatrix}
\quad\quad\quad
\textbf{y} =
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_m
\end{bmatrix}
\end{gather}$"

# ╔═╡ 51545508-59b9-11eb-3f94-59e743bfd1a8
md" where,
- The matrix $X$ is the data matrix. $X$ is a $(m\times (n+1)) matrix$. $m$ is the number of data points and $n$ is the n is the degree of the plyonomial. The choice of this polynomial becomes clear in the next step.
- The vector $\textbf{w}$ represents the weights vector. $w_0$ the bias. $w_j$ is the weight corresponding to $x^j$.
- The vector $\textbf{y}$ represents the results vector. $y_i$ is the result of $x_i$ on applying the function $f$"

# ╔═╡ 70b0cea0-59ba-11eb-2ae9-439069c94467
md"$\begin{gather}
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_m
\end{bmatrix}
=
\begin{bmatrix}
1 & x_1 & \ldots & x_1^n\\
\vdots & & \ddots &\vdots \\
1 & x_m & & x_m^n
\end{bmatrix}
\odot
\begin{bmatrix}
w_o \\
w_1 \\
\vdots\\
w_n
\end{bmatrix}
\end{gather}$"

# ╔═╡ ea14bfa4-59ba-11eb-35a1-6548f64bed7e
md"$\begin{gather}
\quad \quad \
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_m
\end{bmatrix}
=
\begin{bmatrix}
w_o + x_1\times w_1 &\ldots  &+x_1^n \times w_n\\
\vdots \\
w_0 + x_m\times w_1  &\ldots &+x_m^n\times w_n
\end{bmatrix}
\end{gather}$"

# ╔═╡ 32bb7c2c-59bc-11eb-0914-d99847d081a3
load(download("https://thehritikjain.files.wordpress.com/2020/05/1_e4_ptjctmaofsrpzczbv-g-1.jpeg?w=935"))

# ╔═╡ d2ea16ee-57b7-11eb-2d27-f3c04343dc59
md"Before we approach the real problem of solving the approximate function, we need a predicting function and loss function to evaluate good our approximation is."

# ╔═╡ fa7fe6b4-59c1-11eb-2486-d90418eaee4a
md"### Predict and loss functions : "

# ╔═╡ e8864a02-59c1-11eb-12ff-c9a4ffde849d
md"#### Predict function "

# ╔═╡ bf75963a-56ac-11eb-2b34-21b724a59f96
function predict(x, ω::Vector)
	n = length(ω)
	x_vec = [x^i for i= 0:n-1]
	return dot(x_vec, ω)
end

# ╔═╡ c8a25b68-59c1-11eb-363f-0772ec4edbf2
md"**Usage Example :**"

# ╔═╡ 10be275c-59c1-11eb-10dc-4becfaf27f8d
begin
	# Try changing the values
	local ω = rand(5) # 4 polynomial weight vector
	local x = 3. # input
	predict(x, ω)
end

# ╔═╡ d928b270-59c1-11eb-09d0-a9c3c878fe2d
md"#### Square Loss function "

# ╔═╡ 0129d404-57b8-11eb-3dde-efa8a687c1d4
function sq_loss(func, X, ω::Vector)
	n = length(ω)
	y = func.(X)
	ŷ = [predict(x, ω) for x in X]
	return sum((y-ŷ).^2)/n
end

# ╔═╡ bc758f7c-59c1-11eb-3f4a-5b987c7c6e0b
md"**Usage Example :**"

# ╔═╡ 55bd943c-59c1-11eb-19d3-a7bba6fcbd30
begin
	local ω = rand(5) # 4 polynomial weight vector
	local X = [3.] # input
	local loss = sq_loss(sin, X, ω)
end

# ╔═╡ 00f8133c-57be-11eb-19d0-33006b7bbb06
md"## Data Generation"

# ╔═╡ 6d407f0e-59ce-11eb-0ed6-33f817c914a3
md"Choose the number of data points : "

# ╔═╡ 4b28aee8-56ac-11eb-3be7-0de4e9cb211b
@bind n Slider(0:10, show_value = true, default=8)

# ╔═╡ e6d69724-57bf-11eb-0337-db41600a8ebc
md"We define a function that we want to approximate. **org\_func** is the original function. We are trying to approximate *sin(x)* here. You can try your own function that you want to approximat inside org\_func and chnages will be reflected in rest of the notebook."

# ╔═╡ 9801c9ca-57bf-11eb-3f30-93055ccd6a0f
function org_func(x)
	return sin(x)
end

# ╔═╡ 2fcaeb98-56b1-11eb-0198-a3a5088dc007
begin
	X = [x for x in range(randn(), 12, length = n)]
	y = org_func.(X) +0.1* randn((n))
	plot(-0.5:0.01:12, org_func.(-0.5:0.01:12), label = L"sin(x)")
	scatter!(X, y, label = "data points", title = "Our Data points")
end

# ╔═╡ bd151b36-56b1-11eb-03d1-094237a3bafa
md"## Approximation"

# ╔═╡ c7c2e3f2-57be-11eb-15c0-dbe72c95d654
md"We will try to solve the problem of approximating the function by two differnt ways:
- Linear Algebra Methods
- Gradient Descent using AD"

# ╔═╡ 377b50e4-57bf-11eb-1741-a73fda781579
md" ### Approximation by Linear Algebra"

# ╔═╡ 50c00e44-59c4-11eb-2726-27bee50e4fdd
md"$\begin{gather}
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_m
\end{bmatrix}
=
\begin{bmatrix}
1 & x_1 & \ldots & x_1^n\\
\vdots & & \ddots &\vdots \\
1 & x_m & & x_m^n
\end{bmatrix}
\odot
\begin{bmatrix}
w_o \\
w_1 \\
\vdots\\
w_n
\end{bmatrix}
\end{gather}$"

# ╔═╡ 586240fe-59c4-11eb-1427-312c47374f8a
md"$
\ \ \ \
\begin{gather}
\begin{bmatrix}
w_o \\
w_1 \\
\vdots\\
w_n
\end{bmatrix}
=
\begin{bmatrix}
1 & x_1 & \ldots & x_1^n\\
\vdots & & \ddots &\vdots \\
1 & x_m & & x_m^n
\end{bmatrix}^{-1}
\odot
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_m
\end{bmatrix}

\end{gather}$"

# ╔═╡ 4f5c808e-57bf-11eb-2a7a-3995d6f8f866
md"This may seem like a very viable method, but try considering the case where $X$ is invertible. It is left as an exercise to the reader."

# ╔═╡ 5dd40e02-57bf-11eb-2764-0f3f2b95517b
md" ### Approximation by Gradient Descent"

# ╔═╡ ea205094-59c9-11eb-0f04-818f051b7d1e
md"We will be using gradient descent algorithm to find $\textbf{w}$ here. We are tired of mathematics by now and we somehow want the differentaition to happen automatically(almost magically). We will be using **Automatic Differentation**(or a more intuitive name : **Algorthimic Differentiation**). While how automatic differintation works is another day's agenda, we can have a little taste of magic of what AD does for us."

# ╔═╡ b13969a4-59ca-11eb-15d8-0f9b5f225fed
md"The gradient of a function $f(x)$ at $x_0$ is defined as $\displaystyle\frac{df(x)}{dx}\Bigg|_{x=x_0}$

$f^\prime(x_0) = \displaystyle\frac{df(x)}{dx}\Bigg|_{x=x_0}$"

# ╔═╡ e4f40500-59cb-11eb-1020-3b345469fb1f
md"The Zygote package gives us a *gradient* API that we can use to find gradients"

# ╔═╡ 88002972-59cc-11eb-3d75-ab7e7bb45954
md"$f(x) = x^2+4x+1$
$f^\prime(x) = 2x+4$
$f^\prime(5) = 2\times 5 + 4 = 14$"

# ╔═╡ 3255da08-59cc-11eb-0bdc-61a05344ad50
gradient(x->x^2+4x+1,  5)

# ╔═╡ cd6a28f8-59cc-11eb-28f8-35c3a04414d8
md"$f(x, y) = x^2+y^2+3xy$
$\displaystyle\frac{\partial f(x,y)}{\partial x} = 2x+3y$
$\displaystyle\frac{\partial f(x,y)}{\partial y} = 2y+3x$"

# ╔═╡ 4506b2d0-59cc-11eb-22e2-5f43abeb3ccd
gradient( (x,y) -> x^2 + y^2 + 3x*y, 2,3)

# ╔═╡ 608cc19c-59cd-11eb-092d-558d322966c6
md"A train step function that returns the updated weights vector $\textbf{w}$"

# ╔═╡ 84566232-57bf-11eb-390b-f33370fbd9a3
function train_step(func, X, ω, η)::typeof(ω)
	grad = gradient(ω->sq_loss(func, X, ω), ω)[1]
	ω -= η * grad
	return ω
end

# ╔═╡ de0aeb62-59cd-11eb-3e5f-a7cb176093c9
md"Choose the degree of polynomial of  $f_{ML}(x)$. You may need to choose different $\eta$ and number of epochs to account for the model change"

# ╔═╡ ae1c0e18-59cd-11eb-2e61-fbe810bd5fa0
@bind n_model Slider(1:1:7, default = 4, show_value = true)

# ╔═╡ 5c85511c-5811-11eb-0f2e-2352309f2261
## Intializing weights
ω = rand(n_model+1)

# ╔═╡ 11337388-57c1-11eb-3418-3979d91e6af2
function train(func, X, ω, η, epochs)::typeof(ω)
	for i = 1:1:epochs
		ω =  train_step(func, X, ω, η)
		loss = sq_loss(func, X, ω)
		if isnan(loss) || isinf(loss) || iszero(loss)
			throw(DomainError("Loss cannot take $(loss) value"))
		end
	end
	return ω
end

# ╔═╡ d722598c-57c2-11eb-3452-531fb6522e34
ω_trained = train(org_func, X, ω, 7e-9, 1e6)

# ╔═╡ 11d38860-59c6-11eb-06e2-b10697f21f1d
md"It is upto the user to decide how much accuracy he wants given the time constraints"

# ╔═╡ 8b6bc65e-57c7-11eb-372f-bd6f700363e1
loss = sq_loss(org_func, X, ω_trained)

# ╔═╡ 328c4584-57c4-11eb-3cf5-a57e033cf574
pred = [predict(x, ω_trained) for x in -0.5:0.01:13]

# ╔═╡ a0dfe352-57c3-11eb-348f-8380ee6b3caf
begin
	plot(-0.5:0.01:13, org_func.(-0.5:0.01:13), label = L"f_{org}(x)")
	scatter!(X, y, label = "Orginal data points")
	ŷ = [predict(x, ω_trained) for x in X]
	scatter!(X, ŷ, label = "Predicted data points")
	plot!(-0.5:0.01:13, pred, label = L"f_{ML}(x)")
end

# ╔═╡ Cell order:
# ╟─45f1331e-59b1-11eb-2bc0-c138b1421b55
# ╟─fc1a9e20-56aa-11eb-29b1-3be8043c6ede
# ╟─60a23536-59d0-11eb-368a-5588b71a6093
# ╟─832e2c5a-59ce-11eb-32d7-85b05440f919
# ╟─e35400cc-56ad-11eb-2794-592020c6ac7e
# ╟─761fc420-57af-11eb-24b3-9728b06c482e
# ╟─021f6506-57b1-11eb-2179-7b6f153fce2f
# ╟─df0e42cc-59be-11eb-367e-cf25ab5f105f
# ╟─2730398c-59bf-11eb-060e-276c0138395e
# ╟─e50c3cc4-59b7-11eb-1479-8bba33d251fd
# ╟─51545508-59b9-11eb-3f94-59e743bfd1a8
# ╟─70b0cea0-59ba-11eb-2ae9-439069c94467
# ╟─ea14bfa4-59ba-11eb-35a1-6548f64bed7e
# ╟─32bb7c2c-59bc-11eb-0914-d99847d081a3
# ╟─d2ea16ee-57b7-11eb-2d27-f3c04343dc59
# ╟─fa7fe6b4-59c1-11eb-2486-d90418eaee4a
# ╟─e8864a02-59c1-11eb-12ff-c9a4ffde849d
# ╠═bf75963a-56ac-11eb-2b34-21b724a59f96
# ╟─c8a25b68-59c1-11eb-363f-0772ec4edbf2
# ╠═10be275c-59c1-11eb-10dc-4becfaf27f8d
# ╟─d928b270-59c1-11eb-09d0-a9c3c878fe2d
# ╠═0129d404-57b8-11eb-3dde-efa8a687c1d4
# ╟─bc758f7c-59c1-11eb-3f4a-5b987c7c6e0b
# ╠═55bd943c-59c1-11eb-19d3-a7bba6fcbd30
# ╟─00f8133c-57be-11eb-19d0-33006b7bbb06
# ╟─6d407f0e-59ce-11eb-0ed6-33f817c914a3
# ╟─4b28aee8-56ac-11eb-3be7-0de4e9cb211b
# ╟─e6d69724-57bf-11eb-0337-db41600a8ebc
# ╠═9801c9ca-57bf-11eb-3f30-93055ccd6a0f
# ╠═2fcaeb98-56b1-11eb-0198-a3a5088dc007
# ╟─bd151b36-56b1-11eb-03d1-094237a3bafa
# ╟─c7c2e3f2-57be-11eb-15c0-dbe72c95d654
# ╟─377b50e4-57bf-11eb-1741-a73fda781579
# ╟─50c00e44-59c4-11eb-2726-27bee50e4fdd
# ╟─586240fe-59c4-11eb-1427-312c47374f8a
# ╟─4f5c808e-57bf-11eb-2a7a-3995d6f8f866
# ╟─5dd40e02-57bf-11eb-2764-0f3f2b95517b
# ╟─ea205094-59c9-11eb-0f04-818f051b7d1e
# ╟─b13969a4-59ca-11eb-15d8-0f9b5f225fed
# ╟─e4f40500-59cb-11eb-1020-3b345469fb1f
# ╟─88002972-59cc-11eb-3d75-ab7e7bb45954
# ╠═3255da08-59cc-11eb-0bdc-61a05344ad50
# ╟─cd6a28f8-59cc-11eb-28f8-35c3a04414d8
# ╠═4506b2d0-59cc-11eb-22e2-5f43abeb3ccd
# ╟─608cc19c-59cd-11eb-092d-558d322966c6
# ╠═84566232-57bf-11eb-390b-f33370fbd9a3
# ╟─de0aeb62-59cd-11eb-3e5f-a7cb176093c9
# ╟─ae1c0e18-59cd-11eb-2e61-fbe810bd5fa0
# ╠═5c85511c-5811-11eb-0f2e-2352309f2261
# ╠═11337388-57c1-11eb-3418-3979d91e6af2
# ╠═d722598c-57c2-11eb-3452-531fb6522e34
# ╟─11d38860-59c6-11eb-06e2-b10697f21f1d
# ╠═8b6bc65e-57c7-11eb-372f-bd6f700363e1
# ╠═328c4584-57c4-11eb-3cf5-a57e033cf574
# ╟─a0dfe352-57c3-11eb-348f-8380ee6b3caf

max_issues_repo_name: chaitak-gorai/Study-Group

max_forks_repo_name: chaitak-gorai/Study-Group


text:
module TestWithSequential

using FoldsThreads
using FoldsThreads.Implementations: SpawnAllEx
using Folds.Testing: test_with_sequential

executors = [
    NondeterministicEx(),
    NondeterministicEx(ntasks = 2),
    NondeterministicEx(basesize = 3),
    NondeterministicEx(basesize = 3, ntasks = 2),
    TaskPoolEx(),
    TaskPoolEx(ntasks = 2),
    TaskPoolEx(basesize = 3),
    TaskPoolEx(basesize = 3, ntasks = 2),
    WorkStealingEx(),
    WorkStealingEx(basesize = 3),
    DepthFirstEx(),
    DepthFirstEx(basesize = 3),
    SpawnAllEx(),
    SpawnAllEx(basesize = 3),
    NonThreadedEx(),
    NonThreadedEx(basesize = 3),
]

if Threads.nthreads() > 1
    append!(
        executors,
        [
            TaskPoolEx(background = true, basesize = 3),
            TaskPoolEx(background = true, basesize = 3, ntasks = 2),
        ],
    )
end

test_with_sequential(executors)

end


text:
function save_noise!(W::VirtualBrownianTree)

end

function interpolate!(W::VirtualBrownianTree,u,p,t)
  ts,timeseries,timeseries2 = W.t,W.W,W.Z
  sign(W.dt)*t > sign(W.dt)*ts[end] && error("Bridge cannot extrapolate past the final timepoint. Build a longer VBT to cover the integration.")
  sign(W.dt)*t < sign(W.dt)*ts[1] && error("Bridge cannot extrapolate before the first timepoint. Build a longer VBT to cover the integration.")
  tdir = sign(ts[end]-ts[1])

  if t isa Union{Rational,Integer}
    @inbounds i = searchsortedfirst(ts,t,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  else
    @inbounds i = searchsortedfirst(ts,t-tdir*10eps(typeof(t)),rev=tdir<0)
  end

  @inbounds if (t isa Union{Rational,Integer} && ts[i] == t) || (isapprox(t, ts[i]; atol = 100eps(typeof(t)), rtol = 100eps(t)))
    val1 = timeseries[i]
    W.dZ != nothing ? val2 = timeseries2[i] : val2 = nothing
  elseif ts[i-1] == t # Can happen if it's the first value!
    val1 = timeseries[i-1]
    W.dZ != nothing ? val2 = timeseries2[i-1] : val2 = nothing
  else
    if W.dZ != nothing
      val1, val2 = search_VBT(t,W.seeds[i-1],ts[i-1],ts[i],timeseries[i-1],timeseries[i],
        timeseries2[i-1],timeseries2[i], W, W.rng)
    else
      val1, val2 = search_VBT(t,W.seeds[i-1],ts[i-1],ts[i],timeseries[i-1],timeseries[i],
        nothing,nothing, W, W.rng)
    end
  end
  val1,val2
end

function interpolate!(out1,out2,W::VirtualBrownianTree,u,p,t)
  ts,timeseries,timeseries2 = W.t,W.W,W.Z
  sign(W.dt)*t > sign(W.dt)*ts[end] && error("Bridge cannot extrapolate past the final timepoint. Build a longer VBT to cover the integration.")
  sign(W.dt)*t < sign(W.dt)*ts[1] && error("Bridge cannot extrapolate before the first timepoint. Build a longer VBT to cover the integration.")
  tdir = sign(ts[end]-ts[1])

  if t isa Union{Rational,Integer}
    @inbounds i = searchsortedfirst(ts,t,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  else
    @inbounds i = searchsortedfirst(ts,t-tdir*10eps(typeof(t)),rev=tdir<0)
  end

  @inbounds if (t isa Union{Rational,Integer} && ts[i] == t) || (isapprox(t, ts[i]; atol = 100eps(typeof(t)), rtol = 100eps(t)))
    copyto!(out1,timeseries[i])
    W.dZ != nothing && copyto!(out2,timeseries2[i])
  elseif ts[i-1] == t # Can happen if it's the first value!
    copyto!(out1,timeseries[i-1])
    W.dZ != nothing && copyto!(out2,timeseries2[i-1])
  else
    if W.dZ != nothing
      search_VBT!(out1,out2,t,seeds[i-1],ts[i-1],ts[i],timeseries[i-1],timeseries[i],
        timeseries2[i-1],timeseries2[i], W, W.rng)
    else
      search_VBT!(out1,out2,t,seeds[i-1],ts[i-1],ts[i],timeseries[i-1],timeseries[i],
        nothing,nothing, W, W.rng)
    end
  end
  return nothing
end

function calculate_step!(W::VirtualBrownianTree,dt,u,p)
  t = W.curt+dt
  if typeof(t) <: AbstractFloat && abs(t - W.t[end]) < 100eps(typeof(dt))
    t = W.t[end]
  end
  if isinplace(W)
    interpolate!(W.dW,W.dZ,W,u,p,t)
    W.dW .-= W.curW
    if W.dZ != nothing
      W.dZ .-= W.curZ
    end
  else
    new_W, new_Z = interpolate!(W,u,p,t)
    W.dW = new_W - W.curW
    if W.dZ != nothing
      W.dZ = new_Z - W.curZ
    end
  end
  W.dt = dt
  return nothing
end

function accept_step!(W::VirtualBrownianTree,dt,u,p,setup_next=true)
  W.step_setup == false && error("Stepped past the defined domain for the VBT.")

  if isinplace(W)
    W.curW .+= W.dW
  else
    W.curW += W.dW
  end
  W.curt += W.dt
  if W.dZ != nothing
    if isinplace(W)
      W.curZ .+= W.dZ
    else
      W.curZ += W.dZ
    end
  end

  W.dt = dt #dtpropose
  if (W.dt isa Union{Rational,Integer})
    if sign(W.dt)*(W.curt + W.dt) > sign(W.dt)*W.t[end]
      setup_next = false
      W.step_setup = false
    end
  else
    if sign(W.dt)*(W.curt + W.dt) > sign(W.dt)*(W.t[end]+sign(W.dt)*10eps(typeof(dt)))
      setup_next = false
      W.step_setup = false
    end
  end

  if setup_next
    calculate_step!(W,dt,u,p)
  end
  return nothing
end

function reject_step!(W::VirtualBrownianTree,dtnew,u,p)
  calculate_step!(W,dtnew,u,p)
  return nothing
end

function setup_next_step!(W::VirtualBrownianTree,u,p)
  calculate_step!(W,W.dt,u,p)
  return nothing
end

# split seeds

# for counter-based RNGs
function split_VBT_seed(rng::Random123.AbstractR123, parent_seed, current_depth, Nt)

  # seed left
  seed_l = convert(typeof(parent_seed),parent_seed - (Nt-1)/2^(current_depth+1))
  # seed right
  seed_r = convert(typeof(parent_seed),parent_seed + (Nt-1)/2^(current_depth+1))

  RandomNumbers.Random123.set_counter!(rng, parent_seed)
  return seed_l, seed_r, parent_seed
end

# create the cache of the VBT with depth tree_depth
function create_VBT_cache(bridge,t0,W0,Z0,tend,Wend,Zend,rng::Random123.AbstractR123,tree_depth,search_depth)
  # total number of cached time steps and W values
  Nt = Int(2^search_depth+1)

  ts = [t0, tend]
  Ws = [W0, Wend]

  dW = Wend-W0 # for input of bridges

  if Z0 != nothing
    Zs = [Z0, Zend]
  else
    Zs = [nothing]
  end

  seeds = [convert(typeof(rng.ctr1), rng.ctr1+(Nt+1)/2)]

  q = 1//2

  for level in 1:Int(tree_depth)

    new_ts = Vector{typeof(t0)}(undef, 0)
    new_Ws = Vector{typeof(W0)}(undef, 0)
    if Z0 != nothing
      new_Zs = Vector{typeof(Z0)}(undef, 0)
    else
      new_Zs = [nothing]
    end
    new_seeds = Vector{typeof(rng.ctr1)}(undef, 0)
    for (i, parent) in enumerate(seeds)
      seed_l, seed_r, seed_v = split_VBT_seed(rng,parent,level,Nt)
      append!(new_seeds,[seed_l, seed_r])

      t0, t1 = ts[i], ts[i+1]
      W0tmp, W1tmp = Ws[i], Ws[i+1]
      t = (t0 + t1) / 2

      h = t1-t0
      #q = (t-t0)/h # == 1//2 defined above

      w = bridge(dW,nothing,W0tmp,W1tmp,q,h,nothing,nothing,nothing,rng)
      if Z0 != nothing
        Z0tmp, Z1tmp = Zs[i], Zs[i+1]
        z = bridge(dW,nothing,Z0tmp,Z1tmp,q,h,nothing,nothing,nothing,rng)
        append!(new_Zs,Z0)
        append!(new_Zs,z)
      end

      append!(new_ts,[t0,t])
      append!(new_Ws,[W0tmp,w])
    end
    push!(new_ts,tend)
    push!(new_Ws,Wend)

    if Z0 != nothing
      push!(new_Zs,Zend)
      Zs = new_Zs
    end

    ts = new_ts
    Ws = new_Ws

    seeds = new_seeds
  end

  return ts, Ws, Zs, seeds
end

function search_VBT(t,seed,t0,t1,W0,W1,Z0,Z1,W::VirtualBrownianTree,rng::Random123.AbstractR123)

  Nt = Int(2^W.search_depth+1)
  depth = Int(W.tree_depth+1)
  seed_l, seed_r, seed_v = split_VBT_seed(rng,seed,depth,Nt)

  tmid = (t0 + t1) / 2
  h = t1-t0
  q = 1//2

  out1 = W.bridge(W.dW,nothing,W0,W1,q,h,nothing,nothing,nothing,rng)

  if Z0 != nothing
    out2 = W.bridge(W.dW,nothing,Z0,Z1,q,h,nothing,nothing,nothing,rng)
  else
    out2 = nothing
  end

  # create a buffer
  W0tmp, W1tmp = copy(W0), copy(W1)
  if Z0 != nothing
    Z0tmp, Z1tmp = copy(Z0), copy(Z1)
  end

  while (abs(t-tmid) > W.atol && depth<W.search_depth)
    depth += 1
    if t < tmid
      t1 = tmid
      W0tmp, W1tmp = W0tmp, out1
      seed_v = seed_l
      if Z0 != nothing
        Z0tmp, Z1tmp = Z0tmp, out2
      end
    else
      t0 = tmid
      W0tmp, W1tmp = out1, W1tmp
      seed_v = seed_r
      if Z0 != nothing
        Z0tmp, Z1tmp = out2, Z1tmp
      end
    end

    seed_l, seed_r, seed_v = split_VBT_seed(rng,seed_v,depth,Nt)
    tmid = (t0 + t1) / 2
    h = t1-t0

    out1 = W.bridge(W.dW,nothing,W0tmp,W1tmp,q,h,nothing,nothing,nothing,rng)
    if Z0 != nothing
      out2 = W.bridge(W.dW,nothing,Z0tmp,Z1tmp,q,h,nothing,nothing,nothing,rng)
    end
  end

  return out1, out2
end

function search_VBT!(out1,out2,t,seed,t0,t1,W0,W1,Z0,Z1,W::VirtualBrownianTree,rng::Random123.AbstractR123)

  Nt = Int(2^W.search_depth+1)
  depth = Int(W.tree_depth+1)
  seed_l, seed_r, seed_v = split_VBT_seed(rng,seed,depth,Nt)

  tmid = (t0 + t1) / 2
  h = t1-t0
  q = 1//2

  # create a buffer
  copyto!(W.W0tmp,W0)
  copyto!(W.W1tmp,W1)
  if Z0 != nothing
    copyto!(W.Z0tmp,Z0)
    copyto!(W.Z1tmp,Z1)
  end

  W.bridge(out1,W.dW,nothing,W.W0tmp,W.W1tmp,q,h,nothing,nothing,nothing,rng)

  if Z0 != nothing
    W.bridge(out2,W.dW,nothing,W.Z0tmp,W.Z1tmp,q,h,nothing,nothing,nothing,rng)
  end

  while (abs(t-tmid) > W.atol && depth<W.search_depth)
    depth += 1
    if t < tmid
      t1 = tmid
      copyto!(W.W1tmp,out1)
      seed_v = seed_l
      if Z0 != nothing
        copyto!(W.Z1tmp,out2)
      end
    else
      t0 = tmid
      copyto!(W.W0tmp,out1)
      seed_v = seed_r
      if Z0 != nothing
        copyto!(W.Z0tmp,out2)
      end
    end

    seed_l, seed_r, seed_v = split_VBT_seed(rng,seed_v,depth,Nt)
    tmid = (t0 + t1) / 2
    h = t1-t0

    W.bridge(out1,W.dW,nothing,W.W0tmp,W.W1tmp,q,h,nothing,nothing,nothing,rng)
    if Z0 != nothing
      W.bridge(out2,W.dW,nothing,W.Z0tmp,W.Z1tmp,q,h,nothing,nothing,nothing,rng)
    end
  end

  return nothing
end


text:
abstract type AbstractMetaDatum{T<:Any} end

struct MetaDatum{T<:Any} <: AbstractMetaDatum{T}
    id::Symbol
    kind::Symbol
    value::T
end

abstract type AbstractMetaData{T} end

struct MetaDict{T} <: AbstractMetaData{T}
    kind::Symbol
    values::Dict{Symbol,MetaDatum{T}}
    links::Set{Symbol}
end

struct MetaDatums{T} <: AbstractMetaData{T}
    metadicts::Dict{Symbol, MetaDict}
    links::Dict{Symbol,Set{Symbol}}}
end

kind(md::AbstractMetaDatum) = md.kind
value(md::AbstractMetaDatum) = md.value
identifier(md::AbstractMetaDatum) = md.id

kind(md::AbstractMetaData) = md.kind
links(md::AbstractMetaData) = md.links

values(md::MetaDict) = values(md.values)
identifiers(md::MetaDict) = keys(md.ids)

values(md::MetaDatums) = Set([values(d.values)... for d in md.metadicts])
identifiers(md::MetaDatums) = Set([keys(d.values)... for d in md.metadicts])


text:
# `UnionFinder{T <: Integer}` is a graph containing a constant number of nodes
# which allows for union-find operations. All nodes are indexed by an integer
# of type `T` which is between 1 an the number of internal nodes.
mutable struct UnionFinder{T <: Integer}
    sizes :: Vector{T}
    parents :: Vector{T}

    # `UnionFinder(nodes)` returns a `UnionFinder` with `nodes` unconnected
    # internal nodes.
    function UnionFinder{T}(nodes :: T) where T
        if nodes <= 0
            throw(ArgumentError("Non-positive nodes, $nodes."))
        end

        uf = new(Vector{T}(undef, Int(nodes)), Vector{T}(undef, Int(nodes)))
        reset!(uf)

        return uf
    end
end

UnionFinder(nodes :: Integer) = UnionFinder{typeof(nodes)}(nodes)

# `reset(uf)` disconnects all the nodes within `uf`.
function reset!(uf :: UnionFinder)
    for i in 1:length(uf.parents)
        uf.sizes[i] = 1
        uf.parents[i] = i
    end
end

# `union!(uf, iterator)` iterates through `iterator` which returns integer
# edges, (`u`, `v`), and connects them within `uf`. `u` and `v` must be valid
# node indices for `uf`.
function union!(uf :: UnionFinder{T}, iterator) where T <: Integer
    for (u, v) in iterator
        union!(uf, u, v)
    end
end

# `union!(uf, us, vs)` connects nodes within `uf` which are bridged by
# the edges (`us[i]`, `vs[i]`). All values in `us` and `vs` must be valid node
# indices for `uf` and `us` and `vs` must be the same length.
function union!(uf :: UnionFinder{T},
                us :: Vector{T}, vs :: Vector{T}) where T <: Integer
    if length(us) != length(vs)
        throw(ArgumentError("us and vs not of the same length."))
    end

    for i in 1:length(us)
        union!(uf, us[i], vs[i])
    end
end

# `union!(uf, edges)` conncts all nodes within `uf` which are bridged by an
# edge within `edges`. Both vertices for each edge must be valid node indices
# into `uf`.
function union!(uf :: UnionFinder{T},
                edges :: Vector{Tuple{T,T}}) where T <: Integer
    for (u, v) in edges
        union!(uf, u, v)
    end
end

# `union!(uf, node1, node2)` connects the nodes within `uf` with indices
# `node1` and `node2`. `node1` and `node2` must be valid indices into `uf`.
function union!(uf :: UnionFinder{T}, node1 :: T, node2 :: T) where T <: Integer
    root1 = find!(uf, node1)
    root2 = find!(uf, node2)

    # TODO: Test whether using rank or using size is better for performance.
    if root1 == root2
        return
    elseif uf.sizes[root1] < uf.sizes[root2]
        uf.parents[root1] = root2
        uf.sizes[root2] += uf.sizes[root1]
    else
        uf.parents[root2] = root1
        uf.sizes[root1] += uf.sizes[root2]
    end
end

# `find!(uf, node)` returns the group ID of `node`. `node` must be a valid
# index into `uf`.
function find!(uf :: UnionFinder{T}, node :: T) where T <: Integer
    if node > length(uf.parents) || node <= 0
        throw(BoundsError())
    end
    if uf.parents[node] != uf.parents[uf.parents[node]]
        compress!(uf, node)
    end
    return uf.parents[node]
end

# `compress(uf, node)` compresses the internal parental node tree so that
# all nodes between `node` and the root of its group will point directly to the
# root. `node` must be a valid index into `uf`.
#
# Not publicly exported.
function compress!(uf :: UnionFinder{T}, node :: T) where T <: Integer
    child = node
    parent = uf.parents[child]

    while child != parent
        child = parent
        parent = uf.parents[child]
    end
    root = child

    child = node
    parent = uf.parents[child]
    uf.parents[child] = root

    while child != parent
        child = parent
        parent = uf.parents[child]
        uf.parents[child] = root
    end
end

# `length(uf)` returns the number of nodes in `uf`.
function Base.length(uf :: UnionFinder)
    return length(uf.parents)
end

# `size!(uf, node)` returns the size of the group containing `node`.
function size!(uf :: UnionFinder{T}, node :: T) where T <: Integer
    return uf.sizes[find!(uf, node)]
end


text:
## Colwise
# TODO: replace with future mapslices specialization, see https://github.com/JuliaArrays/StaticArrays.jl/pull/99
"""
    $(SIGNATURES)

Equivalent to one of

```julia
mapslices(x -> f(a, x), B, dims=1)
mapslices(x -> f(x, b), A, dims=1)
```

but optimized for statically-sized matrices.
"""
function colwise end

colwise(f, a::AbstractVector, B::AbstractMatrix) = mapslices(x -> f(a, x), B, dims=1)
colwise(f, A::AbstractMatrix, b::AbstractVector) = mapslices(x -> f(x, b), A, dims=1)

@inline function colwise(f, a::StaticVector, B::StaticMatrix)
    Sa = Size(a)
    SB = Size(B)
    Sa[1] === SB[1] || throw(DimensionMismatch())
    _colwise(f, Val(SB[2]), a, B)
end

@inline function _colwise(f, ::Val{0}, a::StaticVector, B::StaticMatrix)
    zero(similar_type(B, promote_eltype(a, B)))
end

@inline function _colwise(f, M::Val, a::StaticVector, B::StaticMatrix)
    cols = ntuple(i -> f(a, B[:, i]), M)
    hcat(cols...)
end

@inline function colwise(f, A::StaticMatrix, b::StaticVector)
    SA = Size(A)
    Sb = Size(b)
    SA[1] === Sb[1] || throw(DimensionMismatch())
    _colwise(f, Val(SA[2]), A, b)
end

@inline function _colwise(f, ::Val{0}, A::StaticMatrix, b::StaticVector)
    zero(similar_type(A, promote_eltype(A, b)))
end

@inline function _colwise(f, M::Val, A::StaticMatrix, b::StaticVector)
    cols = ntuple(i -> f(A[:, i], b), M)
    hcat(cols...)
end

@inline function vector_to_skew_symmetric(v::SVector{3, T}) where {T}
    @SMatrix [zero(T) -v[3] v[2];
            v[3] zero(T) -v[1];
            -v[2] v[1] zero(T)]
end
const hat = vector_to_skew_symmetric

@inline function vector_to_skew_symmetric_squared(a::SVector{3})
    a1² = a[1] * a[1]
    a2² = a[2] * a[2]
    a3² = a[3] * a[3]
    b11 = -a2² - a3²
    b12 = a[1] * a[2]
    b13 = a[1] * a[3]
    b22 = -a1² - a3²
    b23 = a[2] * a[3]
    b33 = -a1² - a2²
    @SMatrix [b11 b12 b13;
            b12 b22 b23;
            b13 b23 b33]
end
const hat_squared = vector_to_skew_symmetric_squared

# The 'Bortz equation'.
# Bortz, John E. "A new mathematical formulation for strapdown inertial navigation."
# IEEE transactions on aerospace and electronic systems 1 (1971): 61-66.
#
# Or, interpreted in a Lie setting:
# d/dt(exp(ϕ(t))) = ̂(dexp(ϕ(t)) * ϕ̇(t)) * exp(ϕ(t)) (where dexp is the 'right trivialized' tangent of the exponential map)
# ̂(dexp(ϕ(t)) * ϕ̇(t)) = d/dt(exp(ϕ(t))) * exp(ϕ(t))⁻¹  (hat form of angular velocity in world frame)
#                      = ̂(exp(ϕ(t)) ω)      (with ω angular velocity in body frame)
# ϕ̇(t) = dexp⁻¹(ϕ(t)) * exp(ϕ(t)) * ω
function rotation_vector_rate(rotation_vector::AbstractVector{T}, angular_velocity_in_body::AbstractVector{T}) where {T}
    ϕ = rotation_vector
    ω = angular_velocity_in_body
    @boundscheck length(ϕ) == 3 || error("ϕ has wrong length")
    @boundscheck length(ω) == 3 || error("ω has wrong length")
    θ = norm(ϕ)
    ϕ̇ = ω
    if θ > eps(θ)
        s, c = sincos(θ)
        ϕ̇ += (ϕ × ω) / 2 + 1 / θ^2 * (1 - (θ * s) / (2 * (1 - c))) * ϕ × (ϕ × ω)
    end
    ϕ̇
end

@inline function transform_spatial_motion(angular::SVector{3}, linear::SVector{3}, rot::R, trans::SVector{3}) where {R <: Rotation{3}}
    angular = rot * angular
    linear = rot * linear + trans × angular
    angular, linear
end

@inline function mul_inertia(J, c, m, ω, v)
    angular = J * ω + c × v
    linear = m * v - c × ω
    angular, linear
end

# also known as 'spatial motion cross product'
@inline function se3_commutator(xω, xv, yω, yv)
    angular = xω × yω
    linear = xω × yv + xv × yω
    angular, linear
end

function quaternion_derivative end
function spquat_derivative end
function angular_velocity_in_body end

@inline function velocity_jacobian(::typeof(quaternion_derivative), q::Quat)
    (@SMatrix [
        -q.x -q.y -q.z;
        q.w -q.z  q.y;
        q.z  q.w -q.x;
        -q.y  q.x  q.w]) / 2
end

@inline function velocity_jacobian(::typeof(spquat_derivative), q::SPQuat)
    quat = Quat(q)
    dQuat_dW = velocity_jacobian(quaternion_derivative, quat)
    dSPQuat_dQuat = Rotations.jacobian(SPQuat, quat)
    dSPQuat_dQuat * dQuat_dW
end

@inline function velocity_jacobian(::typeof(angular_velocity_in_body), q::Quat)
    2 * @SMatrix [
    -q.x  q.w  q.z -q.y;
    -q.y -q.z  q.w  q.x;
    -q.z  q.y -q.x  q.w]
end

@inline function velocity_jacobian(::typeof(angular_velocity_in_body), q::SPQuat)
    quat = Quat(q)
    dW_dQuat = velocity_jacobian(angular_velocity_in_body, quat)
    dQuat_dSPQuat = Rotations.jacobian(Quat, q)
    dW_dQuat * dQuat_dSPQuat
end

@inline function quaternion_derivative(q::Quat, angular_velocity_in_body::AbstractVector)
    @boundscheck length(angular_velocity_in_body) == 3 || error("size mismatch")
    velocity_jacobian(quaternion_derivative, q) * angular_velocity_in_body
end

@inline function spquat_derivative(q::SPQuat, angular_velocity_in_body::AbstractVector)
    @boundscheck length(angular_velocity_in_body) == 3 || error("size mismatch")
    velocity_jacobian(spquat_derivative, q) * angular_velocity_in_body
end

@inline function angular_velocity_in_body(q::Quat, quat_derivative::AbstractVector)
    @boundscheck length(quat_derivative) == 4 || error("size mismatch")
    velocity_jacobian(angular_velocity_in_body, q) * quat_derivative
end

@inline function angular_velocity_in_body(q::SPQuat, spq_derivative::AbstractVector)
    @boundscheck length(spq_derivative) == 3 || error("size mismatch")
    velocity_jacobian(angular_velocity_in_body, q) * spq_derivative
end

function linearized_rodrigues_vec(r::RotMatrix) # TODO: consider moving to Rotations
    x = (r[3, 2] - r[2, 3]) / 2
    y = (r[1, 3] - r[3, 1]) / 2
    z = (r[2, 1] - r[1, 2]) / 2
    RodriguesVec(x, y, z)
end


text:
abstract type AbstractUpperTriangle{T} <: AbstractMatrix{T} end
# abstract type AbstractLowerTriangle{T} <: AbstractMatrix{T} end
abstract type AbstractSymmetric{T} <: AbstractMatrix{T} end

struct UpperTriangle2{T} <: AbstractUpperTriangle{T}
    data::SVector{3,T}
end
struct UpperTriangle3{T} <: AbstractUpperTriangle{T}
    data::SVector{6,T}
end
UpperTriangle2(A...) = UpperTriangle2(SVector{3}(A...))
UpperTriangle3(A...) = UpperTriangle3(SVector{6}(A...))
# struct LowerTriangle3{T} <: AbstractLowerTriangle{T}
#     data::SVector{6,T}
# end
# LowerTriangle3(A...) = LowerTriangle3(SVector{6}(A...))

struct SymmetricM2{T} <: AbstractSymmetric{T}
    data::SVector{3,T}
end
struct SymmetricM3{T} <: AbstractSymmetric{T}
    data::SVector{6,T}
end
SymmetricM2(A...) = SymmetricM2(SVector{3}(A...))
SymmetricM3(A...) = SymmetricM3(SVector{6}(A...))
@inline Base.getindex(A::Union{<: AbstractUpperTriangle, AbstractSymmetric}, i) = @inbounds A.data[i]
@inline function Base.getindex(A::AbstractUpperTriangle, i, j)
    @boundscheck begin
        i1, j1 = minmax(i,j)
        (i1 < 1 || j1 > 3) && throwboundserror()
    end
    i > j && return 0.0
    @inbounds A.data[i + (j-1)*j ÷ 2]
end
function LinearAlgebra.logdet(U::UpperTriangle3)
    log(U[1]) + log(U[3]) + log(U[6])
end
@inline function Base.getindex(A::AbstractSymmetric, i, j)
    i, j = minmax(i, j)
    @inbounds A.data[i + (j-1)*j ÷ 2]
end
Base.size(::Union{UpperTriangle2,SymmetricM2}) = (2,2)
Base.size(::Union{UpperTriangle3,SymmetricM3}) = (3,3)
Base.length(::Union{UpperTriangle2,SymmetricM2}) = 4
Base.length(::Union{UpperTriangle3,SymmetricM3}) = 9
Base.Array(S::SymmetricM2{T}) where T = T[S[1] S[2]; S[2] S[3]]
Base.Array(U::UpperTriangle2{T}) where T = T[U[1] U[2] ; 0 U[3]]
Base.Array(S::SymmetricM3{T}) where T = T[S[1] S[2] S[4]; S[2] S[3] S[5]; S[4] S[5] S[6]]
Base.Array(U::UpperTriangle3{T}) where T = T[U[1] U[2] U[4]; 0 U[3] U[5]; 0 0 U[6]]

@inline function Base.:+(A1::T, A2::T) where {T <: Union{AbstractUpperTriangle,AbstractSymmetric}}
    T( A1.data + A2.data )
end

@inline function Base.:*(A::UpperTriangle2, B::UpperTriangle2)
    @fastmath @inbounds UpperTriangle2(
        A[1]*B[1],
        A[1]*B[2] + A[2]*B[3],
        A[3]*B[3]
    )
end
@inline function Base.:*(A::UpperTriangle3, B::UpperTriangle3)
    @fastmath @inbounds UpperTriangle3(
        A[1]*B[1],
        A[1]*B[2] + A[2]*B[3],
        A[3]*B[3],
        A[1]*B[4] + A[2]*B[5] + A[4]*B[6],
        A[3]*B[5] + A[5]*B[6],
        A[6]*B[6]
    )
end
# @inline function Base.:*(A::LowerTriangle3, B::LowerTriangle3)
#     @fastmath @inbounds LowerTriangle3(
#         A[1]*B[1],
#         A[1]*B[2] + A[2]*B[3],
#         A[3]*B[3],
#         A[1]*B[4] + A[2]*B[5] + A[4]*B[6],
#         A[3]*B[5] + A[5]*B[6],
#         A[6]*B[6]
#     )
# end
@inline function Base.:*(a::UpperTriangle2, b::SVector{2})
    @fastmath @inbounds SVector(
        a[1]*b[1] + a[2]*b[2],
        a[3]*b[2]
    )
end
@inline function Base.:*(a::UpperTriangle3, b::SVector{3})
    @fastmath @inbounds SVector(
        a[1]*b[1] + a[2]*b[2] + a[4]*b[3],
        a[3]*b[2] + a[5]*b[3],
        a[6]*b[3]
    )
end
@inline function Base.:*(a::Adjoint{T,SVector{3,T}}, S::SymmetricM3{T}) where T
    @fastmath @inbounds SVector(
        S[1]*a[1] + S[2]*a[2] + S[4]*a[3],
        S[2]*a[1] + S[3]*a[2] + S[5]*a[3],
        S[4]*a[1] + S[5]*a[2] + S[6]*a[3]
    )'
end
@inline function Base.:*(S::SymmetricM3, a::SVector{3})
    @fastmath @inbounds SVector(
        S[1]*a[1] + S[2]*a[2] + S[4]*a[3],
        S[2]*a[1] + S[3]*a[2] + S[5]*a[3],
        S[4]*a[1] + S[5]*a[2] + S[6]*a[3]
    )
end
@inline function quadform(z::SVector{3}, Σ::SymmetricM3)
    @inbounds @fastmath begin
         a = z[1]^2*Σ[1] + 2z[1]*z[2]*Σ[2] + 2z[1]*z[3]*Σ[4]
         b = z[2]^2*Σ[3] + 2z[2]*z[3]*Σ[5] + z[3]^2*Σ[6]
         a + b
     end
end
@inline function quadform(U::UpperTriangle3, Σ::SymmetricM3) # U * Σ * U'
    @fastmath begin
        @inbounds U11, U12, U22, U13, U23, U33 = U[1], U[2], U[3], U[4], U[5], U[6]
        @inbounds Σ11, Σ12, Σ22, Σ13, Σ23, Σ33 = Σ[1], Σ[2], Σ[3], Σ[4], Σ[5], Σ[6]

        UΣ12 = (U11*Σ12 + U12*Σ22 + U13*Σ23)
        UΣ13 = (U11*Σ13 + U12*Σ23 + U13*Σ33)
        S11 = (U11*Σ11 + U12*Σ12 + U13*Σ13) * U11 + UΣ12 * U12 + UΣ13 * U13
        S12 = UΣ12 * U22 + UΣ13 * U23
        S13 = UΣ13 * U33
        UΣ23 = U22*Σ23 + U23*Σ33
        S22 = (U22*Σ22 + U23*Σ23) * U22 + UΣ23 * U23
        S23 = UΣ23*U33
        S33 = U33*Σ33*U33
        SymmetricM3(
            S11, S12, S22, S13, S23, S33
        )
    end
end
@inline function Base.:+(a::U, b::Number) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a.data + b)
end
@inline function Base.:+(a::Number, b::U) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a + b.data)
end
@inline function Base.:-(a::U, b::Number) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a.data - b)
end
@inline function Base.:-(a::Number, b::U) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a - b.data)
end
@inline function Base.:*(a::U, b::Number) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a.data * b)
end
@inline function Base.:*(a::Number, b::U) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a * b.data)
end
@inline function Base.:/(a::U, b::Number) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a.data / b)
end
@inline function Base.:/(a::Number, b::U) where {U <: Union{AbstractSymmetric,AbstractUpperTriangle}}
    U(a / b.data)
end
@inline function xxt(A::UpperTriangle2)
    @inbounds SymmetricM2(
        A[1]^2 + A[2]^2,
        A[2]*A[3],
        A[3]^2
    )
end
@inline function xxt(A::UpperTriangle3)
    @inbounds SymmetricM3(
        A[1]^2 + A[2]^2 + A[4]^2,
        A[2]*A[3] + A[4]*A[5],
        A[3]^2 + A[5]^2,
        A[4]*A[6],
        A[5]*A[6],
        A[6]^2
    )
end
@inline function xtx(A::UpperTriangle2)
    @inbounds SymmetricM2(
        A[1]^2,
        A[1]*A[2],
        A[2]^2 + A[3]^2
    )
end
@inline function xtx(A::UpperTriangle3)
    @inbounds SymmetricM3(
        A[1]^2,
        A[1]*A[2],
        A[2]^2 + A[3]^2,
        A[1]*A[4],
        A[2]*A[4] + A[3]*A[5],
        A[4]^2 + A[5]^2 + A[6]^2
    )
end

function Base.inv(U::UpperTriangle2)
    @fastmath @inbounds begin
        t11 = 1 / U[1]
        t22 = 1 / U[3]
        t12 = - U[2] * t11 * t22
    end
    UpperTriangle2( t11, t12, t22 )
end
function Base.inv(U::UpperTriangle3)
    @fastmath @inbounds begin
        t11 = 1 / U[1]
        t22 = 1 / U[3]
        t33 = 1 / U[6]
        t12 = - U[2] * t11 * t22
        t13 = - (U[4] * t11 + U[5] * t12) * t33
        t23 = - U[5] * t22 *t33
    end
    UpperTriangle3( t11, t12, t22, t13, t23, t33 )
end

function revchol(S::SymmetricM2{T}) where {T}
    @fastmath @inbounds begin
        ru22 = finite_sqrt(S[3], T(1e-4))
        ru12 = S[2] / ru22
        ru11 = finite_sqrt(S[1] - ru12^2, T(1e-4))
    end
    UpperTriangle2( ru11, ru12, ru22 )
end
function revchol(S::SymmetricM3{T}) where {T}
    @fastmath @inbounds begin
        ru33 = finite_sqrt(S[6], T(1e-4))
        ru13 = S[4] / ru33
        ru23 = S[5] / ru33
        ru22 = finite_sqrt(S[3] - ru23^2, T(1e-4))
        ru12 = (S[2] - ru13*ru23) / ru22
        ru11 = finite_sqrt(S[1] - ru12^2 - ru13^2, T(1e-4))
    end
    UpperTriangle3( ru11, ru12, ru22, ru13, ru23, ru33 )
end
function chol(S::SymmetricM3{T}) where {T}
    @fastmath @inbounds begin
        U11 = finite_sqrt(S[1], T(1e-4))
        U12 = S[2] / U11
        U13 = S[4] / U11
        U22 = finite_sqrt(S[3] - U12*U12, T(1e-4))
        U23 = (S[5] - U12*U13) / U22
        U33 = finite_sqrt(S[6] - U13*U13 - U23*U23, T(1e-4))
        UpperTriangle3(
            U11, U12, U22, U13, U23, U33
        )
    end
end
function chol(S::SymmetricM2{T}) where {T}
    @fastmath @inbounds begin
        U11 = finite_sqrt(S[1], T(1e-4))
        U12 = S[2] / U11
        U22 = finite_sqrt(S[3] - U12*U12, T(1e-4))
        UpperTriangle2(
            U11, U12, U22
        )
    end
end

@generated function Base.:*(A::SMatrix{M,N,T}, B::SMatrix{N,P,T}) where {M,N,P,T}
    outtup = Vector{Expr}(undef, M*P)
    i = 0
    for p ∈ 1:P, m ∈ 1:M
        i += 1
        outtup[i] = :(@inbounds $(Symbol(:C_, p))[$m].value )
    end
    quote
        $(Expr(:meta, :inline))
        A_col = $(Expr(:tuple, [:(@inbounds Core.VecElement(A[$m,1])) for m ∈ 1:M]...))
        Base.Cartesian.@nexprs $P p -> C_p = SIMDPirates.extract_data(SIMDPirates.vmul(A_col, @inbounds B[1,p]))
        @inbounds for n ∈ 2:$N
            A_col = $(Expr(:tuple, [:(@inbounds Core.VecElement(A[$m,n])) for m ∈ 1:M]...))
            Base.Cartesian.@nexprs $P p -> C_p = SIMDPirates.vfma(A_col, (@inbounds B[n,p]), C_p)
        end
        SMatrix{$M,$P,$T}(
            $(Expr(:tuple, outtup...))
        )
    end

end


text:
struct Xref

  n :: Int64
  x :: Matrix{Float64}

end

function compute_xcart!(x,xref,xcart)

  icm = 0
  iang = structures.n*3
  icart = 0
  for i in 1:structures.n
    icm = icm + 1
    iang = iang + 1
    for j in 1:xref[i].n
      icart = icart + 1
      xcart[icart,1:3] = compcart( x[icm:icm+2], xref[i].x[j,1:3],
                                   x[iang], x[iang+1], x[iang+2]  )
    end
  end

end

dist2(x,y) = (x[1]-y[1])^2 + (x[2]-y[2])^2 + (x[3]-y[3])^2

function f_naive(xcart,xref,tol2)

  n = length(xref)
  f = 0.
  for i in 1:xref.n-1
    for j in i+1:xref.n
      for ix in 1:xref[i].n
        for jx in 1:xref[j].n
          d2 = dist2(xref[i].x[ix,1:3], xref.x[jx,1:3])
          if d2 < tol2
            f = f + ( tol2 - d2 )^2
          end
        end
      end
    end
  end

  return f
end

max_stars_repo_name: mcubeg/JPleXL

max_issues_repo_name: mcubeg/JPleXL

max_forks_repo_name: mcubeg/JPleXL


text:
import Base: @kwdef

# structs defining terms of a composable model for the independent-particle
# susceptibility χ0. The struct define a call operator, which does some setup
# and returns an `apply!(δρ, δV, α=1)` function to form `δρ .+= α * χ0 * δV`.

abstract type χ0Model end

@doc raw"""
Represents the LDOS-based ``χ_0`` model
```math
χ_0(r, r') = (-D_\text{loc}(r) δ(r, r') + D_\text{loc}(r) D_\text{loc}(r') / D)
```
where ``D_\text{loc}`` is the local density of states and ``D`` the density of states.
For details see Herbst, Levitt 2020 arXiv:2009.01665
"""
@kwdef struct LdosModel <: χ0Model
    adjust_temperature = IncreaseMixingTemperature()
end
function (χ0::LdosModel)(basis; eigenvalues, ψ, εF, kwargs...)
    n_spin = basis.model.n_spin_components

    # Catch cases that will yield no contribution
    temperature = χ0.adjust_temperature(basis.model.temperature; kwargs...)
    iszero(temperature) && return nothing
    ldos = compute_ldos(εF, basis, eigenvalues, ψ; temperature)
    if maximum(abs, ldos) < eps(eltype(basis))
        return nothing
    end

    tdos = sum(sum, ldos) * basis.dvol  # Integrate LDOS to form total DOS
    function apply!(δρ, δV, α=1)
        δεF = dot(ldos, δV) .* basis.dvol
        δρ .+= α .* (-ldos .* δV .+ ldos .* δεF ./ tdos)
    end
end

@doc raw"""
A localised dielectric model for ``χ_0``:
```math
\sqrt{L(x)} \text{IFFT} \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)} \text{FFT} \sqrt{L(x)}
```
where ``C_0 = 1 - ε_r``, `L(r)` is a real-space localization function
and otherwise the same conventions are used as in [`DielectricMixing`](@ref).
"""
@kwdef struct DielectricModel <: χ0Model
    εr::Real  = 10
    kTF::Real = 0.8
    localization::Function = identity
end
function (χ0::DielectricModel)(basis; kwargs...)
    T   = eltype(basis)
    εr  = T(χ0.εr)
    kTF = T(χ0.kTF)
    C0  = 1 - εr
    iszero(C0) && return nothing  # Will yield no contribution

    Gsq = [sum(abs2, G) for G in G_vectors_cart(basis)]
    apply_sqrtL = identity
    if χ0.localization != identity
        sqrtL = sqrt.(χ0.localization.(r_vectors(basis)))
        apply_sqrtL = x -> sqrtL .* x
    end

    # TODO simplify apply_sqrtL
    function apply!(δρ, δV, α=1)
        loc_δV = r_to_G(basis, apply_sqrtL(δV))
        dielectric_loc_δV = @. C0 * kTF^2 * Gsq / 4T(π) / (kTF^2 - C0 * Gsq) * loc_δV
        δρ .+= α .* apply_sqrtL(G_to_r(basis, dielectric_loc_δV))
        δρ
    end
end

"""
Full χ0 application, optionally dropping terms or disabling Sternheimer.
All keyword arguments passed to [`apply_χ0`](@ref).
"""
struct Applyχ0Model <: χ0Model
    kwargs_apply_χ0
end
Applyχ0Model(; kwargs_apply_χ0...) = Applyχ0Model(kwargs_apply_χ0)
function (χ0::Applyχ0Model)(basis; ham, eigenvalues, ψ, εF, n_ep_extra, kwargs...)
    # self_consistent_field uses a few extra bands, which are not converged by the eigensolver
    # For the χ0 application, bands need to be separated between converged bands
    # and extra bands, which is what we do here
    ψ_cvg   = [@view ψk[:, 1:end-n_ep_extra]     for ψk in ψ]
    ψ_extra = [@view ψk[:, end-n_ep_extra+1:end] for ψk in ψ]
    ε_cvg   = [εk[1:end-n_ep_extra]              for εk in eigenvalues]
    ε_extra = [εk[end-n_ep_extra+1:end]          for εk in eigenvalues]

    function apply!(δρ, δV, α=1)
        χ0δV = apply_χ0(ham, ψ_cvg, εF, ε_cvg, δV; ψ_extra, ε_extra,
                        χ0.kwargs_apply_χ0...)
        δρ .+= α .* χ0δV
    end
end


text:
#abstract type TR1DModel{T} end

abstract type CubicModel{T} <: TR1DModel{T} end

abstract type CubicModel2Pt{T} <: CubicModel{T} end

"Cubic model solution"
function solveModel(m:: CubicModel{T}, t :: T) where T
    disc = m.H^2 - T(2)*m.C*m.g
    d = disc < 0 ? -sign(m.C)*Inf : (-m.H + sqrt(disc)) / m.C
    return d
end

"Cubic model evaluation"
function m(M:: CubicModel{T}, t :: T) where T
    # Horner rule
    return M.f + t*(M.g + t*( T(0.5)*M.H + t*( T(1.0/6.0)*M.C)))
end

"Cubic model derivative evaluation"
function dm(M:: CubicModel{T}, t :: T) where T
    return M.g  +  t*(M.H  +  t*T(0.5)*M.C)
end

"Cubic model predicted reduction"
function Δm(M:: CubicModel{T}, t :: T) where T
    return M.g * t  +  T(0.5) * M.H * t^2  + T(1.0/6.0)*M.C*t^3
end

"Cubic model predicted reduction given two points"
function Δm(M:: CubicModel{T}, t1 :: T, t2 :: T) where T
    c1 = M.g * t1  +  T(0.5) * M.H * t1^2  + T(1.0/6.0)*M.C*t1^3
    c2 = M.g * t2  +  T(0.5) * M.H * t2^2  + T(1.0/6.0)*M.C*t2^3
    return c1 - c2
end

"Taylor cubic polynomial model"
mutable struct Taylor3{T} <: CubicModel{T}
    f :: T
    g :: T
    H :: T
    C :: T
end
Taylor3{T}() where T = Taylor3{T}(0,0,0,0)

"2 points interpolation cubic polynomial model"
mutable struct TwoPoints3{T} <: CubicModel2Pt{T}
    f :: T
    g :: T
    H :: T
    C :: T
end
TwoPoints3{T}() where T = TwoPoints3{T}(0,0,0,0)

function initModel(M:: Taylor3{T}, ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H = hess(ϕ, t)
    M.C = tens(ϕ,t)
    return M
end

function updateModel!(M:: Taylor3{T}, ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g  = dϕt
    M.H = hess(ϕ, t)
    M.C = tens(ϕ,t)
    return M
end

"Secant cubic model"
mutable struct secant3{T} <: CubicModel2Pt{T}
    f :: T
    g :: T
    H :: T
    C :: T
    Q :: T  # to test plain secant
end
secant3{T}() where T = secant3{T}(0,0,0,0,0)

function initModel(M:: CubicModel2Pt{T}, ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H, M.C = initHC_s!(ϕ, t, ϕt, dϕt)
    return M
end

function updateModel!(M:: CubicModel2Pt{T}, ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H, M.C = updateHC_s!(ϕ, t, ϕt, dϕt)
    return M
end

function initModel(M:: secant3{T}, ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H, M.C, M.Q = initHC_s!(ϕ, t, ϕt, dϕt)
    return M
end

function updateModel!(M:: secant3{T}, ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H, M.C, M.Q = updateHC_s!(ϕ, t, ϕt, dϕt)
    return M
end

let (tpred, ϕpred, dϕpred) = (NaN, NaN, NaN)

    global function initHC_s!(ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
        (tpred, ϕpred, dϕpred) = (t, ϕt, dϕt)
        return T(1), T(0), T(1)
    end

    global function updateHC_s!(ϕ :: LineModel, t :: T, ϕt :: T, dϕt :: T) where T
        # estimate quadratic and cubic coefficients from the hermite interpolant
        a = coefHIN3([t;tpred], [ϕt dϕt; ϕpred dϕpred])
        s = t - tpred
        y = dϕt - dϕpred
        Q = y/s

        tpred, ϕpred, dϕpred = t, ϕt, dϕt

        return T(2.0)*a[3], T(6.0)*a[4], Q
    end

    "Two consecutive secant steps on the cubic model"
    global function solveModel(M:: secant3{T}, t :: T) where T
        # the "Taylor" cubic polynomial fits with the function and its derivative at t and tpred
        # Thus, secant iteration may be performed on the cubic model without knowledge of ϕ.
        # first secant iteration
        #s = t - tpred
        dmt = M.g  #  dm(M, t)
        #dmpred = dm(M, tpred)
        #y = dmt - dmpred
        dS = - (dmt / M.Q)

        # second secant iteration
        t2 = t + dS
        dm2 = dm(M, t2)
        y2 = dm2 - dmt
        d2 = - (dm2 * dS / y2)
        dS2 = dS + d2

        #return dS # for testing, should be equivalent to secant
        #return dS2 # for testing, should be equivalent to secant

        # try the cubic's root
        disc = M.H^2 - 2.0*M.C*M.g
        dC = disc < 0 ? -sign(M.C)*Inf : (-M.H + sqrt(disc)) / M.C
        return dC
    end
end


text:
"""
    Acrobot
"""

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

    m1   # mass link 1
    J1    # inertia link 1
    l1    # length link 1
    lc1   # length to COM link 1

    m2    # mass link 2
    J2    # inertia link 2
    l2    # length link 2
    lc2   # length to COM link 2

    g     # gravity

    b1    # joint friction
    b2
end

function M(model::Acrobot, x, w)
    a = (model.J1 + model.J2 + model.m2 * model.l1 * model.l1
         + 2.0 * model.m2 * model.l1 * model.lc2 * cos(x[2]))

    b = model.J2 + model.m2 * model.l1 * model.lc2 * cos(x[2])

    c = model.J2

    @SMatrix [a b;
              b c]
end

function τ(model::Acrobot, x, w)
    a = (-1.0 * model.m1 * model.g * model.lc1 * sin(x[1])
         - model.m2 * model.g * (model.l1 * sin(x[1])
         + model.lc2 * sin(x[1] + x[2])))

    b = -1.0 * model.m2 * model.g * model.lc2 * sin(x[1] + x[2])

    @SVector [a,
              b]
end

function C(model::Acrobot, x, w)
    a = -2.0 * model.m2 * model.l1 * model.lc2 * sin(x[2]) * x[4]
    b = -1.0 * model.m2 * model.l1 * model.lc2 * sin(x[2]) * x[4]
    c = model.m2 * model.l1 * model.lc2 * sin(x[2]) * x[3]
    d = 0.0

    @SMatrix [a b;
              c d]
end

function B(model::Acrobot, x, w)
    @SMatrix [0.0;
              1.0]
end

function f(model::Acrobot, x, u, w)
    q = view(x, 1:2)
    v = view(x, 3:4)
    qdd = M(model, q, w) \ (-1.0 * C(model, x, w) * v
            + τ(model, q, w) + B(model, q, w) * u[1:1] - [model.b1; model.b2] .* v)
    @SVector [x[3],
              x[4],
              qdd[1],
              qdd[2]]
end

function k_mid(model::Acrobot, x)
    @SVector [model.l1 * sin(x[1]),
              -1.0 * model.l1 * cos(x[1])]
end

function k_ee(model::Acrobot, x)
    @SVector [model.l1 * sin(x[1]) + model.l2 * sin(x[1] + x[2]),
              -1.0 * model.l1 * cos(x[1]) - model.l2 * cos(x[1] + x[2])]
end

n, m, d = 4, 1, 4
model = Acrobot{Midpoint, FixedTime}(4, 1, 4, 1.0, 0.33, 1.0, 0.5, 1.0, 0.33, 1.0, 0.5, 9.81, 0.1, 0.1)

# visualization
function visualize!(vis, model::Acrobot, x;
        color=RGBA(0.0, 0.0, 0.0, 1.0),
        r = 0.1, Δt = 0.1)

    default_background!(vis)

    i = 1
    l1 = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l1),
        convert(Float32, 0.025))
    setobject!(vis["l1$i"], l1, MeshPhongMaterial(color = color))
    l2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l2),
        convert(Float32, 0.025))
    setobject!(vis["l2$i"], l2, MeshPhongMaterial(color = color))

    setobject!(vis["elbow$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 1.0)))
    setobject!(vis["ee$i"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 1.0)))

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

    T = length(x)
    for t = 1:T

        MeshCat.atframe(anim,t) do
            p_mid = [k_mid(model, x[t])[1], 0.0, k_mid(model, x[t])[2]]
            p_ee = [k_ee(model, x[t])[1], 0.0, k_ee(model, x[t])[2]]

            settransform!(vis["l1$i"], cable_transform(zeros(3), p_mid))
            settransform!(vis["l2$i"], cable_transform(p_mid, p_ee))

            settransform!(vis["elbow$i"], Translation(p_mid))
            settransform!(vis["ee$i"], Translation(p_ee))
        end
    end

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

    MeshCat.setanimation!(vis, anim)
end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# Maximum adjacency visit / traversal

#################################################
#
#  Maximum adjacency visit
#
#################################################

struct MaximumAdjacency <: AbstractGraphVisitAlgorithm
end

abstract type AbstractMASVisitor <: AbstractGraphVisitor end

function maximum_adjacency_visit_impl!{T}(
    g::AbstractGraph,                               # the graph
    pq::DataStructures.PriorityQueue{Int, T},         # priority queue
    visitor::AbstractMASVisitor,                      # the visitor
    colormap::Vector{Int})                            # traversal status

    while !isempty(pq)
        u = DataStructures.dequeue!(pq)
        discover_vertex!(visitor, u)
        for v in out_neighbors(g, u)
            examine_neighbor!(visitor, u, v, 0, 0, 0)

            if haskey(pq,v)
                ed = visitor.distmx[u, v]
                pq[v] += ed
            end
        end
        close_vertex!(visitor, u)
    end

end

function traverse_graph!(
    g::AbstractGraph,
    T::DataType,
    alg::MaximumAdjacency,
    s::Integer,
    visitor::AbstractMASVisitor,
    colormap::Vector{Int})

    pq = DataStructures.PriorityQueue(Int,T,Base.Order.Reverse)

    # Set number of visited neighbors for all vertices to 0
    for v in vertices(g)
        pq[v] = zero(T)
    end

    @assert haskey(pq,s)
    @assert nv(g) >= 2

    #Give the starting vertex high priority
    pq[s] = one(T)

    #start traversing the graph
    maximum_adjacency_visit_impl!(g, pq, visitor, colormap)
end

#################################################
#
#  Visitors
#
#################################################

#################################################
#
#  Minimum Cut Visitor
#
#################################################

mutable struct MinCutVisitor{T, U<:Integer} <: AbstractMASVisitor
    graph::AbstractGraph
    parities::BitVector
    colormap::Vector{Int}
    bestweight::T
    cutweight::T
    visited::Int
    distmx::AbstractMatrix{T}
    vertices::Vector{U}
end

function MinCutVisitor{T}(g::AbstractGraph, distmx::AbstractMatrix{T})
    U = eltype(g)
    n = nv(g)
    parities = falses(n)
    return MinCutVisitor(
        g,
        falses(n),
        zeros(Int,n),
        typemax(T),
        zero(T),
        zero(Int),
        distmx,
        Vector{U}()
    )
end

function discover_vertex!(vis::MinCutVisitor, v::Integer)
    vis.parities[v] = false
    vis.colormap[v] = 1
    push!(vis.vertices,v)
    return true
end

function examine_neighbor!(vis::MinCutVisitor, u::Integer, v::Integer, ucolor::Int, vcolor::Int, ecolor::Int)
    ew = vis.distmx[u, v]

    # if the target of e is already marked then decrease cutweight
    # otherwise, increase it

    if vis.colormap[v] != vcolor # here vcolor is 0
        vis.cutweight -= ew
    else
        vis.cutweight += ew
    end
    return true
end

function close_vertex!(vis::MinCutVisitor, v::Integer)
    vis.colormap[v] = 2
    vis.visited += 1

    if vis.cutweight < vis.bestweight && vis.visited < nv(vis.graph)
        vis.bestweight = vis.cutweight
        for u in vertices(vis.graph)
            vis.parities[u] = ( vis.colormap[u] == 2)
        end
    end
    return true
end

#################################################
#
#  MAS Visitor
#
#################################################

struct MASVisitor{T, U<:Integer} <: AbstractMASVisitor
    io::IO
    vertices::Vector{U}
    distmx::AbstractMatrix{T}
    log::Bool
end

function discover_vertex!{T}(visitor::MASVisitor{T}, v::Integer)
    push!(visitor.vertices,v)
    visitor.log ? println(visitor.io, "discover vertex: $v") : nothing
    return true
end

function examine_neighbor!(visitor::MASVisitor, u::Integer, v::Integer, ucolor::Int, vcolor::Int, ecolor::Int)
    visitor.log ? println(visitor.io, " -- examine neighbor from $u to $v") : nothing
    return true
end

function close_vertex!(visitor::MASVisitor, v::Integer)
    visitor.log ? println(visitor.io, "close vertex: $v") : nothing
    return true
end

#################################################
#
#  Minimum Cut
#
#################################################

"""
    mincut(g, distmx=DefaultDistance())

Return a tuple `(parity, bestcut)`, where `parity` is a vector of integer
values that determines the partition in `g` (1 or 2) and `bestcut` is the
weight of the cut that makes this partition. An optional `distmx` matrix may
be specified; if omitted, edge distances are assumed to be 1.
"""
function mincut{T}(
    g::AbstractGraph,
    distmx::AbstractMatrix{T}=DefaultDistance()
)
    visitor = MinCutVisitor(g, distmx)
    colormap = zeros(Int, nv(g))
    traverse_graph!(g, T, MaximumAdjacency(), 1, visitor, colormap)
    return(visitor.parities + 1, visitor.bestweight)
end

"""
    maximum_adjacency_visit(g[, distmx][, log][, io])

Return the vertices in `g` traversed by maximum adjacency search. An optional
`distmx` matrix may be specified; if omitted, edge distances are assumed to
be 1. If `log` (default `false`) is `true`, visitor events will be printed to
`io`, which defaults to `STDOUT`; otherwise, no event information will be
displayed.
"""
function maximum_adjacency_visit{T}(
    g::AbstractGraph,
    distmx::AbstractMatrix{T},
    log::Bool,
    io::IO
)
    visitor = MASVisitor(io, Vector{Int}(), distmx, log)
    traverse_graph!(g, T, MaximumAdjacency(), 1, visitor, zeros(Int, nv(g)))
    return visitor.vertices
end

maximum_adjacency_visit(g::AbstractGraph) = maximum_adjacency_visit(
    g,
    DefaultDistance(),
    false,
    STDOUT
)

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
"""
     add_singleton(x::AbstractArray, ::Val{dim}) where {dim}

Add an additional dimension `dim` of size 1 to array `x`.
"""
function add_singleton(x::AbstractArray, ::Val{dim}) where {dim}
    shape = ntuple(max(ndims(x) + 1, dim)) do i
        return i < dim ? size(x, i) : (i > dim ? size(x, i - 1) : 1)
    end
    return reshape(x, shape)
end

"""
    dot_matwise(x::AbstractArray, y::AbstractArray)

Compute the inner product of all matrices in `x` and `y`.

At least one of `x` and `y` has to be a matrix.
"""
dot_matwise(x::AbstractMatrix, y::AbstractMatrix) = dot(x, y)
function dot_matwise(x::AbstractArray, y::AbstractMatrix)
    xmat = reshape(x, size(x, 1) * size(x, 2), :)
    return reshape(reshape(y, 1, :) * xmat, size(x)[3:end])
end
dot_matwise(x::AbstractMatrix, y::AbstractArray) = dot_matwise(y, x)

function dot_vecwise(x::AbstractMatrix, y::AbstractMatrix)
    return [
        dot(u, v) for (u, v) in
        zip((view(x, :, i) for i in axes(x, 2)), (view(y, :, i) for i in axes(y, 2)))
    ]
end

dot_vecwise(x::AbstractMatrix, y::AbstractVector) = x' * y

dot_vecwise(x::AbstractVector, y::AbstractMatrix) = dot_vecwise(y, x)
dot_vecwise(x::AbstractVector, y::AbstractVector) = dot(x, y)

"""
    checksize(μ::AbstractVecOrMat, ν::AbstractVecOrMat, C::AbstractMatrix)

Check that dimensions of source and target marginals `μ` and `ν` are consistent with cost
matrix `C`.
"""
function checksize(μ::AbstractVecOrMat, ν::AbstractVecOrMat, C::AbstractMatrix)
    size(C) == (size(μ, 1), size(ν, 1)) || throw(
        DimensionMismatch("cost matrix `C` must be of size `(size(μ, 1), size(ν, 1))`")
    )
    return nothing
end

"""
    checksize2(x::AbstractVecOrMat, y::AbstractVecOrMat)

Check if arrays `x` and `y` are compatible, then return a tuple of its broadcasted second
dimension.
"""
checksize2(::AbstractVector, ::AbstractVector) = ()
function checksize2(μ::AbstractVecOrMat, ν::AbstractVecOrMat)
    size_μ_2 = size(μ, 2)
    size_ν_2 = size(ν, 2)
    if size_μ_2 > 1 && size_ν_2 > 1 && size_μ_2 != size_ν_2
        throw(DimensionMismatch("size of source and target marginals is not compatible"))
    end
    return (max(size_μ_2, size_ν_2),)
end

"""
     isallapprox(x::AbstractVector)

Check that all entries of `x` are approximately equal
"""
function isallapprox(x::AbstractVecOrMat)
    return all(y -> isapprox(y, x[1]), x[2:end])
end

"""
     checkbalanced(μ::AbstractVecOrMat, ν::AbstractVecOrMat)

Check that source and target marginals `μ` and `ν` are balanced.
"""
function checkbalanced(μ::AbstractVector, ν::AbstractVector)
    sum(μ) ≈ sum(ν) || throw(ArgumentError("source and target marginals are not balanced"))
    return nothing
end
function checkbalanced(x::AbstractVecOrMat, y::AbstractVecOrMat)
    all(isapprox.(sum(x; dims=1), sum(y; dims=1))) ||
        throw(ArgumentError("source and target marginals are not balanced"))
    return nothing
end
function checkbalanced(x::AbstractMatrix)
    isallapprox(sum(x; dims=1)) ||
        throw(ArgumentError("source and target marginals are not balanced"))
    return nothing
end

"""
    A_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector)

Compute the matrix-vector product `Ab` and write the result to `c`.
"""
A_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector) = mul!(c, A, b)

"""
    A_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

Compute the matrix-matrix product `AB` and write the result to `C`.
"""
A_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix) = mul!(C, A, B)

"""
    A_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)

Compute the matrix-vector products `A[:, :, i] B[:, i]` and write them to `C[:, i]`.
"""
function A_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)
    return NNlib.batched_mul!(add_singleton(C, Val(2)), A, add_singleton(B, Val(2)))
end

"""
    At_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector)

Compute the matrix-vector product `transpose(A)b` and write the result to `c`.
"""
function At_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector)
    return mul!(c, transpose(A), b)
end

"""
    At_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

Compute the matrix-matrix product `transpose(A)B` and write the result to `C`.
"""
function At_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
    return mul!(C, transpose(A), B)
end

"""
    At_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)

Compute the matrix-vector products `transpose(A[:, :, i]) B[:, i]` and write them to
`C[:, i]`.
"""
function At_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)
    return NNlib.batched_mul!(add_singleton(C, Val(1)), add_singleton(B, Val(1)), A)
end


text:
#=
1 - download the zip file
2 - extract the zip file
3 - move inc/ and lib/ to compile/
4 - move libODriveLite.1.4.3.a to compile/
5 - compile
=#

@linux ? nothing : error("OpenDriveManagerLite only supports Linux at this time")

if isfile("odrMgrLite.1.4.3.zip")
    rm("odrMgrLite.1.4.3.zip")
end
if isdir("compile/lib")
    rm("compile/lib", recursive=true)
end
if isdir("compile/inc")
    rm("compile/inc", recursive=true)
end
if isfile("compile/libODriveLite.1.4.3.a")
    rm("compile/libODriveLite.1.4.3.a")
end

run(`wget --user-agent="Chrome" http://www.vires.com/download/OpenDRIVE/odrMgrLite.1.4.3.zip`)
run(`unzip odrMgrLite.1.4.3.zip`)
if isfile("odrMgrLite.1.4.3.zip")
    rm("odrMgrLite.1.4.3.zip")
end

run(`mv V1.4.3lite/inc compile/inc`)
run(`mv V1.4.3lite/lib compile/lib`)
run(`cp compile/lib/Linux/lib64/libODriveLite.1.4.3.a compile/`)

if isdir("V1.4.3lite")
    rm("V1.4.3lite", recursive=true)
end

cd("compile") do
    run(`bash -x compile.sh`)
end

max_stars_repo_name: sisl/OpenDriveManagerLite

max_issues_repo_name: sisl/OpenDriveManagerLite

max_forks_repo_name: sisl/OpenDriveManagerLite


text:
"""
    big_m_reformulation!(constr::ConstraintRef, bin_var, M, i, j, k)

Perform Big-M reformulation on a linear or quadratic constraint at index k of constraint j in disjunct i.

    big_m_reformulation!(constr::NonlinearConstraintRef, bin_var, M, i, j, k)

Perform Big-M reformulaiton on a nonlinear constraint at index k of constraint j in disjunct i.

    big_m_reformulation!(constr::AbstractArray{<:ConstraintRef}, bin_var, M, i, j, k)

Perform Big-M reformulation on a constraint at index k of constraint j in disjunct i.
"""
function big_m_reformulation!(constr::ConstraintRef, bin_var, M, i, j, k)
    M = get_reform_param(M, i, j, k; constr)
    bin_var_ref = constr.model[bin_var][i]
    add_to_function_constant(constr, -M)
    set_normalized_coefficient(constr, bin_var_ref , M)
end
function big_m_reformulation!(constr::NonlinearConstraintRef, bin_var, M, i, j, k)
    M = get_reform_param(M, i, j, k; constr)
    #create symbolic variables (using Symbolics.jl)
    for var_ref in constr.model[:gdp_variable_refs]
        var_sym = Symbol(var_ref)
        eval(:(Symbolics.@variables($var_sym)[1]))
    end
    bin_var_sym = Symbol("$bin_var[$i]")
    λ = Num(Symbolics.Sym{Float64}(bin_var_sym))

    #parse constr
    op, lhs, rhs = parse_constraint(constr)
    replace_Symvars!(lhs, constr.model) #convert JuMP variables into Symbolic variables
    gx = eval(lhs) #convert the LHS of the constraint into a Symbolic expression
    gx = gx - M*(1-λ) #add bigM

    #update constraint
    replace_constraint(constr, gx, op, rhs)
end
big_m_reformulation!(constr::AbstractArray{<:ConstraintRef}, bin_var, M, i, j, k) =
    big_m_reformulation(constr[k], bin_var, M, i, j, k)


text:
using PrototypesCriticisms
using Clustering
using Distances
using KernelFunctions
using LinearAlgebra
using Statistics
using Test

@testset verbose=true "PrototypesCriticisms.jl" begin
    @testset "sqmmd (mmd²)" begin
        X = rand(5, 10)
        k = RBFKernel()

        @test sqmmd(X, X, k) ≈ 0

        Y = X
        XX = kernelmatrix(k, X, obsdim=2)
        XY = kernelmatrix(k, X, Y, obsdim=2)
        YY = kernelmatrix(k, Y, obsdim=2)
        @test sqmmd(XX, XY, YY) ≈ 0

        @test mmd²(X, X, k) == sqmmd(X, X, k)
        @test mmd²(XX, XY, YY) == sqmmd(XX, XY, YY)
    end

    @testset "witness" begin
        X = zeros(5, 10)
        Y = ones(5, 10)
        x = X[:, 1]
        y = Y[:, 1]
        k = RBFKernel()
        @test witness(x, X, X, k) ≈ 0
        @test witness(x, X, Y, k) > 0
        @test witness(y, X, Y, k) < 0
    end

    @testset "prototypes" begin
        n = 50
        X = rand(5, n)
        k = RBFKernel()

        @test length(prototypes(X, k, 0)) == 0
        @test Set(prototypes(X, k, n)) == Set(1:n)

        K = kernelmatrix(k, X, obsdim=2)
        @test length(prototypes(K, 0)) == 0
        @test Set(prototypes(K, n)) == Set(1:n)

        k = 2
        c = kmedoids(pairwise(Euclidean(), X, dims=2), k)
        protoids = prototypes(c)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == k
        for i = 1:k
            @test length(protoids[i]) == 1
        end
        protoids = prototypes(c, 2)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == k
        for i = 1:k
            @test length(protoids[i]) == 2
        end

        k = 2
        c = kmeans(X, k)
        protoids = prototypes(c)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == k
        for i = 1:k
            @test length(protoids[i]) == 1
        end
        protoids = prototypes(c, 2)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == k
        for i = 1:k
            @test length(protoids[i]) == 2
        end

        k = 2
        c = fuzzy_cmeans(X, k, 2)
        protoids = prototypes(c)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == k
        for i = 1:k
            @test length(protoids[i]) == 1
        end
        protoids = prototypes(c, 2)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == k
        for i = 1:k
            @test length(protoids[i]) == 2
        end

        S = -pairwise(Euclidean(), X, dims=2)
        S = S - diagm(0 => diag(S)) + median(S) * I
        c = affinityprop(S)
        protoids = prototypes(c, X)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == nclusters(c)
        for i = 1:nclusters(c)
            @test length(protoids[i]) == 1
        end
        protoids = prototypes(c, X, 2)
        @test typeof(protoids) == Vector{Vector{Int}}
        @test length(protoids) == nclusters(c)
        for i = 1:nclusters(c)
            @test length(protoids[i]) == 2
        end
    end

    @testset "criticisms" begin
        n = 50
        X = rand(5, n)
        k = RBFKernel()

        @test length(criticisms(X, k, [1, 2], 0)) == 0
        @test Set(criticisms(X, k, [1, 2], n - 2)) == Set(3:n)

        K = kernelmatrix(k, X, obsdim=2)
        @test length(criticisms(K, [1, 2], 0)) == 0
        @test Set(criticisms(K, [1, 2], n - 2)) == Set(3:n)

        k = 2
        c = kmedoids(pairwise(Euclidean(), X, dims=2), k)
        critids = criticisms(c)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == k
        for i = 1:k
            @test length(critids[i]) == 1
        end
        critids = criticisms(c, 2)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == k
        for i = 1:k
            @test length(critids[i]) == 2
        end

        k = 2
        c = kmeans(X, k)
        critids = criticisms(c)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == k
        for i = 1:k
            @test length(critids[i]) == 1
        end
        critids = criticisms(c, 2)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == k
        for i = 1:k
            @test length(critids[i]) == 2
        end

        k = 2
        c = fuzzy_cmeans(X, k, 2)
        critids = criticisms(c)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == k
        for i = 1:k
            @test length(critids[i]) == 1
        end
        critids = criticisms(c, 2)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == k
        for i = 1:k
            @test length(critids[i]) == 2
        end

        S = -pairwise(Euclidean(), X, dims=2)
        S = S - diagm(0 => diag(S)) + median(S) * I
        c = affinityprop(S)
        critids = criticisms(c, X)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == nclusters(c)
        for i = 1:nclusters(c)
            @test length(critids[i]) == 1
        end
        critids = prototypes(c, X, 2)
        @test typeof(critids) == Vector{Vector{Int}}
        @test length(critids) == nclusters(c)
        for i = 1:nclusters(c)
            @test length(critids[i]) == 2
        end
    end
end


text:
macro UInt()
  Symbol("UInt", __BITS)
end

macro Int()
  Symbol("Int", __BITS)
end

macro u(value)
  esc(:(reinterpret(@UInt,$value)))
end

macro s(value)
  esc(:(reinterpret(@Int,$value)))
end

macro signbit()
  one(@UInt) << (__BITS - 1)
end

macro invertbit()
  one(@UInt) << (__BITS - 2)
end

macro mask(N)
  :(((one(@UInt) << N) - one(@UInt)) << ((sizeof(@UInt) * 8) - N))
end

doc"""
  `SigmoidNumbers.@breakdown(S::Symbol, [:arithmetic | :numeric])`

  takes a small sigmoid number referenced by a symbol and creates in the containing
  scope the following symbols:

  s_S::@Int         -signed int value
  u_S::@UInt        -unsigned, int of the absolute value.
  S_sgn::Bool       -true if negative, false if positive.
  S_reg::Int        -regime value for the sigmoid number.
  S_exp::Int        -signed integer representation of the exponent.
  S_inv::Bool       -true if S_exp is negative, false if positive
  S_frc::@UInt      -left-shifted fraction value, in a UInt whose size matches
                     library resolution.
  S_ubt::Bool       -true if S has a ubit.  Always false for non-ubit-type

  in the "arithmetic" mode, the exponent is set to -1 in the case of [0..1], and the
  fraction is flushed left to have the "leading one" present, aka the "hidden bit"
  in the representation by default.  Default mode is numeric.

  NB:  arithmetic mode is only useful when (ES == 0), and this macro will insert
  an error if that is not the case.

  This does not handle the zero or infinity cases, which should always be dealt with separately.
"""
macro breakdown(args...)

  goodsymbols = [:arithmetic, :numeric]

  #argument list handling.
  if length(args) == 0
    throw(ArgumentError("@breakdown macro requires an argument"))
  elseif length(args) == 1
    numeric = true
  elseif length(args) == 2
    (args[2] in goodsymbols) || throw(ArgumentError("option $(args[2]) unknown to @breakdown macro"))
    numeric = (args[2] == :numeric)
  else
    throw(ArgumentError("too many arguments to @breakdown macro"))
  end
  x = args[1] #reassign the x value to the first symbol.

  #throw an error if we don't have a simple symbol for the @breakdown macro.
  (typeof(x) == Symbol) || throw(ArgumentError("@breakdown macro can only operate on symbols."))

  #a small chart of the things we'll generate.
  s_x   = Symbol("s_", x)
  u_x   = Symbol("u_", x)
  x_sgn = Symbol(x, "_sgn")
  x_inv = Symbol(x, "_inv")
  x_reg = Symbol(x, "_reg")
  x_exp = Symbol(x, "_exp")
  x_frc = Symbol(x, "_frc")
  x_ubt = Symbol(x, "_ubt")

  precode = quote
    #generate the signed integer representation for easy sign checking.
    $s_x = @s($x)
    #evaluate the sign of the value.
    $x_sgn::Bool = ($s_x < 0)

    #convert to an unsigned integer representing the absolute value of the integer
    #representation.
    $u_x = @u($x_sgn ? -$s_x : $s_x)

    #probe the second position to decide if it's less than one or not.
    $x_inv = ($u_x & @invertbit) == 0
  end

  #codegen for the exponent and fraction code.
  if numeric
    esc(quote
      #import the common precode
      $precode

      #in the normal case, calculate the exponent by counting the UNARY portion of
      #the number.  This is leading zeros for less than one, and leading ones for
      #greater than one.  Since greater than one starts at exponent zero, subtract
      #one in that case.

      $x_reg::Int = (($x_inv) ? leading_zeros($u_x) : leading_ones($u_x | @signbit)) - 1 - !$x_inv

      if (mode == :ubit)
        ($x_reg = min($x_reg, N - 3))
        const __UB_mask = ((@signbit) >> (N - 1))
      else
        const __UB_mask = zero(@UInt)
      end

      $x_exp::Int = $x_reg * ($x_inv ? -1 : 1) * (1 << ES)

      if (ES != 0)
        const __ES_mask = (1 << ES) - 1
        #it's OK if this is negative, yo.
        __ES_shift = __BITS - ES - 3 - $x_reg + $x_inv

        $x_exp += @s(((($u_x & ~__UB_mask) >> __ES_shift) & __ES_mask))
      end
      #mask out the ubit, then flush the fraction all the way to the left.  This
      #also means we don't have to do any masking operations for higher order bits.
      $x_ubt = ($u_x & __UB_mask != 0)
      $x_frc = ($u_x & ~__UB_mask) << ($x_reg + ($x_inv ? 2 : 3) + ES)
    end)
  else
    esc(quote
      #insert a check precluding the use of arithmetic mode when ES is not zero.
      (ES != 0) && throw(AssertionError("arithmetic mode used when ES == $ES is not zero."))

      $precode
      #the arithmetic case for the exponent is similar, except we truncate at -1.
      $x_exp::Int = ($x_inv ? 0 : (leading_ones($u_x | @signbit) - 1) - 1)

      #the shift scales linerarly with the "truncated exponent".  We need to graft
      #the "hidden one" in the case that the number is >= 1, aka not "subnormal."
      $x_frc = ($u_x << ($x_exp + 2 - $x_inv)) | (!$x_inv * @signbit)
    end)
  end
end

#for test purposes

function __numeric_breakdown_wrapper{N, ES, mode}(x::Sigmoid{N, ES, mode})
  @breakdown x numeric
  (x_sgn, x_inv, x_exp, x_frc, x_ubt)
end

function __arithmetic_breakdown_wrapper{N, ES, mode}(x::Sigmoid{N, ES, mode})
  @breakdown x arithmetic
  (x_sgn, x_inv, x_exp, x_frc, x_ubt)
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
function init_values(x0::Array{Float64,2}, size_total_x::Int, size_total_y::Int,
                    temp1_init::Float64, temp2_init::Float64)

   #Setup temp1_init on borders
   for i = 1:size_total_x - 1
     x0[i,1] = temp1_init
     x0[i,size_total_y] = temp1_init
   end

   for j = 1:size_total_y
     x0[1,j] = temp1_init
     x0[size_total_x, j] = temp1_init
   end

   for i = 2:size_total_x-1
     x0[i,2] = temp1_init
     x0[i,size_total_y-1] = temp1_init
   end

   for j = 2:size_total_y-1
     x0[2,j] = temp1_init
     x0[size_total_x-1,j] = temp1_init
   end

   #Setup temp2_init inside
   for i = 3:size_total_x-2
     for j = 3:size_total_y-2
       x0[i,j] = temp2_init
     end
   end

   return
end

function neighbors(my_id::Int, nproc::Int, nx_domains::Int, ny_domains::Int)
    #find all processes that are my neighbors
    id_pos = Array{Int,2}(undef, nx_domains, ny_domains)
    for id = 0:nproc-1
        n_row = (id+1) % nx_domains > 0 ? (id+1) % nx_domains : nx_domains
        n_col = ceil(Int, (id + 1) / nx_domains)
        if (id == my_id)
            global my_row = n_row
            global my_col = n_col
        end
        id_pos[n_row, n_col] = id
    end

    neighbor_N = my_row + 1 <= nx_domains ? my_row + 1 : -1
    neighbor_S = my_row - 1 > 0 ? my_row - 1 : -1
    neighbor_E = my_col + 1 <= ny_domains ? my_col + 1 : -1
    neighbor_W = my_col - 1 > 0 ? my_col - 1 : -1

    neighbors = Dict{String,Int}()
    neighbors["N"] = neighbor_N >= 0 ? id_pos[neighbor_N, my_col] : -1
    neighbors["S"] = neighbor_S >= 0 ? id_pos[neighbor_S, my_col] : -1
    neighbors["E"] = neighbor_E >= 0 ? id_pos[my_row, neighbor_E] : -1
    neighbors["W"] = neighbor_W >= 0 ? id_pos[my_row, neighbor_W] : -1

    return neighbors
end

function process_coordinates!(xs::Array{Int}, ys::Array{Int}, xe::Array{Int},
    ye::Array{Int}, xcell::Int, ycell::Int, nx_domains::Int, ny_domains::Int, nproc::Int)

    #Computation of starting ys,ye on (Ox) standard axis
    #for the first column of global domain,
    #Convention x(i,j) with i row and j column
    ys[1:nx_domains] .= 3
    #Here, ye[1:x_domains] = 3 + ycell-1
    ye[1:nx_domains] = ys[1:nx_domains] .+ ycell .- 1

    #Computation of ys,ye on (Ox) standard axis
    #for all other cells of global domain
    for i = 1:ny_domains-1
        ys[i*nx_domains+1:(i+1)*nx_domains] = ys[(i-1)*nx_domains+1:i*nx_domains] .+ ycell .+ 2
        ye[i*nx_domains+1:(i+1)*nx_domains] = ys[i*nx_domains+1:(i+1)*nx_domains] .+ ycell .- 1
    end

    #Computation of starting xs,xe on (Oy) standard axis
    #for the first row of global domain,
    #Convention x(i,j) with i row and j column
    for i = 1:ny_domains
        xs[(i-1)*nx_domains+1] = 3
        #Here, xe(i*x_domains) = 2+xcell-1
        xe[(i-1)*nx_domains+1] = xs[(i-1)*nx_domains+1] + xcell - 1
    end

    #Computation of xs,xe on (Oy) standard axis
    #for all other cells of global domain
    for i = 1:ny_domains
        for j = 2:nx_domains
            xs[(i-1)*nx_domains+j] = xs[(i-1)*nx_domains+(j-1)] + xcell + 2
            xe[(i-1)*nx_domains+j] = xs[(i-1)*nx_domains+j] + xcell - 1
        end
    end

    return xs, xe, ys, ye
end

max_stars_repo_name: ornl-training/julia-basics

max_issues_repo_name: ornl-training/julia-basics

max_forks_repo_name: ornl-training/julia-basics


text:
using HMMs
using Brobdingnag
using Distributions
using Random

Random.seed!(1)
x = rand( 10, 5 )
name = tempname()
file = open( name, "w" )
HMMs.writearray( file, x )
close(file)
file = open( name, "r" )
y = HMMs.readarray( file, Matrix{Float64} )
close(file)
rm(name)
@assert( x == y )

N = 2
T = 1000
Random.seed!(1)
hmm1 = HMMs.rand( HMMs.HMM{N,Normal,BigFloat,Float64} )
y1 = rand( hmm1, T );
HMMs.setobservations( hmm1, y1 );
alpha = HMMs.forwardprobabilities( hmm1 );
beta = HMMs.backwardprobabilities( hmm1 );

@assert(maximum(abs.(diff([sum(alpha[i,:].*beta[i,:]) for i in 1:T]))) < 1e-8)

pi1 = HMMs.stationary( hmm1 )
@assert( maximum(abs.(pi1 * hmm1.transitionprobabilities - pi1)) < 1e-8 )

gamma = HMMs.conditionalstateprobabilities( hmm1 );
xi = HMMs.conditionaljointstateprobabilities( hmm1 );
@assert( maximum(abs.([maximum(abs.(sum(xi[i,:,:],dims=2) - gamma[i,:])) for i in 1:T-1])) < 1e-8 )
@assert( maximum(abs.([maximum(abs.(sum(xi[i,:,:],dims=1)' - gamma[i+1,:])) for i in 1:T-1])) < 1e-8 )

hmm2 = deepcopy( hmm1 );
@time HMMs.emstep( hmm1, hmm2 )
@assert( maximum(abs.( sum(hmm2.transitionprobabilities,dims=2) .- 1 )) < 1e-8 )

Random.seed!(2)
hmm3 = HMMs.rand( HMMs.HMM{2,Normal,BigFloat,Float64} )
y2 = rand( hmm1, 100000 );
HMMs.setobservations( hmm3, y2 );
@time HMMs.em( hmm3, debug=2 )

error = HMMs.permutederror( hmm1, hmm3 )
@assert( error.transitionprobabilities < 1e-2 )
@assert( error.means < 1e-2 )
@assert( error.stds < 1e-2 )

# check derivatives; need to do this before checking statistical convergence
# since MLE has 0 likelihood derivative
function ddffd( hmm, parameter, f, df; delta = 1e-6, relative=false )
    original = parameter[1]

    parameter[1] += delta
    HMMs.clear( hmm )
    dpos = copy( f( hmm ) )
    parameter[1] = original

    parameter[1] -= delta
    HMMs.clear( hmm )
    dneg = copy( f( hmm ) )
    parameter[1] = original

    fd = (dpos - dneg)/(2*delta)
    dims = length(size(fd))
    outtype = dims == 0 ? Float64 : Array{Float64,dims}
    result = convert(outtype,df .- fd)
    result = abs.(relative ? result ./ (1 .+ abs.(convert(outtype,df))) : result)
    return result
end

function testfd( hmm, parameter, f, df; delta = 1e-6, epsilon = 1e-4, relative=false )
    diffs = ddffd( hmm, parameter, f, df, delta=delta, relative=relative )
    (max, index) = findmax( diffs )
    location = join( Tuple(index), ", " )
    println( "Max is $max (out of $epsilon) at ($location)" )
    @assert( max < epsilon )
end

Random.seed!(2)
hmm4 = rand( HMMs.HMM{2,Normal,Brob,Float64} )
HMMs.setobservations( hmm4, y2 );

b = copy( HMMs.probabilities( hmm4 ) );
(T,m) = size(b)
dlogb = copy( HMMs.dlogprobabilities( hmm4 ) );
d2logb = copy( HMMs.d2logprobabilities( hmm4 ) );
d2b = copy( HMMs.d2probabilities( hmm4 ) );
alpha = copy( HMMs.forwardprobabilities( hmm4 ) );
dalpha = copy( HMMs.dforwardprobabilities( hmm4 ) );
d2alpha = copy( HMMs.d2forwardprobabilities( hmm4 ) );
l = copy( HMMs.likelihood( hmm4 ) );
dl = copy( HMMs.dlikelihood( hmm4 ) );
d2l = copy( HMMs.d2likelihood( hmm4 ) );
d2logl = copy( HMMs.d2loglikelihood( hmm4 ) );

f1 = m -> log.(HMMs.probabilities(m))
f2 = m -> permutedims(HMMs.probabilities(m) .* permutedims(HMMs.dlogprobabilities(m),[3,2,1]), [3,2,1])
f3 = m -> (HMMs.dlikelihood(m)'./HMMs.likelihood(m))'
# note the constraint to add to 1 is handled elsewhere
for index in 1:m*(m+2)
    if index <= m^2
        (i,j) = divrem(index - 1, m) .+ (1,1)

        parameter = view( hmm4.transitionprobabilities, i, j )
        println( "\n\nTesting transition probability ($i,$j)" )
    elseif index <= m*(m+1)
        i = index - m^2
        parameter = view( hmm4.stateparameters, 1, i )
        println( "\n\nTesting mean $i" )
    else
        i = index - m*(m+1)
        parameter = view( hmm4.stateparameters, 2, i )
        println( "\n\nTesting std $i" )
    end

    print( "dlogb: " )
    testfd( hmm4, parameter, f1, dlogb[index,:,:]' )
    print( "d2logb: " )
    testfd( hmm4, parameter, HMMs.dlogprobabilities, d2logb[index,:,:,:] )
    print( "d2b: " )
    testfd( hmm4, parameter, f2, d2b[index,:,:,:] )

    print( "dalpha: " )
    testfd( hmm4, parameter, HMMs.forwardprobabilities, dalpha[index,:,:]', epsilon=1e-2, relative=true )
    print( "d2alpha: " )
    testfd( hmm4, parameter, HMMs.dforwardprobabilities, d2alpha[index,:,:,:], epsilon=1e-2, relative=true )

    print( "dl: " )
    testfd( hmm4, parameter, HMMs.likelihood, dl[index,:], epsilon=1e-2 )
    print( "d2l: " )
    testfd( hmm4, parameter, HMMs.dlikelihood, d2l[index,:,:], epsilon=1e-2 )
    print( "d2logl: " )
    testfd( hmm4, parameter, f3, d2logl[index,:,:], epsilon=1e-2, delta=1e-4, relative=true )
end

@time HMMs.em( hmm4, debug=2 )
error = HMMs.permutederror( hmm1, hmm4 )
@assert( error.transitionprobabilities < 1e-2 )
@assert( error.means < 1e-2 )
@assert( error.stds < 1e-2 )

convert( Matrix{Float64}, HMMs.sandwich( hmm4 ) )

file = open( name, "w" )
write( file, hmm4 )
close(file)
file = open( name, "r" )
hmm4a = read( file, typeof(hmm4) )
close(file)
rm(name)
fields = [:initialprobabilities, :transitionprobabilities, :stateparameters]
@assert( all([==( getfield.( [hmm4,hmm4a], field )... ) for field in fields]) )

Random.seed!(1)
hmm5 = HMMs.rand( HMMs.HMM{3,Normal,Brob,Float64} )
y3 = rand( hmm5, 100000 );
Random.seed!(2)
hmm6 = HMMs.rand( HMMs.HMM{3,Normal,Brob,Float64} )
HMMs.setobservations( hmm6, y3 );
@time HMMs.em( hmm6, debug=2, maxiterations=100 )

Random.seed!(1)
hmm10 = HMMs.rand( HMMs.HMM{2,Laplace,Brob,Float64} )
y10 = rand( hmm10, 10_000 )

Random.seed!(2)
hmm11 = HMMs.rand( HMMs.HMM{2,Laplace,Brob,Float64} )
HMMs.setobservations( hmm11, y10 )
HMMs.em( hmm11, debug=2 )

HMMs.reorder!( hmm10 )
HMMs.reorder!( hmm11 )

Random.seed!(1)
hmm12 = HMMs.rand( HMMs.HMM{3,Laplace,Brob,Float64} )

Random.seed!(1)
hmm13 = HMMs.rand( HMMs.HMM{2,HMMs.GenTDist,Brob,Float64} )

y13 = rand( hmm13, 10_000 );

Random.seed!(2)
hmm14 = HMMs.rand( HMMs.HMM{2,HMMs.GenTDist,Brob,Float64} )
HMMs.setobservations( hmm14, y13 );

b = copy( HMMs.probabilities( hmm14 ) );
(T,m) = size(b)
dlogb = copy( HMMs.dlogprobabilities( hmm14 ) );
d2logb = copy( HMMs.d2logprobabilities( hmm14 ) );
d2b = copy( HMMs.d2probabilities( hmm14 ) );
alpha = copy( HMMs.forwardprobabilities( hmm14 ) );
dalpha = copy( HMMs.dforwardprobabilities( hmm14 ) );
d2alpha = copy( HMMs.d2forwardprobabilities( hmm14 ) );
l = copy( HMMs.likelihood( hmm14 ) );
dl = copy( HMMs.dlikelihood( hmm14 ) );
d2l = copy( HMMs.d2likelihood( hmm14 ) );
d2logl = copy( HMMs.d2loglikelihood( hmm14 ) );

for index in 1:m*(m+3)
    if index <= m^2
        (i,j) = divrem(index - 1, m) .+ (1,1)

        parameter = view( hmm14.transitionprobabilities, i, j )
        println( "\n\nTesting transition probability ($i,$j)" )
    elseif index <= m*(m+1)
        i = index - m^2
        parameter = view( hmm14.stateparameters, 1, i )
        println( "\n\nTesting mean $i" )
    elseif index <= m*(m+2)
        i = index - m*(m+1)
        parameter = view( hmm14.stateparameters, 2, i )
        println( "\n\nTesting std $i" )
    else
        i = index - m*(m+2)
        parameter = view( hmm14.stateparameters, 3, i )
        println( "\n\nTesting nu $i" )
    end

    print( "dlogb: " )
    testfd( hmm14, parameter, f1, dlogb[index,:,:]' )
    print( "d2logb: " )
    testfd( hmm14, parameter, HMMs.dlogprobabilities, d2logb[index,:,:,:] )
    print( "d2b: " )
    testfd( hmm14, parameter, f2, d2b[index,:,:,:] )

    print( "dalpha: " )
    testfd( hmm14, parameter, HMMs.forwardprobabilities, dalpha[index,:,:]', epsilon=1e-2, relative=true )
    print( "d2alpha: " )
    testfd( hmm14, parameter, HMMs.dforwardprobabilities, d2alpha[index,:,:,:], epsilon=1e-2, relative=true )

    print( "dl: " )
    testfd( hmm14, parameter, HMMs.likelihood, dl[index,:], epsilon=1e-2 )
    print( "d2l: " )
    testfd( hmm14, parameter, HMMs.dlikelihood, d2l[index,:,:], epsilon=1e-2 )
    print( "d2logl: " )
    testfd( hmm14, parameter, f3, d2logl[index,:,:], epsilon=1e-2, delta=1e-4, relative=true )
end


text:
using GridPotentialFlow
using Plots
Δx = 0.01
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);
c = Lx/2 # chord length
α = -π/3 # angle of attack
plate = Plate(c,4*cellsize(g));
Tr = RigidTransform((0.0,0.0),α)
Tr(plate)
Δs = dlengthmid(plate);

σLE_list = [0.0,0.05,0.1];
pfb_list = [PotentialFlowBody(plate,edges=[1,length(plate)],σ=[SuctionParameter(σLE),SuctionParameter(0.0)]) for σLE in σLE_list];

Δt = 2e-2
vLE = Vortex(plate.x[1]+3Δt*plate.len*cos(plate.α+π/2),plate.y[1]+3Δt*plate.len*sin(plate.α+π/2),0.0);
vTE = Vortex(plate.x[end]+3Δt*cos(α+π/2),plate.y[end]+3Δt*sin(α+π/2),0.0);

function createsheddedvortices(plate,oldvortices)

    vLE = Vortex(2/3*plate.x[1]+1/3*oldvortices[end-1].x,2/3*plate.y[1]+1/3*oldvortices[end-1].y,0.0)
    vTE = Vortex(2/3*plate.x[end]+1/3*oldvortices[end].x,2/3*plate.y[end]+1/3*oldvortices[end].y,0.0)

    return vLE, vTE
end

model_list = [VortexModel(g,bodies=[pfb],vortices=[vLE,vTE],U∞=(1.0,0.0)) for pfb in pfb_list];

T = 0:Δt:0.2
for t in T
    for vm in model_list
        X = getvortexpositions(vm)
        Ẋ = vortexvelocities!(vm)
        X .= X .+ Ẋ*Δt
        setvortexpositions!(vm, X)
        vLEnew, vTEnew = createsheddedvortices(plate,vm.vortices[end-1:end])
        pushvortices!(vm,vLEnew,vTEnew)
    end
end

colors = [:red,:blue,:green];
plot(plate,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")
for i in 1:length(model_list)
    plot!(model_list[i].vortices.x[4:2:end],model_list[i].vortices.y[4:2:end],color=colors[i],marker=:circle,markersize=2)
    plot!(model_list[i].vortices.x[3:2:end],model_list[i].vortices.y[3:2:end],color=colors[i],marker=:circle,markersize=2)
    scatter!(model_list[i].vortices.x[1:2],model_list[i].vortices.y[1:2],color=colors[i],marker=:circle,markersize=2,label="σLE=$(σLE_list[i])")
end
plot!()

plot(xlabel="body point index",ylabel="f̃")
for i in 1:length(model_list)
    sol = solve(model_list[i]);
    plot!(sol.f./model_list[i].system.f₀,linecolor=colors[i],label="σLE=$(σLE_list[i])")
end
plot!()

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl