text:
#
# Solver for Lambert's problem
#
# References:
# [1] David, A. "Vallado. Fundamentals of Astrodynamics and Applications." (2013).
#

"""
Solves Lambert's problem through the use of univeral variables.
"""
function lambert_universal(r̅₁::AbstractVector, r̅₂::AbstractVector, μ::Number, Δt::Number;
                           trajectory=:short, tolerance=1e-12, max_iter=25)

    # Specify short way, or long way trajectory
    if trajectory == :short
        tₘ = 1
    elseif trajectory == :long
        tₘ = -1
    else
        throw(ArgumentError("`trajectory` must be set to `:short` or `:long`"))
    end

    r₁ = norm(r̅₁)
    r₂ = norm(r̅₂)

    cosΔν = (r̅₁⋅r̅₂) / (r₁*r₂)
    Δν    = asin(u"rad", tₘ * √(1 - (cosΔν)^2))

    A = upreferred(tₘ * √(r₂*r₁ * (1 + cosΔν)))

    if A ≈ 0
        throw(ErrorException("Can't calculate the orbit."))
    end

    ψₙ = 0.0
    C₂ = 1/2
    C₃ = 1/6

    ψ₊ = 4π^2
    ψ₋ = -4π
    yₙ = r₁ + r₂ + (A * (ψₙ*C₃ - 1) / √(C₂))

    Δtₙ = Δt + 1u"s"
    iter = 0

    while (iter < max_iter) &&
          (abs(Δtₙ - Δt) > (tolerance * oneunit(Δt))) || (A > 0 *oneunit(A) && yₙ < 0 * oneunit(yₙ))

        yₙ = r₁ + r₂ + (A * (ψₙ*C₃ - 1) / √(C₂))
        χₙ = √(yₙ / C₂)
        Δtₙ = (χₙ^3 * C₃ + A*√(yₙ)) / √(μ)

        if Δtₙ < Δt
            ψ₋ = ψₙ
        else
            ψ₊ = ψₙ
        end

        ψₙ = (ψ₊ + ψ₋) / 2
        if ψₙ > tolerance
            C₂ = (1 - cos(√(ψₙ))) /  ψₙ
            C₃ = (√(ψₙ) - sin(√(ψₙ))) / √(ψₙ^3)
        elseif ψₙ < -tolerance
            C₂ = (1 - cosh(√(-ψₙ))) / ψₙ
            C₃ = (sinh(√(-ψₙ)) - √(-ψₙ)) / √((-ψₙ)^3)
        else
            C₂ = 1.0 / 2.0
            C₃ = 1.0 / 6.0
        end

        iter += 1

    end

    f = 1 - yₙ/r₁
    ġ = 1 - yₙ/r₂
    g = A * √(yₙ/μ)

    v̅₁ = upreferred.((r̅₂ .- (f .* r̅₁)) ./ g)
    v̅₂ = upreferred.(((ġ .* r̅₂) .- r̅₁) ./ g)

    return v̅₁, v̅₂

end

"""
The following code was converted to Julia, from a [GitHub repository](https://github.com/rodyo/FEX-Lambert)
that hosts a MATLAB implementation. At the time of writing, this respository has a BSD license. I'm providing the copyright notice here, as instructed by the license text.

```
Copyright (c) 2018, Rody Oldenhuis
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of this project.
```
"""
function σmax(y)

    an = [
        4.000000000000000e-001,
        2.142857142857143e-001,
        4.629629629629630e-002,
        6.628787878787879e-003,
        7.211538461538461e-004,
        6.365740740740740e-005,
        4.741479925303455e-006,
        3.059406328320802e-007,
        1.742836409255060e-008,
        8.892477331109578e-010,
        4.110111531986532e-011,
        1.736709384841458e-012,
        6.759767240041426e-014,
        2.439123386614026e-015,
        8.203411614538007e-017,
        2.583771576869575e-018,
        7.652331327976716e-020,
        2.138860629743989e-021,
        5.659959451165552e-023,
        1.422104833817366e-024,
        3.401398483272306e-026,
        7.762544304774155e-028,
        1.693916882090479e-029,
        3.541295006766860e-031,
        7.105336187804402e-033
    ]

    powers  = y.^(1:25)
    σ       = 4/3 * powers ⋅ an
    ∂σ_∂x   = ((1:25) .* vcat(1, powers[1:24])) ⋅ an
    ∂²σ_∂x² = ((1:25) .* (0:24) .* vcat(1/y, 1, powers[1:23])) ⋅ an
    ∂³σ_∂x³ = ((1:25) .* (0:24) .* (-1:23) .*
                vcat(1/y/y, 1/y, 1, powers[1:22])) ⋅ an

    return σ, ∂σ_∂x, ∂²σ_∂x², ∂³σ_∂x³

end;

"""
The following code was converted to Julia, from a [GitHub repository](https://github.com/rodyo/FEX-Lambert)
that hosts the MATLAB implementation. At the time of writing, the respository has a BSD license. I'm providing the copyright notice here, as instructed by the license text.

```
Copyright (c) 2018, Rody Oldenhuis
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of this project.
```

"""
function LancasterBlanchard(x, q, m)

    # Validate input
    if x < -one(x)
        @warn "`x` provided implies negative eccentricity. Correcting..."
        x = abs(x) - 2 * one(x)
    elseif x == -one(x)
        @warn "`x` provided is invalid. Setting `x` to the next float..."
        x = nextfloat(x)
    end

    # Parameter E
    E = x*x - one(x)

    if x == 1 # parabolic input ⟹ analytical solution

        T   = 4/3 * (1 - q^3)
        ∂T  = 4/5 * (q^5 - 1)
        ∂²T = ∂T + 120/70 * (1 - q^7)
        ∂³T = 3 * (∂²T - ∂T) + 2400 / 1080 * (q^9 - 1)

    elseif abs(x - one(x)) < 1e-2 # almost parabolic ⟹ use series

        # Series expansion for T & associated partials
        σ₁, ∂σ_∂x₁, ∂²σ_∂x₂₁, ∂³σ_∂x₃₁ = σmax(-E)
        σ₂, ∂σ_∂x₂, ∂²σ_∂x₂₂, ∂³σ_∂x₃₂ = σmax(-E * q * q)

        T   = σ₁ - q^3 * σ₂
        ∂T  = 2 * x * (q^5 * ∂σ_∂x₂ - ∂σ_∂x₁)
        ∂²T = ∂T/x + 4 * x^2 * (∂²σ_∂x₂₁ - q^7 * ∂²σ_∂x₂₂)
        ∂³T = 3 * (∂²T - ∂T/x) / x + 8 * x * x * (q^9 * ∂³σ_∂x₃₂ - ∂³σ_∂x₃₁)

    else

        y = √(abs(E))
        z = √(1 + q^2 * E)
        f = y * (z - q * x)
        g = x * z - q * E

        if E < zero(E)
            δ = atan(f, g) + m*π
        elseif E == zero(E)
            δ = zero(E)
        else
            δ = log(ℯ, max(0, f+g))
        end

        T   = 2 * (x - q * z - δ / y) / E

        ∂T  = (4 - 4 * q^3 * x/z - 3 * x * T) / E

        ∂²T = (-4 * q^3/z * (1 - q^2 * x^2 / z^2) - 3*T - 3 * x * ∂T) / E

        ∂³T = (4 * q^3 / z^2 *
                ((1 - q^2 * x^2 / z^2) + 2 * q^2 * x/z^2 * (z - x)) -
                    8*∂T - 7*x*∂²T) / E

    end

    return T, ∂T, ∂²T, ∂³T

end

"""
The following code was converted to Julia, from a [GitHub repository](https://github.com/rodyo/FEX-Lambert)
that hosts a MATLAB implementation. At the time of writing, this respository has a BSD license. I'm providing the copyright notice here, as instructed by the license text.

```
Copyright (c) 2018, Rody Oldenhuis
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of this project.
```
"""
function minmax_distances(r̲₁, r̲₂, r₁, r₂, δₜ, a, v̲₁, v̲₂, m, μ)

    min_distance = min(r₁, r₂)
    max_distance = max(r₁, r₂)

    longway = abs(δₜ) > π # check if longway trajectory was selected

    # Find eccentricity vector using triple product identity
    e̲ = ((v̲₁⋅v̲₁) * r̲₁ - (v̲₁ ⋅ r̲₁) * v̲₁) / μ - r̲₁ / r₁

    # Scalar eccentricity
    e = norm(e̲)

    # Pericenter / apocenter
    pericenter = a * (1 - e)
    apocenter  = Inf # For parabolic / hyperbolic case
    if e < one(e)
        apocenter = a * (1 + e) # elliptical case
    end

    # We have the eccentricity vector, so we know exactly
    # where the pericenter is. Given δₜ and that fact
    # to cross-check if the trajectory goes past it

    if m > zero(m) # always elliptical, both apses traversed
        min_distance = pericenter
        max_distance = apocenter
    else # more complicated

        pm1 = sign(r₁ * r₁ * (e̲⋅v̲₁) - (r̲₁⋅e̲) * (r̲₁⋅v̲₁))
        pm2 = sign(r₂ * r₂ * (e̲⋅v̲₂) - (r̲₂⋅e̲) * (r̲₂⋅v̲₂))

        # Make sure θ₁, θ₂ ∈ (-1, 1)
        θ₁ = pm1 * acos(max(-1, min(1, (r̲₁./r₁) ⋅ (e̲./e))))
        θ₂ = pm2 * acos(max(-1, min(1, (r̲₂./r₂) ⋅ (e̲./e))))

        if θ₁+θ₂ < zero(θ₁)

            if abs(θ₁) + abs(θ₁) == δₜ   # pericenter was passed
                min_distance = pericenter
            else
                min_distance = apocenter # apocenter was passed
            end

        elseif longway
            min_distance = pericenter
            if e < one(e)
                max_distance = apocenter
            end
        end

    end

    return min_distance, max_distance

end

"""
The following code was converted to Julia, from a [GitHub repository](https://github.com/rodyo/FEX-Lambert)
that hosts a MATLAB implementation. At the time of writing, this respository has a BSD license. I'm providing the copyright notice here, as instructed by the license text.

```
Copyright (c) 2018, Rody Oldenhuis
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of this project.
```
"""
function lambert_lancaster_blanchard(
        r̲₁::AbstractVector, r̲₂::AbstractVector,
        μ::Number, Δt::Number;
        revolutions = 0,
        branch=:left,
        trajectory=:short,
        tolerance=1e-12,
        max_iter=25,
        output_extrema=Val{false})

    m = revolutions

    if output_extrema == Val{true}
        error_output = [NaN, NaN, NaN] * u"km/s", [NaN, NaN, NaN] * u"km/s", (NaN * u"km", NaN * u"km")
    elseif output_extrema == Val{false}
        error_output = [NaN, NaN, NaN] * u"km/s", [NaN, NaN, NaN] * u"km/s"
    else
        throw(ArgumentError("Keyword argument `output_extrema` must be set to Val{true} or Val{false}."))
    end

    if trajectory ∉ (:short, :long)
        throw(ArgumentError("Must specify :long or :short for trajectory kwarg."))
    end

    if Δt < zero(Δt)
        throw(ArgumentError(string(
            "Time of flight must be non-negative!",
            "Use the `trajectory` kwarg to specify",
            "longway or shortway trajectories.")))
    end

    if m < zero(m)
        throw(ArgumentError(string(
            "Number of revolutions must be non-negative!",
            "Use the `branch` kwarg to specify",
            "left or righht branch solutions.")))
    end

    # Collect unit vectors and magnitudes
    # of provided position vectors
    r₁ = norm(r̲₁)
    r₂ = norm(r̲₂)
    r̂₁ = normalize(r̲₁)
    r̂₂ = normalize(r̲₂)

    # Unit position vector orthogonal to
    # position vectors
    r̂ₒ = normalize(r̲₁ × r̲₂)

    # Tangential-direction unit vectors
    t̂₁ = r̂ₒ × r̂₁
    t̂₂ = r̂ₒ × r̂₂

    # Turn angle
    δₜ = acos(max(-1, min(1, (r̲₁⋅r̲₂)/r₁/r₂)))

    # If longway, account for final velocity
    # direction by flipping δₜ sign
    if trajectory == :long
        δₜ -= 2π
    end

    # Constants
    c = √(r₁^2 + r₂^2 - 2*r₁*r₂*cos(δₜ))
    s = (r₁ + r₂ + c) / 2
    T = √(8μ / s^3) * Δt
    q = √(r₁ * r₂) / s * cos(δₜ/2)

    # Initial values
    T₀  = first(LancasterBlanchard(0, q, m))
    δT  = T₀ - T
    phr = mod(2 * atan(1 - q^2, 2 * q), 2π)

    # Assume failure
    v₁ = [NaN, NaN, NaN]
    v₂ = v₁
    rₑ = (NaN, NaN)

    if m == zero(m) # single revolution

        x₀₁ = T₀ * δT / 4/ T
        if δT > zero(Δt)
            x₀ = x₀₁
        else
            x₀₁ = δT / (4 - δT)
            x₀₂ = -√(-δT / (T + T₀/2))

            # TODO potential bug, see https://github.com/rodyo/FEX-Lambert/issues/1
            W   = x₀₁ + 1.7 * √(2 - phr/π)  # one of these is right! (1)
            # W   = x₀₁ + 1.7 * √(2 - δₜ/π)     # is it this one?     (2)

            if W ≥ zero(W)
                x₀₃ = x₀₁
            else
                x₀₃ = x₀₁ + ((-W)^(1/16) * (x₀₂ - x₀₁)) # one of these is right! (1)
                # x₀₃ = x₀₁ + (-W^(1/16) * (x₀₂ - x₀₁)) # is it this one?        (2)
            end

            λ  = 1 + x₀₃ * (1 + x₀₁)/2 - 0.03 * x₀₃^2 * √(1 + x₀₁)
            x₀ = λ * x₀₃
        end

        # Check if solution can't be found
        if x₀ < -one(x₀)
            @error "Unable to find solution."
            return error_output
        end

    else

        # Minumum ∂T(x)
        xMπ = 4 / (3π * (2m + 1))
        if phr < π
            xM₀ = xMπ * (phr / π)^(1/8)
        elseif phr > π
            xM₀ = xMπ * (2 - (2 - phr/π)^(1/8))
        else
            xM₀ = zero(xMπ)
        end

        # Halley's method
        xM = xM₀
        ∂T = Inf
        iter = 0

        while abs(∂T) > tolerance

            iter += 1

            _, ∂T, ∂²T, ∂³T = LancasterBlanchard(xM, q, m)

            xMp = xM
            xM = xM - 2 * ∂T .* ∂²T ./ (2 * ∂²T.^2 - ∂T .* ∂³T)

            if mod(iter, 7) != 0
                xM = (xMp + xM)/2
            end

            if iter > max_iter
                @error "Unable to find solution."
                return error_output
            end

        end

        # xM should be elliptic: xM ∈ (-1, 1)
        if xM < -one(xM) || xM > one(xM)
            @error "Unable to find solution."
            return error_output
        end

        # Corresponding time
        TM = first(LancasterBlanchard(xM, q, m))

        # Check that T > minimum
        if TM > T
            @error "Unable to find solution."
            return error_output
        end

        # Move onto initial values for second solution!

        # Initial values
        TmTM = T - TM
        T0mTM = T₀ - TM
        _, ∂T, ∂²T, __ = LancasterBlanchard(xM, q, m)

        # First estimate if m > 0
        if branch == :left

            x = √(TmTM / (∂²T/2 + TmTM / (1 - xM)^2))
            W = xM + x
            W = 4 * W / (4 + TmTM) + (1 - W)^2
            x₀  = x * (1 - (1 + m + (δₜ - 1/2)) /
                    (1 + 0.15m) * x * (W/2 + 0.03x * √W)) + xM
            if x₀ > one(x₀)
                @error "Unable to find solution."
                return error_output
            end

        else # Second estimate if m > 0

            if δT > zero(δT)
                x₀ = xM - √(TM / (∂²T/2 - TmTM * (∂²T/2/T0mTM - 1/xM^2)))
            else
                x₀₀ = δT / (4 - δT)
                W   = x₀₀ + 1.7 * √Complex(2 * (1 - phr))

                if real(W) ≥ zero(real(W))
                    x₀₃ = x₀₀
                else
                    x₀₃ = x₀₀ - √((-W)^(1/8)) * (x₀₀ + √(-δT / (1.5*T₀ - δT)))
                end

                W = 4 / (4 - δT)
                λ = (1 + (1 + m + 0.24*(δₜ - 1/2)) /
                        (1 + 0.15m) * x₀₃ * (W/2 - 0.03x₀₃ * √(W)))
                x₀ = x₀₃ * λ
            end

            if real(x₀) < -one(real(x₀))
                @error "Unable to find solution."
                return error_output
            end

        end

    end

    # Finally, find root of Lancaster & Blancard's function

    # Halley's method
    x  = x₀
    Tₓ = Inf
    iter = 0

    while abs(Tₓ) > tolerance

        iter += 1

        Tₓ, ∂T, ∂²T, _ = LancasterBlanchard(x, q, m)

        # Find the root of the difference between Tₓ and
        # required time T
        Tₓ = Tₓ - T

        # New value of x
        xₚ = x
        x  = x - 2 * Tₓ * ∂T ./ (2 * ∂T^2 - Tₓ * ∂²T)

        if mod(iter, 7) != 0
            x = (xₚ + x) / 2
        end

        if iter > max_iter
            @error "Unable to find solution."
            return error_output
        end

    end

    # Calculate terminal velocities

    # Constants
    γ = √(μ * s/2)

    if c == zero(c)
        σ = one(x)
        ρ = zero(x)
        z = abs(x)
    else
        σ = 2 * √(r₁*r₂ / c^2) * sin(δₜ/2)
        ρ = (r₁ - r₂) / c
        z = √(1 + q^2 * (x^2 - 1))
    end

    # Radial component
    vᵣ₁ =  γ * ((q * z - x) - ρ * (q * z + x)) / r₁
    vᵣ₂ = -γ * ((q * z - x) + ρ * (q * z + x)) / r₂
    v̲ᵣ₁ = vᵣ₁ * r̂₁
    v̲ᵣ₂ = vᵣ₂ * r̂₂

    # Tangential component
    vₜ₁ = σ * γ * (z + q*x) / r₁
    vₜ₂ = σ * γ * (z + q*x) / r₂
    v̲ₜ₁  = vₜ₁ * t̂₁
    v̲ₜ₂  = vₜ₂ * t̂₂

    # Cartesian velocity
    v̲₁ = v̲ₜ₁ .+ v̲ᵣ₁
    v̲₂ = v̲ₜ₂ .+ v̲ᵣ₂

    if output_extrema == Val{true}

        # Find min / max distances
        a = s/2 / (1 - x^2)
        rₘ = minmax_distances(r̲₁, r̲₂, r₁, r₂, δₜ, a, v̲₁, v̲₂, m, μ)

        return v̲₁, v̲₂, rₘ

    else

        return v̲₁, v̲₂

    end

end

"""
The following code was converted to Julia, from a [GitHub repository](https://github.com/rodyo/FEX-Lambert)
that hosts a MATLAB implementation. At the time of writing, this respository has a BSD license. I'm providing the copyright notice here, as instructed by the license text.

```
Copyright (c) 2018, Rody Oldenhuis
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of this project.
```
"""
function lambert(r1vec, r2vec, tf, m, muC)

    # original documentation:
    #=
    This routine implements a new algorithm that solves Lambert's problem. The
    algorithm has two major characteristics that makes it favorable to other
    existing ones.

    1) It describes the generic orbit solution of the boundary condition
    problem through the variable X=log(1+cos(alpha/2)). By doing so the
    graph of the time of flight become defined in the entire real axis and
    resembles a straight line. Convergence is granted within few iterations
    for all the possible geometries (except, of course, when the transfer
    angle is zero). When multiple revolutions are considered the variable is
    X=tan(cos(alpha/2)*pi/2).

    2) Once the orbit has been determined in the plane, this routine
    evaluates the velocity vectors at the two points in a way that is not
    singular for the transfer angle approaching to pi (Lagrange coefficient
    based methods are numerically not well suited for this purpose).

    As a result Lambert's problem is solved (with multiple revolutions
    being accounted for) with the same computational effort for all
    possible geometries. The case of near 180 transfers is also solved
    efficiently.

    We note here that even when the transfer angle is exactly equal to pi
    the algorithm does solve the problem in the plane (it finds X), but it
    is not able to evaluate the plane in which the orbit lies. A solution
    to this would be to provide the direction of the plane containing the
    transfer orbit from outside. This has not been implemented in this
    routine since such a direction would depend on which application the
    transfer is going to be used in.

    please report bugs to dario.izzo@esa.int
    =#

    # adjusted documentation:
    #=
    By default, the short-way solution is computed. The long way solution
    may be requested by giving a negative value to the corresponding
    time-of-flight [tf].

    For problems with |m| > 0, there are generally two solutions. By
    default, the right branch solution will be returned. The left branch
    may be requested by giving a negative value to the corresponding
    number of complete revolutions [m].
    =#

    # Authors
    # .-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.-`-.
    # Name       : Dr. Dario Izzo
    # E-mail     : dario.izzo@esa.int
    # Affiliation: ESA / Advanced Concepts Team (ACT)

    # Made more readible and optimized for speed by Rody P.S. Oldenhuis
    # Code available in MGA.M on   http://www.esa.int/gsp/ACT/inf/op/globopt.htm

    # last edited 12/Dec/2009

    # ADJUSTED FOR EML-COMPILATION 24/Dec/2009

    # initial values
    tol = 1e-14;    bad = false;     days = 86400;

    # work with non-dimensional units
    r1 = sqrt(dot(r1vec, r1vec));  r1vec = r1vec/r1;
    V = sqrt(muC/r1);          r2vec = r2vec/r1;
    T = r1/V;                  tf    = tf*days/T; # also transform to seconds

    # relevant geometry parameters (non dimensional)
    mr2vec = sqrt(dot(r2vec, r2vec));
    # make 100# sure it's in (-1 <= dth <= +1)
    dth = acos( max(-1, min(1, (dot(r1vec, r2vec))/mr2vec)) );

    # decide whether to use the left or right branch (for multi-revolution
    # problems), and the long- or short way
    leftbranch = sign(m);   longway = sign(tf);
    m = abs(m);             tf = abs(tf);
    if (longway < 0)
        dth = 2*pi - dth
    end

    # derived quantities
    c      = sqrt(1 + mr2vec^2 - 2*mr2vec*cos(dth)); # non-dimensional chord
    s      = (1 + mr2vec + c)/2;                     # non-dimensional semi-perimeter
    a_min  = s/2;                                    # minimum energy ellipse semi major axis
    Lambda = sqrt(mr2vec)*cos(dth/2)/s;              # lambda parameter (from BATTIN's book)
    crossprd = [r1vec[2]*r2vec[3] - r1vec[3]*r2vec[2],
                r1vec[3]*r2vec[1] - r1vec[1]*r2vec[3], # non-dimensional normal vectors
                r1vec[1]*r2vec[2] - r1vec[2]*r2vec[1]];
    mcr       = sqrt(dot(crossprd, crossprd));           # magnitues thereof
    nrmunit   = crossprd/mcr;                        # unit vector thereof

    # Initial values
    # ---------------------------------------------------------

    # ELMEX requires this variable to be declared OUTSIDE the IF-statement
    logt = log(tf); # avoid re-computing the same value

    # single revolution (1 solution)
    if (m == 0)

        # initial values
        inn1 = -0.5233;      # first initial guess
        inn2 = +0.5233;      # second initial guess
        x1   = log(1 + inn1);# transformed first initial guess
        x2   = log(1 + inn2);# transformed first second guess

        # multiple revolutions (0, 1 or 2 solutions)
        # the returned soltuion depends on the sign of [m]
    else
        # select initial values
        if (leftbranch < 0)
            inn1 = -0.5234; # first initial guess, left branch
            inn2 = -0.2234; # second initial guess, left branch
        else
            inn1 = +0.7234; # first initial guess, right branch
            inn2 = +0.5234; # second initial guess, right branch
        end
        x1 = tan(inn1*pi/2);# transformed first initial guess
        x2 = tan(inn2*pi/2);# transformed first second guess
    end

    # since (inn1, inn2) < 0, initial estimate is always ellipse
    xx   = [inn1, inn2];  aa = a_min ./ (1 .- xx.^2);
    bbeta = longway * 2*asin.(sqrt.((s-c) / 2 ./ aa));
    # make 100.4# sure it's in (-1 <= xx <= +1)
    aalfa = 2*acos.(  max.(-1, min.(1, xx)) );

    # evaluate the time of flight via Lagrange expression
    y12  = @. aa.*sqrt(aa).*((aalfa - sin(aalfa)) - (bbeta-sin(bbeta)) + 2*pi*m);

    # initial estimates for y
    if m == 0
        y1 = log(y12[1]) - logt;
        y2 = log(y12[2]) - logt;
    else
        y1 = y12[1] - tf;
        y2 = y12[2] - tf;
    end

    # Solve for x
    # ---------------------------------------------------------

    # Newton-Raphson iterations
    # NOTE - the number of iterations will go to infinity in case
    # m > 0  and there is no solution. Start the other routine in
    # that case
    err = Inf;  iterations = 0; xnew = 0;
    while (err > tol)
        # increment number of iterations
        iterations = iterations + 1;
        # new x
        xnew = (x1*y2 - y1*x2) / (y2-y1);
        # copy-pasted code (for performance)
        if m == 0
            x = exp(xnew) - 1;
        else
            x = atan(xnew)*2/pi;
        end

        a = a_min/(1 - x^2);
        if (x < 1) # ellipse
            beta = longway * 2*asin(sqrt((s-c)/2/a));
            # make 100.4# sure it's in (-1 <= xx <= +1)
            alfa = 2*acos( max(-1, min(1, x)) );
        else # hyperbola
            alfa = 2*acosh(x);
            beta = longway * 2*asinh(sqrt((s-c)/(-2*a)));
        end
        # evaluate the time of flight via Lagrange expression
        if (a > 0)
            tof = a*sqrt(a)*((alfa - sin(alfa)) - (beta-sin(beta)) + 2*pi*m);
        else
            tof = -a*sqrt(-a)*((sinh(alfa) - alfa) - (sinh(beta) - beta));
        end
        # new value of y
        if m ==0
            ynew = log(tof) - logt;
        else
            ynew = tof - tf;
        end

        # save previous and current values for the next iterarion
        # (prevents getting stuck between two values)
        x1 = x2;  x2 = xnew;
        y1 = y2;  y2 = ynew;
        # update error
        err = abs(x1 - xnew);
        # escape clause
        if (iterations > 15)
            bad = true; break;
        end
    end

    # If the Newton-Raphson scheme failed, try to solve the problem
    # with the other Lambert targeter.
    if bad
        # NOTE: use the original, UN-normalized quantities
        _branch = leftbranch > 0 ? :left : :right
        _traj   = longway > 0 ? :long : :short
        _m = m

        V1, V2, extremal_distances = lambert_lancaster_blanchard(
            r1vec*r1, r2vec*r1, tf*T, muC;
            revolutions = _m,
            branch = _branch,
            trajectory = _traj,
            output_extrema = Val{true});
        return V1, V2, extremal_distances, 1
    end

    # convert converged value of x
    if m==0
        x = exp(xnew) - 1;
    else
        x = atan(xnew)*2/pi;
    end

    #=
    The solution has been evaluated in terms of log(x+1) or tan(x*pi/2), we
    now need the conic. As for transfer angles near to pi the Lagrange-
    coefficients technique goes singular (dg approaches a zero/zero that is
    numerically bad) we here use a different technique for those cases. When
    the transfer angle is exactly equal to pi, then the ih unit vector is not
    determined. The remaining equations, though, are still valid.
    =#

    # Solution for the semi-major axis
    a = a_min/(1-x^2);

    # Calculate psi
    if (x < 1) # ellipse
        beta = longway * 2*asin(sqrt((s-c)/2/a));
        # make 100.4# sure it's in (-1 <= xx <= +1)
        alfa = 2*acos( max(-1, min(1, x)) );
        psi  = (alfa-beta)/2;
        eta2 = 2*a*sin(psi)^2/s;
        eta  = sqrt(eta2);
    else       # hyperbola
        beta = longway * 2*asinh(sqrt((c-s)/2/a));
        alfa = 2*acosh(x);
        psi  = (alfa-beta)/2;
        eta2 = -2*a*sinh(psi)^2/s;
        eta  = sqrt(eta2);
    end

    # unit of the normalized normal vector
    ih = longway * nrmunit;

    # unit vector for normalized [r2vec]
    r2n = r2vec/mr2vec;

    # cross-products
    # don't use cross() (emlmex() would try to compile it, and this way it
    # also does not create any additional overhead)
    crsprd1 = [ih[2]*r1vec[3]-ih[3]*r1vec[2],
            ih[3]*r1vec[1]-ih[1]*r1vec[3],
            ih[1]*r1vec[2]-ih[2]*r1vec[1]];
    crsprd2 = [ih[2]*r2n[3]-ih[3]*r2n[2],
            ih[3]*r2n[1]-ih[1]*r2n[3],
            ih[1]*r2n[2]-ih[2]*r2n[1]];

    # radial and tangential directions for departure velocity
    Vr1 = 1/eta/sqrt(a_min) * (2*Lambda*a_min - Lambda - x*eta);
    Vt1 = sqrt(mr2vec/a_min/eta2 * sin(dth/2)^2);

    # radial and tangential directions for arrival velocity
    Vt2 = Vt1/mr2vec;
    Vr2 = (Vt1 - Vt2)/tan(dth/2) - Vr1;

    # terminal velocities
    V1 = (Vr1*r1vec + Vt1*crsprd1)*V;
    V2 = (Vr2*r2n + Vt2*crsprd2)*V;

    # exitflag
    exitflag = 1; # (success)

    # also compute minimum distance to central body
    # NOTE: use un-transformed vectors again!
    extremal_distances = minmax_distances(r1vec*r1, r2vec*r1, r1, mr2vec*r1, dth, a*r1, V1, V2, m, muC);

    return V1, V2, extremal_distances, exitflag

end

"""
Wrapper for Unitful inputs.
"""
function lambert(r1::AbstractVector{<:Unitful.Length}, r2::AbstractVector{<:Unitful.Length}, tf::Unitful.Time, m::Integer, mu::MassParameter)
    v1, v2, ex, fl = lambert(
        ustrip.(u"km", r1),
        ustrip.(u"km", r2),
        ustrip.(u"d", tf),
        m,
        ustrip.(u"km^3/s^2", mu)
    )

    return v1 .* u"km/s", v2 .* u"km/s", ex .* u"km", fl
end


text:
using FixedPointNumbers
using Base.Test
using Compat
import Compat.String

@test reinterpret(N0f8, 0xa2).i  === 0xa2
@test reinterpret(N6f10, 0x1fa2).i === 0x1fa2
@test reinterpret(N4f12, 0x1fa2).i === 0x1fa2
@test reinterpret(N2f14, 0x1fa2).i === 0x1fa2
@test reinterpret(N0f16, 0x1fa2).i === 0x1fa2

@test reinterpret(reinterpret(N0f8, 0xa2))    === 0xa2
@test reinterpret(reinterpret(N6f10, 0x00a2)) === 0x00a2
@test reinterpret(reinterpret(N4f12, 0x00a2)) === 0x00a2
@test reinterpret(reinterpret(N2f14, 0x00a2)) === 0x00a2
@test reinterpret(reinterpret(N0f16, 0x00a2)) === 0x00a2

@test 0.635N0f8   == N0f8(0.635)
@test 0.635N6f10 == N6f10(0.635)
@test 0.635N4f12 == N4f12(0.635)
@test 0.635N2f14 == N2f14(0.635)
@test 0.635N0f16 == N0f16(0.635)

@test N0f8(1.0) == reinterpret(N0f8, 0xff)
@test N0f8(0.5) == reinterpret(N0f8, 0x80)
@test N2f14(1.0) == reinterpret(N2f14, 0x3fff)
v = @compat N4f12.([2])
@test v == N4f12[reinterpret(N4f12, 0x1ffe)]
@test isa(v, Vector{N4f12})

UF2 = (Normed{UInt32,16}, Normed{UInt64,3}, Normed{UInt64,51}, Normed{UInt128,7}, Normed{UInt128,51})

for T in (FixedPointNumbers.UF..., UF2...)
    @test zero(T) == 0
    @test one(T) == 1
    @test one(T) * one(T) == one(T)
    @test typemin(T) == 0
    @test realmin(T) == 0
    @test eps(zero(T)) == eps(typemax(T))
    @test sizeof(T) == sizeof(FixedPointNumbers.rawtype(T))
end
@test typemax(N0f8) == 1
@test typemax(N6f10) == typemax(UInt16)//(2^10-1)
@test typemax(N4f12) == typemax(UInt16)//(2^12-1)
@test typemax(N2f14) == typemax(UInt16)//(2^14-1)
@test typemax(N0f16) == 1
@test typemax(N6f10) == typemax(UInt16) // (2^10-1)
@test typemax(N4f12) == typemax(UInt16) // (2^12-1)
@test typemax(N2f14) == typemax(UInt16) // (2^14-1)
@test typemax(Normed{UInt32,16}) == typemax(UInt32) // (2^16-1)
@test typemax(Normed{UInt64,3}) == typemax(UInt64) // (2^3-1)
@test typemax(Normed{UInt128,7}) == typemax(UInt128) // (2^7-1)
@test typemax(Normed{UInt128,100}) == typemax(UInt128) // (UInt128(2)^100-1)

# TODO: change back to InexactError when it allows message strings
@test_throws ArgumentError N0f8(2)
@test_throws ArgumentError N0f8(255)
@test_throws ArgumentError N0f8(0xff)
@test_throws ArgumentError N0f16(2)
@test_throws ArgumentError N0f16(0xff)
@test_throws ArgumentError N0f16(0xffff)
@test_throws ArgumentError convert(N0f8,  typemax(N6f10))
@test_throws ArgumentError convert(N0f16, typemax(N6f10))
@test_throws ArgumentError convert(Normed{UInt128,100}, 10^9)
@test_throws ArgumentError convert(Normed{UInt128,100}, 10.0^9)

x = N0f8(0.5)
@test convert(N0f8, x) === x
@test isfinite(x) == true
@test isnan(x) == false
@test isinf(x) == false

@test convert(N0f8,  1.1/typemax(UInt8)) == eps(N0f8)
@test convert(N6f10, 1.1/typemax(UInt16)*64) == eps(N6f10)
@test convert(N4f12, 1.1/typemax(UInt16)*16) == eps(N4f12)
@test convert(N2f14, 1.1/typemax(UInt16)*4)  == eps(N2f14)
@test convert(N0f16, 1.1/typemax(UInt16))    == eps(N0f16)
@test convert(Normed{UInt32,16}, 1.1/typemax(UInt32)*2^16) == eps(Normed{UInt32,16})
@test convert(Normed{UInt64,3},  1.1/typemax(UInt64)*UInt64(2)^61)  == eps(Normed{UInt64,3})
@test convert(Normed{UInt128,7}, 1.1/typemax(UInt128)*UInt128(2)^121) == eps(Normed{UInt128,7})

@test convert(N0f8,  1.1f0/typemax(UInt8)) == eps(N0f8)

@test convert(Float64, eps(N0f8)) == 1/typemax(UInt8)
@test convert(Float32, eps(N0f8)) == 1.0f0/typemax(UInt8)
@test convert(BigFloat, eps(N0f8)) == BigFloat(1)/typemax(UInt8)
for T in (FixedPointNumbers.UF..., UF2...)
    @test convert(Bool, zero(T)) == false
    @test convert(Bool, one(T))  == true
    @test convert(Bool, convert(T, 0.2)) == true
    @test convert(Int, one(T)) == 1
    @test convert(Rational, one(T)) == 1
end
@test convert(Rational, convert(N0f8, 0.5)) == 0x80//0xff
@test convert(N0f16, one(N0f8)) === one(N0f16)
@test convert(N0f16, N0f8(0.5)).i === 0x8080
@test convert(Normed{UInt16,7}, Normed{UInt8,7}(0.504)) === Normed{UInt16,7}(0.504)

@test  N0f8(0.2) % N0f8  === N0f8(0.2)
@test N2f14(1.2) % N0f16 === N0f16(0.20002)
@test N2f14(1.2) % N0f8  === N0f8(0.196)

for i = 0.0:0.1:1.0
    @test i % N0f8 === N0f8(i)
end
@test ( 1.5 % N0f8).i == round(Int,  1.5*255) % UInt8
@test (-0.3 % N0f8).i == round(Int, -0.3*255) % UInt8

for i = 0.0:0.1:64.0
    @test i % N6f10 === N6f10(i)
end
@test (65.2 % N6f10).i == round(Int, 65.2*1023) % UInt16
@test (-0.3 % N6f10).i == round(Int, -0.3*1023) % UInt16

@test 1 % N0f8 == 1
@test 2 % N0f8 == N0f8(0.996)

x = N0f8(0b01010001, 0)
@test ~x == N0f8(0b10101110, 0)
@test -x == reinterpret(N0f8, 0xaf)

@test isa(float(one(Normed{UInt8,7})),   Float32)
@test isa(float(one(Normed{UInt32,18})), Float64)
@test isa(float(one(Normed{UInt32,25})), Float64)

for T in (FixedPointNumbers.UF..., UF2...)
    x = T(0x10,0)
    y = T(0x25,0)
    fx = float(x)
    fy = float(y)
    @test y > x
    @test y != x
    @test typeof(x+y) == T
    @test typeof((x+y)-y) == T
    @test typeof(x*y) == T
    @test typeof(x/y) == T
    @test (x+y) ≈ T(0x35,0)
    @test ((x+y)-x) ≈ fy
    @test ((x-y)+y) ≈ fx
    @test (x*y)  ≈ convert(T, fx*fy)
    @test (x/y)  ≈ convert(T, fx/fy)
    @test (x^2) ≈ convert(T, fx^2)
    @test (x^2.1f0) ≈ fx^2.1f0
    @test (x^2.1) ≈ convert(Float64, x)^2.1
end

function testtrunc{T}(inc::T)
    incf = convert(Float64, inc)
    tm = reinterpret(typemax(T))/reinterpret(one(T))
    x = zero(T)
    for i = 0 : min(1e6, reinterpret(typemax(T))-1)
        xf = incf*i
        try
            @test typeof(trunc(x)) == T
            @test trunc(x) == trunc(xf)
            @test typeof(round(x)) == T
            @test round(x) == round(xf)
            cxf = ceil(xf)
            if cxf < tm
                @test typeof(ceil(x)) == T
                @test ceil(x) == ceil(xf)
            end
            @test typeof(floor(x)) == T
            @test floor(x) == floor(xf)
            @test trunc(Int,x) == trunc(Int,xf)
            @test round(Int,x) == round(Int,xf)
            @test floor(Int,x) == floor(Int,xf)
            if cxf < tm
                @test ceil(Int,x) == ceil(Int,xf)
            end
        catch err
            println("Failed on x = ", x, ", xf = ", xf)
            rethrow(err)
        end
        x = convert(T, x+inc)
    end
end

for T in (FixedPointNumbers.UF..., UF2...)
    testtrunc(eps(T))
end

function testapprox{T}(::Type{T})
    for x = typemin(T):eps(T):typemax(T)-eps(T)
        y = x+eps(T)
        @test x ≈ y
        @test y ≈ x
        @test !(x ≈ y+eps(T))
    end
end
for T in FixedPointNumbers.UF
    testapprox(T)
end

@test !(N0f8(0.5) < N0f8(0.5))
@test N0f8(0.5) <= N0f8(0.5)

@test div(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == fld(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == 8
@test div(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == fld(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == 7
@test Base.fld1(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == 8
@test Base.fld1(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == 8
@test mod(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == rem(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == 0
@test mod(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == rem(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == reinterpret(N0f8, 0x01)
@test mod1(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == reinterpret(N0f8, 0x02)
@test mod1(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == reinterpret(N0f8, 0x01)

r = reinterpret(N0f8, 0x01):reinterpret(N0f8, 0x01):reinterpret(N0f8, convert(UInt8, 48))
@test length(r) == 48

counter = 0
for x in N0f8(0):eps(N0f8):N0f8(1)
    counter += 1
end
@test counter == 256

# Promotion within Normed
@test @inferred(promote(N0f8(0.2), N0f8(0.8))) ===
    (N0f8(0.2), N0f8(0.8))
@test @inferred(promote(Normed{UInt16,3}(0.2), Normed{UInt8,3}(0.86))) ===
    (Normed{UInt16,3}(0.2), Normed{UInt16,3}(0.86))
@test @inferred(promote(Normed{UInt8,7}(0.197), Normed{UInt8,4}(0.8))) ===
    (Normed{UInt16,7}(0.197), Normed{UInt16,7}(0.8))

@test Normed{UInt16,16}(1)   == Normed{UInt8,8}(1)
@test Normed{UInt16,16}(0.2) == Normed{UInt8,8}(0.2)
@test Normed{UInt16,8}(1)   == Normed{UInt8,8}(1)
@test Normed{UInt16,8}(0.2) == Normed{UInt8,8}(0.2)
@test Normed{UInt16,16}(1)   == Normed{UInt8,6}(1)
@test Normed{UInt16,16}(0.20635) == Normed{UInt8,6}(0.20635)
@test Normed{UInt16,4}(1)   == Normed{UInt8,6}(1)
@test Normed{UInt16,4}(0.2) == Normed{UInt8,6}(0.2)

@test promote_type(N0f8,Float32,Int) == Float32
@test promote_type(N0f8,Int,Float32) == Float32
@test promote_type(Int,N0f8,Float32) == Float32
@test promote_type(Int,Float32,N0f8) == Float32
@test promote_type(Float32,Int,N0f8) == Float32
@test promote_type(Float32,N0f8,Int) == Float32

# Show
x = reinterpret(N0f8, 0xaa)
iob = IOBuffer()
show(iob, x)
str = String(take!(iob))
@test str == "0.667N0f8"
@test eval(parse(str)) == x

# scaledual
function generic_scale!(C::AbstractArray, X::AbstractArray, s::Number)
    length(C) == length(X) || error("C must be the same length as X")
    for i = 1:length(X)
        @inbounds C[i] = X[i]*s
    end
    C
end

a = rand(UInt8, 10)
rfloat = similar(a, Float32)
rfixed = similar(rfloat)
af8 = reinterpret(N0f8, a)

b = 0.5
bd, eld = scaledual(b, af8[1])
@assert b*a[1] == bd*eld

b, ad = scaledual(0.5, a)
@test b == 0.5
@test ad == a
b, ad = scaledual(0.5, ad)
generic_scale!(rfloat, a, 0.5)
generic_scale!(rfixed, ad, b)
@test rfloat == rfixed

# reductions
a = N0f8[reinterpret(N0f8, 0xff), reinterpret(N0f8, 0xff)]
@test sum(a) == 2.0
@test sum(a, 1) == [2.0]

a = N2f14[3.2, 2.4]
acmp = Float64(a[1])*Float64(a[2])
@test prod(a) == acmp
@test prod(a, 1) == [acmp]

x = N0f8(0.3)
for T in (Float16, Float32, Float64, BigFloat)
    y = convert(T, x)
    @test isa(y, T)
end

for T in (Normed{UInt8,8}, Normed{UInt8,6},
          Normed{UInt16,16}, Normed{UInt16,14},
          Normed{UInt32,32}, Normed{UInt32,30},
          Normed{UInt64,64}, Normed{UInt64,62})
    a = rand(T)
    @test isa(a, T)
    a = rand(T, (3, 5))
    @test isa(a, Array{T,2})
    @test size(a) == (3,5)
end

ext: jl


text:
using SparseArrays

struct InitialParams
    centers::Matrix{Float64}
    covs::Array{CovMat,1};
    assignment::Vector{Int64};

    InitialParams(centers::Array{T, 2} where T <: Real, covs::Float64, assignment::Array{Int64,1}) =
        new(deepcopy(centers), [CovMat([covs 0.; 0. covs]) for i in 1:size(centers, 1)], deepcopy(assignment))

    InitialParams(centers::Array{T, 2} where T <: Real, covs::Array{CovMat, 1}, assignment::Array{Int64,1}) =
        new(deepcopy(centers), deepcopy(covs), deepcopy(assignment))
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct RoleType <: AbstractRoleType
    label::Label
    is_root::Bool
end

# Porting note: Java client calls into RequstBuilder but it really
# has nothing to to with requests... I think it's probably better
# migrating the function here.
function proto(t::AbstractRoleType)
    Proto._Type(
        label = t.label.name,
        encoding = encoding(t)
    )
    # return RoleTypeRequestBuilder.proto_role_type(t.label, encoding(t))
end


text:
# List of available algorithms for AIBECS
algs = [CTKAlg]

@testset "Solvers" begin
    nt = length(T_all)
    n = nt * nb
    @unpack xgeo = p
    x = xgeo * ones(n)
    @testset "$(string(alg))" for alg in algs
        @test alg <: DiffEqBase.AbstractSteadyStateAlgorithm
        @testset "$(nx)x, $(np)p" for nx in 1:2, np in 1:2
            testp = AIBECS.reconstruct(typeof(p), np * vec(p))
            prob = SteadyStateProblem(F, ∇ₓF, nx * x, testp)
            s = solve(prob, alg())
            @test s isa SteadyStateSolution
            @test prob isa SteadyStateProblem
            @test norm(s.u) / norm(F(s.u, testp)) > ustrip(upreferred(1e5u"Myr"))
        end
    end
end


text:
import Base: setindex!, getindex

abstract type Node end

struct Variable <: Node
    val::Gaussian
    messages::Dict{Node, Gaussian}
end

function Variable(; mu::Float64=0.0, sigma::Float64=0.0)
    if sigma == 0
        return Variable(Gaussian(0.0, 0.0), Dict())
    end
    pi = sigma^(-2)
    tau = pi * mu
    Variable(Gaussian(pi, tau), Dict())
end

@inline (*)(a::Variable, b::Variable) = Gaussian(a.val.pi + b.val.pi, a.val.tau + b.val.tau)
@inline (*)(a::Variable, b::Gaussian) = Gaussian(a.val.pi + b.pi, a.val.tau + b.tau)
@inline (*)(a::Gaussian, b::Variable) = Gaussian(a.pi + b.val.pi, a.tau + b.val.tau)
@inline (/)(a::Variable, b::Variable) = Gaussian(a.val.pi - b.val.pi, a.val.tau - b.val.tau)
@inline (/)(a::Variable, b::Gaussian) = Gaussian(a.val.pi - b.pi, a.val.tau - b.tau)
@inline (/)(a::Gaussian, b::Variable) = Gaussian(a.pi - b.val.pi, a.tau - b.val.tau)

@inline function mu(r::Variable)
    if r.val.pi != 0
        return r.val.tau / r.val.pi
    else
        return 0
    end
end

@inline function sigma(r::Variable)
    if r.val.pi != 0
        return(sqrt(1/r.val.pi))
    else
        return Inf
    end
end

function setindex!(v::Variable, message, factor::Node)
    v.messages[factor] = message
end
function getindex(v::Variable, factor::Node)
    v.messages[factor]
end

function delta(v::Variable, w::Gaussian)
    pi_delta = abs(v.val.pi - w.pi)
    if pi_delta == Inf
        return 0.
    end
    return max(abs(v.val.tau - w.tau), sqrt(pi_delta))
end

function delta(v::Variable, w::Variable)
    delta(v, w.val)
end

function set(v::Variable, w::Gaussian)
    _delta = delta(v, w)
    v.val.pi, v.val.tau = w.pi, w.tau
    return _delta
end

function set(v::Variable, w::Variable)
    set(v, w.val)
end

function update_message(v::Variable, factor::Node; pi=0, tau=0, message::Union{Nothing, Gaussian}=nothing)
    if message === nothing
        message = Gaussian(pi, tau)
    end
    old_message, v[factor] = v[factor], message
    set(v, v / old_message * message)
end

function update_value(v::Variable, factor::Node; pi=0, tau=0, value::Union{Nothing, Gaussian}=nothing)
    if value === nothing
        value = Gaussian(pi, tau)
    end
    old_message = v[factor]
    v[factor] = value * old_message / v
    return set(v, value)
end

function show(io::IO, v::Variable)
    print(io, "< Variable $(v.val) with $(length(v.messages)) connections >")
end

function show(io::IO, v::Node)
    print(io, "< $(typeof(v))  with $(length(v.vars)) connections >")
end

struct Factor <: Node
    vars::Vector{Variable}
    function Factor(variables::Vector{Variable})
        r = new(variables)
        for var in variables
            var[r] = Gaussian()
        end
        return r
    end
end

@inline down(f::Factor) = 0.0
@inline up(f::Factor) = 0.0

function var(f::Node)
    @assert length(f.vars) == 1
    return f.vars[1]
end

struct PriorFactor <: Node
    vars::Vector{Variable}
    val::Gaussian
    dynamic::Float64
    function PriorFactor(var::Variable, val::GaussianItem, dynamic::Float64=0.0)
        r = new([var], Gaussian(val.pi, val.tau), dynamic)
        var[r] = Gaussian()
        return r
    end
end

function down(f::PriorFactor)
    sig = sqrt(sigma(f.val)^2 + f.dynamic^2)
    value = Gaussian(; mu=mu(f.val), sigma=sig)
    #println("update_value from down(::PriorFactor): $(var(f)), $f; value=$value")
    return update_value(var(f), f; value=value)
end

struct LikelihoodFactor <: Node
    vars::Vector{Variable}
    mean::Variable
    value::Variable
    variance::Float64
    function LikelihoodFactor(mean_var::Variable, value_var::Variable, variance::Float64)
        r = new([mean_var, value_var], mean_var, value_var, variance)
        mean_var[r] = Gaussian()
        value_var[r] = Gaussian()
        return r
    end
end

@inline function calc_a(f::LikelihoodFactor, g::Gaussian)
    return 1 ./ (1. + f.variance * g.pi)
end

@inline function calc_a(f::LikelihoodFactor, var::Variable)
    return 1. / (1. + f.variance * var.val.pi)
end

function down(f::LikelihoodFactor)
    # update value
    msg = f.mean / f.mean[f]
    a = calc_a(f, msg)
    #println("update_message from down(f::LikelihoodFactor): $(f.value), $f, pi = $(a * msg.pi), tau = $(a * msg.tau)")
    return update_message(f.value, f; pi = a * msg.pi, tau = a * msg.tau)
end

function up(f::LikelihoodFactor)
    # update mean
    msg = f.value / f.value[f]
    a = calc_a(f, msg)
    #println("update_message from up(f::LikelihoodFactor): $(f.mean), $f, pi = $(a * msg.pi), tau = $(a * msg.tau)")
    return update_message(f.mean, f; pi = a * msg.pi, tau = a * msg.tau)
end

struct SumFactor <: Node
    vars::Vector{Variable}
    sum::Variable
    terms::Vector{Variable}
    coeffs::Vector{Float64}
    function SumFactor(sum_var::Variable, term_vars::Vector{Variable}, coeffs::Vector{T}) where T <: Real
        r = new(vcat([sum_var], term_vars), sum_var, term_vars, coeffs)
        sum_var[r] = Gaussian()
        for var in term_vars
            var[r] = Gaussian()
        end
        return r
    end
end

function update(f::SumFactor, var, vals, msgs, coeffs)
    pi_inv = 0.
    _mu = 0.
    for (val, msg, coeff) in zip(vals, msgs, coeffs)
        div = val / msg
        _mu += coeff * mu(div)
        if pi_inv == Inf
            continue
        end
        if div.pi == 0
            pi_inv = Inf
        else
            pi_inv += coeff ^ 2 / div.pi
        end
    end

    pi = 1. / pi_inv
    tau = pi * _mu
    #println("update_message from update(::SumFactor, ...): $var, $f; pi=$pi, tau=$tau")
    return update_message(var, f; pi=pi, tau=tau)
end

function down(f::SumFactor)
    vals = f.terms
    msgs = [var[f] for var in vals]
    #println("update from down(::SumFactor): $f, $(f.sum), $vals, $msgs, $(f.coeffs)")
    return update(f, f.sum, vals, msgs, f.coeffs)
end

function up(f::SumFactor; index=1)
    coeff = f.coeffs[index]
    coeffs = Array{Float64}(undef, length(f.coeffs))
    for (i, c)  in enumerate(f.coeffs)
        if coeff == 0
            coeffs[i] = 0
        elseif i == index
            coeffs[i] = 1. / coeff
        else
            coeffs[i] = -c / coeff
        end
    end
    vals = copy(f.terms)
    vals[index] = f.sum
    msgs = [var[f] for var in vals]
    #println("update from up(::SumFactor): $f, $(f.terms[index]), $vals, $msgs, $coeffs")
    return update(f, f.terms[index], vals, msgs, coeffs)
end

struct TruncateFactor <: Node
    vars::Vector{Variable}
    v_func::Function
    w_func::Function
    draw_margin::Float64
    function TruncateFactor(var, v_func, w_func, draw_margin)
        r = new([var], v_func, w_func, draw_margin)
        var[r] = Gaussian()
        return r
    end
end

function up(f::TruncateFactor)
    val = var(f)
    msg = var(f)[f]
    div = val / msg
    sqrt_pi = sqrt(div.pi)
    args = (div.tau / sqrt_pi, f.draw_margin * sqrt_pi)
    v = f.v_func(args...)
    w = f.w_func(args...)
    denom = (1. - w)
    pi, tau = div.pi / denom, (div.tau + sqrt_pi * v) / denom
    #println("update_value from up(::TruncateFactor): $val, $f; pi=$pi, tau=$tau")
    return update_value(val, f; pi=pi, tau=tau)
end


text:
struct QGramIterator{S <: Union{AbstractString, AbstractVector}}
	s::S   # Collection
	q::Int # Length of Qgram
end
Base.length(qgram::QGramIterator) = max(length(qgram.s) - qgram.q + 1, 0)

# q-grams of AbstractString
function Base.iterate(qgram::QGramIterator{<: AbstractString},
	state = (1, nextind(qgram.s, 0, qgram.q)))
	istart, iend = state
	iend > ncodeunits(qgram.s) && return nothing
	element = SubString(qgram.s, istart, iend)
	nextstate = nextind(qgram.s, istart), nextind(qgram.s, iend)
	element, nextstate
end
Base.eltype(qgram::QGramIterator{SubString{S}}) where {S} = SubString{S}
Base.eltype(qgram::QGramIterator{S}) where {S <: AbstractString} = SubString{S}
qgrams(s::AbstractString, q::Integer) = QGramIterator(s, q)

#q-grams of AbstractVector
# Alternatively, I could also use partition in IterTools but it creates a vector for each iteration
# so it does not seem to be worth it.
function Base.iterate(qgram::QGramIterator{<: AbstractVector}, state = firstindex(qgram.s))
	state + qgram.q - 1 > lastindex(qgram.s) && return nothing
	view(qgram.s, state:(state + qgram.q - 1)), state + 1
end
Base.eltype(qgram::QGramIterator{<: AbstractVector}) = typeof(first(qgram))
qgrams(s::AbstractVector, q::Integer) = QGramIterator(s, q)
qgrams(s, q::Integer) = QGramIterator(collect(s), q)

@doc """
Return an iterator corresponding to the the q-gram of an iterator.
When the iterator is a String, qgrams are SubStrings.

### Arguments
* `s` iterator
* `q::Integer`: length of q-gram

## Examples
```julia
for x in qgrams("hello", 2)
	println(x)
end
```
"""
qgrams

# For two iterators s1 and s2, that define a length and eltype method,
# this returns an iterator that,
# for each element in s1 ∪ s2, returns (numbers of times it appears in s1, numbers of times it appears in s2)
function _count(s1, s2)
	K = promote_type(eltype(s1), eltype(s2))
	d = Dict{K, Tuple{Int, Int}}()
	sizehint!(d, length(s1) + length(s2))
	# I use a faster way to change a dictionary key
	# see setindex! in https://github.com/JuliaLang/julia/blob/master/base/dict.jl#L380
	for x1 in s1
		index = Base.ht_keyindex2!(d, x1)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x1
			@inbounds d.vals[index] = (d.vals[index][1] + 1, 0)
		else
			@inbounds Base._setindex!(d, (1, 0), x1, -index)
		end
	end
	for x2 in s2
		index = Base.ht_keyindex2!(d, x2)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x2
			@inbounds d.vals[index] = (d.vals[index][1], d.vals[index][2] + 1)
		else
			@inbounds Base._setindex!(d, (0, 1), x2, -index)
		end
	end
	return values(d)
end

# Turn a sequence of qgrams to a count dict for them, i.e. map each
# qgram to the number of times it has been seen.
function countdict(qgrams)
    d = Dict{eltype(qgrams), Int32}()
    for qg in qgrams
        index = Base.ht_keyindex2!(d, qg)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = qg
			@inbounds d.vals[index] = d.vals[index][1] + 1
		else
			@inbounds Base._setindex!(d, 1, qg, -index)
		end
    end
    d
end

abstract type AbstractQGramCounts{Q,K} end
q(qc::AbstractQGramCounts{Q,K}) where {Q,K} = Q
counts(qc::AbstractQGramCounts) = qc.counts
Base.length(qc::AbstractQGramCounts{Q}) where Q = length(qc.counts) + Q - 1
"""
	QGramDict(s, q::Integer = 2)

Creates a QGramDict that pre-calculates (pre-counts) the qgrams
of a string or stream. This enables faster calculation of QGram
distances.

Note that the qgram length must correspond with the q length used
in the distance.

## Examples
```julia
str1, str2 = "my string", "another string"
qd1 = QGramDict(str1, 2)
qd2 = QGramDict(str2, 2)
evaluate(Overlap(2), qd1, qd2)
```
"""
struct QGramDict{Q,K} <: AbstractQGramCounts{Q,K}
    counts::Dict{K,Int}
end
function QGramDict(s::Union{AbstractString, AbstractVector}, q::Integer = 2)
    @assert q >= 1
    qgs = qgrams(s, q)
    QGramDict{q, eltype(qgs)}(countdict(qgs))
end
QGramDict(s, q::Integer = 2) = QGramDict(collect(s), q)

"""
	QGramSortedVector(s, q::Integer = 2)

Creates a QGramSortedVector that pre-calculates (pre-counts) the
qgrams of a string or stream. This enables faster calculation of
QGram distances.

Since qgrams are sorted in lexicographic order QGram distances can be
calculated even faster than when using a QGramDict. However, the
sorting means that updating the counts after creation is less
efficient. However, for most use cases QGramSortedVector is preferred
over a QgramDict.

Note that the qgram length must correspond with the q length used
in the distance.

## Examples
```julia
str1, str2 = "my string", "another string"
qs1 = QGramSortedVector(str1, 2)
qs2 = QGramSortedVector(str2, 2)
evaluate(Jaccard(2), qs1, qs2)
```
"""
struct QGramSortedVector{Q,K} <: AbstractQGramCounts{Q,K}
    counts::Vector{Pair{K,Int}}
end
function QGramSortedVector(s::Union{AbstractString, AbstractVector}, q::Integer = 2)
    @assert q >= 1
    qgs = qgrams(s, q)
    countpairs = collect(countdict(qgs))
    sort!(countpairs, by = first)
    QGramSortedVector{q, eltype(qgs)}(countpairs)
end
QGramSortedVector(s, q::Integer = 2) = QGramSortedVector(collect(s), q)

# To implement the distances we will count qgram matches
# between strings or pre-calculated AbstractQgramCounts objects.
# The abstract type defines different fallback versions which can be
# specialied by subtypes for best performance.
abstract type AbstractQGramMatchCounter end
@inline countleft!(c::AbstractQGramMatchCounter, qg, n1::Integer) = countleft!(c, n1)
@inline countright!(c::AbstractQGramMatchCounter, qg, n2::Integer) = countright!(c, n2)
@inline countboth!(c::AbstractQGramMatchCounter, qg, n1::Integer, n2::Integer) =
	countboth!(c, n1, n2)
@inline function countboth!(c::AbstractQGramMatchCounter, n1::Integer, n2::Integer)
	countleft!(c, n1)
	countright!(c, n2)
	countshared!(c, n1, n2)
end
@inline countshared!(c::AbstractQGramMatchCounter, qg, n1::Integer, n2::Integer) = countshared!(c, n1, n2)

# Subtypes must implement these methods:
@inline countleft!(c::AbstractQGramMatchCounter, n1::Integer) =
	error("countleft! not implemented for $(typeof(c))")
@inline countright!(c::AbstractQGramMatchCounter, n2::Integer) =
	error("countright! not implemented for $(typeof(c))")

# Subtypes either must overwrite countboth! from above (so it not uses countshared!) or implement:
@inline countshared!(c::AbstractQGramMatchCounter, n1::Integer, n2::Integer) =
	error("countshared! not implemented for $(typeof(c))")

function countmatches!(mc::AbstractQGramMatchCounter, d1::Vector{Pair{K,I}}, d2::Vector{Pair{K,I}}) where {K,I<:Integer}
    i1 = i2 = 1
    while i1 <= length(d1) || i2 <= length(d2)
        if i2 > length(d2)
			for i in i1:length(d1)
				@inbounds countleft!(mc, d1[i][1], d1[i][2])
            end
            return
        elseif i1 > length(d1)
			for i in i2:length(d2)
				@inbounds countright!(mc, d2[i][1], d2[i][2])
            end
            return
        end
        @inbounds k1, n1 = d1[i1]
        @inbounds k2, n2 = d2[i2]
        cmpval = Base.cmp(k1, k2)
		if cmpval == -1 # k1 < k2
			countleft!(mc, k1, n1)
            i1 += 1
        elseif cmpval == +1 # k2 < k1
			countright!(mc, k2, n2)
            i2 += 1
		else
			countboth!(mc, k1, n1, n2)
            i1 += 1
            i2 += 1
        end
    end
end

function countmatches!(mc::AbstractQGramMatchCounter, d1::Dict{K,I}, d2::Dict{K,I}) where {K,I<:Integer}
    for (k1, c1) in d1
        index = Base.ht_keyindex2!(d2, k1)
		if index > 0
			countboth!(mc, k1, c1, d2.vals[index])
		else
			countleft!(mc, k1, c1)
        end
    end
    for (k2, c2) in d2
        index = Base.ht_keyindex2!(d1, k2)
		if index <= 0
			countright!(mc, k2, c2)
        end
    end
end

abstract type QGramDistance <: SemiMetric end

function (dist::QGramDistance)(s1, s2)
	((s1 === missing) | (s2 === missing)) && return missing
	counter = newcounter(dist)
	for (n1, n2) in _count(qgrams(s1, dist.q), qgrams(s2, dist.q))
		countboth!(counter, n1, n2)
	end
	calculate(dist, counter)
end

function (dist::QGramDistance)(qc1::QC, qc2::QC) where {QC<:AbstractQGramCounts}
    @assert dist.q == q(qc1)
	@assert dist.q == q(qc2)
	counter = newcounter(dist)
	countmatches!(counter, counts(qc1), counts(qc2))
    calculate(dist, counter)
end

"""
	QGram(q::Int)

Creates a QGram distance.

The distance corresponds to

``||v(s1, q) - v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q,
that contains the number of times a q-gram appears for the string s
"""
struct QGram <: QGramDistance
	q::Int
end

mutable struct SingleCounter{T, QD<:QGramDistance} <: AbstractQGramMatchCounter
	n::T
end

newcounter(d::QGram) = SingleCounter{Int, QGram}(0)

@inline countleft!(c::SingleCounter{Int, QGram}, n1::Integer) = c.n += n1 # n1 === abs(n1 - 0)
@inline countright!(c::SingleCounter{Int, QGram}, n2::Integer) = c.n += n2 # n2 === abs(0 - n2)
@inline countboth!(c::SingleCounter{Int, QGram}, n1::Integer, n2::Integer) = c.n += abs(n1 - n2)

calculate(dist::QGram, c::SingleCounter{Int, QGram}) = c.n

"""
	Cosine(q::Int)

Creates a Cosine distance.

The distance corresponds to

`` 1 - v(s1, q).v(s2, q)  / ||v(s1, q)|| * ||v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q,
that contains the  number of times a q-gram appears for the string s
"""
struct Cosine <: QGramDistance
	q::Int
end

mutable struct ThreeCounters{T, QD<:QGramDistance} <: AbstractQGramMatchCounter
	left::T
	right::T
	shared::T
end

newcounter(d::Cosine) = ThreeCounters{Int, Cosine}(0, 0, 0)

@inline countleft!(c::ThreeCounters{Int, Cosine}, n1::Integer) = c.left += n1^2
@inline countright!(c::ThreeCounters{Int, Cosine}, n2::Integer) = c.right += n2^2
@inline countshared!(c::ThreeCounters{Int, Cosine}, n1::Integer, n2::Integer) = c.shared += n1 * n2

calculate(d::Cosine, c::ThreeCounters{Int, Cosine}) =
	1.0 - c.shared / (sqrt(c.left) * sqrt(c.right))

"""
	Jaccard(q::Int)

Creates a Jaccard distance.

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)| / |Q(s1, q) ∪ Q(s2, q))|``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Jaccard <: QGramDistance
	q::Int
end

calculate(d::Jaccard, c::ThreeCounters{Int, Jaccard}) =
	1.0 - c.shared / (c.left + c.right - c.shared)

"""
	SorensenDice(q::Int)

Creates a SorensenDice distance.

The distance corresponds to

``1 - 2 * |Q(s1, q) ∩ Q(s2, q)|  / (|Q(s1, q)| + |Q(s2, q))|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct SorensenDice <: QGramDistance
	q::Int
end

calculate(d::SorensenDice, c::ThreeCounters{Int, SorensenDice}) =
	1.0 - 2.0 * c.shared / (c.left + c.right)

"""
	Overlap(q::Int)

Creates a Overlap distance.

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)|  / min(|Q(s1, q)|, |Q(s2, q)|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Overlap <: QGramDistance
	q::Int
end

const IntersectionDist = Union{Jaccard, SorensenDice, Overlap}
newcounter(d::IntersectionDist) = ThreeCounters{Int, typeof(d)}(0, 0, 0)

@inline countleft!(c::ThreeCounters{Int, QD}, n1::Integer) where {QD<:IntersectionDist} =
	c.left += (n1 > 0)
@inline countright!(c::ThreeCounters{Int, QD}, n2::Integer) where {QD<:IntersectionDist} =
	c.right += (n2 > 0)
@inline countshared!(c::ThreeCounters{Int, QD}, n1::Integer, n2::Integer) where {QD<:IntersectionDist} =
	c.shared += (n1 > 0) & (n2 > 0)

calculate(d::Overlap, c::ThreeCounters{Int, Overlap}) =
	1.0 - c.shared / min(c.left, c.right)

"""
	MorisitaOverlap(q::Int)

Creates a MorisitaOverlap distance, a general, statistical measure of
dispersion which can also be used on dictionaries such as created
from q-grams. See https://en.wikipedia.org/wiki/Morisita%27s_overlap_index
This is more fine-grained than many of the other QGramDistances since
it is based on the counts per q-gram rather than only which q-grams are
in the strings.

The distance corresponds to

``(2 * sum(m(s1) .* m(s2)) / (sum(m(s1).^2)*M(s2)/M(s1) + sum(m(s2).^2)*M(s1)/M(s2))``

where ``m(s)`` is the vector of q-gram counts for string ``s`` and ``M(s)`` is the
sum of those counts.
"""
struct MorisitaOverlap <: QGramDistance
	q::Int
end

mutable struct FiveCounters{T, QD<:QGramDistance} <: AbstractQGramMatchCounter
	leftsum::T    # sum(m(s1))
	rightsum::T   # sum(m(s2))
	leftsq::T     # sum(m(s1).^2)
	rightsq::T    # sum(m(s2).^2)
	shared::T     # sum(m(s1) .* m(s2))
end

newcounter(d::MorisitaOverlap) = FiveCounters{Int, MorisitaOverlap}(0, 0, 0, 0, 0)

@inline function countleft!(c::FiveCounters{Int, MorisitaOverlap}, n1::Integer)
	c.leftsum += n1
	c.leftsq += (n1^2)
end

@inline function countright!(c::FiveCounters{Int, MorisitaOverlap}, n2::Integer)
	c.rightsum += n2
	c.rightsq += (n2^2)
end

@inline countshared!(c::FiveCounters{Int, MorisitaOverlap}, n1::Integer, n2::Integer) =
	c.shared += (n1 * n2)

calculate(d::MorisitaOverlap, c::FiveCounters{Int, MorisitaOverlap}) =
	1.0 - ((2 * c.shared) / (c.leftsq*c.rightsum/c.leftsum + c.rightsq*c.leftsum/c.rightsum))

"""
	NMD(q::Int)

Creates a NMD (Normalized Multiset Distance) as introduced by Besiris and
Zigouris 2013. The goal with this distance is to behave similarly to a normalized
compression distance without having to do any actual compression (and thus being
faster to compute).

The distance corresponds to

``(sum(max.(m(s1), m(s2)) - min(M(s1), M(s2))) / max(M(s1), M(s2))``

where ``m(s)`` is the vector of q-gram counts for string ``s`` and ``M(s)`` is the
sum of those counts.

For details see:
https://www.sciencedirect.com/science/article/pii/S1047320313001417
"""
struct NMD <: QGramDistance
	q::Int
end

newcounter(d::NMD) = ThreeCounters{Int, NMD}(0, 0, 0)

@inline function countleft!(c::ThreeCounters{Int, NMD}, n1::Integer)
	c.left += n1
	c.shared += n1 # max(n1, 0) == n1
end

@inline function countright!(c::ThreeCounters{Int, NMD}, n2::Integer)
	c.right += n2
	c.shared += n2 # max(n2, 0) == n2
end

@inline function countboth!(c::ThreeCounters{Int, NMD}, n1::Integer, n2::Integer)
	c.left += n1
	c.right += n2
	c.shared += max(n1, n2)
end

calculate(d::NMD, c::ThreeCounters{Int, NMD}) =
	(c.shared - min(c.left, c.right)) / max(c.left, c.right)


text:
using JLD
data = readdlm("simdata")
save("CTSV.jld", "simdata", data)


text:
using Conda, Pkg

const condapkgs = [
	"transformers",
	"emoji",
	"sentencepiece"
]

function build()
	Conda.add_channel("conda-forge")
	foreach(Conda.add, condapkgs)
	ENV["PYTHON"] = joinpath(Conda.PYTHONDIR, "python")
	Pkg.build("PyCall")
end

build()


text:
using NNlib
using Base.Test

@testset "NNlib" begin

include("activation.jl")

end


text:
"""
    struct OptimizerWithAttributes
        optimizer_constructor
        params::Vector{Pair{AbstractOptimizerAttribute,<:Any}}
    end

Object grouping an optimizer constructor and a list of optimizer attributes.
Instances are created with [`instantiate`](@ref).
"""
struct OptimizerWithAttributes
    # Function that takes zero arguments and returns a new optimizer.
    # The type of the function could be
    # * `Function`: a function, or
    # * `Type`: a type, or
    # * `UnionAll`: a type with missing parameters.
    optimizer_constructor::Any
    params::Vector{Pair{AbstractOptimizerAttribute,Any}}
end

_to_param(param::Pair{<:AbstractOptimizerAttribute}) = param
function _to_param(param::Pair{String})
    return RawOptimizerAttribute(param.first) => param.second
end
function _to_param(param::Pair)
    return error(
        "Expected an optimizer attribute or a string, got `$(param.first)` which is a `$(typeof(param.first))`.",
    )
end

"""
    OptimizerWithAttributes(optimizer_constructor, params::Pair...)

Create an [`OptimizerWithAttributes`](@ref) with the parameters `params`.
"""
function OptimizerWithAttributes(
    optimizer_constructor,
    args::Vararg{Pair,N},
) where {N}
    if !applicable(optimizer_constructor)
        error(_INSTANTIATE_NOT_CALLABLE_MESSAGE)
    end
    params =
        Pair{AbstractOptimizerAttribute,Any}[_to_param(arg) for arg in args]
    return OptimizerWithAttributes(optimizer_constructor, params)
end

const _INSTANTIATE_NOT_CALLABLE_MESSAGE =
    "The provided `optimizer_constructor` is invalid. It must be callable with zero " *
    "arguments. For example, \"Ipopt.Optimizer\" or " *
    "\"() -> ECOS.Optimizer()\". It should not be an instantiated optimizer " *
    "like \"Ipopt.Optimizer()\" or \"ECOS.Optimizer()\". " *
    "(Note the difference in parentheses!)"

"""
    _instantiate_and_check(optimizer_constructor)

Create an instance of optimizer by calling `optimizer_constructor`.
Then check that the type returned is an empty [`AbstractOptimizer`](@ref).
"""
function _instantiate_and_check(optimizer_constructor)
    if !applicable(optimizer_constructor)
        error(_INSTANTIATE_NOT_CALLABLE_MESSAGE)
    end
    optimizer = optimizer_constructor()
    if !isa(optimizer, AbstractOptimizer)
        error(
            "The provided `optimizer_constructor` returned an object of type " *
            "$(typeof(optimizer)). Expected a " *
            "MathOptInterface.AbstractOptimizer.",
        )
    end
    if !is_empty(optimizer)
        error(
            "The provided `optimizer_constructor` returned a non-empty optimizer.",
        )
    end
    return optimizer
end

"""
    _instantiate_and_check(optimizer_constructor::OptimizerWithAttributes)

Create an instance of optimizer represented by [`OptimizerWithAttributes`](@ref).
Then check that the type returned is an empty [`AbstractOptimizer`](@ref).
"""
function _instantiate_and_check(optimizer_constructor::OptimizerWithAttributes)
    optimizer =
        _instantiate_and_check(optimizer_constructor.optimizer_constructor)
    for param in optimizer_constructor.params
        set(optimizer, param.first, param.second)
    end
    return optimizer
end

"""
    instantiate(
        optimizer_constructor,
        with_bridge_type::Union{Nothing, Type} = nothing,
    )

Creates an instance of optimizer by either:
 * calling `optimizer_constructor.optimizer_constructor()` and setting the
   parameters in `optimizer_constructor.params` if `optimizer_constructor` is a
   [`OptimizerWithAttributes`](@ref)
*  calling `optimizer_constructor()` if `optimizer_constructor` is callable.

If `with_bridge_type` is not `nothing`, it enables all the bridges defined in
the MathOptInterface.Bridges submodule with coefficient type `with_bridge_type`.

If the optimizer created by `optimizer_constructor` does not support loading the
problem incrementally (see [`supports_incremental_interface`](@ref)), then a
[`Utilities.CachingOptimizer`](@ref) is added to store a cache of the bridged
model.
"""
function instantiate(
    optimizer_constructor;
    with_bridge_type::Union{Nothing,Type} = nothing,
)
    optimizer = _instantiate_and_check(optimizer_constructor)
    if with_bridge_type === nothing
        return optimizer
    end
    if !supports_incremental_interface(optimizer)
        cache = default_cache(optimizer, with_bridge_type)
        optimizer = Utilities.CachingOptimizer(cache, optimizer)
    end
    return Bridges.full_bridge_optimizer(optimizer, with_bridge_type)
end

"""
    default_cache(optimizer::ModelLike, ::Type{T}) where {T}

Return a new instance of the default model type to be used as cache for
`optimizer` in a [`Utilities.CachingOptimizer`](@ref) for holding constraints
of coefficient type `T`. By default, this returns
`Utilities.UniversalFallback(Utilities.Model{T}())`. If copying from a instance
of a given model type is faster for `optimizer` then a new method returning
an instance of this model type should be defined.
"""
function default_cache(::ModelLike, ::Type{T}) where {T}
    return Utilities.UniversalFallback(Utilities.Model{T}())
end


text:
using DelaunayTriangulations.GeometryPredicates:
    two_sum_tail,
    two_sum,
    fast_two_sum_tail,
    fast_two_sum,
    two_diff_tail,
    two_diff,
    two_prod_tail,
    two_prod,
    grow_expansion,
    expansion_sum,
    expansion_diff,
    expansion_scale,
    approximate

@testset "Adaptive Arithmetic" begin
    @testset "Two Sum" begin
        @testset "Float64" begin
            @test two_sum_tail(0.1, 0.2, 0.1 + 0.2) < 0.0
            @test two_sum_tail(0.2, 0.2, 0.2 + 0.2) == 0.0
            @test two_sum_tail(0.4, 0.3, 0.4 + 0.3) > 0.0
            @test two_sum_tail(-0.1, -0.2, -0.1 - 0.2) > 0.0
            @test two_sum(0.1, 0.2) == (two_sum_tail(0.1, 0.2, 0.1 + 0.2), 0.1 + 0.2)
            @test two_sum(0.2, 0.2) == (two_sum_tail(0.2, 0.2, 0.2 + 0.2), 0.2 + 0.2)
            @test two_sum(0.4, 0.3) == (two_sum_tail(0.4, 0.3, 0.4 + 0.3), 0.4 + 0.3)
            @test two_sum(-0.1, -0.2) == (two_sum_tail(-0.1, -0.2, -0.1 - 0.2), -0.1 - 0.2)
        end
        @testset "Float32" begin
            @test two_sum_tail(Float32(0.1), Float32(0.2), Float32(0.1) + Float32(0.2)) < Float32(0.0)
            @test two_sum_tail(Float32(0.2), Float32(0.2), Float32(0.2) + Float32(0.2)) == Float32(0.0)
            @test two_sum_tail(Float32(0.3), Float32(0.2), Float32(0.3) + Float32(0.2)) > Float32(0.0)
            @test two_sum_tail(-Float32(0.1), -Float32(0.2), -Float32(0.1) - Float32(0.2)) > Float32(0.0)
            @test two_sum(Float32(0.1), Float32(0.2)) == (two_sum_tail(Float32(0.1), Float32(0.2), Float32(0.1) + Float32(0.2)), Float32(0.1) + Float32(0.2))
            @test two_sum(Float32(0.2), Float32(0.2)) == (two_sum_tail(Float32(0.2), Float32(0.2), Float32(0.2) + Float32(0.2)), Float32(0.2) + Float32(0.2))
            @test two_sum(Float32(0.3), Float32(0.2)) == (two_sum_tail(Float32(0.3), Float32(0.2), Float32(0.3) + Float32(0.2)), Float32(0.3) + Float32(0.2))
            @test two_sum(-Float32(0.1), -Float32(0.2)) == (two_sum_tail(-Float32(0.1), -Float32(0.2), -Float32(0.1) - Float32(0.2)), -Float32(0.1) - Float32(0.2))
        end
        @testset "Float16" begin
            @test two_sum_tail(Float16(0.1), Float16(0.2), Float16(0.1) + Float16(0.2)) > Float16(0.0)
            @test two_sum_tail(Float16(0.2), Float16(0.2), Float16(0.2) + Float16(0.2)) == Float16(0.0)
            @test two_sum_tail(Float16(0.5), Float16(0.2), Float16(0.5) + Float16(0.2)) < Float16(0.0)
            @test two_sum_tail(-Float16(0.1), -Float16(0.2), -Float16(0.1) - Float16(0.2)) < Float16(0.0)
            @test two_sum(Float16(0.1), Float16(0.2)) == (two_sum_tail(Float16(0.1), Float16(0.2), Float16(0.1) + Float16(0.2)), Float16(0.1) + Float16(0.2))
            @test two_sum(Float16(0.2), Float16(0.2)) == (two_sum_tail(Float16(0.2), Float16(0.2), Float16(0.2) + Float16(0.2)), Float16(0.2) + Float16(0.2))
            @test two_sum(Float16(0.5), Float16(0.2)) == (two_sum_tail(Float16(0.5), Float16(0.2), Float16(0.5) + Float16(0.2)), Float16(0.5) + Float16(0.2))
            @test two_sum(-Float16(0.1), -Float16(0.2)) == (two_sum_tail(-Float16(0.1), -Float16(0.2), -Float16(0.1) - Float16(0.2)), -Float16(0.1) - Float16(0.2))
        end
    end
    @testset "Fast Two Sum" begin
        @testset "Float64" begin
            @test fast_two_sum_tail(0.2, 0.1, 0.2 + 0.1) < 0.0
            @test fast_two_sum_tail(0.4, 0.3, 0.4 + 0.3) > 0.0
            @test fast_two_sum_tail(-0.1, -0.2, -0.1 - 0.2) > 0.0
            @test fast_two_sum(0.2, 0.1) == (fast_two_sum_tail(0.2, 0.1, 0.2 + 0.1), 0.2 + 0.1)
            @test fast_two_sum(0.4, 0.3) == (fast_two_sum_tail(0.4, 0.3, 0.4 + 0.3), 0.4 + 0.3)
            @test fast_two_sum(-0.1, -0.2) == (fast_two_sum_tail(-0.1, -0.2, -0.1 - 0.2), -0.1 - 0.2)
        end
        @testset "Float32" begin
            @test fast_two_sum_tail(Float32(0.2), Float32(0.1), Float32(0.2) + Float32(0.1)) < Float32(0.0)
            @test fast_two_sum_tail(Float32(0.3), Float32(0.2), Float32(0.3) + Float32(0.2)) > Float32(0.0)
            @test fast_two_sum_tail(-Float32(0.1), -Float32(0.2), -Float32(0.1) - Float32(0.2)) > Float32(0.0)
            @test fast_two_sum(Float32(0.2), Float32(0.1)) == (fast_two_sum_tail(Float32(0.2), Float32(0.1), Float32(0.2) + Float32(0.1)), Float32(0.2) + Float32(0.1))
            @test fast_two_sum(Float32(0.3), Float32(0.2)) == (fast_two_sum_tail(Float32(0.3), Float32(0.2), Float32(0.3) + Float32(0.2)), Float32(0.3) + Float32(0.2))
            @test fast_two_sum(-Float32(0.1), -Float32(0.2)) == (fast_two_sum_tail(-Float32(0.1), -Float32(0.2), -Float32(0.1) - Float32(0.2)), -Float32(0.1) - Float32(0.2))
        end
        @testset "Float16" begin
            @test fast_two_sum_tail(Float16(0.2), Float16(0.1), Float16(0.2) + Float16(0.1)) > Float16(0.0)
            @test fast_two_sum_tail(Float16(0.5), Float16(0.2), Float16(0.5) + Float16(0.2)) < Float16(0.0)
            @test fast_two_sum_tail(-Float16(0.1), -Float16(0.2), -Float16(0.1) - Float16(0.2)) < Float16(0.0)
            @test fast_two_sum(Float16(0.2), Float16(0.1)) == (fast_two_sum_tail(Float16(0.2), Float16(0.1), Float16(0.2) + Float16(0.1)), Float16(0.2) + Float16(0.1))
            @test fast_two_sum(Float16(0.5), Float16(0.2)) == (fast_two_sum_tail(Float16(0.5), Float16(0.2), Float16(0.5) + Float16(0.2)), Float16(0.5) + Float16(0.2))
            @test fast_two_sum(-Float16(0.1), -Float16(0.2)) == (fast_two_sum_tail(-Float16(0.1), -Float16(0.2), -Float16(0.1) - Float16(0.2)), -Float16(0.1) - Float16(0.2))
        end
    end
    @testset "Two Diff" begin
        @testset "Float64" begin
            @test two_diff_tail(0.3, 0.11, 0.3 - 0.11) < 0.0
            @test two_diff_tail(0.3, 0.1, 0.3 - 0.1) == 0.0
            @test two_diff_tail(0.411, 0.92, 0.411 - 0.92) > 0.0
            @test two_diff(0.3, 0.11) == (two_diff_tail(0.3, 0.11, 0.3 - 0.11), 0.3 - 0.11)
            @test two_diff(0.3, 0.1) == (two_diff_tail(0.3, 0.1, 0.3 - 0.1), 0.3 - 0.1)
            @test two_diff(0.411, 0.92) == (two_diff_tail(0.411, 0.92, 0.411 - 0.92), 0.411 - 0.92)
        end
        @testset "Float32" begin
            @test two_diff_tail(Float32(0.3), Float32(0.1), Float32(0.3) - Float32(0.1)) < Float32(0.0)
            @test two_diff_tail(Float32(0.3), Float32(0.11), Float32(0.3) - Float32(0.11)) == Float32(0.0)
            @test two_diff(Float32(0.3), Float32(0.1)) == (two_diff_tail(Float32(0.3), Float32(0.1), Float32(0.3) - Float32(0.1)), Float32(0.3) - Float32(0.1))
            @test two_diff(Float32(0.3), Float32(0.11)) == (two_diff_tail(Float32(0.3), Float32(0.11), Float32(0.3) - Float32(0.11)), Float32(0.3) - Float32(0.11))
        end
        @testset "Float16" begin
            @test two_diff_tail(Float16(0.11), Float16(0.412), Float16(0.11) - Float16(0.412)) > Float16(0.0)
            @test two_diff_tail(Float16(0.3), Float16(0.11), Float16(0.3) - Float16(0.11)) == Float16(0.0)
            @test two_diff(Float16(0.3), Float16(0.1)) == (two_diff_tail(Float16(0.3), Float16(0.1), Float16(0.3) - Float16(0.1)), Float16(0.3) - Float16(0.1))
            @test two_diff(Float16(0.3), Float16(0.11)) == (two_diff_tail(Float16(0.3), Float16(0.11), Float16(0.3) - Float16(0.11)), Float16(0.3) - Float16(0.11))
        end
    end
    @testset "Two Prod" begin
        @testset "Float64" begin
            @test two_prod_tail(0.3, 0.11, 0.3 * 0.11) < 0.0
            @test two_prod_tail(0.3, 0.1, 0.3 * 0.1) > 0.0
            @test two_prod_tail(2.0, 2.0, 2.0 * 2.0) == 0.0
            @test two_prod(0.3, 0.11) == (two_prod_tail(0.3, 0.11, 0.3 * 0.11), 0.3 * 0.11)
            @test two_prod(0.3, 0.1) == (two_prod_tail(0.3, 0.1, 0.3 * 0.1), 0.3 * 0.1)
            @test two_prod(0.2, 0.2) == (two_prod_tail(0.2, 0.2, 0.2 * 0.2), 0.2 * 0.2)
        end
        @testset "Float32" begin
            @test two_prod_tail(Float32(0.3), Float32(0.11), Float32(0.3) * Float32(0.11)) > Float32(0.0)
            @test two_prod_tail(Float32(0.3), Float32(0.09), Float32(0.3) * Float32(0.09)) < Float32(0.0)
            @test two_prod_tail(Float32(2.0), Float32(2.0), Float32(2.0) * Float32(2.0)) == Float32(0.0)
            @test two_prod(Float32(0.3), Float32(0.11)) == (two_prod_tail(Float32(0.3), Float32(0.11), Float32(0.3) * Float32(0.11)), Float32(0.3) * Float32(0.11))
            @test two_prod(Float32(0.3), Float32(0.09)) == (two_prod_tail(Float32(0.3), Float32(0.09), Float32(0.3) * Float32(0.09)), Float32(0.3) * Float32(0.09))
            @test two_prod(Float32(0.2), Float32(0.2)) == (two_prod_tail(Float32(0.2), Float32(0.2), Float32(0.2) * Float32(0.2)), Float32(0.2) * Float32(0.2))
        end
        @testset "Float16" begin
            @test two_prod_tail(Float16(0.3), Float16(0.11), Float16(0.3) * Float16(0.11)) > Float16(0.0)
            @test two_prod_tail(Float16(0.3), Float16(0.05), Float16(0.3) * Float16(0.05)) < Float16(0.0)
            @test two_prod_tail(Float16(2.0), Float16(2.0), Float16(2.0) * Float16(2.0)) == Float16(0.0)
            @test two_prod(Float16(0.3), Float16(0.11)) == (two_prod_tail(Float16(0.3), Float16(0.11), Float16(0.3) * Float16(0.11)), Float16(0.3) * Float16(0.11))
            @test two_prod(Float16(0.3), Float16(0.05)) == (two_prod_tail(Float16(0.3), Float16(0.05), Float16(0.3) * Float16(0.05)), Float16(0.3) * Float16(0.05))
            @test two_prod(Float16(0.2), Float16(0.2)) == (two_prod_tail(Float16(0.2), Float16(0.2), Float16(0.2) * Float16(0.2)), Float16(0.2) * Float16(0.2))
        end
    end
    @testset "Expansions" begin
        @testset "Grow Expansion" begin
            @test grow_expansion(0.123, 0.456) == two_sum(0.123, 0.456)
            @test grow_expansion(0.42) == 0.42
            e = rand()
            for i in 1:10
                e = grow_expansion(e, rand())
                @test length(e) == i + 1
                for j in 1:i
                    @test e[j] < e[j+1] || e[j] == 0.0 || e[j+1] == 0.0
                end
            end
        end
        @testset "Expansion Sum" begin
            for i in 1:3
                e = grow_expansion(([rand() for _ in 1:i]...,), rand())
                for j in 1:3
                    f = grow_expansion(([rand() for _ in 1:j]...,), rand())
                    g = expansion_sum(e, f)
                    @test length(g) == length(e) + length(f)
                    @test g[end] ≈ e[end] + f[end]
                end
            end
        end
        @testset "Expansion Diff" begin
            e, f = grow_expansion(rand(), rand()), rand()
            g = expansion_diff(e, f)
            @test length(g) == length(e) + length(f)
            @test g[end] ≈ e[end] - f[end]
            e, f = grow_expansion(rand(), rand()), grow_expansion(rand(), rand())
            g = expansion_diff(e, f)
            @test length(g) == length(e) + length(f)
            @test g[end] ≈ e[end] - f[end]
        end
        @testset "Expansion Scale" begin
            for i in 1:10
                e, f = grow_expansion(([rand() for _ in 1:i]...,), rand()), rand()
                g = expansion_scale(e, f)
                @test length(g) == 2 * length(e)
                @test sum(g) ≈ f * sum(e)
            end
        end
        @testset "Approximate" begin
            e = grow_expansion(([rand() for _ in 1:10]...,), rand())
            @test approximate(e) == sum(e)
        end
    end
end

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
################################################################################
#
# Routines for the the evaluation and gradient computation of TWRI
#
################################################################################

export objTWRIdual, objTWRIdual!

### Optim.jl-adapted wrapper

function objTWRIdual!(
	F, Gm, m::Array{R, 2},
	n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R},
	wav::judiVector, dat::judiVector,
	eps::Array{R, 1};
	opt = Options(),
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	if Gm == nothing
		mode = "eval"
	else
		mode = "grad"
	end
	argout = objTWRIdual(m, n, d, o, nothing, wav, dat, eps; opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
	if Gm != nothing
		fval = argout[1]
		Gm .= argout[2]
	else
		fval = argout
	end
	return fval

end

function objTWRIdual!(
	F, Gm, m::Array{R, 2},
	y::judiVector,
	n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R},
	wav::judiVector, dat::judiVector,
	eps::Array{R, 1};
	opt = Options(),
	objfact = 1f0,
	comp_alpha = true, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	if Gm == nothing
		mode = "eval"
	else
		mode = "grad"
	end
	argout = objTWRIdual(m, n, d, o, y, wav, dat, eps; opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = false, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
	if Gm != nothing
		fval = argout[1]
		Gm .= argout[2]
	else
		fval = argout
	end
	return fval

end

### Objective wrapper for array input

function objTWRIdual(
	m::Array{R, 2}, n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R},
	wav::judiVector, dat::judiVector,
	eps::Array{R, 1};
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	return objTWRIdual(m, n, d, o, nothing, wav, dat, eps; opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

function objTWRIdual(
	m::Array{R, 2}, n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R}, y::Union{Nothing, judiVector},
	wav::judiVector, dat::judiVector,
	eps::Array{R, 1};
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	return objTWRIdual(Model(n, d, o, m), y, wav, dat, eps; opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

### Objective wrapper for JUDI data types

function objTWRIdual(
	m::Model, y::Union{Nothing, judiVector},
	wav::judiVector, dat::judiVector,
	eps::Array{R, 1};
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	if y == nothing
		return objTWRIdual_rawinput(m, nothing, wav.geometry, dat.geometry, wav.data, dat.data, eps, 1:length(dat.data); opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
	else
		return objTWRIdual_rawinput(m, y.data, wav.geometry, dat.geometry, wav.data, dat.data, eps, 1:length(dat.data); opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
	end

end

function objTWRIdual(
	m::Model,
	wav::judiVector, dat::judiVector,
	eps::Array{R, 1};
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	return objTWRIdual(m, nothing, wav, dat, eps; opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

### Objective: serial/parallel implementation

## Parallel implementation

function objTWRIdual_rawinput(
	model::Model, y_data,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data, dat_data,
	eps::Array{R, 1},
	src_idx::UnitRange{Int64};
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	# Setup parallelization
	p = default_worker_pool()
	time_modeling_par = remote(objTWRIdual_rawinput)
	time_modeling = retry(time_modeling_par)

	# Initialize output
	nsrc = length(src_idx)
	results = Array{Any}(undef, nsrc)

	# Process shots from source channel asynchronously
	@sync begin
		for j = 1:nsrc

			# Local geometry for current position
			opt_loc = subsample(opt, j)
			src_geom_loc = subsample(src_geom, j)
			rcv_geom_loc = subsample(rcv_geom, j)

			# Selecting variables for current shot index
			src_data_loc = src_data[j]
			if y_data == nothing
				y_data_loc = nothing
			else
				y_data_loc = y_data[j]
			end
			dat_data_loc = dat_data[j]
			if isa(Filter, Array{Array{C, 1}, 1})
				Filter_loc = Filter[j]
			else
				Filter_loc = Filter
			end
			eps_loc = eps[j]
			if dt_comp == nothing
				dt_comp_loc = nothing
			else
				dt_comp_loc = dt_comp[j]
			end

			# Local result
			results[j] = @spawn objTWRIdual_rawinput(model, y_data_loc, src_geom_loc, rcv_geom_loc, src_data_loc, dat_data_loc, eps_loc, opt = opt_loc; mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp_loc, Filter = Filter_loc, gradmprec_fun = gradmprec_fun)

		end
	end

	# Aggregating results
	if mode == "eval"
		argout = fetch(results[1])
		for j = 2:nsrc
			argout += fetch(results[j])
		end
	elseif mode == "grad" && y_data == nothing
		arg = fetch(results[1])
		argout1 = arg[1]
		argout2 = arg[2]
		for j = 2:nsrc
			arg = fetch(results[j])
			argout1 += arg[1]
			argout2 += arg[2]
		end
		argout = (argout1, argout2)
	else
		arg = fetch(results[1])
		argout1 = arg[1]
		argout2 = arg[2]
		argout3 = arg[3]
		for j = 2:nsrc
			arg = fetch(results[j])
			argout1 += arg[1]
			argout2 += arg[2]
			argout3 = [argout3; arg[3]]
		end
		argout = (argout1, argout2, argout3)
	end
	return argout./nsrc

end

## Serial implementation

function objTWRIdual_rawinput(
	model_full::Model, y_data::Union{Nothing, Array{R, 2}},
	src_geom::Geometry, rcv_geom::Geometry,
	src_data::Array{R, 2}, dat_data::Array{R, 2},
	eps::R;
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, R} = nothing,
	Filter::Union{Nothing, Array{C, 1}} = nothing,
	gradmprec_fun = g->g)

	# Setting pre-defined absorbing layer size
	model_full.nb = NB

	# Load full geometry for out-of-core geometry containers
	typeof(rcv_geom) == GeometryOOC && (rcv_geom = Geometry(rcv_geom))
	typeof(src_geom) == GeometryOOC && (src_geom = Geometry(src_geom))
	length(model_full.n) == 3 ? dims = [3, 2, 1] : dims = [2, 1] # model dimensions for Python are (z,y,x) and (z,x)

	# Limit model to area with sources/receivers
	if opt.limit_m == true
		model = deepcopy(model_full)
		model = limit_model_to_receiver_area(src_geom, rcv_geom, model, opt.buffer_size)
	else
		model = model_full
	end

	# Set up Python model structure
	model_py = devito_model(model, "F", 1, opt, 0)

	# Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
	rcv_geom = remove_out_of_bounds_receivers(rcv_geom, model)

	# Call to objective with Julia/Devito interface function
	return objTWRIdual_jldevito(model_py, y_data, model.o, src_geom, rcv_geom, src_data, dat_data, eps; opt = opt, mode = mode, objfact = objfact, comp_alpha = comp_alpha, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

### Julia/devito implementation

function objTWRIdual_jldevito(
	model_py::PyCall.PyObject, y_data::Union{Nothing, Array{R, 2}},
	origin,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data::Array{R, 2}, dat_data::Array{R, 2},
	eps::R;
	opt = Options(),
	mode = "eval",
	objfact = 1f0,
	comp_alpha = true, grad_corr = false, weight_fun_pars = nothing, dt_comp::Union{Nothing, R} = nothing,
	Filter::Union{Nothing, Array{C, 1}} = nothing,
	gradmprec_fun = g->g)

	# Loading python modules for devito implementation of objTWRIdual
	devitomod = load_codegen_TWRIdual()

	# Interpolate input data to computational grid
	if dt_comp == nothing
		dt_comp = R(model_py.critical_dt)
	end
	q_in = time_resample(src_data, src_geom, dt_comp)[1]
	if y_data != nothing
		y_in = time_resample(y_data, rcv_geom, dt_comp)[1]
	else
		y_in = nothing
	end
	dat_in = time_resample(dat_data, rcv_geom, dt_comp)[1]
	nt_comp = size(q_in, 1)
	nt_rcv = Int(trunc(rcv_geom.t[1]/dt_comp+1))

	# Set up coordinates with devito dimensions
	src_coords = setup_grid(src_geom, model_py.shape, origin)
	rcv_coords = setup_grid(rcv_geom, model_py.shape, origin)

	# Computing output
	if y_in == nothing
		y_in_ = nothing
	else
		y_in_ = PyReverseDims(copy(transpose(y_in)))
	end
	if mode == "eval"
		argout = pycall(devitomod.objTWRIdual_devito,
						R,
						model_py, y_in_,
						PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(rcv_coords))),
						PyReverseDims(copy(transpose(q_in))), PyReverseDims(copy(transpose(dat_in))), Filter,
						eps,
						mode,
						objfact,
						comp_alpha, grad_corr, weight_fun_pars, dt_comp, SPACE_ORDER)
	elseif mode == "grad" && y_in != nothing
		argout = pycall(devitomod.objTWRIdual_devito,
						Tuple{R, Array{R, 2}, Array{R, 2}},
						model_py, y_in_,
						PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(rcv_coords))),
						PyReverseDims(copy(transpose(q_in))), PyReverseDims(copy(transpose(dat_in))), Filter,
						eps,
						mode,
						objfact,
						comp_alpha, grad_corr, weight_fun_pars, dt_comp, SPACE_ORDER)
	else
		argout = pycall(devitomod.objTWRIdual_devito,
						Tuple{R, Array{R, 2}},
						model_py, y_in_,
						PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(rcv_coords))),
						PyReverseDims(copy(transpose(q_in))), PyReverseDims(copy(transpose(dat_in))), Filter,
						eps,
						mode,
						objfact,
						comp_alpha, grad_corr, weight_fun_pars, dt_comp, SPACE_ORDER)
	end

	# Output post-processing
	if mode == "grad" && y_in != nothing
		gm = remove_padding(argout[2], model_py.nbpml, true_adjoint = opt.sum_padding) # remove AB layers
		gy = time_resample(argout[3], dt_comp, rcv_geom) # resample
		argout = (argout[1], gradmprec_fun(gm), [gy])
	elseif mode == "grad" && y_in == nothing
		gm = remove_padding(argout[2], model_py.nbpml, true_adjoint = opt.sum_padding) # remove AB layers
		argout = (argout[1], gradmprec_fun(gm))
	end

	return argout

end

# Loader for python devito modeling
function load_codegen_TWRIdual()

	pushfirst!(PyVector(pyimport("sys")."path"), string(pwd(), "/src/"))
	return pyimport("TWRIdualFun_devito")

end


text:
function write_html_table(io::IO, tab; header = Tables.columnnames(first(Tables.rows(tab))), style = (value, i, j) -> "", format = (value, i, j) -> value, include_columns = 1:length(header), standalone = false, css = """
    table, td, th {
        border-collapse: collapse;
        font-family: sans-serif;
    }
    td, th {
        border-bottom: 0;
        background: #fff !important;
        padding: 6px
    }
    tr.header {
        background: #fff !important;
        font-weight: bold;
    }
    """)
    rows = Tables.rows(tab)
    if standalone
        println(io,"""
        <!DOCTYPE html>

        <meta charset="UTF-8">

        """)
    end
    println(io, "
")
    println(io, """
""")
    for (j, h) in enumerate(header)
        j ∈ include_columns || continue
        println(io, """
        $(format(h, 0, j))	""")
    end
    println(io, "
")
    for (i, row) in enumerate(rows)
        print(io, """
""")
        for (j, h) in enumerate(row)
            j ∈ include_columns || continue
            println(io, """
                $(format(h, i, j))	""")
        end
        print(io, "
")
    end
    println(io, "
")
    if standalone
        println(io, """

        """)
    end
end

function html_table(io::IO, ts::TableTestSet; kwargs...)
    tab = Tables.rowtable(ts)
    style = function(value, i, j)
        str = ""

        if j == 1
            str *= "text-align:left;"
            str *= "border-right: solid 2px;"
        else
            str *= "text-align:center;"
        end

        # if it's the header, stop here
        if i == 0
            return str
        end

        if j == 2 && tab[i][j] > 0
            str *= "color:green;"
        elseif j ∈ (3,4) && tab[i][j] > 0
            str *= "color:red;"
        elseif j == 6 && tab[i][j] > 0
            str *= "color:blue;"
        end

        return str
    end

    format = function(value, i, j)
        if i > 0 && j == 1
            depth = tab[i][7]
            value = repeat(" ", 2*depth) * value
        end
        return value
    end

    write_html_table(io::IO, ts; style = style, format = format, include_columns = 1:6, kwargs...)
end