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

using Markdown
using InteractiveUtils

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

# ╔═╡ d5dc9294-34a7-11eb-09b9-e36684398036
using Plots

# ╔═╡ f75649f2-33e9-11eb-0a65-0f48aa7209a6
md"## Heat equation with Neumann boundary conditions

The domain is given by the one-dimensional rod $\Omega = (0,L)$ and the time by $t \in (0,T).$ The rod has the physical properties
- thermal conductivity $\lambda$,
- specific heat capacity $c$,
- mass density $\rho$.

The diffusivity constant is given by $\alpha = \frac{\lambda}{c \rho}$. The heat equation is noted as

$\dot{\vartheta}(t,x) = \alpha \frac{\partial^2}{\partial x^2} \vartheta(t,x)$

for $(t,x) \in (0,T) \times \Omega$ with initial condition

$\vartheta(0,x) = \vartheta_{0}(x)$

for $x \in \overline{\Omega}$ and boundary conditions

$\left.\frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = 0} \cdot \vec{n} ~=~ \left. \frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = L} \cdot \vec{n} ~=~ 0$

with outer [normal vector](https://en.wikipedia.org/wiki/Normal_(geometry)) $\vec{n}$ on the left or right boundary. Here, the normal vector $\vec{n} = -1$ on the left boundary at $x=0$, and  $\vec{n} = 1$ on the right boundary at $x=L$."

# ╔═╡ 038013f8-34a6-11eb-37b1-69150fede3db
md"#### Code: implementing physical constants"

# ╔═╡ 0366d38c-34a6-11eb-2cfe-05e6f8edf17b
# Length of rod
L = 0.5

# ╔═╡ 03494950-34a6-11eb-32f5-7b485d899c42
# Final time
Tf = 10000.0

# ╔═╡ 032aaeb6-34a6-11eb-2742-b12ab7682f17
# Thermal conductivity
λ = 45.0

# ╔═╡ 03002010-34a6-11eb-2b2f-9dd02826187b
# Specific heat capacity
c = 480.0

# ╔═╡ 02e3e238-34a6-11eb-2d89-310870da991b
# Mass density
ρ = 7800.0

# ╔═╡ 44a6fa0c-34a6-11eb-2917-779a9aa42543
# Thermal diffusivity
α = λ/(c*ρ)

# ╔═╡ b1a58804-34a5-11eb-3d49-25b877f17b75
md"### Separation of variables

It is assumed that $\quad \vartheta(t,x) = f(t) ~ g(x) \quad$ holds. Applying [separation of variables](https://en.wikipedia.org/wiki/Separation_of_variables) approach into the heat equation, one yields

$\dot{\vartheta}(t,x) = \dot{f}(t) ~ g(x) =  \alpha ~ f(t) ~ \frac{\partial^2}{\partial x^2} g(x).$

Separating the functions results in

$\frac{\dot{f}(t)}{f(t)} = \alpha ~ \frac{g''(x)}{g(x)} = p.$

The ordinary differential equation $\quad \dot{f}(t) = p ~ f(t) \quad$ is solved by $\quad f(t) = C_{f} ~ \exp\left(p~t\right)$. And the second order ODE $\quad g''(x) = q ~ g(x) \quad$ with $q = \frac{p}{\alpha}$ is solved by

$g(x) ~=~ A_{g} ~ \sin(\sqrt{-q} x) + B_{g} \cos(\sqrt{-q} x).$

### Eigenvalues and Eigenvectors

Applying the initial condition at the boundaries leads to

$\left.\frac{\partial}{\partial x} \vartheta(0, x) \right\rvert_{x = 0} ~=~ f(0) ~ g'(0) = C_{f} ~ \sqrt{-q} ~ A_{g} ~=~ 0$

implying $A_{g}=0$, and

$\left. \frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = L} ~=~ f(0) g'(L) = C_{f} ~ \sqrt{-q} ~ (-B_{g}) ~  \sin(\sqrt{-q} x) ~=~ 0.$

Thus, the Eigenvalues are calculated by $\quad \sqrt{-q}_{i}~\sin\left( \sqrt{-q}_{i} L \right) = 0 \quad$ or equivalent $\quad q_{i} = 0$ for all $i \in \mathbb{N}$ or  $\quad q_{i} =  \left(\frac{i ~ \pi}{L}\right)^{2}$.

Variable $B_{i} = B_{g,i}$ as coefficient of spatial solution

$g_{i}(x) = B_{g,i} \cos\left( \sqrt{-q}_{i} x \right)$

is introduced and used to find the Eigenvectors or Eigenfunctions

$\Phi_{i}(x) ~=~  B_{i} ~ \cos\left( \sqrt{-q_{i}} x \right) ~=~  B_{i}  \cos\left( i ~ \pi ~ \frac{x}{L} \right)$

The Eigenvectors form an [orthonormal basis](https://en.wikipedia.org/wiki/Orthonormal_basis) and thus

$\langle \Phi_{i}, \Phi_{j} \rangle ~=~ \int\limits_{0}^{L} B_{i}  \cos\left( i ~ \pi ~ \frac{x}{L} \right) ~ B_{j}  \cos\left( j ~ \pi ~ \frac{x}{L} \right) \operatorname{d}x ~=~ \delta_{i,j}$

with $\delta_{i,j} ~=~
	\begin{cases}
	1 \quad \text{for} \quad i = j,\\
	0 \quad \text{otherwise}.
	\end{cases}$

In the case $i = j$ one finds with

$\langle \Phi_{i}, \Phi_{i} \rangle = B_{i}^2 \int\limits_{0}^{L}  \cos\left( i ~ \pi ~ \frac{x}{L} \right)^2 \operatorname{d}x = 1$

the coefficients $B_{i} ~=~ \sqrt{\frac{2}{L}}$.

In the case $i \neq j$ the equality

$\langle \Phi_{i}, \Phi_{j} \rangle ~=~ B_{i} B_{j}  \int\limits_{0}^{L} \cos\left( i ~ \pi ~ \frac{x}{L} \right) ~  \cos\left( j ~ \pi ~ \frac{x}{L} \right) \operatorname{d}x = 0$

holds because $\int\limits_{0}^{L} \cos\left( i ~ \pi ~ \frac{x}{L} \right) ~  \cos\left( j ~ \pi ~ \frac{x}{L} \right) \operatorname{d}x = 0.$

Consequently, the Eigenvectors (or Eigenfunctions) are given by

$\Phi_{i}(x) = \sqrt{\frac{2}{L}} ~ \cos\left( i \pi \frac{x}{L}\right).$"

# ╔═╡ 137ff0a4-34a7-11eb-14cf-d5f7c5ad29f7
md"#### Code: Eigenvalues and Eigenvectors"

# ╔═╡ a1008ff4-34a9-11eb-2b19-270ae9e857ef
# Maximum index i
max_order = 10

# ╔═╡ 1365b150-34a7-11eb-3c6e-39d36506722f
# Eigenvalues
let
	i = 0 : 1 : max_order
	q = -1 * (i*pi/L).^2
	scatter(q, label="q", xaxis="Index i", yaxis="Eigenvalue")
end

# ╔═╡ 13426eaa-34a7-11eb-0ee6-f3282b7dcd20
@bind istep html"<input type='range' min='1' max='10' step='1'>"

# ╔═╡ 14490a74-34a8-11eb-0129-ed4370a3dc51
# Number of spatial elements
Nx = 2000

# ╔═╡ 142cef62-34a8-11eb-2df7-898cfde39332
# Finite discretization
Δx = L/(Nx-1)

# ╔═╡ 13ed2e20-34a8-11eb-233e-fd36ac4ceeca
# Discretized x positions
xspan = 0 : Δx : L

# ╔═╡ b1528564-34a5-11eb-20ec-3d3c2933cb36
let
	qi = - (istep*pi/L)^2
	Ai = sqrt(2/L)
	Φi(x) = Ai * cos( sqrt(-qi)*x)
	plot(xspan, Φi.(xspan), xlabel="Position x", ylabel="Eigenvector Phi", legend=false)
end

# ╔═╡ 13cfca10-34a8-11eb-140e-47b3a1b55a56

# ╔═╡ 0412e44a-34a8-11eb-3c66-e18f3e9f8789

# ╔═╡ 8d3b4652-34a5-11eb-1f85-cb954b678ea0
md"### Initial Conditions

The preliminary solution of the heat equation is noted by

$\vartheta(t,x) = \eta_{0} + \sum\limits_{i = 1}^{\infty} T_{i}(t) ~ \Phi_{i}(x) ~=~ \eta_{0} + \sum\limits_{i = 1}^{\infty} C_{f,i} ~ \exp\left( \alpha ~ q_{i} ~ t \right) \sqrt{\frac{2}{L}} \cos\left( i \pi \frac{x}{L} \right)$

**Remark:** the addional coefficient $\eta_{0}$ describes the behaviour of $T_{i}(t) ~ \Phi_{i}(x)$ at index $i=0$.

The coefficients $C_{f,i}$ are found with the initial conditions

$\vartheta(0,x) = \eta_{0} + \sum\limits_{i = 1}^{\infty} C_{i} ~ \Phi_{i}(x) ~=~ \vartheta_{0}(x)$

Both sides of the last equation are multiplied by the Eigenvector $\Phi_{i}$ and the inner product is computed elementwise to yield

$C_{i} \langle \Phi_{i}, \Phi_{i} \rangle ~=~ \langle \Phi_{i}, \vartheta_{0} -  \eta_{0}\rangle$

or equivalent

$C_{i} = \frac{\langle \Phi_{i}, \vartheta_{0}  -  \eta_{0}\rangle}{\langle \Phi_{i}, \Phi_{i} \rangle} = \langle \Phi_{i}, \vartheta_{0}  -  \eta_{0} \rangle$

because of the orthonormal property $\langle \Phi_{i}, \Phi_{i} \rangle = 1$.

Now, an example initial data is assumed with $\quad \vartheta_{0}(x) := m \left[ L~x - x^2 \right] \quad$ with  $m > 0$.

This initial data leads to the coefficients

$C_{i} ~=~ \sqrt{\frac{2}{L}}  ~ (-1) ~ m ~ L ~ \left(\frac{L}{i \pi}\right)^2 \left[\left(-1\right)^{i} + 1 \right] ~=~
\begin{cases}
-2 ~ m ~ L ~ \sqrt{\frac{2}{L}} ~ \left(\frac{L}{i \pi}\right)^2 \quad &\text{if} \quad i \text{ is even} \\
 0 \quad &\text{if} \quad i \text{ is odd}
\end{cases}$"

# ╔═╡ 3ba830e4-34a9-11eb-002b-4952d88fbaad
md"#### Code: Coefficients"

# ╔═╡ 482641dc-34a9-11eb-1794-792fd102aeab
# Amplification (initial data)
@bind amp html"<input type='range' min='0' max='100' step='0.1'>"

# ╔═╡ 48426022-34a9-11eb-2baa-8dfb566d9cce
let
	max_order = 10
	i = 1 : 1 : max_order
	Cf(i) = sqrt(2/L)*amp*(L/(i*pi))^3 * (1 - (-1)^i)
	scatter(i, Cf.(i), label="Cf", xaxis="i")
end

# ╔═╡ 3b902486-34a9-11eb-1529-45521466a6c8

# ╔═╡ 8d573ab0-34a5-11eb-1f0f-59bd9cb87092
md"###  Final solution

Finally, the solution of the heat equation with Dirichlet boundary conditions is given by

$\vartheta(t,x) = \eta_{0} - 2 m ~ \frac{L^2}{\pi^2} \sum\limits_{i = 1}^{\infty} \frac{1}{i^2} \left[\left(-1\right)^{i} + 1 \right] ~ \exp\left(- \alpha \left(\frac{i \pi }{L}\right)^2 ~ t \right) ~  \cos\left( i \pi \frac{x}{L}\right)$

and with odd index $i = 2 k$ one yields

$\vartheta(t,x) = \eta_{0} - m ~ \frac{L^2}{\pi^2}  ~ \sum\limits_{k = 1}^{\infty} \frac{1}{k^2} ~ \exp\left(- \alpha ~ 4 ~ \left(\frac{k \pi }{L}\right)^2 ~ t \right) ~  \cos\left( 2 ~k ~ \pi \frac{x}{L}\right)$"

# ╔═╡ 8ce0b9a8-34a5-11eb-16e2-b1858f9d89c4
md"#### Code: spatial and temporal approximation"

# ╔═╡ 8cc10270-34a5-11eb-0958-b566b8f488f9
function temporal_approx(elem, t, params)
    k = elem
    α = params[1]
    L = params[2]
    return exp( -4*α *(k*π/L)^2 *t)
end

# ╔═╡ 8ca3d774-34a5-11eb-1915-b31f19aad0d5
function spatial_approx(elem, x, params)
    k = elem
    L = params[2]
    return cos( 2*k* π*(x/L))
end

# ╔═╡ 8d1cb34a-34a5-11eb-0a91-cf59727f69b7
md"#### Coefficient $\eta_{0}$

The value of offset $\eta_{0}$ is found by computing the temperature at $(t,x)=(0,0)$ with

$\vartheta(t,x) = \eta_{0} - m ~ \frac{L^2}{\pi^2}  ~ \sum\limits_{k = 1}^{\infty} \frac{1}{k^2} = \vartheta_{0}(0).$

The sum of $\frac{1}{k^2}$ corresponds with the [Riemannian Zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at position 2  and can be used in

$\eta_{0} = \vartheta_{0}(0) + m ~ \frac{L^2}{\pi^2}  ~ \sum\limits_{k = 1}^{\infty} \frac{1}{k^2} = \vartheta_{0}(0) + m ~ \frac{L^2}{\pi^2}  ~ \zeta(2).$

According to the [Basel problem](https://en.wikipedia.org/wiki/Basel_problem), the Zeta function $\zeta(2) = \frac{\pi^2}{6}$ and consequently the coeffient is found with

$\eta_{0} \approx  \vartheta_{0}(0) + m ~ \frac{L^2}{6}.$"

# ╔═╡ 8d00d274-34a5-11eb-23e4-9b9a6999abe9
# Cofficient η₀
η₀ = amp * L^2/6

# ╔═╡ 8c81adb4-34a5-11eb-0221-2f9093992726
md"#### Code: Heat equation and solution"

# ╔═╡ 8c62f89c-34a5-11eb-3136-35c0b89ca4ad
function heat_eq(t,x, order, params)
    α = params[1]
    L = params[2]
    m = params[3]

    sum = 0;
    for i = 1 : order
        sum = sum + (1/i^2) * temporal_approx(i, t, params) * spatial_approx(i, x, params)
    end

    η = m * L^2/6;

    return η  -  m*L^2/(π^2) * sum
end

# ╔═╡ 8c46f208-34a5-11eb-3c16-37149eb6a2c7
# Parameter
params = [α, L, amp]

# ╔═╡ 8beb67c8-34a5-11eb-2d1c-9ba0b4dc71d9
# Number of time steps
Nt = 1001

# ╔═╡ 8bc809ac-34a5-11eb-3fe5-2756444c0e52
# Temporal discretization
Δt = Tf / (Nt - 1)

# ╔═╡ 8bad1f72-34a5-11eb-11b9-fddc7e984050
tspan = 0 : Δt : Tf

# ╔═╡ 8c2ae3c8-34a5-11eb-0a68-950ae952a525
# Solution array
θsol = zeros(length(xspan), length(tspan));

# ╔═╡ 8c0c6054-34a5-11eb-06e1-8335ef9544fb
for i = 1 : length(tspan)
    for k = 1 : length(xspan)
        θsol[k,i] = heat_eq(tspan[i], xspan[k], max_order, params)
    end
end

# ╔═╡ 8b92737c-34a5-11eb-28b1-3dc2e0aa2a64
md"#### Visualisation: spatial distribution"

# ╔═╡ 8b73e0b8-34a5-11eb-048c-892c233c0b05
@bind tstep html"<input type='range' min='1' max='1001' step='1'>"

# ╔═╡ 8b58e428-34a5-11eb-14a8-ed421bee29ba
Time = tspan[tstep]

# ╔═╡ 8b3e3880-34a5-11eb-3dbe-09e2003a6803
plot(xspan, θsol[:,tstep], xlabel="Position x", ylabel="Temperature", legend=false)

# ╔═╡ 5b4d530a-34ab-11eb-3bfb-7985e1617652
md"#### Visualisation: evolution in time"

# ╔═╡ 5ae93d48-34ab-11eb-0444-fd8ebb0108a3
@bind xstep html"<input type='range' min='1' max='2000' step='1'>"

# ╔═╡ ea54dbdc-33fb-11eb-3955-eb0e29d4e3a2
Position = xspan[xstep]

# ╔═╡ 69b29af4-34ab-11eb-26d3-a58d6f015fa1
plot(tspan, θsol[xstep,:], xlabel="Time t", ylabel="Temperature", legend=false)

# ╔═╡ Cell order:
# ╠═f75649f2-33e9-11eb-0a65-0f48aa7209a6
# ╠═038013f8-34a6-11eb-37b1-69150fede3db
# ╠═0366d38c-34a6-11eb-2cfe-05e6f8edf17b
# ╠═03494950-34a6-11eb-32f5-7b485d899c42
# ╠═032aaeb6-34a6-11eb-2742-b12ab7682f17
# ╠═03002010-34a6-11eb-2b2f-9dd02826187b
# ╠═02e3e238-34a6-11eb-2d89-310870da991b
# ╠═44a6fa0c-34a6-11eb-2917-779a9aa42543
# ╟─b1a58804-34a5-11eb-3d49-25b877f17b75
# ╟─137ff0a4-34a7-11eb-14cf-d5f7c5ad29f7
# ╠═d5dc9294-34a7-11eb-09b9-e36684398036
# ╠═a1008ff4-34a9-11eb-2b19-270ae9e857ef
# ╠═1365b150-34a7-11eb-3c6e-39d36506722f
# ╠═13426eaa-34a7-11eb-0ee6-f3282b7dcd20
# ╠═14490a74-34a8-11eb-0129-ed4370a3dc51
# ╠═142cef62-34a8-11eb-2df7-898cfde39332
# ╠═13ed2e20-34a8-11eb-233e-fd36ac4ceeca
# ╠═b1528564-34a5-11eb-20ec-3d3c2933cb36
# ╠═13cfca10-34a8-11eb-140e-47b3a1b55a56
# ╠═0412e44a-34a8-11eb-3c66-e18f3e9f8789
# ╟─8d3b4652-34a5-11eb-1f85-cb954b678ea0
# ╠═3ba830e4-34a9-11eb-002b-4952d88fbaad
# ╠═482641dc-34a9-11eb-1794-792fd102aeab
# ╠═48426022-34a9-11eb-2baa-8dfb566d9cce
# ╠═3b902486-34a9-11eb-1529-45521466a6c8
# ╠═8d573ab0-34a5-11eb-1f0f-59bd9cb87092
# ╠═8ce0b9a8-34a5-11eb-16e2-b1858f9d89c4
# ╠═8cc10270-34a5-11eb-0958-b566b8f488f9
# ╠═8ca3d774-34a5-11eb-1915-b31f19aad0d5
# ╠═8d1cb34a-34a5-11eb-0a91-cf59727f69b7
# ╠═8d00d274-34a5-11eb-23e4-9b9a6999abe9
# ╠═8c81adb4-34a5-11eb-0221-2f9093992726
# ╠═8c62f89c-34a5-11eb-3136-35c0b89ca4ad
# ╠═8c46f208-34a5-11eb-3c16-37149eb6a2c7
# ╠═8beb67c8-34a5-11eb-2d1c-9ba0b4dc71d9
# ╠═8bc809ac-34a5-11eb-3fe5-2756444c0e52
# ╠═8bad1f72-34a5-11eb-11b9-fddc7e984050
# ╠═8c2ae3c8-34a5-11eb-0a68-950ae952a525
# ╠═8c0c6054-34a5-11eb-06e1-8335ef9544fb
# ╠═8b92737c-34a5-11eb-28b1-3dc2e0aa2a64
# ╠═8b73e0b8-34a5-11eb-048c-892c233c0b05
# ╠═8b58e428-34a5-11eb-14a8-ed421bee29ba
# ╠═8b3e3880-34a5-11eb-3dbe-09e2003a6803
# ╠═5b4d530a-34ab-11eb-3bfb-7985e1617652
# ╠═5ae93d48-34ab-11eb-0444-fd8ebb0108a3
# ╠═ea54dbdc-33fb-11eb-3955-eb0e29d4e3a2
# ╠═69b29af4-34ab-11eb-26d3-a58d6f015fa1


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Tests for deprecated functionality.
#
# These can't be run with --depwarn=error, so currently require special
# treatment when run inside the test system.

using Test
using Logging

using Base: remove_linenums!

module DeprecationTests # to test @deprecate
    f() = true

    # test the Symbol path of @deprecate
    @deprecate f1 f
    @deprecate f2 f false # test that f2 is not exported

    # test the Expr path of @deprecate
    @deprecate f3() f()
    @deprecate f4() f() false # test that f4 is not exported
    @deprecate f5(x::T) where T f()

    # test deprecation of a constructor
    struct A{T} end
    @deprecate A{T}(x::S) where {T, S} f()

    # test that @deprecate_moved can be overridden by an import
    Base.@deprecate_moved foo1234 "Foo"
    Base.@deprecate_moved bar "Bar" false
end # module
module Foo1234
    export foo1234
    foo1234(x) = x+1
end

# issue #21972
struct T21972
    @noinline function T21972()
        Base.depwarn("something", :T21972)
        new()
    end
end

@testset "@deprecate" begin
    using .DeprecationTests
    using .Foo1234
    @test foo1234(3) == 4
    @test_throws ErrorException DeprecationTests.bar(3)

    # 22845
    ex = :(module M22845; import ..DeprecationTests: bar;
                          bar(x::Number) = x + 3; end)
    @test_warn "importing deprecated binding" eval(ex)
    @test @test_nowarn(DeprecationTests.bar(4)) == 7

    # enable when issue #22043 is fixed
    # @test @test_warn "f1 is deprecated, use f instead." f1()
    # @test @test_nowarn f1()

    # @test_throws UndefVarError f2() # not exported
    # @test @test_warn "f2 is deprecated, use f instead." DeprecationTests.f2()
    # @test @test_nowarn DeprecationTests.f2()

    # @test @test_warn "f3() is deprecated, use f() instead." f3()
    # @test @test_nowarn f3()

    # @test_throws UndefVarError f4() # not exported
    # @test @test_warn "f4() is deprecated, use f() instead." DeprecationTests.f4()
    # @test @test_nowarn DeprecationTests.f4()

    # @test @test_warn "f5(x::T) where T is deprecated, use f() instead." f5(1)
    # @test @test_nowarn f5(1)

    # @test @test_warn "A{T}(x::S) where {T, S} is deprecated, use f() instead." A{Int}(1.)
    # @test @test_nowarn A{Int}(1.)

    # issue #21972
    @noinline function f21972()
        T21972()
    end
    @test_deprecated "something" f21972()
end

f24658() = depwarn24658()

depwarn24658() = Base.firstcaller(backtrace(), :_func_not_found_)

@testset "firstcaller" begin
    # issue #24658
    @test eval(:(if true; f24658(); end)) == (Ptr{Cvoid}(0),StackTraces.UNKNOWN)
end

# issue #25130
f25130() = Base.depwarn("f25130 message", :f25130)
# The following test is for the depwarn behavior of expressions evaluated at
# top-level, so we can't use the usual `collect_test_logs()` / `with_logger()`
testlogger = Test.TestLogger()
prev_logger = global_logger(testlogger)
# Each call at top level should be distinct. This won't be true if they're
# attributed to internal C frames (including generic dispatch machinery)
f25130()
f25130()
testlogs = testlogger.logs
@test length(testlogs) == 2
@test testlogs[1].id != testlogs[2].id
@test testlogs[1].kwargs[:caller].func == Symbol("top-level scope")
@test all(l.message == "f25130 message" for l in testlogs)
global_logger(prev_logger)

#-------------------------------------------------------------------------------
# BEGIN 0.7 deprecations

@testset "parser syntax deprecations" begin
    # Test empty logs for meta.parse depwarn argument.
    @test_logs Meta.parse("1.+2", depwarn=false)

    # #19089
    @test (@test_deprecated Meta.parse("1.+2")) == :(1 .+ 2)

    # #16356
    @test (@test_deprecated Meta.parse("0xapi")) == :(0xa * pi)

    # #22523 #22712
    @test (@test_deprecated Meta.parse("a?b:c"))    == :(a ? b : c)
    @test (@test_deprecated Meta.parse("a ?b:c"))   == :(a ? b : c)
    @test (@test_deprecated Meta.parse("a ? b:c"))  == :(a ? b : c)
    @test (@test_deprecated Meta.parse("a ? b :c")) == :(a ? b : c)
    @test (@test_deprecated Meta.parse("?")) == Symbol("?")

    # #13079
    @test (@test_deprecated Meta.parse("1<<2*3")) == :(1<<(2*3))

    # ([#19157], [#20418]).
    @test remove_linenums!(@test_deprecated Meta.parse("immutable A; end")) ==
          remove_linenums!(:(struct A; end))
    @test remove_linenums!(@test_deprecated Meta.parse("type A; end")) ==
          remove_linenums!(:(mutable struct A; end))

    # #19987
    @test remove_linenums!(@test_deprecated Meta.parse("try ; catch f() ; end")) ==
          remove_linenums!(:(try ; catch; f() ; end))

    # #15524
    # @test (@test_deprecated Meta.parse("for a=b f() end")) == :(for a=b; f() end)
    @test_broken length(Test.collect_test_logs(()->Meta.parse("for a=b f() end"))[1]) > 0

    # #23076
    @test (@test_deprecated Meta.parse("[a,b;]")) == :([a;b])

    # #24452
    @test (@test_deprecated Meta.parse("(a...)")) == :((a...,))
end

@testset "lowering syntax deprecations" begin
    # #16295
    @test_deprecated Meta.lower(@__MODULE__, :(A.(:+)(a,b) = 1))

    # #11310
    @test_deprecated r"parametric method syntax" Meta.lower(@__MODULE__, :(f{T}(x::T) = 1))

    # #17623
    @test_deprecated r"Deprecated syntax `function .+(...)`" Meta.lower(@__MODULE__, :(function .+(a,b) ; end))

    # #21774 (more uniform let expressions)
    @test_deprecated Meta.lower(@__MODULE__, Expr(:let, :a))
    @test_deprecated Meta.lower(@__MODULE__, Expr(:let, :a, :(a=1), :(b=1)))

    # #23157 (Expression heads for types renamed)
    @test_deprecated Meta.lower(@__MODULE__, Expr(:type, true, :A, Expr(:block)))
    @test_deprecated Meta.lower(@__MODULE__, Expr(:bitstype, 32, :A))

    # #15032
    @test_deprecated Meta.lower(@__MODULE__, :(a.(b) = 1))

    # #5332
    @test_deprecated Meta.lower(@__MODULE__, :(a.'))

    # #19324
    @test_deprecated r"implicit assignment to global" eval(
           :(module M19324
                 x=1
                 for i=1:10
                     x += i
                 end
             end))

    # #24221
    @test_deprecated r"underscores as an rvalue" Meta.lower(@__MODULE__, :(a=_))

    # #22314
    @test_deprecated r"Use of final value of loop variable `i`.*is deprecated. In the future the variable will be local to the loop instead." Meta.lower(@__MODULE__, :(
        function f()
            i=0
            for i=1:10
            end
            i
        end))
    @test_deprecated r"Loop variable `i` overwrites a variable in an enclosing scope" eval(:(
        module M22314
            i=10
            for i=1:10
            end
        end))

    # #6080
    @test_deprecated r"Syntax `&argument`.*is deprecated" Meta.lower(@__MODULE__, :(ccall(:a, Cvoid, (Cint,), &x)))

    @test_logs eval(:(module DotEqualsDep
        a=[1,2]
        a.=3
        0
        end))
    @test_logs include_string(@__MODULE__, """
        a=[1,2]
        a.=3
        0""")
    @test_deprecated include_string(@__MODULE__, """
        a=[1,2]
        a.=3""")
end

module LogTest
    function bar(io)
        info(io,"barinfo")
        warn(io,"barwarn")
        Base.display_error(io,"barerror",backtrace())
    end
    function pooh(io)
        info(io,"poohinfo")
        warn(io,"poohwarn")
        Base.display_error(io,"pooherror",backtrace())
    end
end
function foo(io)
    info(io,"fooinfo")
    warn(io,"foowarn")
    Base.display_error(io,"fooerror",backtrace())
end

# Silence the flurry of depwarns for now.
with_logger(NullLogger()) do

@testset "Deprecated logging" begin

# Test info
@test contains(sprint(info, "test"), "INFO:")
@test contains(sprint(info, "test"), "INFO: test")
@test contains(sprint(info, "test ", 1, 2, 3), "INFO: test 123")
@test contains(sprint(io->info(io,"test", prefix="MYINFO: ")), "MYINFO: test")

# Test warn
@test contains(sprint(Base.warn_once, "test"), "WARNING: test")
@test isempty(sprint(Base.warn_once, "test"))

@test contains(sprint(warn), "WARNING:")
@test contains(sprint(warn, "test"), "WARNING: test")
@test contains(sprint(warn, "test ", 1, 2, 3), "WARNING: test 123")
@test contains(sprint(io->warn(io, "test", prefix="MYWARNING: ")), "MYWARNING: test")
@test contains(sprint(io->warn(io, "testonce", once=true)), "WARNING: testonce")
@test isempty(sprint(io->warn(io, "testonce", once=true)))
@test !isempty(sprint(io->warn(io, "testonce", once=true, key=hash("testonce",hash("testanother")))))
let bt = backtrace()
    ws = split(chomp(sprint(io->warn(io, "test", bt = bt))), '\n')
    bs = split(chomp(sprint(Base.show_backtrace, bt)), '\n')
    @test contains(ws[1],"WARNING: test")
    for (l,b) in zip(ws[2:end],bs[2:end])
        @test contains(l, b)
    end
end

# PR #16213
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest, :bar;  kind=:info)
@test all(contains.(sprint(LogTest.bar), ["WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest;  kind=:info)
@test all(contains.(sprint(LogTest.bar), ["WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull;  kind=:info)
@test all(contains.(sprint(LogTest.bar), ["WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(kind=:info)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest, :bar;  kind=:warn)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest;  kind=:warn)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull;  kind=:warn)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "ERROR: \"fooerror\""]))

logging(kind=:warn)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest, :bar;  kind=:error)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn"]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest;  kind=:error)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn"]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn"]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull;  kind=:error)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn"]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn"]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn"]))

logging(kind=:error)
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest, :bar)
@test sprint(LogTest.bar) == ""
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull, LogTest)
@test sprint(LogTest.bar) == ""
@test sprint(LogTest.pooh) == ""
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

logging(devnull)
@test sprint(LogTest.bar) == ""
@test sprint(LogTest.pooh) == ""
@test sprint(foo) == ""

logging()
@test all(contains.(sprint(LogTest.bar), ["INFO: barinfo", "WARNING: barwarn", "ERROR: \"barerror\""]))
@test all(contains.(sprint(LogTest.pooh), ["INFO: poohinfo", "WARNING: poohwarn", "ERROR: \"pooherror\""]))
@test all(contains.(sprint(foo), ["INFO: fooinfo", "WARNING: foowarn", "ERROR: \"fooerror\""]))

end # @testset

end

# END 0.7 deprecations

ext: jl

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
auth_href(::Provider{BitbucketAPI}) = "$(ROUTES[:BITBUCKET])/auth"

callback_url(::Provider{BitbucketAPI}) = """$(CONFIG["server_url"])/bitbucket/callback"""

function bitbucket(r::HTTP.Request)
    op = split(subtarget(:BITBUCKET, r), r"[/?]")[1]
    route(Route{BitbucketAPI, Symbol(op)}, r)
end

route(F::Type, service::AbstractString, r::HTTP.Request) = route(Route{F, Symbol(service)}, r)

function route(::Type{Route{F, S}}, r::HTTP.Request) where {F, S}
    HTTP.Response(503; body="Invalid request")
end

# OAuth Step 2: Redirect to provider.
function route(::Type{Route{BitbucketAPI, :auth}}, r::HTTP.Request)
    # If the user has already authenticated, skip.
    state = getcookie(r, "state")
    if !isempty(state) && haskey(USERS, state)
        USERS[state].forge isa BitbucketAPI &&
            return HTTP.Response(307, ["Location" => ROUTES[:SELECT]])
        # continue if the user was logged into a different service
        # authentication will change their service
    end
    provider = PROVIDERS["bitbucket"]
    state = String(rand('a':'z', 32))
    return HTTP.Response(307, [
        "Set-Cookie" => String(HTTP.Cookie("state", state; path="/"), false),
        "Location" => provider.auth_url * "?" * HTTP.escapeuri(Dict(
            :response_type => "code",
            :client_id => provider.client_id,
            :state => state,
        )),
    ])
end

# Step 3: OAuth callback -- sent from the browser via a redirect from Bitbucket
function route(::Type{Route{BitbucketAPI, :callback}}, r::HTTP.Request)
    state = getcookie(r, "state")
    (isempty(state) || state != getquery(r, "state")) && return html(400, "Invalid state")
    provider = PROVIDERS["bitbucket"]
    query = Dict(
        :code => getquery(r, "code"),
        :grant_type => "authorization_code",
    )
    provider.include_state && (query[:state] = state)
    url = URI(URI(provider.token_url); userinfo="$(provider.client_id):$(provider.client_secret)")
    resp = HTTP.post(
        string(url),
        headers=["Accept" => "application/json", "User-Agent" => "Registrator.jl"],
        HTTP.Form(query)
    )
    token = JSON.parse(String(resp.body))["access_token"]
    client = typeof(provider.client)(;
        url=GitForge.base_url(provider.client),
        token=Bitbucket.JWT(token),
        has_rate_limits=GitForge.has_rate_limits(provider.client, identity),
    )
    USERS[state] = User(@gf(get_user(client)), client)

    return HTTP.Response(308, ["Location" => ROUTES[:SELECT]])
end

# /bitbucket/describe -- Describe OAUTH app to the repository hosting service
# this is for things like putting a Julia Registrator panel in your Bitbucket repo page

#function bbroute(::Type{BB{:describe}}, r::HTTP.Request)
#    json(;
#         key = "julia-registrator",
#         name = "Julia Registrator",
#         description = "The Julia package registration service",
#         vendor =
#             (;
#              name = "Julia Computing",
#              url = "https://4609-2a0d-6fc2-4632-7600-a89e-5244-9885-febc.ngrok.io",
#              ),
#         baseUrl = "https://4609-2a0d-6fc2-4632-7600-a89e-5244-9885-febc.ngrok.io",
#         authentication = (; type = "jwt"),
#         lifecycle =
#             (;
#              installed = "/bitbucket/installed",
#              uninstalled = "/bitbucket/uninstalled",
#              ),
#         scopes = ["account", "repository"],
#         contexts = ["account"],
#         modules =
#             (;
#              oauthConsumer = (; clientId="dCmqMQYcdtN2RsSmBLTPttN2CwQv97sk"),
#              repoPages = [
#                  (;
#                   url = "/bitbucket/repopage?repoUuid={repository.uuid}",
#                   name = (; value = "Registrator Panel"),
#                   location = "org.bitbucket.repository.navigation",
#                   key = "registrator-repo-page",
#                   params = (;),
#                   ),
#              ],
#              ),
#         )
#end


text:
using CRC
using Test
using Random
import Libz: crc32

@testset "Consistency between no table and tables" begin
    for _ in 1:10
        data = rand(UInt8, rand(100:200))
        for (name, spec) in ALL
                notable = crc(spec, tables=NoTables)(data)
                singletable = crc(spec, tables=Single)(data)
                multitable = crc(spec, tables=Multiple)(data)
                @test notable == singletable == multitable
        end
    end
end

@testset "Handle string" begin
    crc32 = crc(CRC_32)
    @test crc32("abcxyz") == 0xacc462e9
end

@testset "Handle stream" begin
    crc32 = crc(CRC_32)
    @test crc32(IOBuffer("abcxyz")) == 0xacc462e9
end

@testset "Check against spec: $id" for (id, spec) in ALL
    for tables in (NoTables, Single, Multiple)
        result = crc(spec, tables=tables)(CHECK)
        @test result == spec.check
    end
end

@testset "CRC32 consistency with Libz length $len" for len in (0,1,1000)
    data = rand(UInt8, len)
    check = crc32(data)
    @test crc(CRC_32, tables=Single)(data) == check
    @test crc(CRC_32, tables=Multiple)(data) == check
end


text:
function test_craigmr()
  craigmr_tol = 1.0e-6;

  function test_craigmr(A, b; λ=0.0)
    (nrow, ncol) = size(A);
    (x, y, stats) = craigmr(A, b, λ=λ);
    r = b - A * x;
    Ar = A' * r;
    # if λ > 0
    #   s = r / sqrt(λ);
    #   r = r - sqrt(λ) * s;
    # end
    resid = norm(r) / norm(b);
    Aresid = norm(Ar) / (norm(A) * norm(b));
    @printf("CRAIGMR: residual: %7.1e  least-squares: %7.1e\n", resid, Aresid);
    return (x, y, stats, resid, Aresid)
  end

  # Underdetermined consistent.
  A = rand(10, 25); b = A * ones(25);
  (x, y, stats, resid, Aresid) = test_craigmr(A, b);
  @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  @test(resid <= craigmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x);
  @test(norm(xI - xmin) <= cond(A) * craigmr_tol * xmin_norm);

  # Underdetermined inconsistent.
  A = ones(10, 25); b = rand(10); b[1] = -1.0;
  (x, y, stats, resid) = test_craigmr(A, b);
  @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  @test(stats.inconsistent);
  @test(stats.Aresiduals[end] <= craigmr_tol);

  # Square consistent.
  A = rand(10, 10); b = A * ones(10);
  (x, y, stats, resid) = test_craigmr(A, b);
  @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  @test(resid <= craigmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x);
  @test(norm(xI - xmin) <= cond(A) * craigmr_tol * xmin_norm);

  # Square inconsistent.
  A = ones(10, 10); b = rand(10); b[1] = -1.0;
  (x, y, stats, resid) = test_craigmr(A, b);
  @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  @test(stats.inconsistent);
  @test(stats.Aresiduals[end] <= craigmr_tol);

  # Overdetermined consistent.
  A = rand(25, 10); b = A * ones(10);
  (x, y, stats, resid) = test_craigmr(A, b);
  @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  @test(resid <= craigmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x);
  @test(norm(xI - xmin) <= cond(A) * craigmr_tol * xmin_norm);

  # Overdetermined inconsistent.
  A = ones(5, 3); b = rand(5); b[1] = -1.0;
  (x, y, stats, resid) = test_craigmr(A, b);
  @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  @test(stats.inconsistent);
  @test(stats.Aresiduals[end] <= craigmr_tol);

  # With regularization, all systems are underdetermined and consistent.
  # (x, y, stats, resid) = test_craigmr(A, b, λ=1.0e-3);
  # @test(norm(x - A' * y) <= craigmr_tol * norm(x));
  # @test(resid <= craigmr_tol);
  # @test(stats.solved);
  # (xI, xmin, xmin_norm) = check_min_norm(A, b, x, λ=1.0e-3);
  # @test(norm(xI - xmin) <= cond(A) * craigmr_tol * xmin_norm);

  # Code coverage.
  (x, y, stats) = craigmr(sparse(A), b, λ=1.0e-3);
  show(stats);

  # Test b == 0
  (x, y, stats) = craigmr(A, zeros(size(A,1)), λ=1.0e-3)
  @test x == zeros(size(A,2))
  @test y == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [I; rand(1:10, 2, 3)]
  b = A * ones(Int, 3)
  (x, y, stats) = craigmr(A, b)
  @test stats.solved
end

test_craigmr()


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using Base.Test

mutable struct Dummy <: FEMBase.FieldProblem
end

@testset "add elements to problem" begin
    problem = Problem(Dummy, "test", 2)
    element = Element(Quad4, [1, 2, 3, 4])
    elements = [element]
    add_elements!(problem, elements)
    @test problem.elements[1] == element
end


text:
using ShuffleProofs
import ShuffleProofs: PrimeGenerator, Enc, Dec, shuffle, ElGamal

sk = 123 # Only authorithy who do decryption would know

g = PrimeGenerator(3, 23)
options = [g, g^2, g^3]
pk = g^sk

bbord = (;) # A representation of a buletinboard

######## Step 1. ########

bbord = let
    enc = Enc(pk, g)
    ciphertexts_in = ElGamal([enc(options[rand(1:3)], rand(1:10)) for i in 1:10])
    (; bbord..., ciphertexts_in)
end

######## Step 2. ########

bbord = let
    enc = Enc(pk, g)
    proposition, secret = shuffle(bbord.ciphertexts_in, enc) # reencrypts and sorts output
    ciphertexts_out = proposition.𝐞′
    (; bbord..., ciphertexts_out)
end

####### Step 3. ########

bbord = let
    dec = Dec(sk)
    votes = dec(bbord.ciphertexts_out)
    (; bbord..., votes)
end

######### As last step we can tally the votes ############

for (i, o) in enumerate(options)
    s = 0
    for v in bbord.votes
        if v == o
            s += 1
        end
    end

    println("Choice $i received $s votes")
end


text:
using CImGuiExtensions
using GR

function launch_video()

    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "Demo")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsDark()

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)

    stream_webcam = false
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
    show_another_window = true
    counter = 0

    # setup camera
    f = VideoIO.openvideo(video_file)
    # The OpenGL library expects RGBA UInt8 data layed out in a width x height format.
    texture₀ = ImGui_ImplOpenGL3_CreateImageTexture(f.width, f.height, format = GL_RGB)

    # Capture the first frame so that we can initialize a buffer to store each
    # frame that is read from the webcam.
    imageₙ = read(f)

    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        CImGui.Begin("Image Demo")
        @c CImGui.Checkbox("Plot Data", &stream_webcam)

        if stream_webcam
            # consume the next camera frame
            !eof(f) && read!(f, imageₙ)
            imageₙ′ = unsafe_wrap(Array{UInt8,3}, convert(Ptr{UInt8}, pointer(imageₙ)), (Cint(3), f.width, f.height))
            ImGui_ImplOpenGL3_UpdateImageTexture(texture₀ , imageₙ′, f.width, f.height; format = GL_RGB)
            CImGui.Image(Ptr{Cvoid}(texture₀), (f.width, f.height))
        end

        CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())
        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
    close(f)

end


text:
## Banded LDLᵀ decomposition

@inline bandwidth(A::LDLt{T,Symmetric{T,M}}, k) where {T,M<:BandedMatrix{T}} = bandwidth(A.data, k)

factorize(S::Symmetric{<:Any,<:BandedMatrix}) = ldlt(S)
if VERSION ≥ v"1.4-"
    LinearAlgebra._factorize(S::Symmetric{<:Any,<:BandedMatrix}; check=false) = ldlt(S)
end

function ldlt!(F::Symmetric{T,M}) where {T,M<:BandedMatrix{T}}
    A = F.data
    n = size(A, 1)
    b = bandwidth(F)
    if F.uplo == 'L'
        @inbounds for j = 1:n
            @simd for k = max(1,j-b-1):j-1
                A[j,j] -= abs2(A[j,k])*A[k,k]
            end
            for i = j+1:min(n,j+b)
                @simd for k = max(1,j-b-1):j-1
                    A[i,j] -= A[i,k]*A[j,k]*A[k,k]
                end
                A[i,j] /= A[j,j]
            end
        end
    elseif F.uplo == 'U'
        @inbounds for j = 1:n
            @simd for k = max(1,j-b-1):j-1
                A[j,j] -= abs2(A[k,j])*A[k,k]
            end
            for i = j+1:min(n,j+b)
                @simd for k = max(1,j-b-1):j-1
                    A[j,i] -= A[k,i]*A[k,j]*A[k,k]
                end
                A[j,i] /= A[j,j]
            end
        end
    end
    return LDLt{T,Symmetric{T,M}}(F)
end

function ldlt(A::Symmetric{T,M}) where {T,M<:BandedMatrix{T}}
    S = typeof(zero(T)/one(T))
    uplo = A.uplo == 'U' ? :U : :L
    return S == T ? ldlt!(copy(A)) : ldlt!(Symmetric(BandedMatrix{S}(A), uplo))
end

function ldiv!(S::LDLt{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
    require_one_based_indexing(B)
    n, nrhs = size(B, 1), size(B, 2)
    if size(S, 1) != n
        throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))
    end
    A = S.data
    b = bandwidth(A, 1)
    @inbounds for l = 1:nrhs
        for j = 1:n
            @simd for k = max(1,j-b-1):j-1
                B[j,l] -= A[j,k]*B[k,l]
            end
        end
        @simd for j = 1:n
            B[j,l] /= A[j,j]
        end
        for j = n:-1:1
            @simd for k = j+1:min(n,j+b)
                B[j,l] -= A[k,j]*B[k,l]
            end
        end
    end
    B
end

function logabsdet(F::LDLt{T,Symmetric{T,M}}) where {T,M<:BandedMatrix{T}}
    it = (F.data[i,i] for i in 1:size(F, 1))
    return sum(log∘abs, it), prod(sign, it)
end

function getproperty(F::LDLt{<:Any,<:Symmetric{<:Any,<:AbstractBandedMatrix}}, d::Symbol)
    Fdata = getfield(F, :data)
    if d === :d
        return diag(Fdata)
    elseif d === :D
        return Diagonal(diag(Fdata))
    elseif d === :L
        return UnitLowerTriangular(Fdata)
    elseif d === :Lt
        return UnitUpperTriangular(Fdata)
    else
        return getfield(F, d)
    end
end


text:
include("gemm.jl")
include("symm.jl")
include("hemm.jl")
include("syrk.jl")
include("herk.jl")
include("syr2k.jl")
include("her2k.jl")
include("trmm.jl")
include("trsm.jl")


text:
function zero!(x::Mag)
    ccall(@libarb(mag_zero), Cvoid, (Ref{Mag},), x)
    return x
end

function one!(x::Mag)
    ccall(@libarb(mag_one), Cvoid, (Ref{Mag},), x)
    return x
end

function inf!(x::Mag)
    ccall(@libarb(mag_inf), Cvoid, (Ref{Mag},), x)
    return x
end

function zero!(x::ArbFloat{P}) where {P}
    ccall(@libarb(arf_zero), Cvoid, (Ref{ArbFloat},), x)
    return x
end

function one!(x::ArbFloat{P}) where {P}
    ccall(@libarb(arf_one), Cvoid, (Ref{ArbFloat},), x)
    return x
end

function inf!(x::ArbFloat{P}) where {P}
    ccall(@libarb(arf_pos_inf), Cvoid, (Ref{ArbFloat},), x)
    return x
end

function posinf!(x::ArbFloat{P}) where {P}
    ccall(@libarb(arf_pos_inf), Cvoid, (Ref{ArbFloat},), x)
    return x
end

function neginf!(x::ArbFloat{P}) where {P}
    ccall(@libarb(arf_neg_inf), Cvoid, (Ref{ArbFloat},), x)
    return x
end

function nan!(x::ArbFloat{P}) where {P}
    ccall(@libarb(arf_nan), Cvoid, (Ref{ArbFloat},), x)
    return x
end

function zero!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_zero), Cvoid, (Ref{ArbReal},), x)
    return x
end

function one!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_one), Cvoid, (Ref{ArbReal},), x)
    return x
end

function inf!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_pos_inf), Cvoid, (Ref{ArbReal},), x)
    return x
end

function posinf!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_pos_inf), Cvoid, (Ref{ArbReal},), x)
    return x
end

function neginf!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_neg_inf), Cvoid, (Ref{ArbReal},), x)
    return x
end

function pminf!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_zero_pm_inf), Cvoid, (Ref{ArbReal},), x)
    return x
end

function indeterminate!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_indeterminate), Cvoid, (Ref{ArbReal},), x)
    return x
end

function nan!(x::ArbReal{P}) where {P}
    ccall(@libarb(arb_nan), Cvoid, (Ref{ArbReal},), x)
    return x
end

function zero!(x::ArbComplex{P}) where {P}
    ccall(@libarb(acb_zero), Cvoid, (Ref{ArbComplex},), x)
    return x
end

function one!(x::ArbComplex{P}) where {P}
    ccall(@libarb(acb_one), Cvoid, (Ref{ArbComplex},), x)
    return x
end

function inf!(x::ArbComplex{P}) where {P}
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    rea = inf!(rea)
    ima = inf!(ima)
    z = ArbComplex{P}(rea, ima)
    return z
end

function posinf!(x::ArbComplex{P}) where {P}
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    rea = posinf!(rea)
    ima = zero!(ima)
    z = ArbComplex{P}(rea, ima)
    return z
end

function neginf!(x::ArbComplex{P}) where {P}
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    rea = neginf!(rea)
    ima = zero!(ima)
    z = ArbComplex{P}(rea, ima)
    return z
end

function pminf!(x::ArbComplex{P}) where {P}
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    rea = pminf!(rea)
    ima = pminf!(ima)
    z = ArbComplex{P}(rea, ima)
    return z
end

function indeterminate!(x::ArbComplex{P}) where {P}
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    rea = indeterminate!(rea)
    ima = indeterminate!(ima)
    z = ArbComplex{P}(rea, ima)
    return z
end

function nan!(x::ArbComplex{P}) where {P}
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    rea = nan!(rea)
    ima = nan!(ima)
    z = ArbComplex{P}(rea, ima)
    return z
end

zero(::Type{Mag}) = Mag()
zero(::Type{ArbFloat{P}}) where {P} = ArbFloat{P}()
zero(::Type{ArbReal{P}}) where {P} = ArbReal{P}()
zero(::Type{ArbComplex{P}}) where {P} = ArbComplex{P}()
zero(x::Mag) = zero(typeof(x))
zero(x::ArbFloat{P}) where {P} = zero(typeof(x))
zero(x::ArbReal{P}) where {P} = zero(typeof(x))
zero(x::ArbComplex{P}) where {P} = zero(typeof(x))
zero(::Type{ArbFloat}) = zero(ArbFloat{DEFAULT_PRECISION.x})
zero(::Type{ArbReal}) = zero(ArbReal{DEFAULT_PRECISION.x})
zero(::Type{ArbComplex}) = zero(ArbComplex{DEFAULT_PRECISION.x})

function one(::Type{Mag})
    z = Mag()
    one!(z)
    return z
end
function one(::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    one!(z)
    return z
end
function one(::Type{ArbReal{P}}) where {P}
    z = ArbReal{P}()
    one!(z)
    return z
end
function one(::Type{ArbComplex{P}}) where {P}
    z = ArbComplex{P}()
    one!(z)
    return z
end

one(x::Mag) = zero(typeof(x))
one(x::ArbFloat{P}) where {P} = one(typeof(x))
one(x::ArbReal{P}) where {P} = one(typeof(x))
one(x::ArbComplex{P}) where {P} = one(typeof(x))
one(::Type{ArbFloat}) = one(ArbFloat{DEFAULT_PRECISION.x})
one(::Type{ArbReal}) = one(ArbReal{DEFAULT_PRECISION.x})
one(::Type{ArbComplex}) = one(ArbComplex{DEFAULT_PRECISION.x})

function inf(::Type{Mag})
    z = Mag()
    inf!(z)
    return z
end
function inf(::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    inf!(z)
    return z
end
function inf(::Type{ArbReal{P}}) where {P}
    z = ArbReal{P}()
    inf!(z)
    return z
end
function inf(::Type{ArbComplex{P}}) where {P}
    z = ArbComplex{P}()
    inf!(z)
    return z
end

inf(x::Mag) = inf(typeof(x))
inf(x::ArbFloat{P}) where {P} = inf(typeof(x))
inf(x::ArbReal{P}) where {P} = inf(typeof(x))
inf(x::ArbComplex{P}) where {P} = inf(typeof(x))
inf(::Type{ArbFloat}) = inf(ArbFloat{DEFAULT_PRECISION.x})
inf(::Type{ArbReal}) = inf(ArbReal{DEFAULT_PRECISION.x})
inf(::Type{ArbComplex}) = inf(ArbComplex{DEFAULT_PRECISION.x})

posinf(::Type{Mag}) = throw(DomainError("use `inf(Mag)`."))
function posinf(::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    posinf!(z)
    return z
end
function posinf(::Type{ArbReal{P}}) where {P}
    z = ArbReal{P}()
    posinf!(z)
    return z
end
function posinf(::Type{ArbComplex{P}}) where {P}
    z = ArbComplex{P}()
    posinf!(z)
    return z
end

posinf(x::Mag) = posinf(typeof(x))
posinf(x::ArbFloat{P}) where {P} = posinf(typeof(x))
posinf(x::ArbReal{P}) where {P} = posinf(typeof(x))
posinf(x::ArbComplex{P}) where {P} = posinf(typeof(x))
posinf(::Type{ArbFloat}) = posinf(ArbFloat{DEFAULT_PRECISION.x})
posinf(::Type{ArbReal}) = posinf(ArbReal{DEFAULT_PRECISION.x})
posinf(::Type{ArbComplex}) = posinf(ArbComplex{DEFAULT_PRECISION.x})

neginf(::Type{Mag}) = throw(DomainError("use `inf(Mag)`."))
function neginf(::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    neginf!(z)
    return z
end
function neginf(::Type{ArbReal{P}}) where {P}
    z = ArbReal{P}()
    neginf!(z)
    return z
end
function neginf(::Type{ArbComplex{P}}) where {P}
    z = ArbComplex{P}()
    neginf!(z)
    return z
end

neginf(x::Mag) = neginf(typeof(x))
neginf(x::ArbFloat{P}) where {P} = neginf(typeof(x))
neginf(x::ArbReal{P}) where {P} = neginf(typeof(x))
neginf(x::ArbComplex{P}) where {P} = neginf(typeof(x))
neginf(::Type{ArbFloat}) = neginf(ArbFloat{DEFAULT_PRECISION.x})
neginf(::Type{ArbReal}) = neginf(ArbReal{DEFAULT_PRECISION.x})
neginf(::Type{ArbComplex}) = neginf(ArbComplex{DEFAULT_PRECISION.x})

allextendedreals(::Type{ArbReal{P}}) where {P} = pminf!(zero(ArbReal{P}))
allextendedreals(::Type{ArbComplex{P}}) where {P} = ArbComplex{P}(allextendedreals(ArbReal{P}), zero(ArbReal{P}))
allextendedimags(::Type{ArbReal{P}}) where {P} = ArbComplex{P}(zero(ArbReal{P}), allextendedreals(ArbReal{P}))
allextendedcomplex(::Type{ArbComplex{P}}) where {P} = ArbComplex{P}(allextendedreals(ArbReal{P}), allextendedreals(ArbReal{P}))
allextendedreals(::Type{ArbReal}) = allextendedreals(ArbReal{DEFAULT_PRECISION.x})
allextendedreals(::Type{ArbComplex}) = ArbComplex{DEFAULT_PRECISION.x}(allextendedreals(ArbReal), zero(ArbReal))
allextendedimags(::Type{ArbComplex}) = ArbComplex{DEFAULT_PRECISION.x}(zero(ArbReal), allextendedreals(ArbReal))

NaN(::Type{Mag}) = throw(DomainError("nan(Mag) does not exist"))
function NaN(::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    nan!(z)
    return z
end
function NaN(::Type{ArbReal{P}}) where {P}
    z = ArbReal{P}()
    nan!(z)
    return z
end
NaN(::Type{ArbComplex{P}}) where {P} = ArbComplex{P}(NaN(ArbReal{P}), NaN(ArbReal{P}))

NaN(x::Mag) = NaN(typeof(x))
NaN(x::ArbFloat{P}) where {P} = NaN(typeof(x))
NaN(x::ArbReal{P}) where {P} = NaN(typeof(x))
NaN(x::ArbComplex{P}) where {P} = NaN(typeof(x))
NaN(::Type{ArbFloat}) = NaN(ArbFloat{DEFAULT_PRECISION.x})
NaN(::Type{ArbReal}) = NaN(ArbReal{DEFAULT_PRECISION.x})
NaN(::Type{ArbComplex}) = NaN(ArbComplex{DEFAULT_PRECISION.x})

function isspecial(x::Mag)
    0 != ccall(@libarb(mag_is_special), Cint, (Ref{Mag},), x)
end

function isnormal(x::Mag)
    0 == ccall(@libarb(mag_is_special), Cint, (Ref{Mag},), x)
end

function isfinite(x::Mag)
    0 != ccall(@libarb(mag_is_finite), Cint, (Ref{Mag},), x)
end

function isnonfinite(x::Mag)
    0 == ccall(@libarb(mag_is_finite), Cint, (Ref{Mag},), x)
end

function iszero(x::Mag)
    0 != ccall(@libarb(mag_is_zero), Cint, (Ref{Mag},), x)
end

function isnonzero(x::Mag)
    0 == ccall(@libarb(mag_is_zero), Cint, (Ref{Mag},), x)
end

function isinf(x::Mag)
    0 != ccall(@libarb(mag_is_inf), Cint, (Ref{Mag},), x)
end

function isposinf(x::Mag)
    0 != ccall(@libarb(mag_is_inf), Cint, (Ref{Mag},), x)
end

function isneginf(x::Mag)
    false
end

function isspecial(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_special), Cint, (Ref{ArbFloat},), x)
end

function isnormal(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_normal), Cint, (Ref{ArbFloat},), x)
end

function isfinite(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_finite), Cint, (Ref{ArbFloat},), x)
end

function isnonfinite(x::ArbFloat{P}) where {P}
    0 == ccall(@libarb(arf_is_finite), Cint, (Ref{ArbFloat},), x)
end

function iszero(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_zero), Cint, (Ref{ArbFloat},), x)
end

function isone(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_one), Cint, (Ref{ArbFloat},), x)
end

function isnonzero(x::ArbFloat{P}) where {P}
    0 == ccall(@libarb(arf_is_zero), Cint, (Ref{ArbFloat},), x)
end

function isinf(x::ArbFloat{P}) where {P}
    z = abs(x)
    0 != ccall(@libarb(arf_is_inf), Cint, (Ref{ArbFloat},), x)
end

function isposinf(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_pos_inf), Cint, (Ref{ArbFloat},), x)
end

function isneginf(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_neg_inf), Cint, (Ref{ArbFloat},), x)
end

function isnan(x::ArbFloat{P}) where {P}
    0 != ccall(@libarb(arf_is_nan), Cint, (Ref{ArbFloat},), x)
end

function isspecial(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_special), Cint, (Ref{ArbReal},), x)
end

function isnormal(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_normal), Cint, (Ref{ArbReal},), x)
end

function isfinite(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_finite), Cint, (Ref{ArbReal},), x)
end

function isnonfinite(x::ArbReal{P}) where {P}
    0 == ccall(@libarb(arb_is_finite), Cint, (Ref{ArbReal},), x)
end

function iszero(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_zero), Cint, (Ref{ArbReal},), x)
end

function isone(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_one), Cint, (Ref{ArbReal},), x)
end

function isnonzero(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_nonzero), Cint, (Ref{ArbReal},), x)
end

function isexact(x::ArbReal{P}) where {P}
    0 != ccall(@libarb(arb_is_exact), Cint, (Ref{ArbReal},), x)
end

function isinexact(x::ArbReal{P}) where {P}
    0 == ccall(@libarb(arb_is_exact), Cint, (Ref{ArbReal},), x)
end

function isinf(x::ArbReal{P}) where {P}
    isinf(ArbFloat{P}(x))
end

function isposinf(x::ArbReal{P}) where {P}
    isposinf(ArbFloat{P}(x))
end

function isneginf(x::ArbReal{P}) where {P}
    isneginf(ArbFloat{P}(x))
end

function isnan(x::ArbReal{P}) where {P}
    isnan(ArbFloat{P}(x))
end

function iszero(x::ArbComplex{P}) where {P}
    0 != ccall(@libarb(acb_is_zero), Cint, (Ref{ArbComplex},), x)
end

function isnonzero(x::ArbComplex{P}) where {P}
    0 == ccall(@libarb(acb_is_zero), Cint, (Ref{ArbComplex},), x)
end

function isone(x::ArbComplex{P}) where {P}
    0 != ccall(@libarb(acb_is_one), Cint, (Ref{ArbComplex},), x)
end

function isexact(x::ArbComplex{P}) where {P}
    0 != ccall(@libarb(acb_is_exact), Cint, (Ref{ArbComplex},), x)
end

function isinexact(x::ArbComplex{P}) where {P}
    0 == ccall(@libarb(acb_is_exact), Cint, (Ref{ArbComplex},), x)
end

function isfinite(x::ArbComplex{P}) where {P}
    0 != ccall(@libarb(arb_is_finite), Cint, (Ref{ArbComplex},), x)
end

function isnonfinite(x::ArbComplex{P}) where {P}
    0 == ccall(@libarb(arb_is_finite), Cint, (Ref{ArbComplex},), x)
end

function isinf(x::ArbComplex{P}) where {P}
    !isfinite(x) && !isnan(x)
end

function isposinf(x::ArbComplex{P}) where {P}
    isposinf(real(x))
end

function isneginf(x::ArbComplex{P}) where {P}
    isneginf(real(x))
end

function isnan(x::ArbComplex{P}) where {P}
    isnan(real(x)) || isnan(imag(x))
end

typemax(::Type{ArbFloat{P}}) where {P} = posinf(ArbFloat{P})
typemin(::Type{ArbFloat{P}}) where {P} = neginf(ArbFloat{P})

floatmax(::Type{ArbFloat{P}}) where {P} = pow(ArbFloat{P}(P), square(ArbFloat{P}(P)))
floatmin(::Type{ArbFloat{P}}) where {P} = inv(floatmax(ArbFloat{P}))

typemax(::Type{ArbReal{P}}) where {P} = posinf(ArbReal{P})
typemin(::Type{ArbReal{P}}) where {P} = neginf(ArbReal{P})

floatmax(::Type{ArbReal{P}}) where {P} = ArbReal{P}(floatmax(ArbFloat{P}))
floatmin(::Type{ArbReal{P}}) where {P} = inv(floatmax(x))


text:
@testset "α Filter" begin

	@testset "1D α Filter" begin
		# Input Data
		V = [ 1. 3. 2. 6. 7. 8. 12. 10. ]
		DT = AlphaShapes.delaunayTriangulation(V)

		# Expected Output
		VV = [[1],[2],[3],[4],[5],[6],[7],[8]]
		EV = [[1,3],[2,3],[2,4],[4,5],[5,6],[6,8],[7,8]]

		# Evaluation
		filter = AlphaShapes.alphaFilter(V, DT)

		@test length(unique(values(filter))) == 4
		@test sort(unique(values(filter))) == [0.0, 0.5, 1.0, 1.5]
		@test length(unique(keys(filter))) == 15
		@test sort([v for v in unique(keys(filter)) if length(v) == 2]) == EV
	end

	@testset "2D α Filter" begin
		# Input Data
		V = [
			0.0 1.0 0.0 2.0;
			0.0 0.0 1.0 2.0
		]

		# Expected Output
		VV = [[1],[2],[3],[4]]
		EV = [[1, 2],[1, 3],[2, 3],[2, 4],[3, 4]]
		FV = [[1, 2, 3],[2, 3, 4]]

		# Evaluation
		filter = AlphaShapes.alphaFilter(V)

		@test length(unique(values(filter))) == 5
		@test isapprox(sort(unique(values(filter))),
			[0.0, 0.5, 0.7071, 1.1180, 1.1785], atol=1e-4
		)
		@test length(unique(keys(filter))) == 11
		@test sort([v for v in unique(keys(filter)) if length(v) == 1]) == VV
		@test sort([v for v in unique(keys(filter)) if length(v) == 2]) == EV
		@test sort([v for v in unique(keys(filter)) if length(v) == 3]) == FV
	end

	# @testset "3D α Filter" begin
	# 	V = [
	# 		0.0 1.0 0.0 0.0 1.0 0.0;
	# 		0.0 0.0 1.0 0.0 0.0 1.0;
	# 		0.0 0.0 0.0 1.0 1.0 2.0
	# 	]
	#
	# 	# Expected Output
	# 	VV = [[1],[2],[3],[4],[5],[6]]
	# 	CV = [[1, 2, 3, 4], [2, 3, 4, 5], [3, 4, 5, 6]]
	# 	# Evaluation
	# 	filter = AlphaShapes.alphaFilter(V, digits=4)
	#
	# 	@test length(unique(values(filter))) == 6
	# 	@test sort(unique(values(filter))) == [0.0, 0.5, 0.7071, 0.866, 1.0, 1.118]
	# 	@test length(unique(keys(filter))) == 31
	# 	@test sort([v for v in unique(keys(filter)) if length(v) == 1]) == VV
	# 	@test sort([v for v in unique(keys(filter)) if length(v) == 4]) == CV
	# end

end

@testset "α Simplices Evaluation" begin

	@testset "1D α Simplex" begin
		# Input Data
		V = [
			0.0 1.0 2.0 7.0 5.0 3.0 14.0 11.0
		]
		filtration = AlphaShapes.alphaFilter(V)

		# α = 0.0
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 0.0)
		@test α_simplices[1] == [[1],[2],[3],[4],[5],[6],[7],[8]]
		@test isempty(α_simplices[2])
		# α = 0.5
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 0.5)
		@test α_simplices[2] == [[1,2],[2,3],[3,6]]
		# α = 1.0
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 1.0)
		@test α_simplices[2] == [[1, 2], [2, 3], [3, 6], [4, 5], [5, 6]]
		# α = 1.5
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 1.5)
		@test α_simplices[2] == [[1, 2], [2, 3], [3, 6], [4, 5], [5, 6], [7, 8]]

	end

	@testset "2D α Simplex" begin
		# Input Data
		V = [
			0.0 1.0 0.0 2.0;
			0.0 0.0 1.0 0.3
		]
		filtration = AlphaShapes.alphaFilter(V)

		# α = 0.0
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 0.0)
		@test α_simplices[1] == [[1],[2],[3],[4]]
		@test isempty(α_simplices[2])
		@test isempty(α_simplices[3])
		# α = 0.5
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 0.5)
		@test α_simplices[2] == [[1,2],[1,3]]
		@test isempty(α_simplices[3])
		# α = 1.0
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 1.0)
		@test α_simplices[2] == [[1,2],[1,3],[2,3],[2,4]]
		@test α_simplices[3] == [[1,2,3]]
		# α = 1.5
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 1.5)
		@test α_simplices[2] == [[1,2],[1,3],[2,3],[2,4],[3,4]]
		@test α_simplices[3] == [[1,2,3],[2,3,4]]

	end

	@testset "3D α Simplex" begin
		# Input Data
		V = [
			0.0  1.0  0.0  0.0  2.0;
			0.0  0.0  1.0  0.0  2.0;
			0.0  0.0  0.0  1.0  2.0
		]
		filtration = AlphaShapes.alphaFilter(V)

		# α = 0.0
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 0.0)
		@test α_simplices[1] == [[1],[2],[3],[4],[5]]
		@test isempty(α_simplices[2])
		@test isempty(α_simplices[3])
		@test isempty(α_simplices[4])
		# α = 0.5
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 0.5)
		@test α_simplices[2] == [[1,2],[1,3],[1,4]]
		@test isempty(α_simplices[3])
		@test isempty(α_simplices[4])
		# α = 1.0
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 1.0)
		@test α_simplices[2] == [[1,2],[1,3],[1,4],[2,3],[2,4],[3,4]]
		@test α_simplices[3] == [[1,2,3],[1,2,4],[1,3,4],[2,3,4]]
		@test α_simplices[4] == [[1,2,3,4]]

		# α = 1.5
		α_simplices = AlphaShapes.alphaSimplex(V, filtration, 1.6)
		@test α_simplices[4] == [[1,2,3,4],[2,3,4,5]]
	end

end


text:
module DateDataFrames

using Base: padding
export DateDataFrame, resample_nearest, resample, Nearest,
        timestamp, values, names, rename, rename!, lead, moving,
        ⊕, ⊖, ⊗, ⊘

using DataFrames
import DataFrames: names, rename, rename!
using Dates
using Statistics
using RecipesBase

include("date_dataframes.jl")
include("broadcasting.jl")
include("resample.jl")
include("misc.jl")
include("infix_ops.jl")
include("recipes.jl")

end


text:
# # Tutorial 5
#
# ## Mesh sizes, holes in volumes
#
import gmsh
#
gmsh.initialize()
#
gmsh.model.add("t5")
#
lcar1 = .1
lcar2 = .0005
lcar3 = .055
# If we wanted to change these mesh sizes globally (without changing the above
# definitions), we could give a global scaling factor for all mesh sizes with
# e.g.
#
# `gmsh.option.setNumber("Mesh.MeshSizeFactor", 0.1);`
# See Tutorial 10 for more information about mesh sizes.
#
# We proceed by defining some elementary entities describing a truncated cube:
gmsh.model.geo.addPoint(0.5,0.5,0.5, lcar2, 1)
gmsh.model.geo.addPoint(0.5,0.5,0, lcar1, 2)
gmsh.model.geo.addPoint(0,0.5,0.5, lcar1, 3)
gmsh.model.geo.addPoint(0,0,0.5, lcar1, 4)
gmsh.model.geo.addPoint(0.5,0,0.5, lcar1, 5)
gmsh.model.geo.addPoint(0.5,0,0, lcar1, 6)
gmsh.model.geo.addPoint(0,0.5,0, lcar1, 7)
gmsh.model.geo.addPoint(0,1,0, lcar1, 8)
gmsh.model.geo.addPoint(1,1,0, lcar1, 9)
gmsh.model.geo.addPoint(0,0,1, lcar1, 10)
gmsh.model.geo.addPoint(0,1,1, lcar1, 11)
gmsh.model.geo.addPoint(1,1,1, lcar1, 12)
gmsh.model.geo.addPoint(1,0,1, lcar1, 13)
gmsh.model.geo.addPoint(1,0,0, lcar1, 14)
#
gmsh.model.geo.addLine(8,9, 1);   gmsh.model.geo.addLine(9,12, 2)
gmsh.model.geo.addLine(12,11, 3); gmsh.model.geo.addLine(11,8, 4)
gmsh.model.geo.addLine(9,14, 5);  gmsh.model.geo.addLine(14,13, 6)
gmsh.model.geo.addLine(13,12, 7); gmsh.model.geo.addLine(11,10, 8)
gmsh.model.geo.addLine(10,13, 9); gmsh.model.geo.addLine(10,4, 10)
gmsh.model.geo.addLine(4,5, 11);  gmsh.model.geo.addLine(5,6, 12)
gmsh.model.geo.addLine(6,2, 13);  gmsh.model.geo.addLine(2,1, 14)
gmsh.model.geo.addLine(1,3, 15);  gmsh.model.geo.addLine(3,7, 16)
gmsh.model.geo.addLine(7,2, 17);  gmsh.model.geo.addLine(3,4, 18)
gmsh.model.geo.addLine(5,1, 19);  gmsh.model.geo.addLine(7,8, 20)
gmsh.model.geo.addLine(6,14, 21);
#
gmsh.model.geo.addCurveLoop([-11,-19,-15,-18], 22)
gmsh.model.geo.addPlaneSurface([22], 23)
gmsh.model.geo.addCurveLoop([16,17,14,15], 24)
gmsh.model.geo.addPlaneSurface([24], 25)
gmsh.model.geo.addCurveLoop([-17,20,1,5,-21,13], 26)
gmsh.model.geo.addPlaneSurface([26], 27)
gmsh.model.geo.addCurveLoop([-4,-1,-2,-3], 28)
gmsh.model.geo.addPlaneSurface([28], 29)
gmsh.model.geo.addCurveLoop([-7,2,-5,-6], 30)
gmsh.model.geo.addPlaneSurface([30], 31)
gmsh.model.geo.addCurveLoop([6,-9,10,11,12,21], 32)
gmsh.model.geo.addPlaneSurface([32], 33)
gmsh.model.geo.addCurveLoop([7,3,8,9], 34)
gmsh.model.geo.addPlaneSurface([34], 35)
gmsh.model.geo.addCurveLoop([-10,18,-16,-20,4,-8], 36)
gmsh.model.geo.addPlaneSurface([36], 37)
gmsh.model.geo.addCurveLoop([-14,-13,-12,19], 38)
gmsh.model.geo.addPlaneSurface([38], 39)
#
shells = []
#
sl = gmsh.model.geo.addSurfaceLoop([35,31,29,37,33,23,39,25,27])
push!(shells, sl)
#
# This function will create a spherical hole in a volume. We don't specify
# tags manually, and let the functions return them automatically:
function cheeseHole(x, y, z, r, lc, shells)
    p1 = gmsh.model.geo.addPoint(x,  y,  z,   lc)
    p2 = gmsh.model.geo.addPoint(x+r,y,  z,   lc)
    p3 = gmsh.model.geo.addPoint(x,  y+r,z,   lc)
    p4 = gmsh.model.geo.addPoint(x,  y,  z+r, lc)
    p5 = gmsh.model.geo.addPoint(x-r,y,  z,   lc)
    p6 = gmsh.model.geo.addPoint(x,  y-r,z,   lc)
    p7 = gmsh.model.geo.addPoint(x,  y,  z-r, lc)

    c1 = gmsh.model.geo.addCircleArc(p2,p1,p7)
    c2 = gmsh.model.geo.addCircleArc(p7,p1,p5)
    c3 = gmsh.model.geo.addCircleArc(p5,p1,p4)
    c4 = gmsh.model.geo.addCircleArc(p4,p1,p2)
    c5 = gmsh.model.geo.addCircleArc(p2,p1,p3)
    c6 = gmsh.model.geo.addCircleArc(p3,p1,p5)
    c7 = gmsh.model.geo.addCircleArc(p5,p1,p6)
    c8 = gmsh.model.geo.addCircleArc(p6,p1,p2)
    c9 = gmsh.model.geo.addCircleArc(p7,p1,p3)
    c10 = gmsh.model.geo.addCircleArc(p3,p1,p4)
    c11 = gmsh.model.geo.addCircleArc(p4,p1,p6)
    c12 = gmsh.model.geo.addCircleArc(p6,p1,p7)

    l1 = gmsh.model.geo.addCurveLoop([c5,c10,c4])
    l2 = gmsh.model.geo.addCurveLoop([c9,-c5,c1])
    l3 = gmsh.model.geo.addCurveLoop([c12,-c8,-c1])
    l4 = gmsh.model.geo.addCurveLoop([c8,-c4,c11])
    l5 = gmsh.model.geo.addCurveLoop([-c10,c6,c3])
    l6 = gmsh.model.geo.addCurveLoop([-c11,-c3,c7])
    l7 = gmsh.model.geo.addCurveLoop([-c2,-c7,-c12])
    l8 = gmsh.model.geo.addCurveLoop([-c6,-c9,c2])

    s1 = gmsh.model.geo.addSurfaceFilling([l1])
    s2 = gmsh.model.geo.addSurfaceFilling([l2])
    s3 = gmsh.model.geo.addSurfaceFilling([l3])
    s4 = gmsh.model.geo.addSurfaceFilling([l4])
    s5 = gmsh.model.geo.addSurfaceFilling([l5])
    s6 = gmsh.model.geo.addSurfaceFilling([l6])
    s7 = gmsh.model.geo.addSurfaceFilling([l7])
    s8 = gmsh.model.geo.addSurfaceFilling([l8])

    sl = gmsh.model.geo.addSurfaceLoop([s1, s2, s3, s4, s5, s6, s7, s8])
    v = gmsh.model.geo.addVolume([sl])
    push!(shells, sl)
    return v
end
#
# We create five holes in the cube:
x = 0
y = 0.75; z = 0; r = 0.09
for t in 1:5
    global x, z
    x += 0.166
    z += 0.166
    v = cheeseHole(x, y, z, r, lcar3, shells)
    gmsh.model.geo.addPhysicalGroup(3, [v], t)
end
#
# The volume of the cube, without the 5 holes, is defined by 6 surface loops:
# the first surface loop defines the exterior surface; the surface loops other
# than the first one define holes:
gmsh.model.geo.addVolume(shells, 186);
#
gmsh.model.geo.synchronize()
#
# We finally define a physical volume for the elements discretizing the cube,
# without the holes (for which physical groups were already defined in the
# `cheeseHole()` function):
gmsh.model.addPhysicalGroup(3, [186], 10);
#
# To generate a curvilinear mesh and optimize it to produce provably valid
# curved elements (see A. Johnen, J.-F. Remacle and C. Geuzaine. Geometric
# validity of curvilinear finite elements. Journal of Computational Physics
# 233, pp. 359-372, 2013; and T. Toulorge, C. Geuzaine, J.-F. Remacle,
# J. Lambrechts. Robust untangling of curvilinear meshes. Journal of
# Computational Physics 254, pp. 8-26, 2013), you can uncomment the following
# lines:
#
# `gmsh.option.setNumber("Mesh.ElementOrder", 2)`
# `gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)`
gmsh.model.mesh.generate(3)
gmsh.write("t5.msh")
#
gmsh.finalize()


text:
heapleft(i) = 2i
heapright(i) = 2i + 1
heapparent(i) = div(i,2)

function heapify!(v, order)
    # https://en.wikipedia.org/wiki/Binary_heap#Building_a_heap
    n = length(v)
    for i = heapparent(n):-1:1
        bubble!((v,n), v[i], i, order, Val(:down))
    end
    return v
end

function bubble!((v,n), vi, i, order)
    if n == 0; return i; end
    if lt(order, vi, v[i])
        return bubble!((v,n), vi, i, order, Val(:up))
    else
        return bubble!((v,n), vi, i, order, Val(:down))
    end
end

function bubble!((v,n), vi, i, order, dir)
    if n == 0; return i; end
    i = walk(
        Bubbler!(v,vi,order,dir),
        n, i, dir
    )
    v[i] = vi
    return i
end

function Bubbler!(v, vi, order, ::Val{:up})
    function UpBubblerF(i,np,nc)
        if np == 0; return :break; end
        vp = v[heapparent(i)]
        if lt(order, vi, vp)
            v[i] = vp
            return :continue
        end
        return :break
    end
    return UpBubblerF
end

function Bubbler!(v, vi, order, ::Val{:down})
    function DownBubblerF(i,np,nc)
        if nc == 0; return :break; end
        if nc == 1
            if lt(order, v[heapleft(i)], vi)
                v[i] = v[heapleft(i)]
                return :continue
            end
            return :break
        end
        vl,vr = v[heapleft(i)], v[heapright(i)]
        d,vj = if lt(order, vl, vr)
            :left, vl
        else
            :right, vr
        end
        if lt(order, vj, vi)
            v[i] = vj
            return d
        end
        return :break
    end
    return DownBubblerF
end


text:
using Distributions
using Base.LinAlg.BLAS

### Base MF definitions
#### Naive mean field
function mag_vis_naive(rbm::RBM, m_hid::Mat{Float64}) ## to be constrained to being only Bernoulli
    buf = gemm('T', 'N', rbm.W, m_hid) .+ rbm.vbias
    return logsig(buf)
end
# Defining a method with same arguments as other mean field approxiamtions
mag_vis_naive(rbm::RBM, m_vis::Mat{Float64}, m_hid::Mat{Float64})=mag_vis_naive(rbm, m_hid)

function mag_hid_naive(rbm::RBM, m_vis::Mat{Float64})
    buf = gemm('N', 'N', rbm.W, m_vis) .+ rbm.hbias
    return logsig(buf)
end

mag_hid_naive(rbm::RBM, m_vis::Mat{Float64}, m_hid::Mat{Float64})=mag_hid_naive(rbm, m_vis)

#### Second order development
function mag_vis_tap2(rbm::RBM, m_vis::Mat{Float64}, m_hid::Mat{Float64}) ## to be constrained to being only Bernoulli
    buf = gemm('T', 'N', rbm.W, m_hid) .+ rbm.vbias

    dv = 0.5-m_vis

    ggg =  gemm('T', 'N', rbm.W2, m_hid-abs2(m_hid))

    second_order = gemm('T', 'N', rbm.W2, m_hid-abs2(m_hid)).*(0.5-m_vis)
    axpy!(1.0, second_order, buf)
    return logsig(buf)
end

function mag_hid_tap2(rbm::RBM, m_vis::Mat{Float64}, m_hid::Mat{Float64})
    buf = gemm('N', 'N', rbm.W, m_vis) .+ rbm.hbias
    second_order = gemm('N', 'N', rbm.W2, m_vis-abs2(m_vis)).*(0.5-m_hid)
    axpy!(1.0, second_order, buf)
    return logsig(buf)
end

#### Third order development

function mag_vis_tap3(rbm::RBM, m_vis::Mat{Float64}, m_hid::Mat{Float64}) ## to be constrained to being only Bernoulli
    buf = gemm('T', 'N', rbm.W, m_hid) .+ rbm.vbias
    second_order = gemm('T', 'N', rbm.W2, m_hid-abs2(m_hid)).*(0.5-m_vis)
    third_order = gemm('T', 'N', rbm.W3, abs2(m_hid).*(1.-m_hid)).*(1/3-2*(m_vis-abs2(m_vis)))
    axpy!(1.0, second_order, buf)
    axpy!(1.0, third_order, buf)
    return logsig(buf)
end

function mag_hid_tap3(rbm::RBM, m_vis::Mat{Float64}, m_hid::Mat{Float64})
    buf = gemm('N', 'N', rbm.W, m_vis) .+ rbm.hbias
    second_order = gemm('N', 'N', rbm.W2, m_vis-abs2(m_vis)).*(0.5-m_hid)
    third_order = gemm('N', 'N', rbm.W3, abs2(m_vis).*(1.-m_vis)).*(1/3-2*(m_hid-abs2(m_hid)))
    axpy!(1.0, second_order, buf)
    axpy!(1.0, third_order, buf)
    return logsig(buf)
end

function equilibrate(rbm::RBM, vis_init::Mat{Float64}, hid_init::Mat{Float64}; iterations=3, approx="tap2", damp=0.5)
    # Redefine names for clarity
    m_vis = copy(vis_init)
    m_hid = copy(hid_init)

    # Set the proper iteration based on the approximation type
    if approx == "naive"
        mag_vis = mag_vis_naive
        mag_hid = mag_hid_naive
    elseif approx == "tap3"
        mag_vis = mag_vis_tap3
        mag_hid = mag_hid_tap3
    else
        mag_vis = mag_vis_tap2
        mag_hid = mag_hid_tap2
    end

    # Take the desired number of steps
    for i=1:iterations
       m_vis = 0.5 * mag_vis(rbm, m_vis, m_hid) + 0.5 * m_vis
       m_hid = 0.5 * mag_hid(rbm, m_vis, m_hid) + 0.5 * m_hid

    end

    return m_vis, m_hid
end

##-----------------------------------------------------------------------------##
### To be deprecated
function iter_mag(rbm::RBM, vis::Mat{Float64}; n_times=3, approx="tap2")
    v_pos = vis
    h_pos = ProbHidCondOnVis(rbm, v_pos)

    if approx == "naive"
        mag_vis = mag_vis_naive
        mag_hid = mag_hid_naive
    elseif approx == "tap3"
        mag_vis = mag_vis_tap3
        mag_hid = mag_hid_tap3
    else
        mag_vis = mag_vis_tap2
        mag_hid = mag_hid_tap2
    end

    m_vis = 0.5 * mag_vis(rbm, vis, h_pos) + 0.5 * vis
    m_hid = 0.5 * mag_hid(rbm, m_vis, h_pos) + 0.5 * h_pos
    for i=1:2

       m_vis = 0.5 * mag_vis(rbm, m_vis, m_hid) + 0.5 * m_vis
       m_hid = 0.5 * mag_hid(rbm, m_vis, m_hid) + 0.5 * m_hid
    end

    return v_pos, h_pos, m_vis, m_hid
end

function iter_mag_persist!(rbm::RBM, vis::Mat{Float64}; n_times=3, approx="tap2")
    v_pos = vis
    h_pos = ProbHidCondOnVis(rbm, v_pos)

    if approx == "naive"
        mag_vis = mag_vis_naive
        mag_hid = mag_hid_naive
    elseif approx == "tap3"
        mag_vis = mag_vis_tap3
        mag_hid = mag_hid_tap3
    else
        mag_vis = mag_vis_tap2
        mag_hid = mag_hid_tap2
    end

    m_vis = rbm.persistent_chain_vis
    m_hid = rbm.persistent_chain_hid

    for i=1:n_times
       m_vis = 0.5 * mag_vis(rbm, m_vis, m_hid) + 0.5 * m_vis
       m_hid = 0.5 * mag_hid(rbm, m_vis, m_hid) + 0.5 * m_hid
    end

    rbm.persistent_chain_vis = m_vis
    rbm.persistent_chain_hid = m_hid

    return v_pos, h_pos, m_vis, m_hid
end


text:
module Displaz
using Compat
using FixedSizeArrays

# Quick and nasty matlab-like plotting interface
export plot3d, plot3d!, clearplot

# Convert julia array into a type name and type appropriate for putting in the
# ply header
ply_type_convert(a::AbstractArray{UInt8})    = ("uint8",   a)
ply_type_convert(a::AbstractArray{UInt16})   = ("uint16",  a)
ply_type_convert(a::AbstractArray{UInt32})   = ("uint32",  a)
ply_type_convert(a::AbstractArray{Int8})     = ("int8",    a)
ply_type_convert(a::AbstractArray{Int16})    = ("int16",   a)
ply_type_convert(a::AbstractArray{Int32})    = ("int32",   a)
ply_type_convert(a::AbstractArray{Float32})  = ("float32", a)
ply_type_convert(a::AbstractArray{Float64})  = ("float64", a)
# Generic cases - actually do a conversion
ply_type_convert{T<:Unsigned}(a::AbstractArray{T}) = ("uint32",  map(UInt32,a))
ply_type_convert{T<:Integer }(a::AbstractArray{T}) = ("int32",   map(Int32,a))
ply_type_convert{T<:Real    }(a::AbstractArray{T}) = ("float64", map(Float64,a))

const array_semantic = 0
const vector_semantic = 1
const color_semantic = 2

# get ply header property name for given field index and semantic
function ply_property_name(semantic, idx)
    if semantic == array_semantic
        string(idx-1)
    elseif semantic == vector_semantic
        ("x", "y", "z", "w")[idx]
    elseif semantic == color_semantic
        ("r", "g", "b")[idx]
    end
end

# Write a set of points to displaz-native ply format
function write_ply_points(filename, nvertices, fields)
    converted_fields = [ply_type_convert(value) for (_,__,value) in fields]
    open(filename, "w") do fid
        write(fid, "ply\n")
        write(fid, "format binary_little_endian 1.0\n")
        write(fid, "comment Displaz native\n")
        for ((name,semantic,_), (typename, value)) in zip(fields, converted_fields)
            n = size(value,2)
            assert(n == nvertices || n == 1)
            write(fid, "element vertex_$name $n\n")
            for i = 1:size(value,1)
                propname = ply_property_name(semantic, i)
                write(fid, "property $typename $propname\n")
            end
        end
        write(fid, "end_header\n")
        for (_,value) in converted_fields
            write(fid, value)
        end
    end
end

function write_ply_lines(filename, position, color, linebreak)
    nvalidvertices = size(position,2)

    # Create and write to ply file
    fid = open(filename, "w")
    write(fid,
        """
        ply
        format binary_little_endian 1.0
        element vertex $nvalidvertices
        property double x
        property double y
        property double z
        element color $nvalidvertices
        property float r
        property float g
        property float b
        element edge $(length(linebreak))
        property list int int vertex_index
        end_header
        """
    )

    write(fid,convert(Array{Float64,2},position))
    write(fid,color)

    realstart = 0
    linelen = []
    range = []
    for i = 1:size(linebreak,1)
        if i != size(linebreak,1)
            linelen = linebreak[i+1] - linebreak[i]
            range = realstart:realstart + linelen-1
        else
            linelen = size(position,2) - linebreak[i] + 1
            range = realstart:realstart + linelen - 1
        end
        write(fid,Int32(linelen))
        write(fid,UnitRange{Int32}(range))
        realstart = realstart + linelen
    end
    close(fid)
end

#const standard_elements = [:position  => (vector_semantic,3),
#                           :color     => (color_semantic,3),
#                           :marksize  => (array_semantic,1),
#                           :markshape => (array_semantic,1)]

const _color_names = @compat Dict('r' => [1.0, 0,   0],
                                  'g' => [0.0, 0.8, 0],
                                  'b' => [0.0, 0,   0.8],
                                  'c' => [0.0, 1,   1],
                                  'm' => [1.0, 0,   1],
                                  'y' => [1.0, 1,   0],
                                  'k' => [0.0, 0,   0],
                                  'w' => [1.0, 1,   1])

const _shape_ids = @compat Dict('.' => 0,
                                's' => 1,
                                'o' => 2,
                                'x' => 3,
                                '+' => 4,
                                '-' => '-')

interpret_color(color) = color
interpret_color(s::AbstractString) = length(s) == 1 ? interpret_color(s[1]) : error("Unknown color abbreviation $s")
interpret_color(c::Char) = _color_names[c]
interpret_shape(markershape) = markershape
interpret_shape(c::Char) = [_shape_ids[c]]
interpret_shape(s::Vector{Char}) = Int[_shape_ids[c] for c in s]
interpret_shape(s::AbstractString) = Int[_shape_ids[c] for c in s]
interpret_linebreak(nvertices, linebreak) = linebreak
interpret_linebreak(nvertices, i::Integer) = i == 1 ? [1] : 1:i:nvertices

# Multiple figure window support
# TODO: Consider how the API below relates to Plots.jl and its tendency to
# create a lot of new figure windows rather than clearing existing ones.
type DisplazWindow
    name::AbstractString
end

_current_figure = DisplazWindow("default")
"Get handle to current figure window"
function current()
    _current_figure
end

"Get figure window by name may be new or already existing"
function figure(name::AbstractString)
    global _current_figure
    _current_figure = DisplazWindow(name)
end
"Get figure window with given id"
figure(id::Integer) = figure("Figure $id")

_figure_id = 1
"Create next incrementally named figure window, counting automatically from \"Figure 1\""
function newfigure()
    global _figure_id
    id = _figure_id
    _figure_id += 1
    figure(id)
end

"""
Add 3D points or lines to the current plot.

```
  plot3d([plotobj,] position; attr1=value1, ...)
```

The `position` array should be a 3xN array of vertex positions. (This is
consistent with treating the components of a single point as a column vector.
It also has the same memory layout as an array of short `Vector3` instances
from ImmutableArrays, for example).

The `plotobj` argument is optional and determines which plot window to send the
data to.  If it's not used the data will be sent to the plot window returned
by `current()`.

TODO: Revisit the nasty decision of the shape of position again - the above
choice is somewhat inconsistent with supplying markersize / markershape as a
column vector :-(  Can we have a set of consistent broadcasting rules for this?
It seems like the case of a 3x3 matrix will always be ambiguous if we try
to guess what the user wants.

### Data set attributes

The following attributes can be attached to a dataset on each call to `plot3d`:

  * `label` - A string labeling the data set

### Vertex attributes

Each point may have a set of vertex attributes attached to control the visual
representation and tag the point for inspection. The following are supported by
the default shader:

  * `color`       - A 3xN array of vertex colors
  * `markersize`  - Vertex marker size
  * `markershape` - Vector of vertex marker shapes.  Shape can be represented
                    either by a Char or a numeric identifier between 0 and 4:

```
                    sphere - '.' 0    square - 's' 1
                    circle - 'o' 2    times  - 'x' 3
                    cross  - '+' 4
```

### Plotting points

To plot 10000 random points with distance from the origin determining the
color, and random marker shapes:
```
  P = randn(3,10000)
  c = 0.5./sumabs2(P,1) .* [1,0,0]
  plot3d(P, color=c, markershape=rand(1:4,10000))
```

### Plotting lines

To plot a piecewise linear curve between 10000 random vertices
```
  plot3d(randn(3,10000), markershape="-")
```

When plotting lines, the `linebreak` keyword argument can be used to break the
position array into multiple line segments.  Each index in the line break array
is the initial index of a line segment.
"""
function plot3d(plotobj::DisplazWindow, position; color=[1,1,1], markersize=[0.1], markershape=[0],
                label=nothing, linebreak=[1], _overwrite_label=false)
    nvertices = size(position, 2)
    color = interpret_color(color)
    markershape = interpret_shape(markershape)
    linebreak = interpret_linebreak(nvertices, linebreak)
    size(position, 1) == 3 || error("position must be a 3xN array")
    size(color, 1)    == 3 || error("color must be a 3xN array")
    # FIXME in displaz itself.  No repmat waste should be required.
    if size(color,2) == 1
        color = repmat(color, 1, nvertices)
    end
    # Ensure all fields are floats for now, to avoid surprising scaling in the
    # shader
    color = map(Float32,color)
    markersize = map(Float32,markersize)
    size(color,2) == nvertices || error("color must have same number of rows as position array")
    filename = tempname()*".ply"
    seriestype = "Points"
    if '-' in markershape # Plot lines
        seriestype = "Line"
        write_ply_lines(filename, position, color, linebreak)
    else # Plot points
        if length(markersize) == 1
            markersize = repmat(markersize, 1, nvertices)
        end
        if length(markershape) == 1
            markershape = repmat(markershape, 1, nvertices)
        end
        write_ply_points(filename, nvertices, (
                         (:position, vector_semantic, position),
                         (:color, color_semantic, color),
                         # FIXME: shape of markersize??
                         (:markersize, array_semantic, vec(markersize)'),
                         (:markershape, array_semantic, vec(markershape)'),
                         ))
    end
    if label === nothing
        label = "$seriestype [$nvertices vertices]"
    end
    addopt = _overwrite_label ? [] : "-add"
    run(`displaz -script $addopt -server $(plotobj.name) -label $label -shader generic_points.glsl -rmtemp $filename`)
    nothing
end

# Interop with FixedSizeArrays.
plot3d{V<:FixedVector{3}}(plotobj::DisplazWindow, position::Vector{V}; kwargs...) =
    plot3d(plotobj, destructure(position); kwargs...)

"""
Overwrite points or lines with the same label on the 3D plot

See plot3d for documentation
"""
function plot3d!(plotobj::DisplazWindow, position; kwargs...)
    plot3d(plotobj, position; _overwrite_label=true, kwargs...)
end

# Plot to current window
plot3d!(position; kwargs...) = plot3d!(current(), position; kwargs...)
plot3d(position; kwargs...)  = plot3d(current(), position; kwargs...)

"""
    clearplot([plotobj::DisplazWindow=current()], [pattern1, ...])

Clear all or a subset of datasets in a plot window.

If not specified, `plotobj` is the current plot window.  If no patterns are
supplied, clears all data sets.  If one or more patterns are given, the dataset
labels matching those patterns will be removed.  Patterns shoudl be specified
using unix shell glob pattern syntax.
"""
function clearplot(plotobj::DisplazWindow, patterns...)
    unload_args = isempty(patterns) ? ["-clear"] : ["-unload", patterns...]
    run(`displaz -script -server $(plotobj.name) $unload_args`)
    nothing
end
clearplot(patterns...) = clearplot(current(), patterns...)

end


text:
__precompile__()

module Xpress

    using Libdl

    # Load in `deps.jl`, complaining if it does not exist
    const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")

    if !isfile(depsjl_path)
        error("XPRESS cannot be loaded. Please run Pkg.build(\"Xpress\").")
    end

    include(depsjl_path)

    ### imports

    import Base.show, Base.copy

    # ### exports
    # export

    # ## xprs_env
    # free_env,

    # ## xprs_params
    # getparam, setparam!, setparams!,

    # ## xprs_model
    # set_sense!, update_model!, reset_model!,
    # read_model, write_model,
    # copy,

    # ## xprs_attrs
    # model_name, model_sense, model_type,
    # num_vars, num_constrs, num_sos, num_qconstrs,
    # num_cnzs, num_qnzs, num_qcnzs,
    # is_qp, is_qcp, is_mip,

    # lowerbounds, upperbounds, objcoeffs, set_objcoeffs!,

    # ## xprs_vars
    # XPRS_CONTINUOUS, XPRS_BINARY, XPRS_INTEGER,
    # add_var!, add_vars!, add_cvar!, add_cvars!,
    # add_bvar!, add_bvars!, add_ivar!, add_ivars!,
    # del_vars!,

    # ## xprs_constrs
    # add_constr!, add_constrs!, add_constrs_t!,
    # add_rangeconstr!, add_rangeconstrs!, add_rangeconstrs_t!,
    # get_constrmatrix, add_sos!, del_constrs!, chg_coeffs!,

    # ## xprs_quad
    # add_qpterms!, add_qconstr!,

    # ## higher level
    # xpress_model,

    # ## xprs_solve
    # optimize, computeIIS, get_solution,
    # get_status, OptimInfo, get_optiminfo, get_objval

    ### include source files
    module Lib

        import ..Xpress
        const libxprs = Xpress.libxprs

        include("ctypes.jl")
        include("common.jl")
        include("lib.jl")

        # include("custom.jl")

    end

    include("utils.jl")

    include("helper.jl")

    include("api.jl")

    # # license checker
    include("license.jl")

    function initialize()
        userlic()
        init() # Call XPRSinit for initialization
        atexit(free) # Call free when process terminates
    end
    # include("xprs_full_defines.jl")

    # include("xprs_common.jl")
    # include("xprs_env.jl")

    # include("xprs_model.jl")
    # include("xprs_params.jl")
    # include("xprs_attrs.jl")

    # include("xprs_vars.jl")

    # include("xprs_constrs.jl")
    # include("xprs_quad.jl")
    # include("xprs_highlevel.jl")

    # include("xprs_solve.jl")
    # include("xprs_callbacks.jl")
    # include("xprs_iis.jl")

    # include("XpressSolverInterface.jl")
    include("MOI/MOI_Wrapper.jl")

end


text:
module Components

using LinearAlgebra
using StaticArrays, ComponentArrays
using UnPack

using Flight.Modeling
using Flight.Plotting
import Flight.Plotting: plots

using Flight.Attitude

using Flight.Dynamics

export HasMass, HasNoMass, get_mass_properties
export HasWrench, HasNoWrench, get_wr_b
export HasAngularMomentum, HasNoAngularMomentum, get_hr_b

############################ MassTrait #################################

abstract type MassTrait end
struct HasMass <: MassTrait end
struct HasNoMass <: MassTrait end

MassTrait(::S) where {S<:System} = error(
    "Please extend Components.MassTrait for $S")

get_mass_properties(sys::System) = get_mass_properties(MassTrait(sys), sys)

get_mass_properties(::HasMass, sys::System) = error(
    "$(typeof(sys)) has the HasMass trait, but no get_mass_properties method is defined for it")

get_mass_properties(::HasNoMass, sys::System) = error(
    "$(typeof(sys)) has the HasNoMass trait")

###################### AngularMomentumTrait ##########################

abstract type AngularMomentumTrait end
struct HasAngularMomentum <: AngularMomentumTrait end
struct HasNoAngularMomentum <: AngularMomentumTrait end

#prevents the trait system from failing silently when wrongly extended
AngularMomentumTrait(::S) where {S<:System} = error(
    "Please extend Components.AngularMomentumTrait for $S")

get_hr_b(sys::System) = get_hr_b(AngularMomentumTrait(sys), sys)

get_hr_b(::HasAngularMomentum, sys::System) = error(
    "$(typeof(sys)) has angular momentum, but no method get_hr_b was defined for it")

get_hr_b(::HasNoAngularMomentum, sys::System) = zeros(SVector{3})

#default implementation for a SystemGroup with the HasAngularMomentum trait, tries
#to sum the angular momentum from its individual components. override as required
@inline @generated function get_hr_b(::HasAngularMomentum, sys::System{D}) where {D<:SystemGroupDescriptor}

    # Core.print("Generated function called")
    ex = Expr(:block)
    push!(ex.args, :(h = SVector(0., 0., 0.))) #allocate
    for label in fieldnames(D)
        push!(ex.args,
            :(h += get_hr_b(sys.subsystems[$(QuoteNode(label))])))
    end
    return ex

end

###################### WrenchTrait ##########################

abstract type WrenchTrait end
struct HasWrench <: WrenchTrait end
struct HasNoWrench <: WrenchTrait end

#prevents the trait system from failing silently when wrongly extended
WrenchTrait(::S) where {S<:System} = error(
    "Please extend Components.WrenchTrait for $S")

get_wr_b(sys::System) = get_wr_b(WrenchTrait(sys), sys)

get_wr_b(::HasWrench, sys::System) = error(
    "$(typeof(sys)) is a Wrench source, but no method get_wr_b was defined for it")

get_wr_b(::HasNoWrench, sys::System) = Wrench()

#default implementation for a SystemGroup with the HasWrench trait, tries
#to sum all the Wrenches from its individual components. override as required
@inline @generated function get_wr_b(::HasWrench, sys::System{D}) where {D<:SystemGroupDescriptor}

    # Core.print("Generated function called")
    ex = Expr(:block)
    push!(ex.args, :(wr = Wrench())) #allocate a zero wrench
    for label in fieldnames(D)
        push!(ex.args,
            :(wr += get_wr_b(sys.subsystems[$(QuoteNode(label))])))
    end
    return ex

end

end #module


text:
module PowerLAPACK

include("lapack.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Priority <: SwaggerModel

    function Priority(;)
        o = new()
        o
    end
end # type Priority

const _property_map_Priority = Dict{Symbol,Symbol}()
const _property_types_Priority = Dict{Symbol,String}()
Base.propertynames(::Type{ Priority }) = collect(keys(_property_map_Priority))
Swagger.property_type(::Type{ Priority }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Priority[name]))}
Swagger.field_name(::Type{ Priority }, property_name::Symbol) =  _property_map_Priority[property_name]

function check_required(o::Priority)
    true
end

function validate_property(::Type{ Priority }, name::Symbol, val)
end


text:
using RecipesBase
using LaTeXStrings

@recipe function f(to::TransferOperatorTriangulationExact)
    seriestype := :heatmap

    xlabel --> "i"
    ylabel --> "j"
    colorbar_title --> L"p(S_i \to im(S_j))"
    xticks --> false
    yticks --> false
    size --> (600, 500)
    to.transfermatrix
end

@recipe function f(to::TransferOperatorTriangulationApprox)
    seriestype := :heatmap

    xlabel --> "i"
    ylabel --> "j"
    colorbar_title --> L"p(S_i \to im(S_j))"
    xticks --> false
    yticks --> false
    size --> (600, 500)

    to.transfermatrix
end

@recipe function f(to::TransferOperatorRectangularBinning)
    seriestype := :heatmap

    xlabel --> "i"
    ylabel --> "j"
    colorbar_title --> L"p(b_i \to im(b_j))"
    xticks --> false
    yticks --> false
    size --> (600, 500)

    to.transfermatrix
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1beta2DaemonSetList <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiAppsV1beta2DaemonSet} } # spec name: items
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ListMeta } # spec name: metadata

    function IoK8sApiAppsV1beta2DaemonSetList(;apiVersion=nothing, items=nothing, kind=nothing, metadata=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2DaemonSetList, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAppsV1beta2DaemonSetList, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiAppsV1beta2DaemonSetList, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAppsV1beta2DaemonSetList, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type IoK8sApiAppsV1beta2DaemonSetList

const _property_map_IoK8sApiAppsV1beta2DaemonSetList = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("items")=>Symbol("items"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_IoK8sApiAppsV1beta2DaemonSetList = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("items")=>"Vector{IoK8sApiAppsV1beta2DaemonSet}", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ListMeta")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2DaemonSetList }) = collect(keys(_property_map_IoK8sApiAppsV1beta2DaemonSetList))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2DaemonSetList }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAppsV1beta2DaemonSetList[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2DaemonSetList }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2DaemonSetList[property_name]

function check_required(o::IoK8sApiAppsV1beta2DaemonSetList)
    (getproperty(o, Symbol("items")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta2DaemonSetList }, name::Symbol, val)
end


text:
# Matern 3/2 ARD covariance function

"""
    Mat32Ard <: MaternARD

ARD Matern 3/2 kernel (covariance)
```math
k(x,x') = σ²(1 + √3|x-x'|/L)\\exp(- √3|x-x'|/L)
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)`` and signal standard deviation ``σ`` where
``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct Mat32Ard <: MaternARD
    "Inverse squared length scale"
    iℓ2::Vector{Float64}
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Mat32Ard(ll::Vector{Float64}, lσ::Float64)

    Create `Mat32Ard` with length scale `exp.(ll)` and signal standard deviation `exp(lσ)`.
    """
    Mat32Ard(ll::Vector{Float64}, lσ::Float64) = new(exp.(-2 .* ll), exp(2 * lσ), [])
end

function set_params!(mat::Mat32Ard, hyp::VecF64)
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat32 kernel only has $(num_params(mat)) parameters"))
    @views @. mat.iℓ2 = exp(-2 * hyp[1:(end-1)])
    mat.σ2 = exp(2 * hyp[end])
end

get_params(mat::Mat32Ard) = [-log.(mat.iℓ2) / 2; log(mat.σ2) / 2]
get_param_names(mat::Mat32Ard) = [get_param_names(mat.iℓ2, :ll); :lσ]
num_params(mat::Mat32Ard) = length(mat.iℓ2) + 1

Statistics.cov(mat::Mat32Ard, r::Float64) =
    (s = √3 * r; mat.σ2 * (1 + s) * exp(-s))

dk_dll(mat::Mat32Ard, r::Float64, wdiffp::Float64) = 3 * mat.σ2 * wdiffp * exp(-√3 * r)


text:
"""
    Sobel-sample x+y in [0,10]x[0,10],
    then minimize it on Section([NaN,10.0]),
    and verify that the minimum is on x,y=(0,10)
    rather than in (0,0)
"""

using Surrogates
using Test

lb = [ 0.0, 0.0, 0.0]
ub = [10.0,10.0,10.0]
x = Surrogates.sample(10,lb,ub,LatinHypercubeSample())
f = x -> x[1]+x[2]+x[3]
y = f.(x)
f([0,0,0]) == 0

f_hat = Kriging(x,y,lb,ub)

f_hat([0,0,0])

isapprox(f([0,0,0]), f_hat([0,0,0]))

""" The global minimum is at (0,0) """

(xy_min, f_hat_min) = surrogate_optimize(
    f,
    DYCORS(), lb, ub,
    f_hat,
    SobolSample())

isapprox(xy_min[1], 0.0, atol=1e-3)

""" The minimum on the (0,10) section is around (0,10) """

section_sampler_z_is_10 = SectionSample(
    [NaN64, NaN64, 10.0],
    Surrogates.UniformSample())

@test [3] == Surrogates.fixed_dimensions(section_sampler_z_is_10)
@test [1,2] == Surrogates.free_dimensions(section_sampler_z_is_10)

Surrogates.sample(5, lb, ub, section_sampler_z_is_10)

(xy_min, f_hat_min) = surrogate_optimize(
    f,
    EI(), lb, ub,
    f_hat,
    section_sampler_z_is_10, maxiters=1000)

isapprox(xy_min[1],  0.0, atol=0.1)
isapprox(xy_min[2],  0.0, atol=0.1)
isapprox(xy_min[3], 10.0, atol=0.1)