text:
#-----------------------------------------------------------------------# Bin2D
mutable struct Bin2D
    x::Float64
    y::Float64
    z::Int
end
distance(a::Bin2D, b::Bin2D) = sqrt((a.x - b.x)^2 + (a.y - b.y)^2)
function _merge!(a::Bin2D, b::Bin2D)
    γ = b.z / (a.z += b.z)
    a.x = smooth(a.x, b.x, γ)
    a.y = smooth(a.y, b.y, γ)
end
nobs(o::Bin2D) = o.z

#-----------------------------------------------------------------------# Hist2D
"""
    KHist2D(b=300)

Approximate scatterplot of `b` centers.  This implementation is slow!  See [`IndexedPartition`](@ref)
and [`KIndexedPartition`](@ref) instead.

# Example

    x = randn(10^4)
    y = x + randn(10^4)
    plot(fit!(KHist2D(), zip(x, y)))
"""
struct KHist2D <: OnlineStat{TwoThings{<:Number, <:Number}}
    value::Vector{Bin2D}
    b::Int
end
KHist2D(b::Int=300) = KHist2D(Bin2D[], b)

Base.@deprecate_binding PlotNN KHist2D

nobs(o::KHist2D) = length(o.value) > 0 ? sum(nobs, o.value) : 0

# Works, but much room to optimize
function _fit!(o::KHist2D, xy)
    newbin = Bin2D(xy..., 1)
    push!(o.value, newbin)
    if length(o.value) > o.b
        i, j = 1, 2
        mindist = Inf
        for _i in 1:o.b
            v1 = o.value[_i]
            for _j in (_i + 1):o.b
                dist = distance(v1, o.value[_j])
                if dist < mindist
                    mindist = dist
                    i = _i
                    j = _j
                end
            end
        end
        _merge!(o.value[i], o.value[j])
        deleteat!(o.value, j)
    end
end

@recipe function f(o::KHist2D; maxsize=10, minsize=1)
    x = [v.x for v in o.value]
    y = [v.y for v in o.value]
    z = [v.z for v in o.value]
    marker_z --> z
    markersize --> z / maximum(nobs, o.value) .* maxsize .+ minsize
    markerstrokewidth --> 0
    seriestype --> :scatter
    seriescolor --> :viridis
    x, y
end


text:
include("utils.jl")

@testset "Grassmann" begin
    @testset "Real" begin
        M = Grassmann(3,2)
        @testset "Basics" begin
            @test repr(M) == "Grassmann(3, 2, ℝ)"
            @test representation_size(M) == (3,2)
            @test manifold_dimension(M) == 2
            @test !is_manifold_point(M,[1., 0., 0., 0.])
            @test !is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], [0., 0., 1., 0.])
            @test_throws DomainError is_manifold_point(M, [2. 0.; 0. 1.; 0. 0.],true)
            @test_throws DomainError is_tangent_vector(M, [2. 0.; 0. 1.; 0. 0.],zeros(3,2),true)
            @test_throws DomainError is_tangent_vector(M, [1. 0.; 0. 1.; 0. 0.],ones(3,2),true)
            @test is_manifold_point(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], true)
            @test_throws DomainError is_manifold_point(M, 1im*[1.0 0.0; 0.0 1.0; 0.0 0.0], true)
            @test is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], zero_tangent_vector(M,[1.0 0.0; 0.0 1.0; 0.0 0.0]), true)
            @test_throws DomainError is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], 1im*zero_tangent_vector(M,[1.0 0.0; 0.0 1.0; 0.0 0.0]), true)
            @test injectivity_radius(M) == π/2
            @test injectivity_radius(M, ExponentialRetraction()) == π/2
            @test injectivity_radius(M, [1.0 0.0; 0.0 1.0; 0.0 0.0]) == π/2
            @test injectivity_radius(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], ExponentialRetraction()) == π/2
        end
        types = [Matrix{Float64}, ]
        TEST_STATIC_SIZED && push!(types, MMatrix{3, 2, Float64})

        TEST_FLOAT32 && push!(types, Matrix{Float32})
        basis_types = (ProjectedOrthonormalBasis(:gram_schmidt),)
        @testset "Type $T" for T in types
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            v = [0.0 0.0; 0.0 0.0; 0.0 1.0]
            y = exp(M,x,v)
            w = [0.0 1.0; -1.0 0.0; 1.0 0.0]
            z = exp(M,x,w)
            pts = convert.(T, [x,y,z])
            test_manifold(
                M,
                pts,
                test_exp_log = true,
                test_injectivity_radius = false,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vee_hat = false,
                retraction_methods = [PolarRetraction(), QRRetraction()],
                inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()],
                #basis_types_vecs = basis_types,
                exp_log_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 10.0,
            )

            @testset "inner/norm" begin
                v1 = inverse_retract(M, pts[1], pts[2], PolarInverseRetraction())
                v2 = inverse_retract(M, pts[1], pts[3], PolarInverseRetraction())

                @test real(inner(M, pts[1], v1, v2)) ≈ real(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v2)) ≈ -imag(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v1)) ≈ 0

                @test norm(M, pts[1], v1) isa Real
                @test norm(M, pts[1], v1) ≈ sqrt(inner(M, pts[1], v1, v1))
            end
        end
    end

    @testset "Complex" begin
        M = Grassmann(3,2,ℂ)
        @testset "Basics" begin
            @test repr(M) == "Grassmann(3, 2, ℂ)"
            @test representation_size(M) == (3,2)
            @test manifold_dimension(M) == 4
            @test !is_manifold_point(M,[1., 0., 0., 0.])
            @test !is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], [0., 0., 1., 0.])
            @test Manifolds.allocation_promotion_function(M,exp!,(1,)) == complex
            @test_throws DomainError is_manifold_point(M, [2. 0.; 0. 1.; 0. 0.],true)
            @test_throws DomainError is_tangent_vector(M, [2. 0.; 0. 1.; 0. 0.],zeros(3,2),true)
            @test_throws DomainError is_tangent_vector(M, [1. 0.; 0. 1.; 0. 0.],ones(3,2),true)
            @test is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], 1im*zero_tangent_vector(M,[1.0 0.0; 0.0 1.0; 0.0 0.0]))
            @test is_manifold_point(M, [1.0 0.0; 0.0 1.0; 0.0 0.0])
            @test injectivity_radius(M) == π/2
        end
        types = [
            Matrix{ComplexF64},
        ]
        @testset "Type $T" for T in types
            x = [0.5+0.5im 0.5+0.5im; 0.5+0.5im -0.5-0.5im; 0.0 0.0]
            v = [0.0 0.0; 0.0 0.0; 0.0 1.0]
            y = exp(M,x,v)
            w = [0.0 1.0; -1.0 0.0; 1.0 0.0]
            z = exp(M,x,w)
            pts = convert.(T, [x,y,z])
            test_manifold(
                M,
                pts,
                test_exp_log = true,
                test_injectivity_radius = false,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vee_hat = false,
                retraction_methods = [PolarRetraction(), QRRetraction()],
                inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()],
                exp_log_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 10.0,
            )

            @testset "inner/norm" begin
                v1 = inverse_retract(M, pts[1], pts[2], PolarInverseRetraction())
                v2 = inverse_retract(M, pts[1], pts[3], PolarInverseRetraction())

                @test real(inner(M, pts[1], v1, v2)) ≈ real(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v2)) ≈ -imag(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v1)) ≈ 0

                @test norm(M, pts[1], v1) isa Real
                @test norm(M, pts[1], v1) ≈ sqrt(inner(M, pts[1], v1, v1))
            end
        end
    end
end

ext: jl


text: minWinWidth = Int32(20)  # basically 0.


text:
# Module written by Ed Scheinerman, ers@jhu.edu
# distributed under terms of the MIT license

module SimpleGraphs
using DataStructures
using SimpleRandom
using Primes # needed for Paley graphs
using Polynomials
using LinearAlgebra
using Statistics
using IterTools
using Random
using AbstractLattices

import AbstractLattices.dist

"""
`AbstractSimpleGraph` is a parent class for `SimpleGraph` and `SimpleDigraph`.
"""
abstract type AbstractSimpleGraph end
export AbstractSimpleGraph

include("simple_core.jl")
include("d_simple_core.jl")

include("cache.jl")
include("simple_ops.jl")
include("simple_constructors.jl")
include("platonic.jl")
include("simple_connect.jl")
include("simple_matrices.jl")
# include("disjoint_sets_helper.jl")
# include("simple_converters.jl")
include("simple_coloring.jl")
include("simple_euler.jl")
include("simple_ham.jl")
include("simple_girth.jl")
include("bisect.jl")
include("trans_orient.jl")
include("interlace.jl")
include("matching_poly.jl")
include("indep_poly.jl")
include("twins.jl")

include("d_simple_constructors.jl")
include("d_simple_matrices.jl")
include("d_dist.jl")
include("d_euler.jl")
include("d_ham.jl")

include("embedding/GraphEmbeddings.jl")

end # module SimpleGraphs


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

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

# ╔═╡ f1295720-40c7-11eb-12be-8bd0509deca0
begin
	using Markdown
	using InteractiveUtils
	using LinearAlgebra
end

# ╔═╡ 0a3822a0-40c8-11eb-342c-a77110d06445
begin
    using PlutoUI
    using Plots
    using Printf
    using Interact
end

# ╔═╡ 142d35f0-4192-11eb-244d-8db5ddc08a84
md"_Compressible Similarity (Interactive), version 1_"

# ╔═╡ 024ab760-40c8-11eb-3343-2f5d5cd7535e
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**

**Computational Aerodynamics and Hypersonics Laboratory**

**Oklahoma State University**

**foz@okstate.edu**
"""

# ╔═╡ 295debb0-40c8-11eb-1d32-a37eb331b81c
function Y1(η,y₁,y₂,y₃,y₄,y₅)
    return y₂
end

# ╔═╡ 32e77ac0-40c8-11eb-3f9a-313f3cb78e79
function Y2(η,y₁,y₂,y₃,y₄,y₅)
    return y₃
end

# ╔═╡ 373ec6f0-40c8-11eb-06c9-df8d801f6678
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

# ╔═╡ 3c949760-40c8-11eb-1ccd-75955be21825
function Y4(η,y₁,y₂,y₃,y₄,y₅)
    return y₅
end

# ╔═╡ 41c04a90-40c8-11eb-23a8-13a2fef3fb74
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

# ╔═╡ 132d8800-40c8-11eb-04f8-350ebed9c8ee
"""
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,

Computational Aerodynamics and Hypersonics Laboratory,

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

# ╔═╡ 730c0e5e-4193-11eb-25d0-8394adde63bb
md"""
-------------------------------------------------------------------------
"""

# ╔═╡ 5a236880-4193-11eb-3d79-d5668e165e21
md"""
Please use the sliders to observe the change with the Mach number and temperature.
"""

# ╔═╡ 5fb1b250-40c8-11eb-10e3-1bff8856c8ae
md"""
M∞:
$(@bind Mtest Slider(0.1:0.1:10; default=1, show_value=true))

"""

# ╔═╡ b3831040-40c8-11eb-20ad-f1c8c2f59e96
md"""
T∞:
$(@bind Ttest Slider(80:20:500; default=300, show_value=true))

"""

# ╔═╡ a838d380-40ca-11eb-22c3-13a7d164bf17
ηᵢ,yᵢ,Uᵢ,Tᵢ,Nᵢ = selfsimilar(Mtest,Ttest);

# ╔═╡ f0cf33be-40c8-11eb-341d-13c3c8b2b99b
begin
	fp = plot(Uᵢ,yᵢ,

		label = "U",
		legend = :topleft,
		ylims = (0,40),
		xlabel = "U",
		ylabel = "y/√(νx/U)",
		linewidth = 2,
		linecolor = :black,
		markershape = :circle,
		markercolor = :red,
		)
	sp = plot(Tᵢ,yᵢ,

		label = "T",
		legend = :topright,
		xlims = (0,18),
		ylims = (0,40),
		xlabel = "T",
		ylabel = "y/√(νx/U)",
		linewidth = 2,
		linecolor = :black,
		markershape = :circle,
		markercolor = :red,
		)
	plot(fp,sp,title = "Similarity (M=$Mtest - T=$Ttest)",layout=2)
end

# ╔═╡ Cell order:
# ╟─142d35f0-4192-11eb-244d-8db5ddc08a84
# ╟─024ab760-40c8-11eb-3343-2f5d5cd7535e
# ╟─f1295720-40c7-11eb-12be-8bd0509deca0
# ╟─0a3822a0-40c8-11eb-342c-a77110d06445
# ╟─132d8800-40c8-11eb-04f8-350ebed9c8ee
# ╟─295debb0-40c8-11eb-1d32-a37eb331b81c
# ╟─32e77ac0-40c8-11eb-3f9a-313f3cb78e79
# ╟─373ec6f0-40c8-11eb-06c9-df8d801f6678
# ╟─3c949760-40c8-11eb-1ccd-75955be21825
# ╟─41c04a90-40c8-11eb-23a8-13a2fef3fb74
# ╟─730c0e5e-4193-11eb-25d0-8394adde63bb
# ╟─5a236880-4193-11eb-3d79-d5668e165e21
# ╟─5fb1b250-40c8-11eb-10e3-1bff8856c8ae
# ╟─b3831040-40c8-11eb-20ad-f1c8c2f59e96
# ╟─a838d380-40ca-11eb-22c3-13a7d164bf17
# ╟─f0cf33be-40c8-11eb-341d-13c3c8b2b99b


text:
import Base.ht_keyindex
function msumby{T,S}(by::Vector{T}, val::Vector{S})::Array{Dict{T, S}}
  res = Array{Dict{T, S}}(Threads.nthreads())
  for j in 1:Threads.nthreads()
    res[j] = Dict{T, S}()
  end
  Threads.@threads for i = 1:length(by)
    # ti = Threads.threadid()
    # @inbounds id = by[i]
    # @inbounds val = val[i]
    # @inbounds dict =  res[ti]
    # index = ht_keyindex(dict, id)
    # if index > 0
    #   @inbounds vw = dict.vals[index]
    #   new_vw = vw + val
    #   @inbounds dict.vals[index] = new_vw
    # else
    #   @inbounds dict[id] = val
    # end
  end
  return res
end

@time res = msumby(id4, v1);

function msumby_test{T,S}(by::Vector{T}, val::Vector{S})
  res = Array{Dict{T, S}}(Threads.nthreads())
  for j in 1:Threads.nthreads()
    res[j] = Dict{T, S}()
  end
  for i = 1:length(by)
    # ti = Threads.threadid()
    # @inbounds id = by[i]
    # @inbounds val = val[i]
    # @inbounds dict =  res[ti]
    # index = ht_keyindex(dict, id)
    # if index > 0
    #   @inbounds vw = dict.vals[index]
    #   new_vw = vw + val
    #   @inbounds dict.vals[index] = new_vw
    # else
    #   @inbounds dict[id] = val
    # end
  end
  return res
end

@time msumby_test(id4, v1)

function nothread_test(v)
  for i = 1:length(v)
    @inbounds v[i] = 0.5
  end
  sum(v)
end

function thread_test(v)
  Threads.@threads for i = 1:length(v)
      @inbounds v[i] = 0.5
  end
  sum(v)
end

v = zeros(Float64, 250_000_000)
@benchmark nothread_test(v)
v = zeros(Float64, 250_000_000)
@benchmark thread_test(v)


text:
using BinDeps

@BinDeps.setup

version = "1.0.2"
url = "https://github.com/kbarbary/sep/archive/v$(version).tar.gz"

libsep = library_dependency("libsep")
println(libsep)
downloadsdir = BinDeps.downloadsdir(libsep)
libdir = BinDeps.libdir(libsep)
srcdir = BinDeps.srcdir(libsep)
if is_apple()
    libfilename = "libsep.dylib"
elseif is_unix()
    libfilename = "libsep.so"
end

# Unix
prefix = joinpath(BinDeps.depsdir(libsep), "usr")
provides(Sources, URI(url), libsep, unpacked_dir="sep-$(version)")
provides(BuildProcess,
         (@build_steps begin
             GetSources(libsep)
             @build_steps begin
                 ChangeDirectory(joinpath(srcdir, "sep-$(version)"))
                 FileRule(joinpath(libdir, libfilename),
                          @build_steps begin
                              `make`
                              `make PREFIX=$(prefix) install`
                          end)
             end
          end), libsep, os = :Unix)

@BinDeps.install Dict(:libsep => :libsep)


text:
import DataFrames

function countdict(x)
    typ = eltype(x)
    cnts = Dict{typ, Int}()
    n = length(x)
    for i = 1:n
        cnts[x[i]] = get(cnts, x[i], 0) + 1
    end
    cnts
end

X = rand(200, 10)
y = vcat(zeros(180), ones(20))

X2, y2 = ClassImbalance.smote(X, y, k = 5, pct_over = 200, pct_under = 200)

cnts = countdict(y2)

Test.@test cnts[0.0] == 80
Test.@test cnts[1.0] == 60

df = DataFrames.DataFrame(hcat(X, y))
X3, y3 = ClassImbalance.smote(df[:, 1:10], df[:, 11], k = 5, pct_over = 300, pct_under = 300)
cnts = countdict(y3)

Test.@test cnts[0.0] == 180
Test.@test cnts[1.0] == 80


text:
struct FileFormat{sym} end

struct File{Format <: FileFormat}
  io :: IO
end

macro format_str(format)
  FileFormat{Symbol(format)}
end

filespec = Tuple{Type{<:FileFormat},Function,String}[]

add_format(format, test::Function, ext) = push!(filespec, (format, test, ext))
function add_format(format, magicbytes::Vector{UInt8}, ext)
  len = length(magicbytes)
  test = io -> read(io, len) == magicbytes
  push!(filespec, (format, test, ext))
end
function add_format(format, magictext::String, ext)
  lines = split(magictext, r"\r?\n")
  test = io -> all(readline(io) == line for line in lines)
  push!(filespec, (format, test, ext))
end

load(path; kwargs...) = open(path, "r") do io
  ext = splitext(path)[2]
  for (format, test, targetext) in filespec
    targetext == ext || continue
    seekstart(io)
    test(io) || continue
    seekstart(io)
    return load(File{format}(io); kwargs...)
  end
end


text:
using LSHFunctions

#========================
Global variables
========================#

const RANDOM_SEED = 0

#========================
Helper types
========================#

# ShiftedSine so that we can quickly generate functions of the form
# f(x) = sin(αx+δ) without having to use closures, which can be relatively more
# time-intensive.
struct ShiftedSine{S <: Real, T <: Real}
    α :: T
    δ :: T
end

ShiftedSine(α::S, δ::T) where {S,T} = ShiftedSine{S,T}(α,δ)

(f::ShiftedSine)(x) = sin(f.α * x + f.δ)

#========================
Helper functions
========================#

mean(x) = sum(x) / length(x)

# Function that draws two random vectors and computes their hashes with
# hashfn. Returns true if the number of collisions of individual hashes
# is within δ of the expected single-hash collision probability, and
# false otherwise.
function test_collision_probability(
        hashfn :: LSHFunctions.LSHFunction,
        δ :: AbstractFloat,
        sampler = () -> randn(4))

    similarity_fn = similarity(hashfn)
    x, y = sampler(), sampler()
    sim = similarity_fn(x,y)

    hx, hy = hashfn(x), hashfn(y)
    prob = collision_probability(hashfn, sim; n_hashes=1)
    coll_freq = mean(hx .== hy)

    prob - δ ≤ coll_freq ≤ prob + δ
end

"""
    create_step_function(N::Integer; step_generator = i -> randn())

Create a random function on the interval `[0,N]` that is piecewise constant over the intervals `[i,i+1]`, `i = 1, ..., N`. These functions have the nice property that you can isomorphically embed them in `R^N` by mapping a function to its length-`N` list of steps.

# Arguments
- `N::Integer`: the number of steps that the function should take.

# Keyword parameters
- `step_generator`: a function that takes a step number `i` as input and returns a value for the function at the `i`th step. Note that the steps start from index `1`.

# Returns
Returns the step function as well as the `N` steps.
"""
function create_step_function(N::Integer; step_generator = i -> randn())
    steps = map(step_generator, 1:N)

    return (x -> _step_function(x, steps), steps)
end

function _step_function(x, steps)
    @assert 0 ≤ x ≤ length(steps)
    current_step  = Int64(ceil(x))
    steps[current_step]
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*csc(x))^-1*tan(x)^4, x) == :(x*a^-1+-1//5*a^-1*tan(x)^5*(1+-1*csc(x))+-1//15*a^-1*(15+-8*csc(x))*tan(x)+1//15*a^-1*tan(x)^3*(5+-4*csc(x)))
@test integrate((a+a*csc(x))^-1*tan(x)^3, x) == :(-1//8*a^-1*(1+sin(x))^-2+1//8*a^-1*(1+-1*sin(x))^-1+3//4*a^-1*(1+sin(x))^-1+5//16*a^-1*log(1+-1*sin(x))+11//16*a^-1*log(1+sin(x)))
@test integrate((a+a*csc(x))^-1*tan(x)^2, x) == :(-1*x*a^-1+-1//3*a^-1*tan(x)^3*(1+-1*csc(x))+1//3*a^-1*(3+-2*csc(x))*tan(x))
@test integrate((a+a*csc(x))^-1*tan(x), x) == :(-3//4*a^-1*log(1+sin(x))+-1//2*a^-1*(1+sin(x))^-1+-1//4*a^-1*log(1+-1*sin(x)))
@test integrate((a+a*csc(x))^-1*cot(x), x) == :(a^-1*log(1+sin(x)))
@test integrate((a+a*csc(x))^-1*cot(x)^2, x) == :(-1*x*a^-1+-1*a^-1*arctanh(cos(x)))
@test integrate((a+a*csc(x))^-1*cot(x)^3, x) == :(-1*a^-1*csc(x)+-1*a^-1*log(sin(x)))
@test integrate((a+a*csc(x))^-1*cot(x)^4, x) == :(x*a^-1+(1/2)*a^-1*arctanh(cos(x))+(1/2)*a^-1*(2+-1*csc(x))*cot(x))
@test integrate((a+a*csc(x))^-1*cot(x)^5, x) == :(a^-1*csc(x)+a^-1*log(sin(x))+(1/2)*a^-1*csc(x)^2+-1//3*a^-1*csc(x)^3)
@test integrate((a+a*csc(x))^-1*cot(x)^6, x) == :(-1*x*a^-1+-3//8*a^-1*arctanh(cos(x))+-1//8*a^-1*(8+-3*csc(x))*cot(x)+1//12*a^-1*cot(x)^3*(4+-3*csc(x)))
@test integrate((a+a*csc(x))^-1*cot(x)^7, x) == :(-1*a^-1*csc(x)^2+-1*a^-1*csc(x)+-1*a^-1*log(sin(x))+-1//5*a^-1*csc(x)^5+1//4*a^-1*csc(x)^4+2//3*a^-1*csc(x)^3)
@test integrate((a+b*csc(x))^-1*tan(x)^5, x) == :((1+-1*csc(x))^-2*(16a+16b)^-1+(1+csc(x))^-2*(-16b+16a)^-1+-1*a^-1*log(sin(x))+-1//16*(a+b)^-3*(8*a^2+15*b^2+21*a*b)*log(1+-1*csc(x))+-1//16*(a+-1b)^-3*(8*a^2+15*b^2+-21*a*b)*log(1+csc(x))+1//16*(1+-1*csc(x))^-1*(a+b)^-2*(5a+7b)+1//16*(1+csc(x))^-1*(a+-1b)^-2*(-7b+5a)+a^-1*b^6*(a^2+-1*b^2)^-3*log(a+b*csc(x)))
@test integrate((a+b*csc(x))^-1*tan(x)^3, x) == :(a^-1*log(sin(x))+-1*(1+-1*csc(x))^-1*(4a+4b)^-1+-1*(1+csc(x))^-1*(-4b+4a)^-1+1//4*(a+b)^-2*(2a+3b)*log(1+-1*csc(x))+1//4*(a+-1b)^-2*(-3b+2a)*log(1+csc(x))+a^-1*b^4*(a^2+-1*b^2)^-2*log(a+b*csc(x)))
@test integrate((a+b*csc(x))^-1*tan(x), x) == :(-1*a^-1*log(sin(x))+-1*(-2b+2a)^-1*log(1+csc(x))+-1*(2a+2b)^-1*log(1+-1*csc(x))+a^-1*b^2*(a^2+-1*b^2)^-1*log(a+b*csc(x)))
@test integrate((a+b*csc(x))^-1*cot(x), x) == :(a^-1*log(a+b*csc(x))+a^-1*log(sin(x)))
@test integrate((a+b*csc(x))^-1*cot(x)^3, x) == :(-1*a^-1*log(sin(x))+-1*b^-1*csc(x)+-1*a^-1*(1+-1*a^2*b^-2)*log(a+b*csc(x)))
@test integrate((a+b*csc(x))^-1*cot(x)^5, x) == :(a^-1*log(sin(x))+-1//3*b^-1*csc(x)^3+(1/2)*a*b^-2*csc(x)^2+-1*b^-3*(a^2+-2*b^2)*csc(x)+a^-1*b^-4*(a^2+-1*b^2)^2*log(a+b*csc(x)))
@test integrate((a+b*csc(x))^-1*cot(x)^7, x) == :(-1*a^-1*log(sin(x))+-1//5*b^-1*csc(x)^5+-1*b^-5*(a^4+3*b^4+-3*a^2*b^2)*csc(x)+-1//3*b^-3*csc(x)^3*(a^2+-3*b^2)+1//4*a*b^-2*csc(x)^4+a^-1*b^-6*(a^2+-1*b^2)^3*log(a+b*csc(x))+(1/2)*a*b^-4*csc(x)^2*(a^2+-3*b^2))
@test integrate((a+b*csc(x))^-1*cot(x)^2, x) == :(-1*x*a^-1+-1*b^-1*arctanh(cos(x))+2*a^-1*b^-1*(a^2+-1*b^2)^(1/2)*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*cot(x)^4, x) == :(x*a^-1+a*b^-2*cot(x)+-1//2*b^-1*cot(x)*csc(x)+-1//2*b^-3*(-3*b^2+2*a^2)*arctanh(cos(x))+2*a^-1*b^-3*(a^2+-1*b^2)^3//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*cot(x)^6, x) == :(-1*x*a^-1+-3//8*b^-1*arctanh(cos(x))+a*b^-2*cot(x)+-1*b^-5*(a^4+3*b^4+-3*a^2*b^2)*arctanh(cos(x))+-3//8*b^-1*cot(x)*csc(x)+-1//2*b^-3*(a^2+-3*b^2)*arctanh(cos(x))+-1//4*b^-1*csc(x)^3*cot(x)+1//3*a*b^-2*cot(x)^3+a*b^-4*(a^2+-3*b^2)*cot(x)+2*a^-1*b^-5*(a^2+-1*b^2)^5//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x)))+-1//2*b^-3*(a^2+-3*b^2)*cot(x)*csc(x))

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.1.4 (d cot)^n (a+b csc)^m.jl


text:
module CoreStarFg

using ..Ahorn, Maple

placements = Maple.CoreStarFg

function Ahorn.canFgBg(effect::Maple.CoreStarFg)
    return true, false
end

end


text:
@polly function kernel_doitgen(A,C4,sum)
    nr,nq,np = size(A)

    for r = 1:nr, q = 1:nq
        for p = 1:np
            sum[p] = zero(eltype(sum))
            for s = 1:np
                sum[p] += A[r,q,s] * C4[s,p]
            end
        end
        for p = 1:np
            A[r,q,p] = sum[p]
        end
    end
end

let
    nq = 140
    nr = 150
    np = 160

    A = zeros(Float32,nr,nq,np)
    sum = zeros(Float32,np)
    C4 = zeros(Float32,np,np)

    for i = 1:nr, j = 1:nq, k = 1:np
        A[i,j,k] = ((i*j + k)%np) / np
    end
    for i = 1:np, j = 1:np
        C4[i,j] = (i*j % np) / np
    end

    SUITE["doitgen"] = @benchmarkable kernel_doitgen(A, C4, sum) setup = (A = copy($A); C4 = copy($C4); sum = copy($sum))
end


text:
module FractionalTransforms

using LinearAlgebra, DSP, ToeplitzMatrices, FFTW

include("frft.jl")
include("frst.jl")
include("frct.jl")

export frft
export freq_shear, time_shear, sinc_interp

export frst
export frct

end


text:
using Turing, Random, Test, LinearAlgebra

include("../test_utils/AllUtils.jl")

@testset "advi.jl" begin
    @turing_testset "advi constructor" begin
        Random.seed!(0)
        N = 500

        s1 = ADVI()
        q = vi(gdemo_default, s1)
        c1 = rand(q, N)
    end
    @numerical_testset "advi inference" begin
        Random.seed!(1)
        N = 500

        alg = ADVI(10, 5000)
        q = vi(gdemo_default, alg)
        samples = transpose(rand(q, N))
        chn = Chains(reshape(samples, size(samples)..., 1), ["s", "m"])

        # TODO: uhmm, seems like a large `eps` here...
        check_gdemo(chn, eps = 0.5)
    end
end


text:
export print, init_output, final_output

###############################################################################
####  Relaxation context
###############################################################################
function print_build_relctx(relax_ctx, pb)
    relaxparams = relax_ctx.relaxparams

    outstream = []
    relaxparams[:opt_outmode]!=1 && push!(outstream, STDOUT)
    relaxparams[:opt_outmode]>0  && push!(outstream, open(relaxparams[:opt_outname], "a"))

    (relaxparams[:opt_outlev] == 0) && return

    # Compute indicators
    nb_densecstrs = 0
    maxdeg_densecstr = Float32[]
    for (cstr, ki_) in relax_ctx.ki
        if relax_ctx.di[cstr] > ki_
            nb_densecstrs += 1
            push!(maxdeg_densecstr, ki_)
        end
    end

    exposet = Set()
    nb_expotot = 0
    degbycstr = Int64[]
    for (cstrname, cstr) in get_constraints(pb)
        push!(degbycstr, max(cstr.p.degree.explvar, cstr.p.degree.conjvar))
        for (expo, λ) in cstr.p
            push!(exposet, expo)
            nb_expotot += 1
        end
    end

    # print
    for outstr in outstream
        if relaxparams[:opt_outlev] ≥ 1
            println(outstr, "\n=== set_relaxation()")

            print(outstr, "Relaxation is ", string(relaxparams[:opt_hierarchykind]))
            relaxparams[:opt_issparse] && print(outstr, ", sparse")
            relaxparams[:opt_multiordered] && print(outstr, ", multiordered")
            println(outstr, ".")

            println(outstr, "Global order is : ", relaxparams[:opt_globalorder], ".")

            print(outstr, "PhaseInvariance ")
            !relaxparams[:opt_sym_phaseinv] && print(outstr, "not ")
            print(outstr, "required and ")
            !relaxparams[:pb_isphaseinv] && print(outstr, "not ")
            println(outstr, "found.")

            println(outstr, "\n=> Problem with:")
            println(outstr, "-> Nb of variables:                        ", length(pb.variables))
            println(outstr, "-> Nb of constraints:                      ", length(pb.constraints))
            println(outstr, "-> Nb of monomials:                        $nb_expotot ($(length(exposet)) different)")
            @printf(outstr, "-> Max total degree by constraint:         %.1f / %.1f (mean/std)\n", mean(degbycstr), std(degbycstr))

            println(outstr, "=> Relaxation characteristics:")
            @printf(outstr, "-> Number of constraints s.t. di > ki:     %i / %i\n", nb_densecstrs, length(pb.constraints))
            @printf(outstr, "-> Max total degree on such constraints:   %.1f / %.1f (mean/std)\n", mean(maxdeg_densecstr), std(maxdeg_densecstr))
            @printf(outstr, "All variables appearing in such constraints will be linked in the sparsity pattern, which will largely densify it.\n")
        end

        (outstr!=STDOUT) && close(outstr)
    end
end

function print(io::IO, relctx::RelaxationContext)
    print(io, "RelaxationContext:\n")
    print(io, "ismultiordered         : $(relctx.ismultiordered)\n")
    print(io, "issparse               : $(relctx.issparse)\n")
    print(io, "symmetries             : $(relctx.symmetries)\n")
    print(io, "hierarchykind          : $(relctx.hierarchykind)\n")
    print(io, "renamevars             : $(relctx.renamevars)\n")
    for cstrname in sort(collect(keys(relctx.di)))
        di = relctx.di[cstrname]
        print(io, "di                     : $cstrname  \t=> $di\n")
    end
    for cstrname in sort(collect(keys(relctx.ki)))
        ki = relctx.ki[cstrname]
        print(io, "ki                     : $cstrname  \t=> $ki\n")
    end
    for cstrname in sort(collect(keys(relctx.cstrtypes)))
        cstrtype = relctx.cstrtypes[cstrname]
        print(io, "bar var types          : $cstrname  \t=> $(string(cstrtype))\n")
    end
end

###############################################################################
####  Moment problem construction
###############################################################################
function print_build_momentrelax(relax_ctx, momentrelaxation, nb_expos)
    relaxparams = relax_ctx.relaxparams

    (relaxparams[:opt_outlev] == 0) && return

    outstream = []
    relaxparams[:opt_outmode]!=1 && push!(outstream, STDOUT)
    relaxparams[:opt_outmode]>0  && push!(outstream, open(relaxparams[:opt_outname], "a"))

    for outstr in outstream
        if relaxparams[:opt_outlev] ≥ 1
            println(outstr, "\n=== build_momentrelaxation(relax_ctx, problem, moment_param::Dict{String, Tuple{Set{String}, Int}}, max_cliques::Dict{String, Set{Variable}})")
            println(outstr, "Compute the moment and localizing matrices associated with the problem constraints and clique decomposition and return a SDPDual object.")

            if relaxparams[:opt_outlev] ≥ 3
                print(outstr, momentrelaxation)
            end

            println(outstr, "Number of moments      : ", nb_expos)
            println(outstr, "Nb exponents coupled   : ", length(momentrelaxation.moments_overlap))

            ## NOTE: which relevant indicators here ?
        end

        (outstr!=STDOUT) && close(outstr)
    end
end

function print(io::IO, momentrelax::SDPDual{T}) where T
    println(io, "Moment Relaxation Problem:")
    println(io, "→ Objective: ")
    momentlen = maximum(x->length(string(x)), keys(momentrelax.objective))
    for moment in sort(collect(keys(momentrelax.objective)))
        coeff = momentrelax.objective[moment]
        print_string(io, string(moment), momentlen)
        println(io, " $coeff")
    end

    println(io, "→ Constraints:")
    for (cstrname, blocname) in sort(collect(keys(momentrelax.constraints)))
        mmtmat = momentrelax.constraints[(cstrname, blocname)]
        println(io, " → $cstrname, $blocname")
        println(io, mmtmat)
    end

    println(io, "→ Moments clique overlap:")
    if length(momentrelax.moments_overlap) > 0
        mmtlength = maximum(x->length(string(x)), keys(momentrelax.moments_overlap))
        for moment in sort(collect(keys(momentrelax.moments_overlap)))
            cliquenames = momentrelax.moments_overlap[moment]
            print(io, " → ")
            print_string(io, string(moment), mmtlength)
            for clique in sort(collect(cliquenames)) print(io, "$clique, ") end
            @printf(io, "\b\b \n")
        end
    else
        print(io, "  None")
    end
end

###############################################################################
####  SOS problem construction
###############################################################################
function print_build_SOSrelax(relax_ctx::RelaxationContext, sosrel::SDPPrimal)
    relaxparams = relax_ctx.relaxparams

    (relaxparams[:opt_outlev] == 0) && return

    outstream = []
    relaxparams[:opt_outmode]!=1 && push!(outstream, STDOUT)
    relaxparams[:opt_outmode]>0  && push!(outstream, open(relaxparams[:opt_outname], "a"))

    ## Compute indicators
    nb_SDPvars = length(Set([(key[2]) for key in keys(sosrel.blocks)]))
    size_SDPvars = length(Set([(key[2], key[3], key[4]) for key in keys(sosrel.blocks)]))
    nb_symvars = length(Set([(key[2]) for key in keys(sosrel.linsym)]))
    size_symvars = length(Set([(key[2], key[3]) for key in keys(sosrel.linsym)]))
    nb_scalvars = length(Set([key[3] for key in keys(sosrel.lin)]))

    momentset = Set{Moment}([key[1] for key in keys(sosrel.blocks)])
    union!(momentset, Set([key[1] for key in keys(sosrel.linsym)]))
    union!(momentset, Set([key[1] for key in keys(sosrel.lin)]))
    union!(momentset, Set([key for key in keys(sosrel.cst)]))
    for moment in momentset
        if product(moment.expl_part, moment.conj_part) == Exponent()
            delete!(momentset, moment)
        end
    end

    for outstr in outstream
        if relaxparams[:opt_outlev] ≥ 1
            println(outstr, "\n=== SOSrelaxation")

            if relaxparams[:opt_outlev] ≥ 2
                warn(outstr, "(C)SDP variables :")
                println(outstr, sort([(key[2]) for key in keys(sosrel.blocks)]))
                warn(outstr, "Symmetric variables :")
                println(outstr, sort([(key[2]) for key in keys(sosrel.linsym)]))
                warn(outstr, "Scalar variables :")
                println(outstr, sort([key[3] for key in keys(sosrel.lin)]))
                warn(outstr, "Constraint keys :")
                println(outstr, sort(collect(momentset)))
            end

            println(outstr, "- nb of (C)SDP matrix vars     : ", nb_SDPvars)
            println(outstr, "- size of (C)SDP matrix vars   : ", size_SDPvars)
            println(outstr, "- nb of (C)Sym matrix vars     : ", nb_symvars)
            println(outstr, "- size of (C)Sym matrix vars   : ", size_symvars)
            println(outstr, "- nb of scalar variables       : ", nb_scalvars)
            println(outstr, "- nb of constraints            : ", length(momentset))

            if relaxparams[:opt_outlev] ≥ 3
                print(outstr, sosrel)
            end
            ## NOTE: which relevant indicators here ?
        end

        (outstr!=STDOUT) && close(outstr)
    end
end

function print(io::IO, sdpinst::SDPPrimal)
    println(io, " -- SDP Blocks:")
    print(io, sdpinst.blocks)
    println(io, " -- linear part:")
    print(io, sdpinst.lin, sdpinst.linsym)
    println(io, " -- const part:")
    print(io, sdpinst.cst)
    println(io, " -- mat var types:")
    for blockname in sort(collect(keys(sdpinst.block_to_vartype)))
        blocktype = sdpinst.block_to_vartype[blockname]
        println(io, "   $blockname  \t $blocktype")
    end
end

function print(io::IO, sdpblocks::Dict{Tuple{Moment, String, Exponent, Exponent}, T}; indentedprint=true) where T
    print_blocksfile(io, sdpblocks; indentedprint=indentedprint, print_header=false)
end

function print(io::IO, sdplin::Dict{Tuple{Moment, Exponent}, T}, sdplinsym::Dict{Tuple{Moment, String, Exponent}, T}; indentedprint=true) where T
    print_linfile(io, sdplin, sdplinsym; indentedprint=indentedprint, print_header=false)
end

function print(io::IO, sdpcst::Dict{Moment, T}; indentedprint=true) where T
    print_cstfile(io, sdpcst; indentedprint=indentedprint, print_header=false)
end

function init_output(relax_ctx::RelaxationContext)
    relaxparams = relax_ctx.relaxparams

    # Create empty log file
    if relaxparams[:opt_outmode] ≥ 1
        isfile(relaxparams[:opt_outname]) && rm(relaxparams[:opt_outname])
        open(relaxparams[:opt_outname], "w") do f
            repo = LibGit2.GitRepo(pwd()); branch = LibGit2.shortname(LibGit2.head(repo))

            println(f, "MomentSOS hierarchy, date: ", String(Dates.format(now(), "mm_dd-HH:MM:SS")))
        end
    end

    # Create minimal csv file
    write_ctxcsv(relax_ctx.relaxparams)
end

function final_output(relax_ctx::RelaxationContext)
    relaxparams = relax_ctx.relaxparams

    # Print CSV file
    write_ctxcsv(relaxparams)
end

function write_ctxcsv(relaxparams)
    if relaxparams[:opt_outcsv] ≥ 1
        isfile(relaxparams[:opt_outcsvname]) && rm(relaxparams[:opt_outcsvname])
        open(relaxparams[:opt_outcsvname], "w") do f
            for key in keys(relaxparams)
                print(f, key, ";")
            end
            println(f)

            for val in values(relaxparams)
                print(f, val, ";")
            end
        end
    end
end


text:
# Truncated normal distribution
"""
    TruncatedNormal(mu, sigma, l, u)

The *truncated normal distribution* is a particularly important one in the family of truncated distributions.
We provide additional support for this type with `TruncatedNormal` which calls `Truncated(Normal(mu, sigma), l, u)`.
Unlike the general case, truncated normal distributions support `mean`, `mode`, `modes`, `var`, `std`, and `entropy`.
"""
TruncatedNormal

@deprecate TruncatedNormal(mu::Real, sigma::Real, a::Real, b::Real) truncated(Normal(mu, sigma), a, b)

### statistics

minimum(d::Truncated{Normal{T},Continuous}) where {T <: Real} = d.lower
maximum(d::Truncated{Normal{T},Continuous}) where {T <: Real} = d.upper

function mode(d::Truncated{Normal{T},Continuous}) where T <: Real
    μ = mean(d.untruncated)
    d.upper < μ ? d.upper :
    d.lower > μ ? d.lower : μ
end

modes(d::Truncated{Normal{T},Continuous}) where {T <: Real} = [mode(d)]

# do not export. Used in mean
# computes mean of standard normal distribution truncated to [a, b]
function _tnmom1(a, b)
    if !(a ≤ b)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return middle(a, b)
    elseif abs(a) > abs(b)
        return -_tnmom1(-b, -a)
    elseif isinf(a) && isinf(b)
        return zero(middle(a, b))
    end
    Δ = (b - a) * middle(a, b)
    if a ≤ 0 ≤ b
        m = √(2/π) * expm1(-Δ) * exp(-a^2 / 2) / (erf(a/√2) - erf(b/√2))
    elseif 0 < a < b
        z = exp(-Δ) * erfcx(b/√2) - erfcx(a/√2)
        iszero(z) && return middle(a, b)
        m = √(2/π) * expm1(-Δ) / z
    end
    return clamp(m, a, b)
end

# do not export. Used in var
# computes 2nd moment of standard normal distribution truncated to [a, b]
function _tnmom2(a::Real, b::Real)
    if !(a ≤ b)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return middle(a, b)^2
    elseif abs(a) > abs(b)
        return _tnmom2(-b, -a)
    elseif isinf(a) && isinf(b)
        return one(middle(a, b))
    elseif isinf(b)
        return 1 + √(2 / π) * a / erfcx(a / √2)
    end

    if a ≤ 0 ≤ b
        ea = √(π/2) * erf(a / √2)
        eb = √(π/2) * erf(b / √2)
        fa = ea - a * exp(-a^2 / 2)
        fb = eb - b * exp(-b^2 / 2)
        m2 = (fb - fa) / (eb - ea)
        return m2
    else # 0 ≤ a ≤ b
        exΔ = exp((a - b)middle(a, b))
        ea = √(π/2) * erfcx(a / √2)
        eb = √(π/2) * erfcx(b / √2)
        fa = ea + a
        fb = eb + b
        m2 = (fa - fb * exΔ) / (ea - eb * exΔ)
        return m2
    end
end

# do not export. Used in var
function _tnvar(a::Real, b::Real)
    if a == b
        return zero(middle(a, b))
    elseif a < b
        m1 = _tnmom1(a, b)
        m2 = √_tnmom2(a, b)
        return (m2 - m1) * (m2 + m1)
    else
        return oftype(middle(a, b), NaN)
    end
end

function mean(d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    if iszero(σ)
        return mode(d)
    else
        a = (d.lower - μ) / σ
        b = (d.upper - μ) / σ
        return μ + _tnmom1(a, b) * σ
    end
end

function var(d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    if iszero(σ)
        return σ
    else
        a = (d.lower - μ) / σ
        b = (d.upper - μ) / σ
        return _tnvar(a, b) * σ^2
    end
end

function entropy(d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    z = d.tp
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    aφa = isinf(a) ? 0.0 : a * normpdf(a)
    bφb = isinf(b) ? 0.0 : b * normpdf(b)
    0.5 * (log2π + 1.) + log(σ * z) + (aφa - bφb) / (2.0 * z)
end

### sampling

## Use specialized sampler, as quantile-based method is inaccurate in
## tail regions of the Normal, issue #343

function rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    if isfinite(μ)
        a = (d.lower - μ) / σ
        b = (d.upper - μ) / σ
        z = randnt(rng, a, b, d.tp)
        return μ + σ * z
    else
        return clamp(μ, d.lower, d.upper)
    end
end

# Rejection sampler based on algorithm from Robert (1995)
#
#  - Available at http://arxiv.org/abs/0907.4010

function randnt(rng::AbstractRNG, lb::AbstractFloat, ub::AbstractFloat, tp::AbstractFloat)
    local r::AbstractFloat
    if tp > 0.3   # has considerable chance of falling in [lb, ub]
        r = randn(rng)
        while r < lb || r > ub
            r = randn(rng)
        end
        return r

    else
        span = ub - lb
        if lb > 0 && span > 2.0 / (lb + sqrt(lb^2 + 4.0)) * exp((lb^2 - lb * sqrt(lb^2 + 4.0)) / 4.0)
            a = (lb + sqrt(lb^2 + 4.0))/2.0
            while true
                r = rand(rng, Exponential(1.0 / a)) + lb
                u = rand(rng)
                if u < exp(-0.5 * (r - a)^2) && r < ub
                    return r
                end
            end
        elseif ub < 0 && ub - lb > 2.0 / (-ub + sqrt(ub^2 + 4.0)) * exp((ub^2 + ub * sqrt(ub^2 + 4.0)) / 4.0)
            a = (-ub + sqrt(ub^2 + 4.0)) / 2.0
            while true
                r = rand(rng, Exponential(1.0 / a)) - ub
                u = rand(rng)
                if u < exp(-0.5 * (r - a)^2) && r < -lb
                    return -r
                end
            end
        else
            while true
                r = lb + rand(rng) * (ub - lb)
                u = rand(rng)
                if lb > 0
                    rho = exp((lb^2 - r^2) * 0.5)
                elseif ub < 0
                    rho = exp((ub^2 - r^2) * 0.5)
                else
                    rho = exp(-r^2 * 0.5)
                end
                if u < rho
                    return r
                end
            end
        end
    end
end


text:
"""
    wasserstein(p1::AbstractParticles,p2::AbstractParticles,p)

Returns the Wasserstein distance (Earth-movers distance) of order `p`, to the `p`th power, between `p1` and `p2`.
I.e., for `p=2`, this returns W₂²
"""
function wasserstein(p1::AbstractParticles,p2::AbstractParticles,p)
    p1 = sort(p1.particles)
    p2 = sort(p2.particles)
    wp = mean(eachindex(p1)) do i
        abs(p1[i]-p2[i])^p
    end
end