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

using Markdown
using InteractiveUtils

# ╔═╡ 49630972-4192-11eb-3dd3-c14371066e1a
begin
	using Markdown
	using InteractiveUtils
	using LinearAlgebra
end

# ╔═╡ 54ef1f40-4192-11eb-3e9d-1b8d62a4d99f
begin
    using PlutoUI
    using Plots
    using Printf
    using Interact
end

# ╔═╡ 3a43bac0-4192-11eb-30c2-e1c845f3e7bb
md"_Compressible Similarity, version 1_"

# ╔═╡ 4072d430-4192-11eb-15ad-d551a4cc6523
md"""

## **Compressible Similarity Solution**

#### **Description:**

This notebook computes the compressible version of Falkner-Skan Similarity Solution coupled with energy equation

		selfsimilar(M∞, T∞, ηmax, N, itermax, ϵProfile, ϵBC)

If the arguments are missing, it will use the default values.

		selfsimilar(M∞=1, T∞=300, ηmax=10, N=50, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

#### **Compressible Similarity Equation**
Boundary-layer velocity and temperature profiles on the flat plate can be projected onto single profile wich is self-similar profile. It can be represented using the ordinary differential equations (ODEs) below:

$$(cf'')'+ff'' =0$$
$$(a_1g'+a_2f'f'')'+fg'=0$$

where

$$f'=\frac{u}{u_e}$$
$$c=\frac{\rho \mu}{\rho_e \mu_e}$$
$$g=\frac{H}{H_e}$$
$$a_1=\frac{c}{\sigma}$$
$$a_2=\frac{(\gamma-1)M^2}{1+(\frac{\gamma-1}{2})M^2}\left(1-\frac{1}{\sigma}\right)c$$

and H is the enthalpy, γ is the ratio of specific heats, M is the edge Mach number, and σ is the Prandtl number. σ and M can be defined as

$$M=\frac{u_e}{\sqrt{\gamma \mathfrak{R}T_e}}$$
$$\sigma=\frac{\mu c_p}{k}$$

In this code, σ is assumed as 0.72. The viscosity μ is a function of T and it is calculated as

\begin{equation}
    μ = c₁\frac{T^{3/2}}{T+c₂}
\end{equation}

c₂ is 110.4 Kelvin. c₁ is disappearing on the nondimensionalizing process. The boundary conditions for the system of ODEs are

$$y=0;  f=f'=0$$
$$y\rightarrow \infty;  f',g \rightarrow 0$$

The resultant equations along with the boundary conditions are solved with the Runge-Kutta Fehlberg scheme with Newton's iteration method for missing boundary condition.

Details of RK Fehlberg:
Numerical Recipes, Cambridge

Details of Similarity solution formulation:
Boundary-Layer Theory, 7ᵗʰ edition, Schlichting

Feel free to ask questions!

*Furkan Oz*

foz@okstate.edu
"""

# ╔═╡ 5b66da70-4192-11eb-3e5a-051e2c64153f
function Y1(η,y₁,y₂,y₃,y₄,y₅)
    return y₂
end

# ╔═╡ 5fff2740-4192-11eb-39b5-11f529ae8fa1
function Y2(η,y₁,y₂,y₃,y₄,y₅)
    return y₃
end

# ╔═╡ 64d910ee-4192-11eb-30ff-b730c8610bd9
function Y3(η,y₁,y₂,y₃,y₄,y₅,cμ,T∞)
    return -y₃*((y₅/(2*(y₄)))-(y₅/(y₄+cμ/T∞)))-y₁*y₃*((y₄+cμ/T∞)/(sqrt(y₄)*(1+cμ/T∞)));
end

# ╔═╡ 6a635ee0-4192-11eb-3d15-79ce34aa522b
function Y4(η,y₁,y₂,y₃,y₄,y₅)
    return y₅
end

# ╔═╡ 6ed84530-4192-11eb-2c32-9fd8d5c0684e
function Y5(η,y₁,y₂,y₃,y₄,y₅,cμ,T∞,Pr,γ,M∞)
    return -y₅^2*((0.5/y₄)-(1/(y₄+cμ/T∞)))-Pr*y₁*y₅/sqrt(y₄)*(y₄+cμ/T∞)/(1+cμ/T∞)-(γ-1)*Pr*M∞^2*y₃^2;
end

# ╔═╡ 564b53e0-4192-11eb-33c2-717cb09752bf
"""
Computes the compressible Similarity Solution
This notebook computes the compressible Similarity Solution

		selfsimilar(M∞, T∞, ηmax, N, itermax, ϵProfile, ϵBC)

If the arguments are missing, it will use the default values.

		selfsimilar(M∞=1, T∞=300, ηmax=10, N=50, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

Furkan Oz,
foz@okstate.edu,

"""
function selfsimilar(M∞=1, T∞=300, ηmax=10, N=50, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)
    Δη = ηmax/N
    Δη²= Δη^2
    γ  = 1.4   # Ratio of specific heat
    cμ = 110.4 # Sutherland law coefficient for [Kelvin]
    Pr = 0.72  # Prandtl Number

    Δ = 1e-7     # Small number for shooting method (Decrease when ϵ is decreased)

    # Initializing the solution vectors
    y₁ = zeros(N+1)   # f
    y₂ = zeros(N+1)   # f'
    y₃ = zeros(N+1)   # f''
    y₄ = zeros(N+1)   # ρ(η)
    y₅ = zeros(N+1)   # ρ(η)'
    η = [(i-1)*Δη for i=1:N+1]

    # Boundary Conditions for Isothermal Case
    y₁[1] = 0
    y₂[1] = 0
    y₅[1] = 0

    α₀ = 0.1 	  # Change that if code does not converge
    β₀ = 3.0      # Change that if code does not converge
    y₃[1] = α₀    # Initial Guess
    y₄[1] = β₀    # Initial Guess

    # We want
    # y2[N+1] = 1
    # y4[N+1] = 1

	# RK Fehlberg Coefficients see reference for details
    b₂ = 0.2; b₃ = 0.3;
    b₄ = 0.6; b₅ = 1.0; b₆ = 7.0/8.0;
    a₂₁= 0.2;
    a₃₁= 3.0/40.0; a₃₂= 9.0/40.0;
    a₄₁= 0.3; a₄₂=-0.9; a₄₃= 6.0/5.0;
    a₅₁=-11.0/54.0;a₅₂= 2.5;  a₅₃=-70.0/27.0; a₅₄= 35.0/27.0;
    a₆₁= 1631.0/55296.0; a₆₂= 175.0/512.0; a₆₃= 575.0/13824.0;
	a₆₄= 44275.0/110592.0; a₆₅= 253.0/4096.0;
    c₁ = 37.0/378.0; c₂ = 0; c₃ = 250.0/621.0; c₄ = 125.0/594.0;
    c₅ = 0; c₆ = 512.0/1771.0;

    iter = 0
    errorProfile = 1.
    errorBC = 1.
    normₒ = 0
    normₙ = 0

    while ϵProfile<=errorProfile && iter<itermax
        y₃[1] = α₀;    # Initial Guess
        y₄[1] = β₀;    # Initial Guess

		# First solution for Newton's iteration
        for i=1:N

            w₁ = Y5(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i],cμ,T∞,Pr,γ,M∞);
            w₂ = Y5(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₂₁*w₁),cμ,T∞,Pr,γ,M∞);
            w₃ = Y5(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₃₁*w₁+a₃₂*w₂),cμ,T∞,Pr,γ,M∞);
            w₄ = Y5(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃),cμ,T∞,Pr,γ,M∞);
            w₅ = Y5(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄),cμ,T∞,Pr,γ,M∞);
            w₆ = Y5(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅),cμ,T∞,Pr,γ,M∞);

            y₅[i+1] = y₅[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y4(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y4(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₂₁*w₁), y₅[i]);
            w₃ = Y4(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₅[i]);
            w₄ = Y4(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₅[i]);
            w₅ = Y4(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₅[i]);
            w₆ = Y4(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₅[i]);

            y₄[i+1] = y₄[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y3(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i],cμ,T∞);
            w₂ = Y3(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i]+Δη*(a₂₁*w₁), y₄[i], y₅[i],cμ,T∞);
            w₃ = Y3(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₄[i], y₅[i],cμ,T∞);
            w₄ = Y3(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₄[i], y₅[i],cμ,T∞);
            w₅ = Y3(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₄[i], y₅[i],cμ,T∞);
            w₆ = Y3(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₄[i], y₅[i],cμ,T∞);

            y₃[i+1] = y₃[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y2(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y2(η[i]+Δη*b₂, y₁[i], y₂[i]+Δη*(a₂₁*w₁), y₃[i], y₄[i], y₅[i]);
            w₃ = Y2(η[i]+Δη*b₃, y₁[i], y₂[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₃[i], y₄[i], y₅[i]);
            w₄ = Y2(η[i]+Δη*b₄, y₁[i], y₂[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₃[i], y₄[i], y₅[i]);
            w₅ = Y2(η[i]+Δη*b₅, y₁[i], y₂[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₃[i], y₄[i], y₅[i]);
            w₆ = Y2(η[i]+Δη*b₆, y₁[i], y₂[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₃[i], y₄[i], y₅[i]);

            y₂[i+1] = y₂[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y1(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y1(η[i]+Δη*b₂, y₁[i]+Δη*(a₂₁*w₁), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₃ = Y1(η[i]+Δη*b₃, y₁[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₄ = Y1(η[i]+Δη*b₄, y₁[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₅ = Y1(η[i]+Δη*b₅, y₁[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₆ = Y1(η[i]+Δη*b₆, y₁[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₂[i], y₃[i], y₄[i], y₅[i]);

            y₁[i+1] = y₁[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

        end

        # Storing the freestream values for Newton's iteration method
        y₂ₒ = y₂[N+1];
        y₄ₒ = y₄[N+1];

        # Small number addition for Newton's iteration method
        y₃[1] = α₀+Δ;  # Initial Guess + Small number
        y₄[1] = β₀;    # Initial Guess

		# Second solution for Newton's iteration
        for i=1:N

            w₁ = Y5(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i],cμ,T∞,Pr,γ,M∞);
            w₂ = Y5(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₂₁*w₁),cμ,T∞,Pr,γ,M∞);
            w₃ = Y5(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₃₁*w₁+a₃₂*w₂),cμ,T∞,Pr,γ,M∞);
            w₄ = Y5(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃),cμ,T∞,Pr,γ,M∞);
            w₅ = Y5(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄),cμ,T∞,Pr,γ,M∞);
            w₆ = Y5(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅),cμ,T∞,Pr,γ,M∞);

            y₅[i+1] = y₅[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y4(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y4(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₂₁*w₁), y₅[i]);
            w₃ = Y4(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₅[i]);
            w₄ = Y4(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₅[i]);
            w₅ = Y4(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₅[i]);
            w₆ = Y4(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₅[i]);

            y₄[i+1] = y₄[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y3(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i],cμ,T∞);
            w₂ = Y3(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i]+Δη*(a₂₁*w₁), y₄[i], y₅[i],cμ,T∞);
            w₃ = Y3(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₄[i], y₅[i],cμ,T∞);
            w₄ = Y3(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₄[i], y₅[i],cμ,T∞);
            w₅ = Y3(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₄[i], y₅[i],cμ,T∞);
            w₆ = Y3(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₄[i], y₅[i],cμ,T∞);

            y₃[i+1] = y₃[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y2(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y2(η[i]+Δη*b₂, y₁[i], y₂[i]+Δη*(a₂₁*w₁), y₃[i], y₄[i], y₅[i]);
            w₃ = Y2(η[i]+Δη*b₃, y₁[i], y₂[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₃[i], y₄[i], y₅[i]);
            w₄ = Y2(η[i]+Δη*b₄, y₁[i], y₂[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₃[i], y₄[i], y₅[i]);
            w₅ = Y2(η[i]+Δη*b₅, y₁[i], y₂[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₃[i], y₄[i], y₅[i]);
            w₆ = Y2(η[i]+Δη*b₆, y₁[i], y₂[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₃[i], y₄[i], y₅[i]);

            y₂[i+1] = y₂[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y1(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y1(η[i]+Δη*b₂, y₁[i]+Δη*(a₂₁*w₁), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₃ = Y1(η[i]+Δη*b₃, y₁[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₄ = Y1(η[i]+Δη*b₄, y₁[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₅ = Y1(η[i]+Δη*b₅, y₁[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₆ = Y1(η[i]+Δη*b₆, y₁[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₂[i], y₃[i], y₄[i], y₅[i]);

            y₁[i+1] = y₁[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

        end

        # Storing the freestream values for Newton's iteration method
        y₂ₙ₁ = y₂[N+1];
        y₄ₙ₁ = y₄[N+1];

        # Small number addition for Newton's iteration method
        y₃[1] = α₀;    # Initial Guess
        y₄[1] = β₀+Δ;  # Initial Guess + Small number

		# Third solution for Newton's iteration
        for i=1:N

            w₁ = Y5(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i],cμ,T∞,Pr,γ,M∞);
            w₂ = Y5(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₂₁*w₁),cμ,T∞,Pr,γ,M∞);
            w₃ = Y5(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₃₁*w₁+a₃₂*w₂),cμ,T∞,Pr,γ,M∞);
            w₄ = Y5(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃),cμ,T∞,Pr,γ,M∞);
            w₅ = Y5(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄),cμ,T∞,Pr,γ,M∞);
            w₆ = Y5(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅),cμ,T∞,Pr,γ,M∞);

            y₅[i+1] = y₅[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y4(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y4(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₂₁*w₁), y₅[i]);
            w₃ = Y4(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₅[i]);
            w₄ = Y4(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₅[i]);
            w₅ = Y4(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₅[i]);
            w₆ = Y4(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i], y₄[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₅[i]);

            y₄[i+1] = y₄[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y3(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i],cμ,T∞);
            w₂ = Y3(η[i]+Δη*b₂, y₁[i], y₂[i], y₃[i]+Δη*(a₂₁*w₁), y₄[i], y₅[i],cμ,T∞);
            w₃ = Y3(η[i]+Δη*b₃, y₁[i], y₂[i], y₃[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₄[i], y₅[i],cμ,T∞);
            w₄ = Y3(η[i]+Δη*b₄, y₁[i], y₂[i], y₃[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₄[i], y₅[i],cμ,T∞);
            w₅ = Y3(η[i]+Δη*b₅, y₁[i], y₂[i], y₃[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₄[i], y₅[i],cμ,T∞);
            w₆ = Y3(η[i]+Δη*b₆, y₁[i], y₂[i], y₃[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₄[i], y₅[i],cμ,T∞);

            y₃[i+1] = y₃[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y2(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y2(η[i]+Δη*b₂, y₁[i], y₂[i]+Δη*(a₂₁*w₁), y₃[i], y₄[i], y₅[i]);
            w₃ = Y2(η[i]+Δη*b₃, y₁[i], y₂[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₃[i], y₄[i], y₅[i]);
            w₄ = Y2(η[i]+Δη*b₄, y₁[i], y₂[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₃[i], y₄[i], y₅[i]);
            w₅ = Y2(η[i]+Δη*b₅, y₁[i], y₂[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₃[i], y₄[i], y₅[i]);
            w₆ = Y2(η[i]+Δη*b₆, y₁[i], y₂[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₃[i], y₄[i], y₅[i]);

            y₂[i+1] = y₂[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

            w₁ = Y1(η[i], y₁[i], y₂[i], y₃[i], y₄[i], y₅[i]);
            w₂ = Y1(η[i]+Δη*b₂, y₁[i]+Δη*(a₂₁*w₁), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₃ = Y1(η[i]+Δη*b₃, y₁[i]+Δη*(a₃₁*w₁+a₃₂*w₂), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₄ = Y1(η[i]+Δη*b₄, y₁[i]+Δη*(a₄₁*w₁+a₄₂*w₂+a₄₃*w₃), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₅ = Y1(η[i]+Δη*b₅, y₁[i]+Δη*(a₅₁*w₁+a₅₂*w₂+a₅₃*w₃+a₅₄*w₄), y₂[i], y₃[i], y₄[i], y₅[i]);
            w₆ = Y1(η[i]+Δη*b₆, y₁[i]+Δη*(a₆₁*w₁+a₆₂*w₂+a₆₃*w₃+a₆₄*w₄+a₆₅*w₅), y₂[i], y₃[i], y₄[i], y₅[i]);

            y₁[i+1] = y₁[i] + Δη*(c₁*w₁+c₂*w₂+c₃*w₃+c₄*w₄+c₅*w₅+c₆*w₆);

        end

        # Storing the freestream values for Newton's iteration method
        y₂ₙ₂ = y₂[N+1];
        y₄ₙ₂ = y₄[N+1];

        # Calculation of the next initial guess with Newton's iteration method
        p₁₁ = (y₂ₙ₁-y₂ₒ)/Δ;
        p₂₁ = (y₄ₙ₁-y₄ₒ)/Δ;
        p₁₂ = (y₂ₙ₂-y₂ₒ)/Δ;
        p₂₂ = (y₄ₙ₂-y₄ₒ)/Δ;
        r₁ = 1-y₂ₒ;
        r₂ = 1-y₄ₒ;
        Δα = (p₂₂*r₁-p₁₂*r₂)/(p₁₁*p₂₂-p₁₂*p₂₁);
        Δβ = (p₁₁*r₂-p₂₁*r₁)/(p₁₁*p₂₂-p₁₂*p₂₁);
        α₀ = α₀ + Δα;
        β₀ = β₀ + Δβ;

        # Profile change between iteration
        normₙ = norm(y₂)
        errorProfile = maximum(abs.(normₙ-normₒ))

		# Convergence of boundary condition
        errorBC = abs(y₂[N+1]-1.)
        iter += 1
        normₒ = normₙ
        @printf("%4.4d %16.6e %16.6e \n", iter, errorProfile, errorBC)

    end

    if errorProfile<=ϵProfile
        println("")
        println("Solution converged!")
        println("The maximum change between consecutive profiles is less than the error criteria ϵProfile=$ϵProfile.")
    end

    if errorBC<=ϵBC
        println("")
        println("Solution for the boundary condition converged!")
        println("The difference between h(N) and h(N+1) is less than the error criteria ϵBC=$ϵBC.")
    end

	# Copying values for logical names
    U = y₂
    T = y₄

    # Integration for η --> y transformation
    y = zeros(N+1);
    for i=2:N+1
       y[i] = y[i-1] + y₄[i]*(η[i]-η[i-1]);
    end
    y = y*sqrt(2);

    return η,y,U,T,N
end

# ╔═╡ 7830b222-4192-11eb-18b6-5548ffd5f535
begin
	η,y,U,T,N = selfsimilar();
	plot(U,η,
	        title = "Similarity (M=1 - T=300K)",
	        label = "U",
	        legend = :topleft,
	        xlabel = "U",
	        ylabel = "η",
	        linewidth = 2,
	        linecolor = :black,
	        markershape = :circle,
	        markercolor = :red,
	    )
end

# ╔═╡ 90bc1cd0-4192-11eb-352b-572716d665fe
plot(T,η,
        title = "Similarity (M=1 - T=300K)",
        label = "T",
        legend = :topleft,
        xlabel = "T",
        ylabel = "η",
        linewidth = 2,
        linecolor = :black,
        markershape = :circle,
        markercolor = :red,
    )

# ╔═╡ 92859230-4192-11eb-0029-9f6dcb343b1e
plot(U,y,
        title = "Similarity Solution (M=1 - T=300K)",
        label = "U",
        legend = :topleft,
        xlabel = "U",
        ylabel = "y/√(νx/U)",
        linewidth = 2,
        linecolor = :black,
        markershape = :circle,
        markercolor = :red,
    )

# ╔═╡ Cell order:
# ╟─3a43bac0-4192-11eb-30c2-e1c845f3e7bb
# ╟─4072d430-4192-11eb-15ad-d551a4cc6523
# ╟─49630972-4192-11eb-3dd3-c14371066e1a
# ╟─54ef1f40-4192-11eb-3e9d-1b8d62a4d99f
# ╟─564b53e0-4192-11eb-33c2-717cb09752bf
# ╟─5b66da70-4192-11eb-3e5a-051e2c64153f
# ╟─5fff2740-4192-11eb-39b5-11f529ae8fa1
# ╟─64d910ee-4192-11eb-30ff-b730c8610bd9
# ╟─6a635ee0-4192-11eb-3d15-79ce34aa522b
# ╟─6ed84530-4192-11eb-2c32-9fd8d5c0684e
# ╟─7830b222-4192-11eb-18b6-5548ffd5f535
# ╟─90bc1cd0-4192-11eb-352b-572716d665fe
# ╟─92859230-4192-11eb-0029-9f6dcb343b1e


text:
using BundleMethod
using SparseArrays
using JuMP
using Ipopt
using Random

const BM = BundleMethod

#=
This example considers:
    minimize \sum_{i=1}^N \sum_{j=1}^n b_i * (x_j - a_ij)^2
    subject to -1 <= x_j <= 1
where objective function is separable such that
    f_i(x) := \sum_{j=1}^n b_i * (x_j - a_ij)^2 =#

# Randomly generate problem data
Random.seed!(1)
n = 3
N = 2
a = rand(N, n)
b = rand(N)

#=
User-defined function:

This function takes a new trial point `y` of dimension `n`
  and returns the array of separable objective function values
  and gradient of the functions. =#
function evaluate_f(y)
    N, n = size(a)
    fvals = zeros(N)
    grads = Dict{Int,SparseVector{Float64}}()
    for i = 1:N
        grad = zeros(n)
        for j = 1:n
            fvals[i] += b[i] * (y[j] - a[i,j])^2
            grad[j] += 2 * b[i] * (y[j] - a[i,j])
        end
        grads[i] = grad
    end
    return fvals, grads
end

# This initializes the trust region bundle method with required arguments.
pm = BM.TrustRegionMethod(n, N, evaluate_f)
BM.set_bundle_tolerance!(pm, 1e-6)

# Set optimization solver to the internal JuMP.Model
model = BM.get_jump_model(pm)
set_optimizer(model, Ipopt.Optimizer)
set_optimizer_attribute(model, "print_level", 0)

# We overwrite the function to have column bounds.
function BM.add_variables!(method::BM.TrustRegionMethod)
    bundle = BM.get_model(method)
    model = BM.get_model(bundle)
    @variable(model, -1 <= x[i=1:bundle.n] <= 1)
    @variable(model, θ[j=1:bundle.N])
end

# This builds the bundle model.
BM.build_bundle_model!(pm)

# This runs the bundle method.
BM.run!(pm)

ext: jl


text:
module InsidePallet

include("./constants.jl")
using JuMP
using .Constants

export cons_inside_pallet

#Constraints ensuring that the only allowed rotations are those of Currence Robotics
function cons_inside_pallet(model)

    #Check not exceeding in X-axis
    @constraint(model,
    xnel_x_not_exceed_length[i=1:length(boxes)],
    (model[:flbx_front_left_bot_x][i] +
    boxes[i][1]*model[:lpx_length_of_box_parallel_to_x][i] +
    boxes[i][2]*model[:wpx_width_of_box_parallel_to_x][i] +
    boxes[i][3]*model[:hpx_height_of_box_parallel_to_x][i] <=
    L)
    )

    #Check not exceeding in Y-axis
    @constraint(model,
    ynew_y_not_exceed_width[i=1:length(boxes)],
    (model[:flby_front_left_bot_y][i] +
    boxes[i][1]*model[:lpy_length_of_box_parallel_to_y][i] +
    boxes[i][2]*model[:wpy_width_of_box_parallel_to_y][i] +
    boxes[i][3]*model[:hpy_height_of_box_parallel_to_y][i] <=
    W)
    )

    #Check not exceeding in Z-axis
    @constraint(model,
    zneh_z_not_exceed_height[i=1:length(boxes)],
    (model[:flbz_front_left_bot_z][i] +
    boxes[i][1]*model[:lpz_length_of_box_parallel_to_z][i] +
    boxes[i][2]*model[:wpz_width_of_box_parallel_to_z][i] +
    boxes[i][3]*model[:hpz_height_of_box_parallel_to_z][i] <=
    H)
    )

    return model
end

end


text:
module SpectralKit

export domain_extrema, roots, augmented_extrema, evaluate, Chebyshev, ChebyshevSemiInf,
    ChebyshevInf, ChebyshevInterval

using ArgCheck: @argcheck
using DocStringExtensions: FUNCTIONNAME, SIGNATURES, TYPEDEF
using Parameters: @unpack

####
#### utilities
####

flip_if(condition, x) = condition ? -x : x

one_like(::Type{T}) where {T} = cos(zero(T))

"""
$(TYPEDEF)

Abstract type for function families.

# Supported interface

For `F::FunctionFamily`, the supported API functions are

- [`domain_extrema`](@ref) for querying the domain,

- [`evaluate`](@ref) for function evaluation,

- [`roots`](@ref) and [`augmented_extrema`](@ref) to obtain collocation points.
"""
abstract type FunctionFamily end

Broadcast.broadcastable(family::FunctionFamily) = Ref(family)

"""
`$(FUNCTIONNAME)(family)`

Return the extrema of the domain of the given `family` as an ordered tuple.

Type is not specified, but guaranteed to be the same for both endpoints.
"""
function domain_extrema end

"""
`$(FUNCTIONNAME)(family, n, x, order)`

Evaluate the `n`th (starting from 0) function in `family` at `x`.

`order` determines the derivatives:

- `Val(d)` returns the `d`th derivatives, starting from `0` (for the function value)

- `Val(0:d)` returns the derivatives up to `d`, starting from the function value, as
   multiple values (ie a tuple).

The implementation is intended to be type stable.

Consequences are undefined for evaluating outside the domain.
"""
function evaluate end

"""
`$(FUNCTIONNAME)([T], family, N)`

Return the roots of the `N`th function in `family`, as a vector of `N` numbers with element
type `T` (default `Float64`).

In the context of collocation, this is also known as the “Gauss-Chebyshev” grid.

Order is monotone, but not guaranteed to be increasing.
"""
roots(family, N::Integer) = roots(Float64, family, N)

"""
`$(FUNCTIONNAME)([T], family, N)`

Return the augmented extrema (extrema + boundary values) of the `N-1`th function in
`family`, as a vector of `N` numbers with element type `T` (default `Float64`).

In the context of collocation, this is also known as the “Gauss-Lobatto” grid.

Order is monotone, but not guaranteed to be increasing.
"""
function augmented_extrema(family::FunctionFamily, N::Integer)
    augmented_extrema(Float64, family, N)
end

####
#### Chebyshev on [-1,1]
####

"""
$(TYPEDEF)

Chebyhev polynomials of the first kind, defined on `[-1,1]`.
"""
struct Chebyshev <: FunctionFamily end

@inline domain_extrema(::Chebyshev) = (-1, 1)

"""
$(SIGNATURES)

Evaluate the `n`th Chebyshev polynomial at `-1` for the given orders.
"""
chebyshev_min(::Type{T}, n::Integer, ::Val{0}) where {T} = flip_if(isodd(n), one_like(T))

"""
$(SIGNATURES)

Evaluate the `n`th Chebyshev polynomial at `1` for the given orders.
"""
chebyshev_max(::Type{T}, n::Integer, ::Val{0}) where {T} = one_like(T)

function chebyshev_min(::Type{T}, n::Integer, ::Val{1}) where {T}
    flip_if(iseven(n), one_like(T) * abs2(n))
end

chebyshev_max(::Type{T}, n::Integer, ::Val{1}) where {T} = one_like(T) * abs2(n)

function chebyshev_min(::Type{T}, n::Integer, ::Val{0:1}) where {T}
    chebyshev_min(T, n, Val(0)), chebyshev_min(T, n, Val(1))
end

function chebyshev_max(::Type{T}, n::Integer, ::Val{0:1}) where {T}
    chebyshev_max(T, n, Val(0)), chebyshev_max(T, n, Val(1))
end

"""
$(SIGNATURES)

Evaluate the `n`th Chebyshev polynomial at `-1 < x < 1` for the given orders.
"""
chebyshev_interior(n::Integer, x::Real, ::Val{0}) = cos(n * acos(x))

function chebyshev_interior(n::Integer, x::Real, ::Val{1})
    t = acos(x)
    n * sin(n * t) / sin(t)
end

function chebyshev_interior(n::Integer, x::Real, ::Val{0:1})
    t = acos(x)
    t′ = 1 / sin(t)
    s, c = sincos(n * t)
    c, s * n * t′
end

function evaluate(family::Chebyshev, n::Integer, x::T, order) where {T <: Real}
    if x == -1
        chebyshev_min(T, n, order)
    elseif x == 1
        chebyshev_max(T, n, order)
    else
        chebyshev_interior(n, x, order)
    end
end

function roots(::Type{T}, ::Chebyshev, N) where {T <: Real}
    cospi.(((2 * N - 1):-2:1) ./ T(2 * N))
end

function augmented_extrema(::Type{T}, family::Chebyshev, N) where {T}
    cospi.(((N-1):-1:0) ./ T(N - 1))
end

####
#### Transformation of the Chebyshev polynomials — generic code
####

"""
$(TYPEDEF)

Define a function family by transforming from the Chebyshev polynomials on `[-1,1]`.

For `family::TransformedChebyshev`, subtypes implement:

- `from_chebyshev(family, x)`, for transforming from `x ∈ [-1,1]` to `y` of the domain,

- `to_chebyshev(family, y, order)` for transforming `y` the domain to `x ∈ [-1,1]`, where
  `order` follows the semantics of [`evaluate`](@ref) and returns `x`, `∂x/∂y`, … as
  requested.

- `domain_extrema(family)` is optional.
"""
abstract type TransformedChebyshev <: FunctionFamily end

function roots(::Type{T}, family::TransformedChebyshev, N) where {T}
    from_chebyshev.(family, roots(T, Chebyshev(), N))
end

function augmented_extrema(::Type{T}, family::TransformedChebyshev, N) where {T}
    from_chebyshev.(family, augmented_extrema(T, Chebyshev(), N))
end

function evaluate(family::TransformedChebyshev, n::Integer, x::Real, order::Val{0})
    evaluate(Chebyshev(), n, to_chebyshev(family, x, order), order)
end

function evaluate(family::TransformedChebyshev, n::Integer, x::Real, ::Val{1})
    x, x′ = to_chebyshev(family, x, Val(0:1))
    t′ = evaluate(Chebyshev(), n, x, Val(1))
    t′ * x′
end

function evaluate(family::TransformedChebyshev, n::Integer, x::Real, ::Val{0:1})
    x, x′ = to_chebyshev(family, x, Val(0:1))
    t, t′ = evaluate(Chebyshev(), n, x, Val(0:1))
    t, t′ * x′
end

function domain_extrema(family::TransformedChebyshev)
    from_chebyshev.(family, domain_extrema(Chebyshev()))
end

####
#### Rational Chebyshev on [A,Inf) or (-Inf,A]
####

"""
$(TYPEDEF)

Chebyshev polynomials transformed to the domain `[A, Inf)` (when `L > 0`) or `(-Inf,A]`
(when `L < 0`) using ``y = A + L * (1 + x) / (1 - x)``.
"""
struct ChebyshevSemiInf{T <: Real} <: TransformedChebyshev
    "The finite endpoint `A`."
    A::T
    "Scale factor `L ≠ 0`."
    L::T
    function ChebyshevSemiInf(A::T, L::T) where {T <: Real}
        @argcheck L ≠ 0
        new{T}(A, L)
    end
end

ChebyshevSemiInf(A::Real, L::Real) = ChebyshevSemiInf(promote(A, L)...)

function domain_extrema(TL::ChebyshevSemiInf)
    @unpack A, L = TL
    if L > 0
        promote(A, Inf)
    else
        promote(-Inf, A)
    end
end

function from_chebyshev(TL::ChebyshevSemiInf, x)
    TL.A + TL.L * (1 + x) / (1 - x)
end

function chebyshev_semiinf_endpoints(y, L, x::T) where {T}
    if (y == Inf && L > 0) || (y == -Inf && L < 0)
        one_like(T)
    else
        x
    end
end

function to_chebyshev(TL::ChebyshevSemiInf, y, ::Val{0})
    @unpack A, L = TL
    z = y - A
    x = (z - L) / (z + L)
    chebyshev_semiinf_endpoints(y, L, x)
end

function to_chebyshev(TL::ChebyshevSemiInf, y, ::Val{1})
    @unpack A, L = TL
    2 * L / abs2(y - A + L)
end

function to_chebyshev(TL::ChebyshevSemiInf, y, ::Val{0:1})
    @unpack A, L = TL
    z = y - A
    num = z - L
    den = z + L
    chebyshev_semiinf_endpoints(y, L, num/den), 2 * L / abs2(den)
end

####
#### Rational Chebyshev on (-Inf,Inf)
####

"""
$(TYPEDEF)

Chebyshev polynomials transformed to the domain `(-Inf, Inf)`.
(when `L < 0`) using ``y = A + L * x / √(L^2 + x^2)``.
"""
struct ChebyshevInf{T <: Real} <: TransformedChebyshev
    "The center `A`."
    A::T
    "Scale factor `L > 0`."
    L::T
    function ChebyshevInf(A::T, L::T) where {T <: Real}
        @argcheck L > 0
        new{T}(A, L)
    end
end

ChebyshevInf(A::Real, L::Real) = ChebyshevInf(promote(A, L)...)

domain_extrema(TB::ChebyshevInf) = (-Inf, Inf)

function from_chebyshev(TB::ChebyshevInf, x)
    TB.A + TB.L * x / √(1 - abs2(x))
end

function chebyshev_inf_endpoints(y, x::T) where {T}
    if y == Inf
        one_like(T)
    elseif y == -Inf
        -one_like(T)
    else
        x
    end
end

function to_chebyshev(TB::ChebyshevInf, y, ::Val{0})
    @unpack A, L = TB
    z = y - A
    x = z / hypot(L, z)
    chebyshev_inf_endpoints(y, x)
end

function to_chebyshev(TB::ChebyshevInf, y, ::Val{1})
    @unpack A, L = TB
    abs2(L) / hypot(L, z)^3
end

function to_chebyshev(TB::ChebyshevInf, y, ::Val{0:1})
    @unpack A, L = TB
    z = y - A
    den = hypot(L, z)
    x = z / den
    x′ = abs2(L) / den^3
    chebyshev_inf_endpoints(y, x), x′
end

####
#### Chebyshev on a finite interval `[a,b]`
####

"""
$(TYPEDEF)

Chebyshev polynomials transformed to the domain `(a, b)`. using ``y = x ⋅ s + m``.

`m` and `s` are calculated and checked by the constructor; `a < b` is enforced.
"""
struct ChebyshevInterval{T <: Real} <: TransformedChebyshev
    "Start of interval `a`."
    a::T
    "End of interval `b`."
    b::T
    "Midpoint `m`."
    m::T
    "Scale `s`."
    s::T
    function ChebyshevInterval(a::T, b::T) where {T <: Real}
        @argcheck isfinite(a) && isfinite(b)
        s = (b - a) / 2
        m = (a + b) / 2
        @argcheck s > 0 "Need `a < b`."
        new{T}(a, b, m, s)
    end
end

ChebyshevInterval(A::Real, L::Real) = ChebyshevInterval(promote(A, L)...)

domain_extrema(TI::ChebyshevInterval) = (TI.a, TI.b)

function from_chebyshev(TI::ChebyshevInterval, x)
    @unpack m, s = TI
    x * s + m
end

function to_chebyshev(TI::ChebyshevInterval, y, ::Val{0})
    @unpack m, s = TI
    (y - m) / s
end

to_chebyshev(TI::ChebyshevInterval, y, ::Val{1}) = 1 / TI.s

function to_chebyshev(TI::ChebyshevInterval, y, ::Val{0:1})
    @unpack m, s = TI
    (y - m) / s, 1 / s
end

end # module


text:
#=
   From Allen Downey
   https://twitter.com/AllenDowney/status/1063460117029535746
   """
   Today's Bayesian problem of the week: Suppose we visit a wild animal preserve where we
   know that the only animals are lions and tigers and bears, but we don't know how
   many of each there are.

   During the tour, we see 3 lions, 2 tigers, and 1 bear. Assuming that every animal had an equal
   chance to appear in our sample, estimate the prevalence of each species.

   What is the probability that the next animal we see is a bear?
   """

   Also see: https://towardsdatascience.com/estimating-probabilities-with-bayesian-modeling-in-python-7144be007815

   Cf ~/blog/lions_tigers_and_bears.blog
      ~/webppl/lions_tigers_and_bears.wppl

    From https://towardsdatascience.com/estimating-probabilities-with-bayesian-modeling-in-python-7144be007815
    """
    Species: lions    Prevalence: 44.44%.
    Species: tigers   Prevalence: 33.33%.
    Species: bears    Prevalence: 22.22%.
    """

=#

using Turing, StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)
include("jl_utils.jl")

@model function lions_tigers_and_bears()
    lion = 1
    tiger = 2
    bear = 3

    # Prior
    # We have a prior of alphas for the Dirichlet distribution.
    # We draw 6 times with the Multinomial distrib with the Dirichlet prior.
    # What is the probability of different combinations of the number of each animal?

    alphas = [1/3, 1/3, 1/3] # Highest prob x: [3, 2, 1]       0.2920000000000001
    # alphas = [3/6, 2/6, 1/6] # Highest prob x: [3, 2, 1]       0.4446000000000001
    # alphas = [2/6, 2/6, 2/6] # Highest prob x: [3, 2, 1]       0.26220000000000004
    # alphas = [1/6, 2/6, 3/6] # Highest prob x: [2, 2, 2]!      0.2613
    # alphas = [7/21, 7/21, 7/21] # Highest prob x: [2, 2, 2]       0.23359999999999997

    # Draw 6 animals
    n = 6
    # Using Multinomial seems to be much slower than using Dirichlet
    x ~ Multinomial(n,alphas)

    # The probabilities to calculate ("aliased" for simplicity)
    probLion  = x[1]
    probTiger = x[2]
    probBear  = x[3]

    n = 7
    o = tzeros(n)
    for i in 1:n
        o[i] ~ Categorical(simplex([probLion,probTiger,probBear]))
    end

    # It shouldn't matter in what order we see the different animals.
    true ~ Dirac(o[1] == lion)
    true ~ Dirac(o[2] == lion)
    true ~ Dirac(o[3] == lion)
    true ~ Dirac(o[4] == tiger)
    true ~ Dirac(o[5] == tiger)
    true ~ Dirac(o[6] == bear)

    # return o[7] == lion
    return o[7] == tiger
    # return o[7] == bear
    # return x

end

model = lions_tigers_and_bears()
num_chains = 4

# chains = sample(model, Prior(), 10_000)

chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, MH(), 10_000)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)
# chains = sample(model, PG(15), 1_000)

# chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 10_000)

# chains = sample(model, IS(), 1_000)

#
display(chains)
show_var_dist_pct(chains,:len,1000)

genq = generated_quantities(model, chains)
show_var_dist_pct(genq,1000)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
mutable struct NestedIntervalControl <: AbstractControl
    isenabled::Bool
end

function (control::NestedIntervalControl)(model::NestedInterval, properties::NestedIntervalDisplayProperties)
    id = get_id(properties)
    caption = get_caption(properties)
    col = get_col(properties)
    rectangle = get_layout(properties)
    # pos = get_pos(rectangle) % TODO Remove this
    totalwidth = get_width(rectangle)
    totalheight = get_height(rectangle)
    padding = get_padding(properties)
    interval = get_interval(model)
    start = get_start(model)
    stop = get_stop(model)
    CImGui.Dummy(ImVec2(0, padding[1]))
    CImGui.Indent(padding[2])
    width = totalwidth - padding[2]
    height = totalheight
    draw_list = CImGui.GetWindowDrawList()
    pos = CImGui.GetCursorScreenPos()
    CImGui.InvisibleButton("###Invisible Nested Interval Button", ImVec2(width, height))
    x₀ = stretch_linearly(start, first(interval), last(interval), pos.x, pos.x + width)
    x₁ = stretch_linearly(stop, first(interval), last(interval), pos.x, pos.x + width)
    io = CImGui.GetIO()
    if CImGui.IsItemHovered()
       mousepos = io.MousePos
       # Relax threshold when dragging to make things smoother
       if (abs(x₀ - mousepos.x) <=  25)
           CImGui.SetMouseCursor(CImGui.ImGuiMouseCursor_ResizeEW)
           if CImGui.IsItemActive()
               i₀ = stretch_linearly(io.MousePos.x, pos.x, pos.x + width, first(interval), last(interval))
               # Round to the nearest multiple of the interval stepsize
               s = step(interval)
               i₀′ = div(i₀ + (s / 2), s) * s
               set_start!(model, i₀′)
           end
       elseif (abs(x₁ - mousepos.x) <= 25)
           CImGui.SetMouseCursor(CImGui.ImGuiMouseCursor_ResizeEW)
           if CImGui.IsItemActive()
               i₁ = stretch_linearly(io.MousePos.x, pos.x, pos.x + width, first(interval), last(interval))
               # Round to the nearest multiple of the interval stepsize
               s = step(interval)
               i₁′ = div(i₁ + (s / 2), s) * s
               set_stop!(model,  i₁′)
           end
       end
       # Make sure that start is always less than or equal to stop
       if get_start(model) > get_stop(model)
           x₀′ = get_start(model)
           set_start!(model, get_stop(model))
           set_stop!(model, x₀′)
       end
    end
    CImGui.Unindent(padding[2])
    CImGui.SetCursorScreenPos(pos.x, pos.y)
    # The plotcontext allows us to plot data that appears underneath the nested interval control.
    plotcontext = get_plotcontext(properties)
    plotcontext()
    CImGui.AddRectFilled(draw_list, ImVec2(x₀, pos.y + 5), ImVec2(x₁,  pos.y + 5 + height), Base.convert(ImU32, ImVec4(0.0, 0, 0.99, 0.1)));
end


text:
"""
Creates the (normalized) half range Chebyshev polynomials of the first kind.

These orthonormal polynomials `(p_k(t))` on [-1,1] are orthogonal with respect to the weight
```
w(t) = [(1-t)(t-m(T))]^{-1/2}
```
with
```
m(T) = 1-2cosec^2(pi/(2T))
```

They appear in the context of Fourier extensions. There, an even function `f(x)` on
`[-1,1]` is approximated with a cosine series in `x` on `[-T,T]`. By the cosine mapping
```
    y = m(x) = cos(pi/T x)
```
this problem is transformed. The Fourier Extension problem on `f` is then equivalent with
the approximation of `f(m^{-1}(y))` with orthogonal polynomials `(q_k(y))` on `[cos(pi/T),1]`

By another mapping the interval `[cos(pi/T),1]` is mapped to `[-1,1]` and the polynomials `(q_k(y))` are mapped
to the polynomials `(p_k(t))`

See Huybrechs 2010, On the Fourier Extension of nonperiodic functions
    Adcock, Huybrechs, Vaquero 2014, On the numerical stability of Fourier Extension

The `indicator_function_nodes` refer to the boundary points of the domain of the function `f(x)` is approximated on.
It is assumed that `indicator_function_nodes` has the form
```
    [-1,b_1,a_2,b_2, ..., a_n, 1 ] with -1<b_1<a_2<b_2<...<a_n<1.
```
By default it is the vector [-1.,1.].

See also `HalfRangeChebyshevIkind`
"""
HalfRangeChebyshevIkind(n::Int, T::ELT, indicator_function_nodes::Vector{ELT}=default_indicator_nodes(ELT); options...) where ELT =
    HalfRangeChebyshev(n, ELT(-1//2), T, indicator_function_nodes; options...)

"""
Creates the (normalized) half range Chebyshev polynomials of the second kind.

These orthonormal polynomials `(p_k(t))` on [-1,1] are orthogonal with respect to the weight
```
w(t) = [(1-t)(t-m(T))]^{1/2}
```
with
```
m(T) = 1-2cosec^2(pi/(2T))
```

They appear in the context of Fourier extensions. There, an odd function `f(x)` on
`[-1,1]` is approximated with a sine series in `x` on `[-T,T]`.

See also `HalfRangeChebyshevIkind`
"""
HalfRangeChebyshevIIkind(n::Int, T::ELT, indicator_function_nodes::Vector{ELT}=default_indicator_nodes(ELT); options...) where ELT =
    HalfRangeChebyshev(n, ELT(1//2), T, indicator_function_nodes; options...)

m_forward = (x,T)->2*(cos(pi/T*x)-cos(pi/T))/(1-cos(pi/T))-1
m_inv = (t,T)->T/pi*acos(cos(pi/T)+(1-cos(pi/T))/2*(t+1))
weight_of_indicator(T,indicator) = t->.5*(indicator(-m_inv(t,T))+indicator(m_inv(t,T)))
default_indicator_nodes(ELT) = [-ELT(1),ELT(1)]

function HalfRangeChebyshev(n::Int, α, T::ELT, indicator_function_nodes::Vector{ELT}; options...) where ELT
    my_quadrature_rule = n->_halfrangechebyshevweights(n, α, T, indicator_function_nodes)
    OrthonormalOPSfromQuadrature(n, my_quadrature_rule, Interval(-one(ELT), one(ELT)); options...)
end

function _halfrangechebyshevweights(n, α::ELT, T::ELT, indicator_function_nodes::Vector{ELT}) where {ELT}
    @assert indicator_function_nodes[1] == -1 && indicator_function_nodes[end] == 1
    @assert reduce(&, indicator_function_nodes[1:end-1] .< indicator_function_nodes[2:end], init=true)
    @assert iseven(length(indicator_function_nodes))

    if α < 0
        C = 2T/pi
    else
        C = T*(1-cos(pi/T))^2/2/pi
    end
    if length(indicator_function_nodes)==2
        nodes, weights = gaussjacobi(n, α, ELT(0))
        Λ = weight_of_indicator(T,x->1)
        modified_weights = weights.*((nodes.-m_forward(T,T)).^α).*Λ.(nodes)*C
        nodes, modified_weights
    else
        indicator = indicator_function(indicator_function_nodes)
        mapped_nodes = sort(unique(map(x->m_forward(x,T), push!(indicator_function_nodes,0))))
        pop!(indicator_function_nodes)
        no_intervals = length(mapped_nodes)-1
        n_interval = n
        N = n_interval*no_intervals
        nodes = zeros(ELT,N)
        weights = zeros(ELT,N)

        for i in 1:no_intervals
            nodes_interval, weights_interval = gaussjacobi(n_interval, ELT(0), ELT(0))
            a = mapped_nodes[i]
            b = mapped_nodes[i+1]
            nodes_interval[:] .= a + (b-a)/2*(nodes_interval+1)
            weights_interval[:] .= weights_interval*(b-a)/2
            Δ = weight_of_indicator(T,indicator)
            weights_interval[:] .= C*weights_interval.*(1-nodes_interval).^α.*(nodes_interval-m_forward(T,T)).^α.*Δ.(nodes_interval)

            nodes[1+(i-1)*n_interval:i*n_interval] .= nodes_interval[:]
            weights[1+(i-1)*n_interval:i*n_interval] .= weights_interval[:]
        end
        nodes, weights
    end
end

indicator_function(nodes) = x-> reduce(|, false, nodes[1:2:end] .<= x .<= nodes[2:2:end])

function WaveOPS(n::Int,omega::ELT; options...) where {ELT}
    my_quadrature_rule = n->_wavePolynomialweight(n, omega)
    BasisFunctions.OrthonormalOPSfromQuadrature(n, my_quadrature_rule, Interval(-one(ELT), one(ELT)), (x->(abs(x)<=1) ? exp(1im*omega*x) : ELT(0)); options...)
end

function _wavePolynomialweight(n, omega::ELT) where ELT
    nodes, weights = gaussjacobi(n,ELT(0),ELT(0))
    weights = weights.*exp.(1im.*omega.*nodes)
    nodes, weights
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 28dec7bc-c846-11eb-2576-0d28a2569db6
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ bcdf3958-9a53-437f-8c0e-529b37e3e833
begin
	Pkg.add(["Plots", "PlutoUI","StatsBase", "DataStructures", "Unitful", "UnitfulRecipes", "PhysicalConstants", "UnitfulEquivalences", "DataFrames", "CSV", "Dates", "TimeZones", "StatsPlots", "ImageFiltering", "OffsetArrays", "FileIO", "StaticArrays", "FastIOBuffers", "TranscodingStreams", "CodecBzip2", "Tokenize", "DataInterpolations", "GR"])
	Pkg.develop(["RobotOSData"])

	using Plots
	gr()
	# plotly()
	using PlutoUI
	using StatsBase
	using DataStructures
	using Unitful
	using PhysicalConstants.CODATA2018: c_0, h, ħ
	using UnitfulEquivalences
	using UnitfulRecipes
	using DataFrames
	using CSV
	using TimeZones
	using Dates
	using StatsPlots
	using ImageFiltering
	using OffsetArrays
	using RobotOSData
	using FileIO
	using StaticArrays
	using DataInterpolations
end

# ╔═╡ c13e43d4-b7d6-4d65-996d-b2568d43b9f4
module ExtraMessages
    module tf2_msgs
        using RobotOSData.Messages # imports Readable, Header
        using RobotOSData.StdMsgs
        using RobotOSData.CommonMsgs # for sensor_msgs
        struct TFMessage <: Readable
	        transforms::Vector{geometry_msgs.TransformStamped}
        end
		struct AprilTagDetection <: Readable
			id::Vector{Int32}
			size::Vector{Float64}
			pose::Vector{geometry_msgs.PoseWithCovarianceStamped}
		end
		struct AprilTagDetectionArray <: Readable
			header::Header
			detections::Vector{AprilTagDetection}
		end
    end
end

# ╔═╡ ebbe541a-7ae5-4b18-b079-239bd733c5fb
bag = load("/home/kris/learning/robotont/bagfiles/bagfiles/new-square-10-tagged.bag", ExtraMessages)

# ╔═╡ a3730cdc-ad7e-4456-85f9-d963f0fde827
odom = bag["/robotont/odom"]

# ╔═╡ 433f3fe3-45c6-43c4-84f1-7b54178ed34a
bag.topic_map

# ╔═╡ 37eced46-67df-4956-b319-f3ef6cb8e827
oxy = [(convert(DateTime, o.data.header.time), o.data.pose.pose.position.y, o.data.pose.pose.position.x) for o = odom]

# ╔═╡ 26ce7dd3-24bb-46ff-8a9f-fea7bc1ad884
doxy = DataFrame(NamedTuple{(:time, :y, :x)}.(oxy))

# ╔═╡ c52dc88c-1d6a-4adf-ac1c-33414ae92b6b
@df doxy plot(:x, [:y])

# ╔═╡ 4efa0158-4098-4976-8590-a956632ef39c
twist = [(convert(DateTime, o.data.header.time), o.data.twist.twist.linear.x, o.data.twist.twist.linear.y, o.data.twist.twist.angular.z) for o = odom];

# ╔═╡ 63f85c0f-aa5e-4dfc-a770-940bc1dc7d29
dtw = DataFrame(NamedTuple{(:time, :x, :y, :phi)}.(twist));

# ╔═╡ ec189794-e5f9-447c-a0ef-18c1e34fbaa6
twist[1][1]

# ╔═╡ 9d500cb5-0344-45a7-b04c-0a3e23b51fbd
begin
	# simple integration
	results = zeros((size(twist)[1],3))
	x = 0
	y = 0
	phi = 0
	prev_time = twist[1][1]
	for (i,(newtime,vx,vy,vphi)) = enumerate(twist)
		dtime = (newtime - prev_time).value / 1000
		prev_time = newtime
		x += vx * dtime
		y += vy * dtime
		phi += vphi * dtime
		results[i, 1] = x
		results[i, 2] = y
		results[i, 3] = phi
	end
end

# ╔═╡ 74355b05-0bff-4d7a-99ea-473ce4768f4b
begin
	plot(results[:, 1], results[:, 2], label="path")
	plot!([results[end,1]], [results[end,2]], seriestype = :scatter, label="end")
	plot!([results[1,1]], [results[1,2]], seriestype = :scatter, label="start")
end

# ╔═╡ 10e08ba1-6e48-4d6c-8949-83e0c8ee82e7
tfs = [tf for msg ∈ bag[RobotOSData.ConnectionFilter(2), time_range] for tf ∈ msg.data.transforms ]

# ╔═╡ 0df647f5-340f-4b8a-9826-4ee55ea66890
countmap([tf.child_frame_id for tf ∈ tfs])

# ╔═╡ dd51603a-96cc-4149-80a0-6a0697b9384f

# ╔═╡ bf83571f-c7f5-45d9-87a9-23c2bcff04fe
zs = [(convert(DateTime, tf.header.time), tf.transform.translation.z) for tf ∈ tfs if tf.child_frame_id ∈ ("tag_5","tag_6")]

# ╔═╡ 8f96e1d2-32e4-46d8-998c-7728f4204fda
df = DataFrame(NamedTuple{(:time, :pos)}.(zs))

# ╔═╡ 79d2d3ed-d65d-44b2-aa6c-6964785a044a
mzs = combine(groupby(df, [:time]), :pos => mean => :mpos)

# ╔═╡ d42aa9df-4991-494e-a30d-24351ede8a0a
@df mzs plot(:time, [:mpos], ticks=:native)

# ╔═╡ fa2f217f-0184-4947-b161-191feba95d42
dist = mzs[!, :mpos] .* u"m" # add units

# ╔═╡ accbb017-a1c4-4122-b2ce-7711a946312a
time = mzs[!, :time] - mzs[1, :time]

# ╔═╡ 02e4c44e-42d7-4064-b522-1384122eee5e
plot(time[2:end], diff(dist) * u"1/s", title="Not smoothed", ylabel="speed")

# ╔═╡ 854e8a1a-3caa-4661-ba85-d0542f1489bb
begin
	smooth_range = -10:10
	n_ = size(smooth_range)[1]
	kernel = OffsetArray(fill(1/n_, n_), smooth_range)
end;

# ╔═╡ 4ae7323d-be6b-45c8-b16f-52c27968efc8
plot(time[2:end], diff(imfilter(dist, kernel)) * u"1/s", title="Smoothed with uniform kernel", ylabel="velocity")

# ╔═╡ c4054e36-b8a9-4743-8b36-bba84f29ec1b
velocity_gauss = diff(imfilter(dist, KernelFactors.gaussian(3))) * u"1/s"

# ╔═╡ faeb3515-b797-4920-b3ff-5672e47b1a51
plot(time[2:end], velocity_gauss, title="Smoothed with gaussian kernel", ylabel="velocity")

# ╔═╡ eb19d479-d778-4f9f-89b7-6825163769aa
acceleration_gauss = diff(imfilter(velocity_gauss, KernelFactors.gaussian(3))) * u"1/s"

# ╔═╡ cfe4a5f3-5d70-467b-9706-bebd17b55305
plot(time[3:end], acceleration_gauss, title="Smoothed with gaussian kernel", ylabel="acceleration")

# ╔═╡ 53699f51-6295-4f3e-8376-0891e4f3afd2

# ╔═╡ 14101c1e-f328-487b-a8a9-570e54fe7aa9
ox = [(convert(DateTime, o.data.header.time), o.data.twist.twist.linear.y, o.data.pose.pose.position.x) for o = odom]

# ╔═╡ 88f25939-55ba-4c26-9976-7f9e733f6966
odf = DataFrame(NamedTuple{(:time, :vel, :pos)}.(ox))

# ╔═╡ 9dfa7cda-1a14-43bf-9b0e-ca2fe8f0dcac
otime = odf[!, :time] - odf[1, :time]

# ╔═╡ 3c548f42-982c-49f7-9768-5f75771f28a7
plot(otime[1:end], odf[!, :vel])

# ╔═╡ 5f003d0a-f6ae-4bb9-8b8e-a1f7be829385
mean(odf[!,:vel])

# ╔═╡ eb52fa86-4c01-4100-91f4-eb94f4ad15d8
begin
	start = mzs[1,:mpos]

	plot(time, mzs[!, :mpos], ticks=:native, label="camera")
	plot!(otime, start .- odf[!, :pos], ticks=:native, label="odom")
end

# ╔═╡ d1929e40-a617-4577-a875-0b7e2f8878e6

# ╔═╡ Cell order:
# ╠═28dec7bc-c846-11eb-2576-0d28a2569db6
# ╠═bcdf3958-9a53-437f-8c0e-529b37e3e833
# ╠═ebbe541a-7ae5-4b18-b079-239bd733c5fb
# ╠═a3730cdc-ad7e-4456-85f9-d963f0fde827
# ╠═433f3fe3-45c6-43c4-84f1-7b54178ed34a
# ╠═37eced46-67df-4956-b319-f3ef6cb8e827
# ╠═26ce7dd3-24bb-46ff-8a9f-fea7bc1ad884
# ╠═c52dc88c-1d6a-4adf-ac1c-33414ae92b6b
# ╠═4efa0158-4098-4976-8590-a956632ef39c
# ╠═63f85c0f-aa5e-4dfc-a770-940bc1dc7d29
# ╠═ec189794-e5f9-447c-a0ef-18c1e34fbaa6
# ╠═9d500cb5-0344-45a7-b04c-0a3e23b51fbd
# ╠═74355b05-0bff-4d7a-99ea-473ce4768f4b
# ╠═10e08ba1-6e48-4d6c-8949-83e0c8ee82e7
# ╠═0df647f5-340f-4b8a-9826-4ee55ea66890
# ╠═c13e43d4-b7d6-4d65-996d-b2568d43b9f4
# ╠═dd51603a-96cc-4149-80a0-6a0697b9384f
# ╠═bf83571f-c7f5-45d9-87a9-23c2bcff04fe
# ╠═8f96e1d2-32e4-46d8-998c-7728f4204fda
# ╠═79d2d3ed-d65d-44b2-aa6c-6964785a044a
# ╠═d42aa9df-4991-494e-a30d-24351ede8a0a
# ╠═fa2f217f-0184-4947-b161-191feba95d42
# ╠═accbb017-a1c4-4122-b2ce-7711a946312a
# ╠═02e4c44e-42d7-4064-b522-1384122eee5e
# ╠═854e8a1a-3caa-4661-ba85-d0542f1489bb
# ╠═4ae7323d-be6b-45c8-b16f-52c27968efc8
# ╠═c4054e36-b8a9-4743-8b36-bba84f29ec1b
# ╠═faeb3515-b797-4920-b3ff-5672e47b1a51
# ╠═eb19d479-d778-4f9f-89b7-6825163769aa
# ╠═cfe4a5f3-5d70-467b-9706-bebd17b55305
# ╠═53699f51-6295-4f3e-8376-0891e4f3afd2
# ╠═14101c1e-f328-487b-a8a9-570e54fe7aa9
# ╠═88f25939-55ba-4c26-9976-7f9e733f6966
# ╠═9dfa7cda-1a14-43bf-9b0e-ca2fe8f0dcac
# ╠═3c548f42-982c-49f7-9768-5f75771f28a7
# ╠═5f003d0a-f6ae-4bb9-8b8e-a1f7be829385
# ╠═eb52fa86-4c01-4100-91f4-eb94f4ad15d8
# ╠═d1929e40-a617-4577-a875-0b7e2f8878e6


text:
module phasespace

export qfunc, wigner

using ..bases, ..states, ..operators, ..operators_dense, ..fock

"""
    qfunc(a, α)
    qfunc(a, x, y)
    qfunc(a, xvec, yvec)

Husimi Q representation ``⟨α|ρ|α⟩/π`` for the given state or operator `a`. The
function can either be evaluated on one point α or on a grid specified by
the vectors `xvec` and `yvec`. Note that conversion from `x` and `y` to `α` is
done via the relation ``α = \\frac{1}{\\sqrt{2}}(x + i y)``.
"""
function qfunc(rho::Operator, alpha::Number)
    b = basis(rho)
    @assert isa(b, FockBasis)
    _qfunc_operator(rho, convert(Complex128, alpha), Ket(b), Ket(b))
end

function qfunc(rho::Operator, xvec::Vector{Float64}, yvec::Vector{Float64})
    b = basis(rho)
    @assert isa(b, FockBasis)
    Nx = length(xvec)
    Ny = length(yvec)
    tmp1 = Ket(b)
    tmp2 = Ket(b)
    result = Matrix{Complex128}(Nx, Ny)
    @inbounds for j=1:Ny, i=1:Nx
        result[i, j] = _qfunc_operator(rho, complex(xvec[i], yvec[j])/sqrt(2), tmp1, tmp2)
    end
    result
end

function qfunc(psi::Ket, alpha::Number)
    b = basis(psi)
    @assert isa(b, FockBasis)
    _alpha = convert(Complex128, alpha)
    _conj_alpha = conj(_alpha)
    N = length(psi.basis)
    s = psi.data[N]/sqrt(N-1)
    @inbounds for i=1:N-2
        s = (psi.data[N-i] + s*_conj_alpha)/sqrt(N-i-1)
    end
    s = psi.data[1] + s*_conj_alpha
    return abs2(s)*exp(-abs2(_alpha))/pi
end

function qfunc(psi::Ket, xvec::Vector{Float64}, yvec::Vector{Float64})
    b = basis(psi)
    @assert isa(b, FockBasis)
    Nx = length(xvec)
    Ny = length(yvec)
    points = Nx*Ny
    N = length(b)::Int
    _conj_alpha = [complex(x, -y)/sqrt(2) for x=xvec, y=yvec]
    q = fill(psi.data[N]/sqrt(N-1), size(_conj_alpha))
    @inbounds for n=1:N-2
        f0_ = 1/sqrt(N-n-1)
        x = psi.data[N-n]
        for i=1:points
            q[i] = (x + q[i]*_conj_alpha[i])*f0_
        end
    end
    result = similar(q, Float64)
    x = psi.data[1]
    @inbounds for i=1:points
        result[i] = abs2(x + q[i]*_conj_alpha[i])*exp(-abs2(_conj_alpha[i]))/pi
    end
    result
end

function qfunc(state::Union{Ket, Operator}, x::Number, y::Number)
    qfunc(state, Complex128(x, y)/sqrt(2))
end

function _qfunc_operator(rho::Operator, alpha::Complex128, tmp1::Ket, tmp2::Ket)
    coherentstate(basis(rho), alpha, tmp1)
    operators.gemv!(complex(1.), rho, tmp1, complex(0.), tmp2)
    a = dot(tmp1.data, tmp2.data)
    return a/pi
end

"""
    wigner(a, α)
    wigner(a, x, y)
    wigner(a, xvec, yvec)

Wigner function for the given state or operator `a`. The
function can either be evaluated on one point α or on a grid specified by
the vectors `xvec` and `yvec`. Note that conversion from `x` and `y` to `α` is
done via the relation ``α = \\frac{1}{\\sqrt{2}}(x + i y)``.
"""
function wigner(rho::DenseOperator, x::Number, y::Number)
    b = basis(rho)
    @assert isa(b, FockBasis)
    N = b.N::Int
    _2α = complex(convert(Float64, x), convert(Float64, y))*sqrt(2)
    abs2_2α = abs2(_2α)
    w = complex(0.)
    coefficient = complex(0.)
    @inbounds for L=N:-1:1
        coefficient = 2*_clenshaw(L, abs2_2α, rho.data)
        w = coefficient + w*_2α/sqrt(L+1)
    end
    coefficient = _clenshaw(0, abs2_2α, rho.data)
    w = coefficient + w*_2α
    exp(-abs2_2α/2)/pi*real(w)
end

function wigner(rho::DenseOperator, xvec::Vector{Float64}, yvec::Vector{Float64})
    b = basis(rho)
    @assert isa(b, FockBasis)
    N = b.N::Int
    _2α = [complex(x, y)*sqrt(2) for x=xvec, y=yvec]
    abs2_2α = abs2.(_2α)
    w = zeros(_2α)
    b0 = similar(_2α)
    b1 = similar(_2α)
    b2 = similar(_2α)
    @inbounds for L=N:-1:1
        _clenshaw_grid(L, rho.data, abs2_2α, _2α, w, b0, b1, b2, 2)
    end
    _clenshaw_grid(0, rho.data, abs2_2α, _2α, w, b0, b1, b2, 1)
    @inbounds for i=eachindex(w)
        abs2_2α[i] = exp(-abs2_2α[i]/2)/pi.*real(w[i])
    end
    abs2_2α
end

wigner(psi::Ket, x, y) = wigner(dm(psi), x, y)
wigner(state, alpha::Number) = wigner(state, real(alpha)*sqrt(2), imag(alpha)*sqrt(2))

function _clenshaw_grid(L::Int, ρ::Matrix{Complex128},
                abs2_2α::Matrix{Float64}, _2α::Matrix{Complex128}, w::Matrix{Complex128},
                b0::Matrix{Complex128}, b1::Matrix{Complex128}, b2::Matrix{Complex128}, scale::Int)
    n = size(ρ, 1)-L-1
    points = length(w)
    if n==0
        f = scale*ρ[1, L+1]
        @inbounds for i=1:points
            w[i] = f + w[i]*_2α[i]/sqrt(L+1)
        end
    elseif n==1
        f1 = 1/sqrt(L+1)
        @inbounds for i=1:points
            w[i] = scale*(ρ[1, L+1] - ρ[2, L+2]*(L+1-abs2_2α[i])*f1) + w[i]*_2α[i]*f1
        end
    else
        f0 = sqrt(float((n+L-1)*(n-1)))
        f1 = sqrt(float((n+L)*n))
        f0_ = 1/f0
        f1_ = 1/f1
        fill!(b1, ρ[n+1, L+n+1])
        @inbounds for i=1:points
            b0[i] = ρ[n, L+n] - (2*n-1+L-abs2_2α[i])*f1_*b1[i]
        end
        @inbounds for k=n-2:-1:1
            b1, b2, b0 = b0, b1, b2
            x = ρ[k+1, L+k+1]
            a1 = -(2*k+1+L)
            a2 = -f0*f1_
            @inbounds for i=1:points
                b0[i] = x + (a1+abs2_2α[i])*f0_*b1[i] + a2*b2[i]
            end
            f1 , f1_ = f0, f0_
            f0 = sqrt((k+L)*k)
            f0_ = 1/f0
        end
        @inbounds for i=1:points
            w[i] = scale*(ρ[1, L+1] - (L+1-abs2_2α[i])*f0_*b0[i] - f0*f1_*b1[i]) + w[i]*_2α[i]*f0_
        end
    end
end

function _clenshaw(L::Int, abs2_2α::Float64, ρ::Matrix{Complex128})
    n = size(ρ, 1)-L-1
    if n==0
        return ρ[1, L+1]
    elseif n==1
        ϕ1 = -(L+1-abs2_2α)/sqrt(L+1)
        return ρ[1, L+1] + ρ[2, L+2]*ϕ1
    else
        f0 = sqrt(float((n+L-1)*(n-1)))
        f1 = sqrt(float((n+L)*n))
        f0_ = 1/f0
        f1_ = 1/f1
        b2 = complex(0.)
        b1 = ρ[n+1, L+n+1]
        b0 = ρ[n, L+n] - (2*n-1+L-abs2_2α)*f1_*b1
        @inbounds for k=n-2:-1:1
            b1, b2 = b0, b1
            b0 = ρ[k+1, L+k+1] - (2*k+1+L-abs2_2α)*f0_*b1 - f0*f1_*b2
            f1, f1_ = f0, f0_
            f0 = sqrt((k+L)*k)
            f0_ = 1/f0
        end
        return ρ[1, L+1] - (L+1-abs2_2α)*f0_*b0 - f0*f1_*b1
    end
end

end #module


text:
#-------------------------------------------------------------------------------
# File: ils_hprc.jl
# Description: ILS for the RENAULT Roadef 2005 challenge
#   inspired by work of
#   Celso C. Ribeiro, Daniel Aloise, Thiago F. Noronha,
#   Caroline Rocha, Sebastián Urrutia
#
# Date: November 25, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

"""
    perturbation_ils_hprc(solution::Solution, instance::Instance,
                          k::Int, critical_cars::Array{Int,1})

Removes `k` cars of `solution` and inserts them back in the sequence.
"""
function perturbation_ils_hprc(solution_init::Solution, instance::Instance,
                               k::Int, critical_cars::Array{Int,1})
    solution = deepcopy(solution_init)
    k = min(k, length(critical_cars))
    remove!(solution, instance, k, critical_cars)
    for i in 1:k
        greedy_add!(solution, instance, instance.nb_cars, 1)
    end
    return solution
end

cost_HPRC(solution::Solution, instance::Instance) = cost(solution, instance, 1)[1]

"""
    local_search_exchange_ils_hprc!(solution::Solution, instance::Instance,
                                    start_time::UInt)

Performs a local search on `solution` using only exchange moves with respect to
`instance`.
"""
function local_search_exchange_ils_hprc!(solution::Solution, instance::Instance,
                                         start_time::UInt)
    while TIME_PART_ILS * TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        overall_delta = 0
        b0 = instance.nb_late_prec_day + 1      #First car of the current production day
        for i in b0:instance.nb_cars
            hprc_current_car = instance.HPRC_keys[solution.sequence[i]]
            best_delta = 0
            new_possible_positions = Array{Int, 1}()
            for j in b0:instance.nb_cars
                if hprc_current_car != instance.HPRC_keys[solution.sequence[j]]
                    delta = cost_move_exchange(solution, i, j, instance,1)[1]
                    if delta < best_delta
                        empty!(new_possible_positions)
                        push!(new_possible_positions, j)
                        best_delta = delta
                    elseif delta == best_delta
                        push!(new_possible_positions, j)
                    end
                end
            end
            if !isempty(new_possible_positions)
                k = rand(new_possible_positions)
                move_exchange!(solution, i, k, instance)
            end
            overall_delta += best_delta
        end
        if overall_delta == 0
            break
        end
    end

    return solution
end

"""
    local_search_insertion_ils_hprc!(solution::Solution, instance::Instance,
                                     start_time::UInt)

Performs a local search on `solution` using only insertion moves with respect to
`instance`.
"""
function local_search_insertion_ils_hprc!(solution::Solution, instance::Instance,
                                          start_time::UInt)
    while TIME_PART_ILS * TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        overall_delta = 0
        b0 = instance.nb_late_prec_day + 1      #First car of the current production day
        for i in b0:instance.nb_cars
            best_delta = 0
            new_possible_positions = Array{Int, 1}()
            costs = cost_move_insertion(solution,i,instance,1)
            for j in b0:instance.nb_cars
                delta = costs[j, 1]
                if delta < best_delta
                    empty!(new_possible_positions)
                    push!(new_possible_positions, j)
                    best_delta = delta
                elseif delta == best_delta
                    push!(new_possible_positions, j)
                end
            end
            if !isempty(new_possible_positions)
                k = rand(new_possible_positions)
                move_insertion!(solution, i, k, instance)
            end
            overall_delta += best_delta
        end
        if overall_delta == 0
            break
        end
    end

    return solution
end

"""
    fast_local_search_exchange_ils_hprc!(solution::Solution, instance::Instance,
                                         critical_cars::Array{Int, 1}, start_time::UInt)

Performs a local search on `solution` using only exchange moves with respect to `instance`
for well chosen cars tagged in `critical_cars`.
"""
function fast_local_search_exchange_ils_hprc!(solution::Solution, instance::Instance,
                                              critical_cars::Array{Int, 1}, start_time::UInt)
    while !isempty(critical_cars) && (TIME_PART_ILS * TIME_LIMIT > (time_ns() - start_time) / 1.0e9)
        b0 = instance.nb_late_prec_day + 1      #First car of the current production day
        position_car_a = rand(critical_cars)
        hprc_current_car = instance.HPRC_keys[solution.sequence[position_car_a]]
        best_delta = 0
        new_possible_positions = Array{Int, 1}()
        for position_car_b in b0:instance.nb_cars
            if hprc_current_car != instance.HPRC_keys[solution.sequence[position_car_b]]
                delta = cost_move_exchange(solution, position_car_a, position_car_b, instance, 1)[1]
                if delta < best_delta
                    empty!(new_possible_positions)
                    push!(new_possible_positions, position_car_b)
                    best_delta = delta
                elseif delta == best_delta
                    push!(new_possible_positions, position_car_b)
                end
            end
        end
        if !isempty(new_possible_positions)
            k = rand(new_possible_positions)
            move_exchange!(solution, position_car_a, k, instance)
            # Update critical cars:
            #TODO
            deleteat!(critical_cars, findfirst(isequal(position_car_a), critical_cars))
        else
            deleteat!(critical_cars, findfirst(isequal(position_car_a), critical_cars))
        end
    end

    return solution
end

"""
    intensification_ils_hprc!(solution::Solution, instance::Instance,
                              start_time::UInt)

Performs a local search on `solution` using insertion moves first then exchange moves with
respect to `instance`.
"""
function intensification_ils_hprc!(solution::Solution, instance::Instance,
                                   start_time::UInt)
    local_search_insertion_ils_hprc!(solution, instance, start_time)
    local_search_exchange_ils_hprc!(solution, instance, start_time)
    return solution
end

"""
    restart_ils_hprc(solution::Solution, instance::Instance)

"""
function restart_ils_hprc(solution_init::Solution, instance::Instance)
    critical_cars = find_critical_cars(solution_init, instance, 1)
    solution = perturbation_ils_hprc(solution_init, instance, NBCAR_DIVERSIFICATION, critical_cars)
    return solution
end

"""
    ILS_HPRC(solution_init::Solution, instance::Instance, start_time::UInt)

Main function of the ILS metaheuristic. Improves the `solution_init` on its first objective
with respect to `instance`.
"""
function ILS_HPRC(solution_init::Solution, instance::Instance, start_time::UInt)
    i = 0                               # Number of iterations since the last improvement
    nb_strong_perturbation = 0                      # Number of restarts done for a solution
    solution = deepcopy(solution_init)
    solution_opt = deepcopy(solution_init)
    lastopt = deepcopy(solution_init)
    cond = 0 #TODO
    while cond < STOPPING_CRITERIA_ILS_HPRC && cost_HPRC(solution_opt, instance) != 0 && (TIME_PART_ILS * TIME_LIMIT > (time_ns() - start_time) / 1.0e9)
        critical_cars = find_critical_cars(solution, instance, 1)
        neighbor = perturbation_ils_hprc(solution, instance, NBCAR_PERTURBATION, critical_cars)
        critical_cars = find_critical_cars(neighbor, instance, 1)
        if length(critical_cars) > (instance.nb_cars * 0.6)
            neighbor = local_search_exchange_ils_hprc!(neighbor, instance, start_time)
        else
            neighbor = fast_local_search_exchange_ils_hprc!(neighbor, instance, critical_cars, start_time)
        end
        if cost_HPRC(solution, instance) <= cost_HPRC(neighbor, instance)
            solution = neighbor
        end
        if i == ALPHA_ILS
            solution = intensification_ils_hprc!(solution, instance, start_time)
        end
        if i == BETA_ILS
            cond = cond + 1
            if cost_HPRC(lastopt, instance) > cost_HPRC(solution_opt, instance)
                lastopt = solution_opt
                cond = 0
            end
            if cost_HPRC(solution, instance) > cost_HPRC(solution_opt, instance)
                solution = solution_opt                   # Restart from solution*
                i = 0
            else
                if nb_strong_perturbation < 3
                    solution = restart_ils_hprc(solution, instance)       # Restart from strong perturbation (50 cars)
                    nb_strong_perturbation = nb_strong_perturbation + 1
                    i = 0
                else
                    solution = greedy(instance)                    # Regenerate an initial solution and restart
                    nb_strong_perturbation = 0
                    i = 0
                end
            end
        end
        if cost_HPRC(solution, instance) < cost_HPRC(solution_opt, instance)            # There is an improvement
            solution_opt = solution
            i = 0                   # So the number of iterations since the last improvement shall return to 0
        else
            i = i + 1
        end
    end
    return solution_opt
end


text:
include("POD_experiment/tspn05.jl")
#include("POD_experiment/ndcc12persp.jl")
include("POD_experiment/FLay02H.jl")
include("basic/gamsworld.jl")

@testset "fp tests" begin

@testset "fp no objective and start value" begin
    println("==================================")
    println("fp No objective and start value")
    println("==================================")
    juniper = DefaultTestSolver(mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper...)
    )

    @variable(m, x, Int, start=3)

    @constraint(m, x >= 0)
    @constraint(m, x <= 5)
    @NLconstraint(m, x^2 >= 17)

    status = solve(m)
    inner = internalmodel(m)
    @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test JuMP.dual_status(m) == MOI.FEASIBLE_POINT
    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.value(x), 5, atol=sol_atol)
    @test Juniper.getnsolutions(inner) == 1
    @test inner.primal_start[1] == 3
end

@testset "FP: no linear" begin
    println("==================================")
    println("FP: no linear")
    println("==================================")

    m = Model()
    @variable(m, x[1:10], Bin)

    special_constr_fct(x,y) = x^2+y^2
    register_args = [:special_constr_fct, 2, special_constr_fct]
    JuMP.register(m, register_args...; autodiff=true)

    # == 1 such that start value of 0 is not optimal (test for issue 189)
    @NLconstraint(m, special_constr_fct(x[1],x[2]) == 1)
    @objective(m, Max, sum(x))

    optimizer = optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:MostInfeasible,
            time_limit = 1,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0),
        )...
    )
    set_optimizer(m, optimizer)

    status = solve(m)

    # should have feasibility_pump be set to true
    @test MOI.get(m, MOI.RawParameter(:feasibility_pump))
    @test Juniper.getnsolutions(internalmodel(m)) >= 1
end

@testset "GLPK Binary bounds #143" begin
    juniper_solver = optimizer_with_attributes(Juniper.Optimizer,
        "log_levels" => [],
        "nl_solver"  => optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0),
        "mip_solver" => optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 0)
    )

    m = Model(juniper_solver)
    @variable(m, z[1:2], Bin)
    @variable(m, x >= 0.5, Bin)
    @constraint(m, sum(z)+x <= 2)
    @objective(m, Max, sum(z)+x)
    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m),2.0,atol=opt_atol)
end

@testset "FP: Integer test" begin
    println("==================================")
    println("FP: Integer Test")
    println("==================================")
    m = Model()

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    # make sure that every solution to lp is solution to nlp
    # FP should definitely find a solution
    for i=1:4
        @constraint(m, x[i] >= i-0.1)
        @constraint(m, x[i] <= i+0.1)
    end
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:MostInfeasible,
            feasibility_pump = true,
            time_limit = 1,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)
    @test Juniper.getnsolutions(internalmodel(m)) >= 1
end

@testset "FP: Integer test2" begin
    println("==================================")
    println("FP: Integer Test 2")
    println("==================================")
    m = Model()

    @variable(m, 0 <= x <= 10, Int)
    @variable(m, y >= 0)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    @objective(m, Min, -3x - y)

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)
    @NLconstraint(m, x^2 >= u*w)

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:MostInfeasible,
            feasibility_pump = true,
            time_limit = 1,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)
    @test Juniper.getnsolutions(internalmodel(m)) >= 1
end

@testset "FP: infeasible cos" begin
    println("==================================")
    println("FP: Infeasible cos")
    println("==================================")
    m = Model()

    @variable(m, 1 <= x <= 5, Int)
    @variable(m, -2 <= y <= 2, Int)

    @objective(m, Min, -x-y)

    @NLconstraint(m, y==2*cos(2*x))

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:MostInfeasible,
            feasibility_pump = true,
            time_limit = 1,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)
    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
    @test Juniper.getnsolutions(internalmodel(m)) == 0
end

@testset "FP: tspn05" begin
    println("==================================")
    println("FP: tspn05")
    println("==================================")

    m = get_tspn05()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            feasibility_pump = true,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m),191.2541,atol=1e0)
end

@testset "FP: FLay02H" begin
    println("==================================")
    println("FP: FLay02H")
    println("==================================")

    # This probably needs a restart in real nlp
    m = get_FLay02H()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            feasibility_pump = true,
            feasibility_pump_time_limit = 10,
            time_limit = 10,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT
    @test Juniper.getnsolutions(internalmodel(m)) >= 1
end

@testset "FP: FLay02H short feasibility_pump_time_limit" begin
    println("==================================")
    println("FP: FLay02H short feasibility_pump_time_limit")
    println("==================================")

    # This probably needs a restart in real nlp
    m = get_FLay02H()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            feasibility_pump = true,
            feasibility_pump_time_limit = 1,
            time_limit = 2,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT
end

@testset "FP: Issue 195: FPump without objective" begin
    println("==================================")
    println("FP: Issue 195: FPump without objective")
    println("==================================")

    ipopt_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "sb" => "yes", "max_iter" => 50000)
    cbc_solver = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    juniper_solver = JuMP.optimizer_with_attributes(Juniper.Optimizer,
        "nl_solver" => ipopt_solver,
        "mip_solver" => cbc_solver, "log_levels" => [])
    m = Model(juniper_solver)
    @variable(m, x[1:3], Int)
    @variable(m, y)
    @NLconstraint(m, x[1]*x[2]*x[3]*y >= 5)
    optimize!(m)

    for i=1:3
        xval = JuMP.value(x[i])
        @test isapprox(round(xval)-xval, 0; atol=sol_atol)
    end
    @test JuMP.objective_value(m) ≈ 0.0
end

end


text:
# ---
# title: 525. Contiguous Array
# id: problem525
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/contiguous-array/description/>
# hidden: true
# ---
#
# Given a binary array, find the maximum length of a contiguous subarray with
# equal number of 0 and 1.
#
# **Example 1:**
#
#
#
#     Input: [0,1]
#     Output: 2
#     Explanation: [0, 1] is the longest contiguous subarray with equal number of 0 and 1.
#
#
# **Example 2:**
#
#
#
#     Input: [0,1,0]
#     Output: 2
#     Explanation: [0, 1] (or [1, 0]) is a longest contiguous subarray with equal number of 0 and 1.
#
#
# **Note:** The length of the given binary array will not exceed 50,000.
#
#
## @lc code=start
using LeetCode

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


text:
module BARON

using MathProgBase
importall MathProgBase.SolverInterface

using Compat

export BaronSolver
immutable BaronSolver <: AbstractMathProgSolver
    options
end
BaronSolver(;kwargs...) = BaronSolver(kwargs)

const baron_exec = ENV["BARON_EXEC"]

type BaronMathProgModel <: AbstractNonlinearModel
    options

    xˡ::Vector{Float64}
    xᵘ::Vector{Float64}
    gˡ::Vector{Float64}
    gᵘ::Vector{Float64}

    nvar::Int
    ncon::Int

    obj
    constrs

    vartypes::Vector{Symbol}

    v_names::Vector{Compat.UTF8String}
    c_names::Vector{Compat.UTF8String}

    sense::Symbol

    x₀::Vector{Float64}

    probfile::Compat.UTF8String
    sumfile::Compat.UTF8String
    resfile::Compat.UTF8String

    objval::Float64
    solution::Vector{Float64}
    status::Symbol

    d::AbstractNLPEvaluator

    function BaronMathProgModel(;options...)
        dir = tempdir()
	    push!(options, (:ResName, joinpath(dir, "res.lst")))
	    push!(options, (:TimName, joinpath(dir, "tim.lst")))
	    push!(options, (:SumName, joinpath(dir, "sum.lst")))
        new(options,
	    zeros(0),
	    zeros(0),
	    zeros(0),
        zeros(0),
	    0,
	    0,
	    :(0),
	    Expr[],
	    Symbol[],
	    Compat.UTF8String[],
	    Compat.UTF8String[],
	    :Min,
	    zeros(0),
        "",
	    "",
	    "",
	    NaN,
	    zeros(0),
        :NotSolved)
    end
end

MathProgBase.NonlinearModel(s::BaronSolver) = BaronMathProgModel(;s.options...)
MathProgBase.LinearQuadraticModel(s::BaronSolver) = MathProgBase.NonlinearToLPQPBridge(MathProgBase.NonlinearModel(s))

verify_support(c) = c

function verify_support(c::Expr)
    if c.head == :comparison
        map(verify_support, c.args)
        return c
    end
    if c.head == :call
        if c.args[1] in (:+, :-, :*, :/, :exp, :log)
            return c
        elseif c.args[1] in (:<=, :>=, :(==))
            map(verify_support, c.args[2:end])
            return c
        elseif c.args[1] == :^
            @assert isa(c.args[2], Real) || isa(c.args[3], Real)
            return c
        else # TODO: do automatic transformation for x^y, |x|
            error("Unsupported expression $c")
        end
    end
    return c
end

function MathProgBase.loadproblem!(m::BaronMathProgModel,
    nvar, ncon, xˡ, xᵘ, gˡ, gᵘ, sense,
    d::MathProgBase.AbstractNLPEvaluator)

    @assert nvar == length(xˡ) == length(xᵘ)
    @assert ncon == length(gˡ) == length(gᵘ)

    m.xˡ, m.xᵘ = xˡ, xᵘ
    m.gˡ, m.gᵘ = gˡ, gᵘ
    m.sense = sense
    m.nvar, m.ncon = nvar, ncon
    m.d = d

    m.v_names = ["x$i" for i in 1:nvar]
    m.c_names = ["e$i" for i in 1:ncon]

    MathProgBase.initialize(d, [:ExprGraph])

    m.obj = verify_support(MathProgBase.obj_expr(d))
    m.vartypes = fill(:Cont, nvar)

    m.constrs = map(1:ncon) do c
        verify_support(MathProgBase.constr_expr(d,c))
    end

    dir = tempdir()
    m.probfile = joinpath(dir, "baron_problem.bar")
    m.sumfile  = joinpath(dir, "sum.lst")
    m.resfile  = joinpath(dir, "res.lst")
    m
end

function MathProgBase.setvartype!(m::BaronMathProgModel, cat::Vector{Symbol})
    @assert all(x-> (x in (:Cont,:Bin,:Int)), cat)
    m.vartypes = cat
end

function print_var_definitions(m, fp, header, condition)
    idx = filter(condition, 1:m.nvar)
    if !isempty(idx)
        println(fp, header, join([m.v_names[i] for i in idx], ", "), ";")
    end
end

# to_str(x::Int) = string(x)
# to_str(x) = (@show(x); string(float(x)))
to_str(x) = string(x)

function to_str(c::Expr)
    if c.head == :comparison
        if length(c.args) == 3
            return join([to_str(c.args[1]), c.args[2], c.args[3]], " ")
        elseif length(c.args) == 5
            return join([c.args[1], c.args[2], to_str(c.args[3]),
                         c.args[4], c.args[5]], " ")
        end
    elseif c.head == :call
        if c.args[1] in (:<=,:>=,:(==))
            if length(c.args) == 3
                return join([to_str(c.args[2]), to_str(c.args[1]), to_str(c.args[3])], " ")
            elseif length(c.args) == 5
                return join([to_str(c.args[1]), to_str(c.args[2]), to_str(c.args[3]), to_str(c.args[4]), to_str(c.args[5])], " ")
            end
        elseif c.args[1] in (:+,:-,:*,:/,:^)
            if all(d->isa(d, Real), c.args[2:end]) # handle unary case
                return string(eval(c))
            elseif c.args[1] == :- && length(c.args) == 2
                return string("(-$(to_str(c.args[2])))")
            else
                return string("(", join([to_str(d) for d in c.args[2:end]], string(c.args[1])), ")")
            end
        elseif c.args[1] in (exp,:log)
            if isa(c.args[2], Real)
                return string(eval(c))
            else
                return string(c.args[1], "( ", to_str(c.args[2]), " )")
            end
        end
    elseif c.head == :ref
        if c.args[1] == :x
            @assert isa(c.args[2], Int)
            return "x$(c.args[2])"
        else
            error("Unrecognized reference expression $c")
        end
    end
end

function write_bar_file(m::BaronMathProgModel)
    fp = open(m.probfile, "w")

    # First: process any options
    println(fp, "OPTIONS{")
    for (opt,setting) in m.options
        if isa(setting, AbstractString) # wrap it in quotes
            println(fp, unescape_string("$opt: $('"')$setting$('"');"))
        else
            println(fp, "$opt: $setting;")
        end
    end
    println(fp, "}")
    println(fp)

    # Next, define variables
    print_var_definitions(m, fp, "BINARY_VARIABLES ",   v->(m.vartypes[v]==:Bin))
    print_var_definitions(m, fp, "INTEGER_VARIABLES ",  v->(m.vartypes[v]==:Int))
    print_var_definitions(m, fp, "POSITIVE_VARIABLES ", v->(m.vartypes[v]==:Cont && m.xˡ[v]==0))
    print_var_definitions(m, fp, "VARIABLE ",           v->(m.vartypes[v]==:Cont && m.xˡ[v]!=0))
    println(fp)

    # Print variable bounds
    if any(c->!isinf(c), m.xˡ)
        println(fp, "LOWER_BOUNDS{")
        for (i,l) in enumerate(m.xˡ)
            if !isinf(l)
                println(fp, "$(m.v_names[i]): $l;")
            end
        end
        println(fp, "}")
        println(fp)
    end
    if any(c->!isinf(c), m.xᵘ)
        println(fp, "UPPER_BOUNDS{")
        for (i,u) in enumerate(m.xᵘ)
            if !isinf(u)
                println(fp, "$(m.v_names[i]): $u;")
            end
        end
        println(fp, "}")
        println(fp)
    end

    # Now let's declare the equations
    if !isempty(m.constrs)
        println(fp, "EQUATIONS ", join(m.c_names, ", "), ";")
        for (i,c) in enumerate(m.constrs)
            str = to_str(c)
            println(fp, "$(m.c_names[i]): $str;")
        end
        println(fp)
    end

    # Now let's do the objective
    print(fp, "OBJ: ")
    print(fp, m.sense == :Min ? "minimize " : "maximize ")
    print(fp, to_str(m.obj))
    println(fp, ";")
    println(fp)

    if !isempty(m.x₀)
        println(fp, "STARTING_POINT{")
        for (i,v) in enumerate(m.x₀)
            println(fp, "$(m.v_names[i]): $v;")
        end
        println(fp, "}")
    end
    close(fp)
end

const user_limits = (
    "Max. allowable nodes in memory reached",
    "Max. allowable BaR iterations reached",
    "Max. allowable CPU time exceeded",
    "Problem is numerically sensitive",
    "Insufficient Memory for Data structures"
)

function read_results(m::BaronMathProgModel)
    # First, we read the summary file to get the solution status
    fp = open(m.sumfile, "r")
    stat = :Undefined
    while true
        line = readline(fp)
        spl = split(chomp(line))
        if !isempty(spl) && spl[1] == "***"
            if spl[2] in user_limits
                stat = :UserLimit
            end
            break
        end
        eof(fp) && error("Reached EOF while searching for termination notice")
    end
    while true
        line = readline(fp)
        spl = split(chomp(line))
        if !isempty(spl) && spl[1:3] == ["Best","solution","found"]
            node = parse(Int,match(r"\d+", line).match)
            if node == -3
                stat = :Infeasible
            end
            break
        elseif spl == ["***", "Problem", "is", "infeasible", "***"]
            stat = :Infeasible
            break
        elseif spl == ["***", "Problem", "is", "unbounded", "***"]
            stat = :Unbounded
            break
        else
            stat = :Optimal
            break
        end
        eof(fp) && error("Reached OEF while looking for node with best solution")
    end
    m.status = stat
    close(fp)

    # Next, we read the results file to get the solution
    x = fill(NaN, m.nvar)
    m.objval = NaN
    if stat == :Optimal
        fp = open(m.resfile, "r")
        while true
            startswith(readline(fp), "The best solution found") && break
            eof(fp) && error()
        end
        readline(fp)
        readline(fp)

        while true
            line = chomp(readline(fp))
            parts = split(line)
            isempty(parts) && break
            mt = match(r"\d+", parts[1])
            mt == nothing && error("Cannot find appropriate variable index from $(parts[1])")
            v_idx = parse(Int,mt.match)
            v_val = float(parts[3])
            x[v_idx] = v_val
        end
        m.solution = x
        line = readline(fp)
        val = match(r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?", chomp(line))
        if val != nothing
            m.objval = float(val.match)
        end
    end
    nothing
end

MathProgBase.setwarmstart!(m::BaronMathProgModel, v::Vector{Float64}) = (m.x₀ = v)

function MathProgBase.optimize!(m::BaronMathProgModel)
    write_bar_file(m)
    run(`$baron_exec $(m.probfile)`)

    read_results(m)
end

MathProgBase.status(m::BaronMathProgModel) = m.status

MathProgBase.numvar(m::BaronMathProgModel) = m.nvar
MathProgBase.getsolution(m::BaronMathProgModel) = m.solution
MathProgBase.getobjval(m::BaronMathProgModel) = m.objval

end


text:
export Problem

"""
    Problem( model, initial, param ; solver, label)

Build an initial-value problem which can then be solved (i.e. integrated in time) through `solve!( problem )`

# Arguments
- `model   :: AbstractModel`,  the system of equation solved.
May be built, e.g., by `WaterWaves(param)`;
- `initial :: InitialData`, the initial data.
May be buit, e.g., by `Init(η,v)` where
`η` is the surface deformation and `v` the derivative of the trace of the velocity potential at the surface;
- `times :: Times` is the time grid, and may be built using the function `Times`.
Alternatively one can simply provide a `NamedTuple` with
    - `T`, the final time of integration
    - `dt`, the timestep
    - optionally, `Ns` the number of computed data or `ns` for storing data every `ns` computation steps (by default, every computed data is stored).

## Optional keyword arguments
- `solver :: TimeSolver`, the solver for time integration (default is explicit Runge-Kutta fourth order solver).
May be built, e.g., by `RK4(model)` or `RK4_naive()`.
- `label   :: String ` is used in future references (e.g. `plot_solution`).
- Information are not printed if `verbose = false` (default is `true`).

"""
mutable struct Problem

    model   :: AbstractModel
    initial :: InitialData
    solver  :: TimeSolver
    times   :: Times
    data    :: Data
    label   :: String

    function Problem(model   :: AbstractModel,
        initial :: InitialData,
        times   :: Times;
        solver = RK4(model)  :: TimeSolver,
        label = nothing
        )

        if isnothing(label)   label = model.label end

        U = model.mapto(initial)
        data  = Data(U)

        #try
        #    mesh  = Mesh(param)
        #catch
        #    η,v,x = model.mapro(U)
        #    mesh = Mesh(x)
        #end

        # A basic check
        try
            step!(solver, model, copy(last(data.U)), 1)
        catch
            @warn "The model and the solver are incompatible. `solve!` will not work."
        end

        new(model, initial, solver, times, data, label)

    end

    function Problem(model   :: AbstractModel,
        initial :: InitialData,
        param   :: NamedTuple;
        solver = RK4(model)  :: TimeSolver,
        label = nothing
        )

        if in(:Ns,keys(param))
        times = Times(param; Ns = param.Ns)
        elseif in(:ns,keys(param))
        times = Times(param; ns = param.ns)
        else
        times = Times(param)
        end

        p = Problem(model, initial, times; solver = solver, label = label)

        new(p.model, p.initial, p.solver, p.times, p.data, p.label)

    end

end

show(io::IO, p::Problem) =
    print(io,"Initial-value problem with...\n\
    ├─Model: $(p.model.label).\n├─Initial data: $(p.initial.label).\n├─Solver: $(p.solver.label).\n\
    └─Grid of times: $(p.times.Nc) computed times on [0, $(p.times.tfin)] (timestep dt=$(p.times.dt)), \
    among which $(p.times.Ns) will be stored.")

#     Discretized with $(p.mesh.N) collocation points on [$(p.mesh.xmin), $(p.mesh.xmax)],\n\

export solve!

"""
    solve!( problem :: Problem; verbose=true )

Solve (i.e. integrate in time) an initial-value problem

The argument `problem` should be of type `Problem`.
It may be buit, e.g., by `Problem(model, initial, param)`

Information are not printed if keyword argument `verbose = false` (default is `true`).

"""
function solve!(problem :: Problem; verbose=true::Bool)

    ci = get(ENV, "CI", nothing) == "true"

    U = copy(last(problem.data.U))

    dt     = problem.times.dt
    solver = problem.solver
    model  = problem.model
    data   = problem.data.U

    if verbose == true
        @info "Now solving the initial-value problem $(problem.label)\n\
            with timestep dt=$(problem.times.dt), final time T=$(problem.times.tfin),\n\
            and N=$(size(U)[1]) collocation points."
    end

    if problem.times.tc == problem.times.ts
        pbar = Progress(problem.times.Ns-1; enabled = !ci)
        for j in 1:problem.times.Ns-1
            step!(solver, model, U, dt)
            push!(data,copy(U))
            next!(pbar)
        end

    elseif length(problem.times.ts) > 25
        pbar = Progress(problem.times.Ns-1; enabled = !ci)
        for j in 1:problem.times.Ns-1
            for l in 1:problem.times.ns[j]
                step!(solver, model, U, dt)
            end
            push!(data,copy(U))
            next!(pbar)
        end
    else
        for j in 1:problem.times.Ns-1
            pbar = Progress(problem.times.ns[j], desc=string("Step ",j,"/",problem.times.Ns-1,"...") ; enabled = !ci)
            for l in 1:problem.times.ns[j]
                step!(solver, model, U, dt)
                next!(pbar)
            end
            push!(data,copy(U))
            println()

        end
    end

    println()

end

using Base.Threads

"""
    solve!( problems; verbose=true )

Solve (i.e. integrate in time) a collection of initial-value problems.

The argument `problems` should be a collection (list, array...) of elements of type `Problem`.

Information are not printed if keyword argument `verbose = false` (default is `true`).
"""
function solve!(problems; verbose=true::Bool)

    ci = get(ENV, "CI", nothing) == "true"

    # Set up
    U=[];nsteps=0;flag = false;
    for problem in problems
        for p in problems # Check whether multi-threading is at risk, using same allocations
            if !(p===problem) && ( p.model===problem.model || p.solver===problem.solver )
                flag = true
            end
        end
        push!(U,copy(last(problem.data.U)))
        nsteps+=problem.times.Ns-1
    end
    if flag
        @warn "`solve!(problems)` uses multi-threading. \
        It is ill-advised to use it with problems sharing the same model or the same solver. \
        Re-define model/solver when building your problems, or solve them independently \
        via `for problem in problems solve!(problem) end`."
    end

    pg=Progress(nsteps;dt=1,enabled = !ci)
    @threads for i in 1:length(problems)
        if verbose == true
            @info "Now solving the initial-value problem $(problems[i].label)\n\
                with timestep dt=$(problems[i].times.dt), final time T=$(problems[i].times.tfin),\n\
                and N=$(size(problems[i].data.U)[1]) collocation points."
        end

        if problems[i].times.tc == problems[i].times.ts

            for j in 1:problems[i].times.Ns-1
                step!(problems[i].solver, problems[i].model, U[i], problems[i].times.dt)
                push!(problems[i].data.U,copy(U[i]))
                next!(pg)
            end

        else
            for j in 1:problems[i].times.Ns-1
                for l in 1:problems[i].times.ns[j]
                    step!(problems[i].solver, problems[i].model, U[i], problems[i].times.dt)
                end
                push!(problems[i].data.U,copy(U[i]))
                next!(pg)
            end
        end
        if verbose == true
            @info "Done solving the problem $(problems[i].label)."
        end

    end

    println()

end

Base.:(≈)(x::NamedTuple{N,T}, y::NamedTuple{N2,T2}) where {N,T,N2,T2} =
    length(N) === length(union(N,N2)) && all(k->getfield(x,k) == getfield(y,k), keys(x))

Base.:(==)(p1::Problem, p2::Problem) =
    p1.model == p2.model &&
    p1.initial == p2.initial &&
    p1.solver == p2.solver &&
    p1.times == p2.times &&
    p1.data == p2.data &&
    p1.label == p2.label