text:
#=
Short characteristic method to solve intensity along rays going through all grid
points.
=#

include("functions.jl")

"""
    function short_characteristics_up(k::Vector{Float64},
                                      S_0::Array{<:UnitsIntensity_λ, 3},
                                      I_0::Matrix{<:UnitsIntensity_λ},
                                      α::Array{<:PerLength, 3},
                                      atmos::Atmosphere;
                                      pt::Bool=false,
                                      n_sweeps::Int=3)

Computes intensity along rays traveling upwards from the bottom to the top
through all grid points in the atmosphere. Initial intensity I_0 = B_λ(T) at the
bottom of the domain.
"""
function short_characteristics_up(k::Vector{Float64},
                                  S_0::Array{<:UnitsIntensity_λ, 3},
                                  I_0::Matrix{<:UnitsIntensity_λ},
                                  α::Array{<:PerLength, 3},
                                  atmos::Atmosphere;
                                  pt::Bool=false,
                                  n_sweeps::Int=3)
    ############################################################################
    # | and - : Grid
    #     x   : Grid points
    #     *   : Ray
    #
    #                  x----------x----------x  z_i+1
    #                  |          |          |
    #                  |          |          |
    #                  |          |          |
    #                  x----------*----------x  z_i
    #                  |          |*         |
    #                  |          | *        |
    #                  |          |  *       |
    #                  x----------x---*------x  z_i-1
    #                                 :
    #                             upwind point
    ############################################################################

    # Allocate array for new intensity
    I = zero(S_0)

    nz = length(atmos.z)
    Δx = atmos.x[2] - atmos.x[1]
    Δy = atmos.y[2] - atmos.y[1]

    # I know that Δx = Δy = constant for all grid points
    # Δxy = atmos.x[2] - atmos.x[1]
    r_x = abs(Δx/k[2])
    r_y = abs(Δy/k[3])

    # find out direction in xy the ray moves (1 for positive, -1 for negative)
    sign_x, sign_y = xy_intersect(k)

    # Boundary condition
    I[1,:,:] = I_0

    # loop upwards through atmosphere
    for idz in 2:nz
        # Grid diffference in z
        Δz = atmos.z[idz] - atmos.z[idz-1]

        # calculate length until ray hits z plane
        r_z = abs(Δz/k[1])

        # This finds which plane the ray intersects with
        plane_cut = argmin([r_z, r_x, r_y])
        if plane_cut == 1
            I[idz,:,:] = xy_up_ray(k, idz, sign_x, sign_y, I[idz-1,:,:], S_0, α, atmos)
            if idz==2 && pt == true
                println("xy_up_ray")
            end
            k::Vector{Float64}, idz::Int, sign_x::Int,
                               sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ},
                               S_0::Array{<:UnitsIntensity_λ, 3},
                               α::Array{<:PerLength, 3}, atmos::Atmosphere
        elseif plane_cut==2
            I[idz,:,:] = yz_up_ray(k, idz, sign_x, sign_y, I[idz-1,:,:], S_0, α, atmos, n_sweeps)
            if idz==2 && pt == true
                println("yz_up_ray")
            end
        elseif plane_cut==3
            I[idz,:,:] = xz_up_ray(k, idz, sign_x, sign_y, I[idz-1,:,:], S_0, α, atmos, n_sweeps)
            if idz==2 && pt == true
                println("xz_up_ray")
            end
        end
    end

    return I
end

"""
    function short_characteristics_down(k::Vector{Float64},
                                        S_0::Array{<:UnitsIntensity_λ, 3},
                                        I_0::Matrix{<:UnitsIntensity_λ},
                                        α::Array{<:PerLength, 3},
                                        atmos::Atmosphere;
                                        pt::Bool=false,
                                        n_sweeps::Int=3)

Computes intensity along rays traveling downwards from the top to the bottom
through all grid points in the atmosphere. Initial intensity I_0 = 0 at the top
of the domain.
"""
function short_characteristics_down(k::Vector{Float64},
                                    S_0::Array{<:UnitsIntensity_λ, 3},
                                    I_0::Matrix{<:UnitsIntensity_λ},
                                    α::Array{<:PerLength, 3},
                                    atmos::Atmosphere;
                                    pt::Bool=false,
                                    n_sweeps::Int=3)
    ############################################################################
    # | and - : Grid
    #     x   : Grid points
    #     *   : Ray
    #
    #                  x----------x----------x  z_i+1
    #                  |          |          *  upwind point
    #                  |          |      *   |
    #                  |          |  *       |
    #                  x----------*----------x  z_i
    #                  |          |          |
    #                  |          |          |
    #                  |          |          |
    #                  x----------x----------x  z_i-1
    #
    ############################################################################

    # Allocate array for new intensity
    I = zero(S_0)

    nz = length(atmos.z)
    Δx = atmos.x[2] - atmos.x[1]
    Δy = atmos.y[2] - atmos.y[1]

    # I know that Δx = Δy = constant for all grid points
    r_x = abs(Δx/k[2])
    r_y = abs(Δy/k[3])

    # find out which plane the upwind part of the ray intersects
    sign_x, sign_y = xy_intersect(k)

    # Boundary condition
    I[end,:,:] = I_0

    # loop downwards through atmosphere
    for idz in nz-1:-1:1
        # Grid diffference in z
        Δz = atmos.z[idz+1] - atmos.z[idz]

        # calculate length until ray hits z plane
        r_z = abs(Δz/k[1])

        # This finds which plane the ray intersects with
        plane_cut = argmin([r_z, r_x, r_y])
        if plane_cut == 1
            I[idz,:,:] = xy_down_ray(k, idz, sign_x, sign_y, I[idz+1,:,:], S_0, α, atmos)
            if idz==nz-1 && pt == true
                println("xy_down_ray")
            end
        elseif plane_cut==2
            I[idz,:,:] = yz_down_ray(k, idz, sign_x, sign_y, I[idz+1,:,:], S_0, α, atmos, n_sweeps)
            if idz==nz-1 && pt == true
                println("yz_down_ray")
            end
        elseif plane_cut==3
            I[idz,:,:] = xz_down_ray(k, idz, sign_x, sign_y, I[idz+1,:,:], S_0, α, atmos, n_sweeps)
            if idz==nz-1 && pt == true
                println("xz_down_ray")
            end
        end
    end

    return I
end

"""
    function xy_up_ray(k::Vector{Float64}, idz::Int, sign_x::Int,
                       sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ},
                       S_0::Array{<:UnitsIntensity_λ, 3},
                       α::Array{<:PerLength, 3}, atmos::Atmosphere)

Ray moving upwards, upwind point intersecting with lower xy plane. Assumes angle
in radians.
"""
function xy_up_ray(k::Vector{Float64}, idz::Int, sign_x::Int,
                   sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ},
                   S_0::Array{<:UnitsIntensity_λ, 3},
                   α::Array{<:PerLength, 3}, atmos::Atmosphere)

    # Allocate space for intensity
    I = zero(I_0)

    nx = length(atmos.x)
    ny = length(atmos.y)

    # Z center and upwind position
    z_centre = atmos.z[idz]
    idz_upwind = idz - 1
    z_upwind = atmos.z[idz_upwind]

    # Length to plane
    Δz = z_upwind-z_centre
    r = abs(Δz/k[1])

    x_increment = r*k[2]
    y_increment = r*k[3]

    α_lower = α[idz_upwind,:,:]

    S_lower = S_0[idz_upwind,:,:]

    for idx in 2:nx-1
        idx_lower = idx - Int((sign_x+1)/2)
        idx_upper = idx_lower + 1

        x_centre = atmos.x[idx]
        x_upwind = x_centre + x_increment

        x_bounds = (atmos.x[idx_lower], atmos.x[idx_upper])

        for idy in 2:ny-1
            # Centre point (c)
            y_centre = atmos.y[idy]

            # upwind coordinate
            y_upwind = y_centre + y_increment

            # Lower corner to interpolate from
            idy_lower = idy - Int((sign_y+1)/2)
            idy_upper = idy_lower + 1

            y_bounds = (atmos.y[idy_lower], atmos.y[idy_upper])

            # Do linear interpolation to find τ and S

            # exinction at centre and upwind point
            α_vals = [α_lower[idx_lower, idy_lower] α_lower[idx_lower, idy_upper]
                      α_lower[idx_upper, idy_lower] α_lower[idx_upper, idy_upper]]

            α_centre = α[idz, idx, idy]
            α_upwind = bilinear(x_upwind, y_upwind, x_bounds, y_bounds, α_vals)

            # Find the Δτ optical path from upwind to grid point
            Δτ_upwind = trapezoidal(r, α_centre, α_upwind)

            # Find source function at grid point and upwind point (intepolate)
            S_vals = [S_lower[idx_lower, idy_lower] S_lower[idx_lower, idy_upper]
                      S_lower[idx_upper, idy_lower] S_lower[idx_upper, idy_upper]]

            S_centre = S_0[idz, idx, idy]
            S_upwind = bilinear(x_upwind, y_upwind, x_bounds, y_bounds, S_vals)

            # Compute weights for formal solution
            α_ijk, β_ijk, expΔτ = linear_weights(Δτ_upwind)
            # println(α_ijk + β_ijk + expΔτ)
            # Interpolate to find intensity at upwind point
            I_vals = [I_0[idx_lower, idy_lower] I_0[idx_lower, idy_upper]
                      I_0[idx_upper, idy_lower] I_0[idx_upper, idy_upper]]

            I_upwind = bilinear(x_upwind, y_upwind, x_bounds, y_bounds, I_vals)

            # Integrate intensity from two-point quadrature
            I[idx, idy] = expΔτ*I_upwind + α_ijk*S_upwind + β_ijk*S_centre
        end

        # Update ghost zones
        I[idx, 1] = I[idx, end-1]
        I[idx, end] = I[idx, 2]
    end
    # Update ghost zones
    I[1,:] = I[end-1,:]
    I[end,:] = I[2,:]
    return I
end

"""
    function xy_down_ray(θ::AbstractFloat, ϕ::AbstractFloat, idz::Int, sign_x::Int,
                         sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                         α::Array{<:PerLength, 3}, atmos::Atmosphere)

Ray moving downwards, upwind point intersecting with upper xy plane. Assumes
angle in radians.
"""
function xy_down_ray(k::Vector{Float64}, idz::Integer, sign_x::Integer,
                     sign_y::Integer, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                     α::Array{<:PerLength, 3}, atmos::Atmosphere)

    # Allocate space for intensity
    I = zero(I_0)

    nx = length(atmos.x)
    ny = length(atmos.y)

    # Z center and upwind position
    z_centre = atmos.z[idz]
    idz_upwind = idz + 1
    z_upwind = atmos.z[idz_upwind]

    # Length to plane
    Δz = z_upwind - z_centre
    r = abs(Δz/k[1])

    x_increment = r*k[2]
    y_increment = r*k[3]

    α_upper = α[idz_upwind,:,:]

    S_upper = S_0[idz_upwind,:,:]

    for idx in 2:nx-1
        idx_lower = idx - Int((sign_x+1)/2)
        idx_upper = idx_lower+1

        x_centre = atmos.x[idx]
        x_upwind = x_centre + x_increment

        x_bounds = (atmos.x[idx_lower], atmos.x[idx_upper])

        for idy in 2:ny-1
            idy_lower = idy - Int((sign_y+1)/2)
            idy_upper = idy_lower+1

            y_centre = atmos.y[idy]
            y_upwind = y_centre + y_increment

            y_bounds = (atmos.y[idy_lower], atmos.y[idy_upper])
            # Do linear interpolation to find τ and S

            # exinction at centre and upwind point
            α_vals = [α_upper[idx_lower, idy_lower] α_upper[idx_lower, idy_upper]
                      α_upper[idx_upper, idy_lower] α_upper[idx_upper, idy_upper]]

            α_centre = α[idz, idx, idy]
            α_upwind = bilinear(x_upwind, y_upwind, x_bounds, y_bounds, α_vals)

            # Find the Δτ optical path from upwind to grid point
            Δτ_upwind = trapezoidal(r, α_centre, α_upwind)

            # Find source function at grid point and upwind point (intepolate)
            S_vals = [S_upper[idx_lower, idy_lower] S_upper[idx_lower, idy_upper]
                      S_upper[idx_upper, idy_lower] S_upper[idx_upper, idy_upper]]

            S_centre = S_0[idz, idx, idy]
            S_upwind = bilinear(x_upwind, y_upwind, x_bounds, y_bounds, S_vals)

            # Compute weights for formal solution
            α_ijk, β_ijk, expΔτ = linear_weights(Δτ_upwind)

            # Interpolate to find intensity at upwind point
            I_vals = [I_0[idx_lower, idy_lower]     I_0[idx_lower, idy_upper]
                      I_0[idx_upper, idy_lower]     I_0[idx_upper, idy_upper]]

            I_upwind = bilinear(x_upwind, y_upwind, x_bounds, y_bounds, I_vals)

            # Integrate intensity from two-point quadrature
            I[idx, idy] = expΔτ*I_upwind + α_ijk*S_upwind + β_ijk*S_centre

        end
        # Update ghost zones
        I[idx, 1] = I[idx, end-1]
        I[idx, end] = I[idx, 2]
    end
    # Update ghost zones
    I[1,:] = I[end-1,:]
    I[end,:] = I[2,:]
    return I
end

"""
    function yz_up_ray(θ::AbstractFloat, ϕ::AbstractFloat, idz::Int, sign_x::Int,
                       sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                       α::Array{<:PerLength, 3}, atmos::Atmosphere)

Ray moving upwards, upwind point intersecting with yz plane. Assumes angle
in radians.
"""
function yz_up_ray(k::Vector{Float64}, idz::Int, sign_x::Int,
                   sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                   α::Array{<:PerLength, 3}, atmos::Atmosphere, n_sweeps::Int)

    nx = length(atmos.x)
    ny = length(atmos.y)

    Δx = atmos.x[2] - atmos.x[1]

    # Loop direction
    start_x, stop_x = range_bounds(sign_x, nx)
    start_y, stop_y = range_bounds(sign_y, ny)

    # Allocate space for intensity
    I = zero(I_0)

    I_upper = zero(I_0[end,:])

    # Z center and interpolation position
    z_centre = atmos.z[idz]
    idz_lower = idz-1

    # calculate the length to upwind position
    r = abs(Δx/k[2])
    z_increment = r*k[1]
    y_increment = r*k[3]

    z_upwind = z_centre + z_increment

    z_bounds = (atmos.z[idz_lower], z_centre)

    α_lower = α[idz_lower,:,:]
    α_upper = α[idz,:,:]

    S_lower = S_0[idz_lower,:,:]
    S_upper = S_0[idz,:,:]

    for sweep in 1:n_sweeps
        for idx in start_x:sign_x:stop_x
            x_centre = atmos.x[idx]
            # x upwind position
            idx_upwind = idx+sign_x
            x_upwind = atmos.x[idx_upwind]
            for idy in start_y:sign_y:stop_y
                # Centre point (c)
                y_centre = atmos.y[idy]

                # Lower corner to interpolate from
                idy_lower = idy-Int((sign_y+1)/2)
                idy_upper = idy_lower+1

                y_upwind = y_centre + y_increment

                y_bounds = (atmos.y[idy_lower], atmos.y[idy_upper])

                # Do linear interpolation to find τ and S

                # exinction at centre and upwind point
                # z is 1st coordinate, changes rows in vals. y is 2nd coordinate
                α_vals = [α_lower[idx_upwind, idy_lower]   α_lower[idx_upwind, idy_upper]
                          α_upper[idx_upwind, idy_lower]   α_upper[idx_upwind, idy_upper]]

                α_centre = α_upper[idx, idy]
                α_upwind = bilinear(z_upwind, y_upwind, z_bounds, y_bounds, α_vals)

                # Find the Δτ optical path from upwind to grid point
                Δτ_upwind = trapezoidal(r, α_centre, α_upwind)

                # Find source function at grid point and upwind point (intepolate)
                S_vals = [S_lower[idx_upwind, idy_lower]    S_lower[idx_upwind, idy_upper]
                          S_upper[idx_upwind, idy_lower]    S_upper[idx_upwind, idy_upper]]

                S_centre = S_upper[idx, idy]
                S_upwind = bilinear(z_upwind, y_upwind, z_bounds, y_bounds, S_vals)

                # Compute weights for formal solution
                α_ijk, β_ijk, expΔτ = linear_weights(Δτ_upwind)

                # Interpolate to find intensity at upwind point
                # I_0 is from the top
                I_vals = [I_0[idx_upwind, idy_lower]    I_0[idx_upwind, idy_upper]
                          I_upper[idy_lower]            I_upper[idy_upper]        ]

                I_upwind = bilinear(z_upwind, y_upwind, z_bounds, y_bounds, I_vals)

                # Integrate intensity from two-point quadrature
                I[idx, idy] = expΔτ*I_upwind + α_ijk*S_upwind + β_ijk*S_centre
            end
            # Update ghost zones
            I[idx, 1] = I[idx, end-1]
            I[idx, end] = I[idx, 2]

            # Update upper interpolation value
            I_upper = I[idx, :]
        end

        # Update ghost zones
        I[1,:] = I[end-1,:]
        I[end,:] = I[2,:]
    end

    return I
end

"""
    function yz_down_ray(θ::AbstractFloat, ϕ::AbstractFloat, idz::Int, sign_x::Int,
                         sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                         α::Array{<:PerLength, 3}, atmos::Atmosphere)

Ray moving downwards, upwind point intersecting with yz plane. Assumes angle
in radians
"""
function yz_down_ray(k::Vector{Float64}, idz::Int, sign_x::Int,
                       sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                       α::Array{<:PerLength, 3}, atmos::Atmosphere, n_sweeps)

    nx = length(atmos.x)
    ny = length(atmos.y)

    Δx = atmos.x[2] - atmos.x[1]

    # Loop direction
    start_x, stop_x = range_bounds(sign_x, nx)
    start_y, stop_y = range_bounds(sign_y, ny)

    # Allocate space for intensity
    I = zero(I_0)

    I_lower = zero(I_0[end,:])

    # Z center and interpolation position
    z_centre = atmos.z[idz]
    idz_upper = idz+1

    # Length to plane
    # Δx = atmos.x[2] - atmos.x[1]

    # calculate the length to upwind position
    r = abs(Δx/k[2])
    z_increment = r*k[1]
    y_increment = r*k[3]

    z_upwind = z_centre + z_increment

    z_bounds = (z_centre, atmos.z[idz_upper])

    α_lower = α[idz,:,:]
    α_upper = α[idz_upper,:,:]

    S_lower = S_0[idz,:,:]
    S_upper = S_0[idz_upper,:,:]

    for sweep in 1:n_sweeps
        for idx in start_x:sign_x:stop_x
            x_centre = atmos.x[idx]
            # X upwind position
            idx_upwind = idx+sign_x
            x_upwind = atmos.x[idx_upwind]

            for idy in start_y:sign_y:stop_y
                # Centre point (c)
                y_centre = atmos.y[idy]

                # Lower corner to interpolate from
                idy_lower = idy-Int((sign_y+1)/2)
                idy_upper = idy_lower+1

                y_upwind = y_centre + y_increment

                y_bounds = (atmos.y[idy_lower], atmos.y[idy_upper])

                # Do linear interpolation to find τ and S

                # exinction at centre and upwind point
                # z is 1st coordinate, changes rows in vals. y is 2nd coordinate
                α_vals = [α_lower[idx_upwind, idy_lower]   α_lower[idx_upwind, idy_upper]
                          α_upper[idx_upwind, idy_lower]   α_upper[idx_upwind, idy_upper]]

                α_centre = α_lower[idx, idy]
                α_upwind = bilinear(z_upwind, y_upwind, z_bounds, y_bounds, α_vals)

                # Find the Δτ optical path from upwind to grid point
                Δτ_upwind = trapezoidal(r, α_centre, α_upwind)

                # Find source function at grid point and upwind point (intepolate)
                S_vals = [S_lower[idx_upwind, idy_lower]    S_lower[idx_upwind, idy_upper]
                          S_upper[idx_upwind, idy_lower]    S_upper[idx_upwind, idy_upper]]

                S_centre = S_lower[idx, idy]
                S_upwind = bilinear(z_upwind, y_upwind, z_bounds, y_bounds, S_vals)

                # Compute weights for formal solution
                α_ijk, β_ijk, expΔτ = linear_weights(Δτ_upwind)

                # Interpolate to find intensity at upwind point
                # I_0 is from the top
                I_vals = [I_lower[idy_lower]         I_lower[idy_upper]
                          I_0[idx_upwind, idy_lower] I_0[idx_upwind, idy_upper]]

                I_upwind = bilinear(z_upwind, y_upwind, z_bounds, y_bounds, I_vals)

                # Integrate intensity from two-point quadrature
                I[idx, idy] = expΔτ*I_upwind + α_ijk*S_upwind + β_ijk*S_centre
            end

            # Update ghost zones
            I[idx, 1] = I[idx, end-1]
            I[idx, end] = I[idx, 2]

            # Update upper interpolate value
            I_lower = I[idx, :]
        end
    end

    # Update ghost zones
    I[1,:] = I[end-1,:]
    I[end,:] = I[2,:]

    return I
end

"""
    function xz_up_ray(θ::AbstractFloat, ϕ::AbstractFloat, idz::Int, sign_x::Int,
                       sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                       α::Array{<:PerLength, 3}, atmos::Atmosphere)

Ray moving upwards, upwind point intersecting with xz plane. Assumes angle
in radians.
"""
function xz_up_ray(k::Vector{Float64}, idz::Int, sign_x::Int,
                   sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                   α::Array{<:PerLength, 3}, atmos::Atmosphere, n_sweeps::Int)

    nx = length(atmos.x)
    ny = length(atmos.y)

    Δy = atmos.y[2] - atmos.y[1]

    # Loop direction
    start_x, stop_x = range_bounds(sign_x, nx)
    start_y, stop_y = range_bounds(sign_y, ny)

    # Allocate space for intensity
    I = zero(I_0)

    I_upper = zero(I_0[end,:])

    # Z center and interpolation position
    z_centre = atmos.z[idz]
    idz_lower = idz-1

    # calculate the length to upwind position
    r = abs(Δy/k[3])
    z_increment = r*k[1]
    x_increment = r*k[2]

    z_upwind = z_centre + z_increment

    z_bounds = (atmos.z[idz_lower], z_centre)

    α_lower = α[idz_lower,:,:]
    α_upper = α[idz,:,:]

    S_lower = S_0[idz_lower,:,:]
    S_upper = S_0[idz,:,:]

    for sweep in 1:n_sweeps
        for idy in start_y:sign_y:stop_y
            y_centre = atmos.y[idy]
            # x upwind position
            idy_upwind = idy+sign_y
            y_upwind = atmos.y[idy_upwind]
            for idx in start_x:sign_x:stop_x
                # Centre point (c)
                x_centre = atmos.x[idx]

                # Lower corner to interpolate from
                idx_lower = idx-Int((sign_x+1)/2)
                idx_upper = idx_lower+1

                x_upwind = x_centre + x_increment

                x_bounds = (atmos.x[idx_lower], atmos.x[idx_upper])

                # Do linear interpolation to find τ and S

                # exinction at centre and upwind point
                # z is 1st coordinate, changes rows in vals. y is 2nd coordinate
                α_vals = [α_lower[idx_lower, idy_upwind]   α_lower[idx_upper, idy_upwind]
                          α_upper[idx_lower, idy_upwind]   α_upper[idx_upper, idy_upwind]]

                α_centre = α_upper[idx, idy]
                α_upwind = bilinear(z_upwind, x_upwind, z_bounds, x_bounds, α_vals)

                # Find the Δτ optical path from upwind to grid point
                Δτ_upwind = trapezoidal(r, α_centre, α_upwind)

                # Find source function at grid point and upwind point (intepolate)
                S_vals = [S_lower[idx_lower, idy_upwind]    S_lower[idx_upper, idy_upwind]
                          S_upper[idx_lower, idy_upwind]    S_upper[idx_upper, idy_upwind]]

                S_centre = S_upper[idx, idy]
                S_upwind = bilinear(z_upwind, x_upwind, z_bounds, x_bounds, S_vals)

                # Compute weights for formal solution
                α_ijk, β_ijk, expΔτ = linear_weights(Δτ_upwind)

                # Interpolate to find intensity at upwind point
                # I_0 is from the top
                I_vals = [I_0[idx_lower, idy_upwind]    I_0[idx_upper, idy_upwind]
                          I_upper[idx_lower]            I_upper[idx_upper]        ]

                I_upwind = bilinear(z_upwind, x_upwind, z_bounds, x_bounds, I_vals)

                # Integrate intensity from two-point quadrature
                I[idx, idy] = expΔτ*I_upwind + α_ijk*S_upwind + β_ijk*S_centre
            end
            # Update ghost zones
            I[1, idy] = I[end-1, idy]
            I[end, idy] = I[2, idy]

            # Update upper interpolate value
            I_upper = I[:, idy]
        end
    end

    # Update ghost zones
    I[:, 1] = I[:, end-1]
    I[:, end] = I[:, 2]

    return I
end

"""
    function xz_down_ray(θ::AbstractFloat, ϕ::AbstractFloat, idz::Int, sign_x::Int,
                           sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                           α::Array{<:PerLength, 3}, atmos::Atmosphere)

Ray moving downwards, upwind point intersecting with xz plane. Assumes angle
in radians.
"""
function xz_down_ray(k::Vector{Float64}, idz::Int, sign_x::Int,
                       sign_y::Int, I_0::Matrix{<:UnitsIntensity_λ}, S_0::Array{<:UnitsIntensity_λ, 3},
                       α::Array{<:PerLength, 3}, atmos::Atmosphere, n_sweeps::Int)

    nx = length(atmos.x)
    ny = length(atmos.y)

    Δy = atmos.y[2] - atmos.y[1]

    # Loop direction
    start_x, stop_x = range_bounds(sign_x, nx)
    start_y, stop_y = range_bounds(sign_y, ny)

    # Allocate space for intensity
    I = zero(I_0)

    I_lower = zero(I_0[end,:])

    # Z center and interpolation position
    z_centre = atmos.z[idz]
    idz_upper = idz+1

    # Length to plane
    # Δx = atmos.x[2] - atmos.x[1]

    # calculate the length to upwind position
    r = abs(Δy/k[3])
    z_increment = r*k[1]
    x_increment = r*k[2]

    z_upwind = z_centre + z_increment

    z_bounds = (z_centre, atmos.z[idz_upper])

    α_lower = α[idz,:,:]
    α_upper = α[idz_upper,:,:]

    S_lower = S_0[idz,:,:]
    S_upper = S_0[idz_upper,:,:]

    for sweep in 1:n_sweeps
        for idy in start_y:sign_y:stop_y
            y_centre = atmos.y[idy]
            # x upwind position
            idy_upwind = idy+sign_y
            y_upwind = atmos.y[idy_upwind]
            for idx in start_x:sign_x:stop_x
                # Centre point (c)
                x_centre = atmos.x[idx]

                # Lower corner to interpolate from
                idx_lower = idx-Int((sign_x+1)/2)
                idx_upper = idx_lower+1

                x_upwind = x_centre + x_increment

                x_bounds = (atmos.x[idx_lower], atmos.x[idx_upper])

                # Do linear interpolation to find τ and S

                # exinction at centre and upwind point
                # z is 1st coordinate, changes rows in vals. y is 2nd coordinate
                α_vals = [α_lower[idx_lower, idy_upwind]   α_lower[idx_upper, idy_upwind]
                          α_upper[idx_lower, idy_upwind]   α_upper[idx_upper, idy_upwind]]

                α_centre = α_upper[idx, idy]
                α_upwind = bilinear(z_upwind, x_upwind, z_bounds, x_bounds, α_vals)

                # Find the Δτ optical path from upwind to grid point
                Δτ_upwind = trapezoidal(r, α_centre, α_upwind)

                # Find source function at grid point and upwind point (intepolate)
                S_vals = [S_lower[idx_lower, idy_upwind]    S_lower[idx_upper, idy_upwind]
                          S_upper[idx_lower, idy_upwind]    S_upper[idx_upper, idy_upwind]]

                S_centre = S_upper[idx, idy]
                S_upwind = bilinear(z_upwind, x_upwind, z_bounds, x_bounds, S_vals)

                # Compute weights for formal solution
                α_ijk, β_ijk, expΔτ = linear_weights(Δτ_upwind)

                # Interpolate to find intensity at upwind point
                # I_0 is from the top
                I_vals = [I_lower[idx_lower]            I_lower[idx_upper]
                          I_0[idx_lower, idy_upwind]    I_0[idx_upper, idy_upwind]]

                I_upwind = bilinear(z_upwind, x_upwind, z_bounds, x_bounds, I_vals)

                # Integrate intensity from two-point quadrature
                I[idx, idy] = expΔτ*I_upwind + α_ijk*S_upwind + β_ijk*S_centre
            end

            # Update ghost zones
            I[1, idy] = I[end-1, idy]
            I[end, idy] = I[2, idy]

            # Update upper interpolate value
            I_lower = I[:, idy]
        end
    end

    # Update ghost zones
    I[:, 1] = I[:, end-1]
    I[:, end] = I[:, 2]

    return I
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
struct POMDP
    γ   # discount factor
    𝒮   # state space
    𝒜   # action space
    𝒪   # observation space
    T   # transition function
    R   # reward function
    O   # observation function
    TRO # sample transition, reward, and observation
end

ext: jl


text:
# rootsah for BigInt and BigFloat
function rootsah(pol::Union{Poly{BigInt},Poly{BigFloat}}, D::Vector{Float64})
    # COMPUTES: the roots of polynomials with all distinct real roots.
    # The computation is forward stable. The program uses SymArrow (arrowhead) companion matrix and
    # corresponding eig routine
    # D are barycentric coordinates - elements od D must interpolate the roots of P,
    # for example
    #              D=roots(polyder(pol))
    # RETURNS: roots E

    # Type is BigFloat
    T = BigFloat

    # τ = [1e2,1e2,1e2,1e2,1e2] or similar is the vector of tolerances for eig
    τ=[1e2,1e2,1e2,1e2,1e2]

    Dm=map(T,D)
    Dm=sort(Dm,rev=true)
    D=sort(D,rev=true)
    p=map(T,[pol[i] for i=0:1:length(pol)-1])
    p=p[end:-1:1]

    # Compute z of the arrowhead
    # First we compute values s=p(D) we use Horner sheme with double the working precision
    setprecision(512)
    n=length(p)-1
    s=Array{BigFloat}(undef,n-1)
    for i=1:n-1
        s[i]=one(T)
    end
    s=p[1]*s
    for i=2:n+1
        r=s.*Dm
        s=r.+p[i]
    end

    # Compute t's
    t=Array{BigFloat}(undef,n-1)
    for j=1:n-1
        h=one(T)
        for i=1:j-1
            g=Dm[j]-Dm[i]
            h=h*g
        end
        for i=j+1:n-1
            g=Dm[j]-Dm[i]
            h=h*g
        end
        t[j]=h
    end
    # Compute alphaD
    αD=p[2]/p[1]
    for i=1:n-1
        αD+=Dm[i]
    end
     αD=-αD
    #  Compute z
    zD=Array{BigFloat}(undef,n-1)
    for i=1:n-1
        zD[i]=sqrt((-s[i])/(t[i]*p[1]))
    end

    z=Array{Float64}(undef,n-1)
    for i=1:n-1
        z[i]=Float64(zD[i])
    end
    α=Float64(αD)
    A=SymArrow(D,z,α,n)
    E=zeros(Float64,n)
    Qout=zeros(Int64,n)
    for k=1:n
        E[k], Qout[k] = eigen( A,zD,αD,k,τ )
    end

    # reset the bigfloat_precision
    setprecision(256)
    E, Qout
end

#------------------------
#--------  Functions

#-------- Inverses

function inv(A::SymArrow{T},zD::Vector{BigFloat}, αD::BigFloat,i::Int64,τ::Vector{Float64}=[1e3, 10*length(A.D)]) where T
    # COMPUTES: inverse of a SymArrow matrix A, inv(A-A.D[i]*I) which is again SymArrow
    # uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # τ=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used

    n=length(A.D)
    D=Array{T}(undef,n)
    z=Array{T}(undef,n)
    wz=1/A.z[i]
    σ=A.D[i]
    for k=1:i-1
        D[k]=one(T)/(A.D[k]-σ)
        z[k]=-A.z[k]*D[k]*wz
    end
    for k=i+1:n
        D[k-1]=one(T)/(A.D[k]-σ)
        z[k-1]=-A.z[k]*D[k-1]*wz
    end
    a=A.a-A.D[i]

    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    nn=n-1
    for k=1:i-1
        D[k]>0.0 ? P+=A.z[k]^2*D[k] : Q+=A.z[k]^2*D[k]
    end
    for k=i:nn
        D[k]>0.0 ? P+=A.z[k+1]^2*D[k] : Q+=A.z[k+1]^2*D[k]
    end
    a<0 ? P-=a : Q-=a

    Kb=(P-Q)/abs(P+Q)
    Kz=maximum(abs,A.z)*abs(wz)
    if Kb<τ[1] ||  Kz<τ[2]
        b=(P+Q)*wz*wz
    else  # recompute in double the working precision
        Qout=1
        AD=map(BigFloat,A.D)
        σd=AD[i]
        Dd=[[AD[k]-σd for k=1:i-1];
            [AD[k]-σd for k=i+1:length(A.D)]]
        wzd=zD[i] #
        ad=αD-σd
        Pd,Qd=map(BigFloat,(0.0,0.0))
        for k=1:i-1
            Dd[k] > 0.0 ? Pd+=zD[k]^2/Dd[k] :
            Qd+=zD[k]^2/Dd[k]
        end
        for k=i:nn
            Dd[k] > 0.0 ? Pd+=zD[k+1]^2/Dd[k] :
            Qd+=zD[k+1]^2/Dd[k]
        end
        ad<0 ?   Pd-=ad : Qd-=ad
        bd=(Pd+Qd)/(wzd*wzd)
        b=Float64(bd)
    end

    if i<A.i
        # SymArrow([D[1:A.i-2],zero(T),D[A.i-1:end]],[z[1:A.i-2],wz,z[A.i-1:end]],b,i),Kb,Kz,Qout
        # This reduces memory allocation (4 vectors) and is 30-40% faster
        for k=n-1:-1:A.i-1
            D[k+1]=D[k]
            z[k+1]=z[k]
        end
        D[A.i-1]=zero(T)
        z[A.i-1]=wz
        SymArrow(D,z,b,i),Kb,Kz,Qout
    else
        # SymArrow([D[1:A.i-1],zero(T),D[A.i:end]],[z[1:A.i-1],wz,z[A.i:end]],b,i+1),Kb,Kz,Qout
        for k=n-1:-1:A.i
            D[k+1]=D[k]
            z[k+1]=z[k]
        end
        D[A.i]=zero(T)
        z[A.i]=wz
        SymArrow(D,z,b,i+1),Kb,Kz,Qout
    end
end # inv

function inv(A::SymArrow{T}, zD::Vector{BigFloat}, αD::BigFloat, σ::Float64, τ::Float64=1.0e3) where T
    # COMPUTES: inverse of the shifted SymArrow A, inv(A-shift*I) which is SymDPR1
    # uses DoubleDouble to compute ρ accurately, if needed.
    # τ is tolerance, usually 1e3,  0.0 forces Double, 1e50 would never use it
    # RETURNS: SymDPR1(D,u,ρ), Kρ, Qout
    # Kρ - condition Kρ, Qout = 1 / 0 - Double was / was not used

    n=length(A.D)
    D=Array{T}(undef,n+1)
    u=Array{T}(undef,n+1)
    for k=1:A.i-1
        D[k]=one(T)/(A.D[k]-σ)
        u[k]=A.z[k]*D[k]
    end
    D[A.i]=zero(T)
    u[A.i]=-one(T)
    for k=A.i:n
        D[k+1]=one(T)/(A.D[k]-σ)
        u[k+1]=A.z[k]*D[k+1]
    end

    # compute ρ and Kρ
    #--- compute the sum in a plain loop
    a=A.a-σ
    P=zero(T)
    Q=zero(T)
    Qout=0
    for k=1:A.i-1
        D[k]>0.0 ? P+=A.z[k]^2*D[k] : Q+=A.z[k]^2*D[k]
    end
    for k=A.i:n
        D[k+1]>0.0 ? P+=A.z[k]^2*D[k+1] : Q+=A.z[k]^2*D[k+1]
    end
    a<0 ? P-=a : Q-=a

    # Condition of ρ
    Kρ=(P-Q)/abs(P+Q)
    if Kρ<τ
        ρ=-1.0/(P+Q)
    else  # recompute in double the working precision
        Qout=1
        Pd,Qd=map(BigFloat,(0.0,0.0))
        AD=map(BigFloat,A.D)
        σd=map(BigFloat,σ)
        ad=αD-σd
        for k=1:A.i-1
            D[k]>0.0 ? Pd+=zD[k]^2/(AD[k]-σd) : Qd+=zD[k]^2/(AD[k]-σd)
        end
        for k=A.i:n
            D[k+1]>0.0 ? Pd+=zD[k]^2/(AD[k]-σd) : Qd+=zD[k]^2/(A.D[k]-σd)
        end
        ad<0 ? Pd-=ad : Qd-=ad
        r=map(BigFloat,1.0)/(Pd+Qd)
        ρ=Float64(-r)
    end

    # returns the following
    SymDPR1(D,u,ρ), Kρ, Qout
end # inv

function inv(A::SymArrow{T}, zD::Vector{BigFloat}, αD::BigFloat, σ::BigFloat) where T
    # COMPUTES: inverse of the shifted SymArrow A, inv(A-σ*I), which is a SymDPR1
    # here σ is BigFloat so it uses double the working precision to compute everything
    # RETURNS: SymDPR1(D1,u1,ρ), Qout
    # Qout = 1 on exit meaning Double was used

    n=length(A.D)
    D=Array{BigFloat}(undef,n+1)
    u=Array{BigFloat}(undef,n+1)
    for k=1:n
        D[k]=map(BigFloat,A.D[k])-σ
    end
    u[1:n]=zD
    a=αD-σ
    i=A.i

    oned=one(BigFloat)
    zerod=zero(BigFloat)
    for k=1:i-1
        u[k]=u[k]/D[k]
        D[k]=oned/D[k]
    end

    for k=i:n
        u[k+1]=u[k]/D[k]
        D[k+1]=oned/D[k+1]
    end
    D[i]=zerod
    u[i]=-oned

    # compute ρ and Kρ
    #--- compute the sum in a plain loop
    P,Q=zerod,zerod
    Qout=1
    for k=1:i-1
        D[k] > 0.0 ? P+=zD[k]*u[k] : Q+=zD[k]*u[k]
    end
    for k=i:n
        D[k+1] > 0.0 ? P+=zD[k]*u[k+1] : Q+=zD[k]*u[k+1]
    end
    a<0.0  ? P-=a : Q-=a
    r=oned/(P+Q)
    ρ=Float64(-r)

    # returns the following
    D₁=Array{T}(undef,n+1)
    u₁=Array{T}(undef,n+1)
    for k=1:n+1
        D₁[k]=Float64(D[k])
    end
    for k=1:n+1
        u₁[k]=Float64(u[k])
    end
    SymDPR1(D1,u1,ρ), Qout
end # inv

function  eigen( A::SymArrow{T},zD::Vector{BigFloat},αD::BigFloat,k::Int64,
    τ::Vector{Float64}=[1e3,10.0*length(A.D),1e3,1e3,1e3]) where T
    # COMPUTES: k-th eigenvalue of an ordered irreducible SymArrow
    # A = [Diagonal(D) z; z' alpha]
    # Specially designed to be used in the polynomial rootfinder rootsah !!!!
    # τ=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: λ, Qout
    # λ - k-th eigenvalue
    # Kb, Kz, Kν, Kρ - condition numbers
    # Qout = 1 / 0 - Double was / was not used

    # If Qout>0, quad precision was used
    # i was the shift index

    # Set the dimension
    n = length(A.D) + 1
    # Set all conditions initially to zero
    Kb,Kz,Kν,Kρ=0.0,0.0,0.0,0.0
    Qout=0
    v=zeros(T,n)
    # Kz is former kappa_nu

    # Determine the shift σ, the shift index i, and whether λ
    # is on the left or the right side of the nearest pole
    # Exterior eigenvalues (k = 1 or k = n):
    if k == 1
        σ,i,side = A.D[1],1,'R'
    elseif k==n
        σ,i,side = A.D[n-1],n-1,'L'
    else
        # Interior eigenvalues (k in (2,...n-1) ):
        # We need to compute Fmiddle carefully
        Dtemp = A.D.-A.D[k]
        middle = Dtemp[k-1]/2.0
        P=zero(Float64)
        Q=zero(Float64)
        P+=sum(A.z[1:k-1].^2 ./(Dtemp[1:k-1].-middle))
        Q+=sum(A.z[k:n-1].^2 ./(Dtemp[k:n-1].-middle))

        P₁=map(BigFloat,P)
        Q₁=map(BigFloat,Q)
        atemp = (αD-map(BigFloat,A.D[k]))-map(BigFloat,middle)
        atemp<0 ? P₁-=atemp : Q₁-=atemp
        # Fmiddle = (atemp-middle)-sum(zD.^2./(Dtemp-middle))
        Fmiddle=-(P₁+Q₁)
        σ,i,side = Fmiddle < zero(BigFloat) ? (A.D[k],k,'R') : (A.D[k-1],k-1,'L')
    end

    # Compute the inverse of the shifted matrix, A_i^(-1), Kb and Kz
    Ainv,Kb,Kz,Qout = inv(A,zD,αD,i)

    # Compute the eigenvalue of the inverse shifted matrix
    ν = bisect( Ainv,side )

    #  ν=fastzero([invD1; 0; invD2], [w1;wz;w2], b, side); # implement later
    #  [ν-ν₁]/ν, [ν-nueva]/ν, pause  # and compare
    if abs(ν)==Inf
        # this is nonstandard
        # Deflation in aheig (ν=Inf)
        v[i]=1.0
        λ=σ
    else
        # standard case, full computation
        # ν₁ is the F- or 1-norm of the inverse of the shifted matrix
        # ν₁0=maximum([sum(abs(Ainv.z))+abs(Ainv.a), maximum(abs(Ainv.D)+abs(Ainv.z))])
        ν₁=0.0
        for l=1:n-1
            ν₁=max(ν₁,abs(Ainv.D[l])+abs(Ainv.z[l]))
        end
        ν₁=max(sum(abs,Ainv.z)+abs(Ainv.a), ν₁)
        Kν=ν₁/abs(ν)
        if Kν<τ[3]
            # Accuracy is fine, compute the eigenvector
            μ = 1.0/ν
            # v=[ A.z./((A.D-σ)-μ);-1.0]
            for l=1:n-1
            	v[l] = A.z[l]/((A.D[l]-σ)-μ)
            end
            v[n]=-1.0
            λ = μ+σ
            normalize!(v)
        else
            # Remedies according to Remark 3 - we shift between original
            # eigenvalues and compute DPR1 matrix
            # 1/ν₁+σ, 1/ν+σ
            # println("Remedy 3 ")
            ν = side=='R' ? abs(ν) : -abs(ν)
            ν₁=-sign(ν)*ν₁
            σ₁=(ν₁+ν)/(2.0*ν*ν₁)+σ
            if count(isequal(σ₁),A.D)>0 # we came back with a pole
                # recompute σᵥ more accurately according with dekker
                setprecision(512)
                σᵥ=(ν₁+ν)/(BigFloat(2.0)*ν*ν₁)+σ
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv,Qout1=inv(A,σᵥ) # Ainv is Float64
                ν₁=bisect(Ainv,side)
                μ₁ = 1.0/ν₁
                D0=map(BigFloat,A.D)-σᵥ
                D1=[D0[l] for l=1:length(D0)]
                if count(isequal(μ₁),D1)>0
                    ind=findall(isequal(μ₁),D1)
                    v=zeros(n)
                    v[ind]=1.0
                else
                    v=[ A.z./(D1-μ₁);-1.0]
                end
                # Shift the eigenvalue back in Double
                lam = BigFloat(1.0)/ν₁+σᵥ
                # Return this
                λ = lam
                normalize!(v)
            else
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv, Kρ,Qout₁=inv(A,σ₁,τ[4]) # Ainv is Float64
                # Compute the eigenvalue by bisect for DPR1
	            # Note: instead of bisection could use dlaed4 (need a wrapper) but
	            # it is not faster. There norm(u)==1
                ν₁= bisect(Ainv,side)
                μ₁=1.0/ν₁
                # standard v
                v=[ A.z./((A.D-σ₁)-μ₁);-1.0]
                # Return this - shift the eigenvalue back and normalize the vector
                λ = μ₁+σ₁
                normalize!(v)
            end
            Qout==1 && (Qout=Qout+2)
        end

        # Remedy according to Remark 1 - we recompute the the eigenvalue
        # near zero from the inverse of the original matrix (a DPR1 matrix).
        if (abs(A.D[i])+abs(1.0/ν))/abs(λ)>τ[5]
            if k==1 && A.D[1]<0.0 || k==n && A.D[n-1]>0.0 || i<n-1 && side=='L' && sign(A.D[i])+sign(A.D[i+1])==0 ||
                side=='R' && sign(A.D[i])+sign(A.D[i-1])==0
                # println("Remedy 1 ")
                # Compute the inverse of the original arrowhead (DPR1)
                Ainv,Kρ,Qout₁ = inv(A,0.0,τ[4]) # Ainv is Float64
                Qout==1 && (Qout=Qout+4)
                if isinf(Ainv.r) || isnan(Ainv.r)
                    λ=0.0
                else
                    # Here we do not need bisection. We compute the Rayleigh
                    # quotient by using already computed vectors which is
                    # componentwise accurate
                    ν₁=sum(v.^2 .*Ainv.D)+Ainv.r*sum(v.*Ainv.u)^2;
                    λ=1.0/ν₁
                end
            end
        end
    end

    # Return this
    Float64(λ), Qout
end # eigen (k)

function rootsWDK(p::Poly{T},x₀::Vector{T1},steps) where {T,T1}
    # Polynomial roots with Weierstrass, Durand, Kerner method
    n=length(x₀)
    x=x₀
    for k=1:steps
        for i=1:n
            x[i]=x[i]-polyval(p,x[i])/(p[end]*prod(x[i].-x[collect([1:i-1;i+1:n])]))
        end
    end
    x
end


text:
"""
    sign_cadlag(x::N)::N where {N<:Real}

This function works like the sign function but is ``1`` for input ``0``.

### Input

- `x` -- real scalar

### Output

``1`` if ``x ≥ 0``, ``-1`` otherwise.

### Notes

This is the sign function right-continuous at zero (see
[càdlàg function](https://en.wikipedia.org/wiki/C%C3%A0dl%C3%A0g)).
It can be used with vector-valued arguments via the dot operator.

### Examples

```jldoctest
julia> LazySets.sign_cadlag.([-0.6, 1.3, 0.0])
3-element Array{Float64,1}:
 -1.0
  1.0
  1.0
```
"""
function sign_cadlag(x::N)::N where {N<:Real}
    return x < zero(x) ? -one(x) : one(x)
end

"""
    substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

### Input

- `substitution` -- substitution (a mapping from an index to a new value)
- `x`            -- vector

### Output

A fresh vector corresponding to `x` after `substitution` was applied.
"""
function substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}
    return substitute!(substitution, copy(x))
end

"""
    substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

### Input

- `substitution` -- substitution (a mapping from an index to a new value)
- `x`            -- vector (modified in this function)

### Output

The same (but see the Notes below) vector `x` but after `substitution` was
applied.

### Notes

The vector `x` is modified in-place if it has type `Vector` or `SparseVector`.
Otherwise, we first create a new `Vector` from it.
"""
function substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}
    return substitute!(Vector(x), substitution)
end

function substitute!(substitution::Dict{Int, T},
                     x::Union{Vector{T}, SparseVector{T}}) where {T}
    for (index, value) in substitution
        x[index] = value
    end
    return x
end

"""
    reseed(rng::AbstractRNG, seed::Union{Int, Nothing})::AbstractRNG

Reset the RNG seed if the seed argument is a number.

### Input

- `rng`  -- random number generator
- `seed` -- seed for reseeding

### Output

The input RNG if the seed is `nothing`, and a reseeded RNG otherwise.
"""
function reseed(rng::AbstractRNG, seed::Union{Int, Nothing})::AbstractRNG
    if seed != nothing
        return Random.seed!(rng, seed)
    end
    return rng
end

"""
    StrictlyIncreasingIndices

Iterator over the vectors of `m` strictly increasing indices from 1 to `n`.

### Fields

- `n` -- size of the index domain
- `m` -- number of indices to choose (resp. length of the vectors)

### Notes

The vectors are modified in-place.

The iterator ranges over ``\\binom{n}{m}`` (`n` choose `m`) possible vectors.

This implementation results in a lexicographic order with the last index growing
first.

### Examples

```jldoctest
julia> for v in LazySets.StrictlyIncreasingIndices(4, 2)
           println(v)
       end
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]
```
"""
struct StrictlyIncreasingIndices
    n::Int
    m::Int

    function StrictlyIncreasingIndices(n::Int, m::Int)
        @assert n >= m > 0 "require n >= m > 0"
        new(n, m)
    end
end

Base.eltype(::Type{StrictlyIncreasingIndices}) = Vector{Int}
Base.length(sii::StrictlyIncreasingIndices) = binomial(sii.n, sii.m)

# initialization
function Base.iterate(sii::StrictlyIncreasingIndices)
    v = [1:sii.m;]
    return (v, v)
end

# normal iteration
function Base.iterate(sii::StrictlyIncreasingIndices, state::AbstractVector{Int})
    v = state
    i = sii.m
    diff = sii.n
    if i == diff
        return nothing
    end
    while v[i] == diff
        i -= 1
        diff -= 1
    end
    # update vector
    v[i] += 1
    for j in i+1:sii.m
        v[j] = v[j-1] + 1
    end
    # detect termination: first index has maximum value
    if i == 1 && v[1] == (sii.n - sii.m + 1)
        return (v, nothing)
    end
    return (v, v)
end

# termination
function Base.iterate(sii::StrictlyIncreasingIndices, state::Nothing)
    return nothing
end

"""
    subtypes(interface, concrete::Bool)

Return the concrete subtypes of a given interface.

### Input

- `interface` -- an abstract type, usually a set interface
- `concrete`  -- if `true`, seek further the inner abstract subtypes of the given
                 interface, otherwise return only the direct subtypes of `interface`

### Output

A list with the subtypes of the abstract type `interface`, sorted alphabetically.

### Examples

Consider the `AbstractPolytope` interface. If we include the abstract subtypes
of this interface,

```jldoctest subtypes
julia> using LazySets: subtypes

julia> subtypes(AbstractPolytope, false)
4-element Array{Any,1}:
 AbstractCentrallySymmetricPolytope
 AbstractPolygon
 HPolytope
 VPolytope
```

We can use this function to obtain the concrete subtypes of
`AbstractCentrallySymmetricPolytope` and `AbstractPolygon` (further until all
concrete types are obtained), using the `concrete` flag:

```jldoctest subtypes
julia> subtypes(AbstractPolytope, true)
14-element Array{Type,1}:
 Ball1
 BallInf
 HPolygon
 HPolygonOpt
 HPolytope
 Hyperrectangle
 Interval
 LineSegment
 Singleton
 SymmetricIntervalHull
 VPolygon
 VPolytope
 ZeroSet
 Zonotope
```
"""
function subtypes(interface, concrete::Bool)

    subtypes_to_test = subtypes(interface)

    # do not seek the concrete subtypes further
    if !concrete
        return sort(subtypes_to_test, by=string)
    end

    result = Vector{Type}()
    i = 0
    while i < length(subtypes_to_test)
        i += 1
        subtype = subtypes_to_test[i]
        new_subtypes = subtypes(subtype)
        if isempty(new_subtypes)
            # base type found
            push!(result, subtype)
        else
            # yet another interface layer
            append!(subtypes_to_test, new_subtypes)
        end
    end
    return sort(result, by=string)
end


text:
"""
    imageaxis(f)

Create a `Makie.Axis` with no interactivity, decorations and aspect distortion
suitable for showing images on.
"""
function imageaxis(f; kwargs...)
    ax = Makie.Axis(f; kwargs...)
    ax.aspect = DataAspect()
    ax.xzoomlock = true
    ax.yzoomlock = true
    ax.xrectzoom = false
    ax.yrectzoom = false
    ax.panbutton = nothing
    ax.xpanlock = true
    ax.ypanlock = true
    ax.bottomspinevisible = false
    ax.leftspinevisible = false
    ax.rightspinevisible = false
    ax.topspinevisible = false
    hidedecorations!(ax)

    return ax
end

@recipe(PlotImage, image) do scene
    Attributes(
        alpha = 1,
        interpolate = false,
    )
end

function Makie.plot!(plot::PlotImage)
    im = plot[:image]
    rim = @lift alphacolor.(copy(rotr90($im)), $(plot.attributes[:alpha]))
    image!(plot, rim; plot.attributes...)
    return plot
end

@recipe(PlotMask, mask, classes) do scene
    Attributes()
end

function Makie.plot!(plot::PlotMask; kwargs...)
    mask = plot[:mask]
    classes = try
        classes = plot[:classes]
    catch
        classes = @lift unique($mask)
    end
    im = @lift maskimage($mask, $classes)
    plotimage!(plot, im; alpha = 1, plot.attributes...)
    return plot
end

function maskimage(mask, classes)
    classtoidx = Dict(class => i for (i, class) in enumerate(classes))
    colors = distinguishable_colors(length(classes), transform=deuteranopic)
    return map(x -> colors[classtoidx[x]], mask)
end

maskimage(mask::AbstractArray{<:Gray{T}}, args...) where T =
    maskimage(reinterpret(T, mask), args...)
maskimage(mask::AbstractArray{<:Normed{T}}, args...) where T =
    maskimage(reinterpret(T, mask), args...)


text:
using BinaryBuilder

name = "JpegTurbo"
version = v"2.1.2"

# Collection of sources required to build Ogg
sources = [
    ArchiveSource("https://github.com/libjpeg-turbo/libjpeg-turbo/archive/$(version).tar.gz",
                  "e7fdc8a255c45bc8fbd9aa11c1a49c23092fcd7379296aeaeb14d3343a3d1bed"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libjpeg-turbo-*/

mkdir build
cd build

cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libjpeg", :libjpeg),
    LibraryProduct("libturbojpeg", :libturbojpeg),
    ExecutableProduct("cjpeg", :cjpeg),
    ExecutableProduct("djpeg", :djpeg),
    ExecutableProduct("jpegtran", :jpegtran),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Avengers; ASSEMBLE!
    HostBuildDependency("YASM_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
module RulesDotProductMarginalsTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_marginalrules
import LinearAlgebra: dot

@testset "marginalrules:DotProduct" begin

     @testset "in1_in2: (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(1.5, 1.0))),
            (input = (m_out = NormalMeanPrecision(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalWeightedMeanPrecision(0.5, 0.5), m_in1 = PointMass(-2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(-2.0), in2 = NormalWeightedMeanPrecision(-0.5, 2.5))),
        ]

        @test_marginalrules [ with_float_conversions = false ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(1.5, 1.0))),
            (input = (m_out = NormalMeanPrecision(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalWeightedMeanPrecision(0.5, 0.5), m_in1 = PointMass(-2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(-2.0), in2 = NormalWeightedMeanPrecision(-0.5, 2.5))),
        ]

        @test_marginalrules [ with_float_conversions = true ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(1.5, 1.0))),
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanPrecision(1.0, 0.5), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = PointMass(-2.0), m_in2 = NormalWeightedMeanPrecision(0.5, 0.5), meta = NoCorrection()), output = (in1 = PointMass(-2.0), in2 = NormalWeightedMeanPrecision(-0.5, 2.5))),
        ]

        @test_marginalrules [ with_float_conversions = false ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(1.5, 1.0))),
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanPrecision(1.0, 0.5), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = PointMass(-2.0), m_in2 = NormalWeightedMeanPrecision(0.5, 0.5), meta = TinyCorrection()), output = (in1 = PointMass(-2.0), in2 = NormalWeightedMeanPrecision(-0.5, 2.5))),
        ]

    end

    @testset "in1_in2: (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(2.0), meta = NoCorrection()), output = (in1 = NormalWeightedMeanPrecision(2.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalMeanPrecision(3.0, 1.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(2.0), meta = NoCorrection()), output = (in1 = NormalWeightedMeanPrecision(6.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalWeightedMeanPrecision(4.0, 1.0), m_in1 = NormalMeanVariance(1.0, 3.0), m_in2 = PointMass(2.0), meta = NoCorrection()), output = (in1 = NormalWeightedMeanPrecision(25/3, 13/3), in2 = PointMass(2.0))),
        ]

        @test_marginalrules [ with_float_conversions = false ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = (in1 = NormalWeightedMeanPrecision(2.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalMeanPrecision(3.0, 1.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = (in1 = NormalWeightedMeanPrecision(6.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalWeightedMeanPrecision(4.0, 1.0), m_in1 = NormalMeanVariance(1.0, 3.0), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = (in1 = NormalWeightedMeanPrecision(25/3, 13/3), in2 = PointMass(2.0))),
        ]

        ##

        @test_marginalrules [ with_float_conversions = true ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(2.0), meta = NoCorrection()), output = (in1 = NormalWeightedMeanPrecision(2.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalMeanVariance(3.0, 1.0), m_in1 = NormalMeanPrecision(1.0, 0.5), m_in2 = PointMass(2.0), meta = NoCorrection()), output = (in1 = NormalWeightedMeanPrecision(6.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalMeanVariance(4.0, 1.0), m_in1 = NormalWeightedMeanPrecision(1.0 / 3.0, 1.0 / 3.0), m_in2 = PointMass(2.0), meta = NoCorrection()), output = (in1 = NormalWeightedMeanPrecision(25/3, 13/3), in2 = PointMass(2.0))),
        ]

        @test_marginalrules [ with_float_conversions = false ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = (in1 = NormalWeightedMeanPrecision(2.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalMeanVariance(3.0, 1.0), m_in1 = NormalMeanPrecision(1.0, 0.5), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = (in1 = NormalWeightedMeanPrecision(6.5, 4.5), in2 = PointMass(2.0))),
            (input = (m_out = NormalMeanVariance(4.0, 1.0), m_in1 = NormalWeightedMeanPrecision(1.0 / 3.0, 1.0 / 3.0), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = (in1 = NormalWeightedMeanPrecision(25/3, 13/3), in2 = PointMass(2.0))),
        ]

        ## ##

        @test_marginalrules [ with_float_conversions = true ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalMeanPrecision(3.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(6.5, 4.5))),
            (input = (m_out = NormalWeightedMeanPrecision(4.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 3.0), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(25/3, 13/3))),
        ]

        @test_marginalrules [ with_float_conversions = false ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalMeanPrecision(3.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(6.5, 4.5))),
            (input = (m_out = NormalWeightedMeanPrecision(4.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 3.0), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(25/3, 13/3))),
        ]

        # ##

        @test_marginalrules [ with_float_conversions = true ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalMeanVariance(3.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanPrecision(1.0, 0.5), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(6.5, 4.5))),
            (input = (m_out = NormalMeanVariance(4.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(1.0 / 3.0, 1.0 / 3.0), meta = NoCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(25/3, 13/3))),
        ]

        @test_marginalrules [ with_float_conversions = false ] typeof(dot)(:in1_in2) [
            (input = (m_out = NormalMeanVariance(1.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(1.0, 2.0), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(2.5, 4.5))),
            (input = (m_out = NormalMeanVariance(3.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanPrecision(1.0, 0.5), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(6.5, 4.5))),
            (input = (m_out = NormalMeanVariance(4.0, 1.0), m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(1.0 / 3.0, 1.0 / 3.0), meta = TinyCorrection()), output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(25/3, 13/3))),
        ]

    end

end

end


text:
file = open("results/test.txt", "a")
i = [1,2,3]
print(file, "MaxSens - test " * string(i) * " - Time: ")
close(file)


text:
#=

# 107: 1D Nonlinear Storage
([source code](SOURCE_URL))

This equation comes from the transformation of the nonlinear diffuision equation.
```math
\partial_t u^\frac{1}{m} -\Delta u = 0
```
in $\Omega=(-1,1)$ with homogeneous Neumann boundary conditions.
We can derive an exact solution from the Barenblatt solution of the previous
example.

=#

module Example107_NonlinearStorage1D
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function barenblatt(x,t,m)
    tx=t^(-1.0/(m+1.0))
    xx=x*tx
    xx=xx*xx
    xx=1- xx*(m-1)/(2.0*m*(m+1));
    if xx<0.0
        xx=0.0
    end
    return tx*xx^(1.0/(m-1.0))
end

function main(;n=20,m=2.0,Plotter=nothing,verbose=false, unknown_storage=:sparse,tend=0.01,tstep=0.0001)

    ## Create a one-dimensional discretization
    h=1.0/convert(Float64,n/2)
    X=collect(-1:h:1)
    grid=VoronoiFVM.Grid(X)

    ## Flux function which describes the flux
    ## between neigboring control volumes
    function flux!(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=u[1,1]-u[1,2]
    end

    ϵ=1.0e-10

    ## Storage term
    ## This needs to be regularized as its derivative
    ## at 0 is infinity
    function storage!(f,u,node)
        f[1]=(ϵ+u[1])^(1.0/m)
    end

    ## Create a physics structure
    physics=VoronoiFVM.Physics(
        flux=flux!,
        storage=storage!)

    ## Create a finite volume system - either
    ## in the dense or  the sparse version.
    ## The difference is in the way the solution object
    ## is stored - as dense or as sparse matrix
    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    ## Add species 1 to region 1
    enable_species!(sys,1,[1])

    ## Create a solution array
    inival=unknowns(sys)
    solution=unknowns(sys)
    t0=0.001

    ## Broadcast the initial value
    inival[1,:].=map(x->barenblatt(x,t0,m)^m,X)
    solution.=inival

    ## Create solver control info
    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    time=t0
    p=GridVisualizer(Plotter=Plotter,layout=(2,1))
    while time<tend
        time=time+tstep
        solve!(solution,inival,sys,control=control,tstep=tstep)
        inival.=solution
        scalarplot!(p[1,1],grid,solution[1,:],title=@sprintf("numerical, t=%.5f",time),clear=true)
        scalarplot!(p[2,1],grid,map(x->barenblatt(x,time,m)^m,grid),title=@sprintf("exact, t=%.4f",time),clear=true)
        sleep(1.0e-5)
#        yield()
    end
    return sum(solution)
end

function test()
    testval=173.84998139534142
    main(unknown_storage=:sparse) ≈ testval && main(unknown_storage=:dense) ≈ testval
end

end


text:
module HipparchusWrapper

using JavaCall

const CalculusFieldElement = @jimport org.hipparchus.CalculusFieldElement
const Class = @jimport java.lang.Class
const Field = @jimport org.hipparchus.Field
const FieldElement = @jimport org.hipparchus.FieldElement
const FieldSinCos = @jimport org.hipparchus.util.FieldSinCos
const Object = @jimport java.lang.Object
const RealFieldElement = @jimport org.hipparchus.RealFieldElement

export CalculusFieldElement
export Field
export FieldElement
export RealFieldElement
export abs
export acos
export acosh
export add
export asin
export asinh
export atan
export atan2
export atanh
export cbrt
export ceil
export copy_sign
export cos
export cosh
export divide
export exp
export expm1
export floor
export get_exponent
export get_field
export get_one
export get_real
export get_runtime_class
export get_zero
export hypot
export is_infinite
export is_na_n
export linear_combination
export log
export log10
export log1p
export multiply
export negate
export new_instance
export pow
export reciprocal
export remainder
export rint
export root_n
export round
export scalb
export signum
export sin
export sin_cos
export sinh
export sqrt
export subtract
export tan
export tanh
export to_degrees
export to_radians

include("AnalysisWrapper/AnalysisWrapper.jl")
include("ComplexWrapper/ComplexWrapper.jl")
include("DfpWrapper/DfpWrapper.jl")
include("DistributionWrapper/DistributionWrapper.jl")
include("FilteringWrapper/FilteringWrapper.jl")
include("FittingWrapper/FittingWrapper.jl")
include("FractionWrapper/FractionWrapper.jl")
include("GeometryWrapper/GeometryWrapper.jl")
include("LinearWrapper/LinearWrapper.jl")
include("OdeWrapper/OdeWrapper.jl")
include("OptimWrapper/OptimWrapper.jl")
include("RandomWrapper/RandomWrapper.jl")
include("SpecialWrapper/SpecialWrapper.jl")
include("StatWrapper/StatWrapper.jl")
include("UtilWrapper/UtilWrapper.jl")
include("calculus_field_element.jl")
include("field.jl")
include("field_element.jl")
include("real_field_element.jl")

end


text:
using LOLTools.GameClient
using .GameClient: lol_replay_server

endpoint = lol_replay_server(host="192.168.0.9")
game = GameClient.replay_game(; endpoint=endpoint)
@info :game game
# particles = GameClient.replay_particles(; endpoint=endpoint)
# @info :particles (keys(particles)[1:3], getfield.(Ref(particles), keys(particles)[1:3]))
playback = GameClient.replay_playback(; endpoint=endpoint)
@info :playback playback
render = GameClient.replay_render(; endpoint=endpoint)
@info :render render
recording = GameClient.replay_recording(; endpoint=endpoint)
@info :recording recording
sequence = GameClient.replay_sequence(; endpoint=endpoint)
@info :sequence sequence
v2 = GameClient.swagger_v2_swagger(; endpoint=endpoint)
v3 = GameClient.swagger_v3_openapi(; endpoint=endpoint)
@info :info v2.info == v3.info == (description = "League of Legends Game Client", title = "LoLClient", version = "1.0.0")


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

##########################################################################################92

"""
    valid_coordinates(img::HDRimage, x::Int, y::Int) :: Bool

Return `true` if (`x`, `y`) are valid coordinates for the
2D matrix of the `HDRimage`, else return `false`.

See also: [`HDRimage`](@ref)
"""
function valid_coordinates(img::HDRimage, x::Int, y::Int)
    x>=0 && y>=0 && x<img.width && y<img.height
end

"""
    pixel_offset(img::HDRimage, x::Int, y::Int) :: Int64

Return the index in the 1D array of the specified pixel (`x`, `y`)
for the given `HDRimage`.
Internally checks also if  (`x`, `y`)  are valid coordinates for `img` through
the `valid_coordinates` function.

See also: [`valid_coordinates`](@ref), [`HDRimage`](@ref)
"""
function pixel_offset(img::HDRimage, x::Int, y::Int)
    @assert valid_coordinates(img, x, y)
    return y*img.width + (x+1)
end

"""
    get_pixel(img::HDRimage, x::Int, y::Int) :: RBG{Float32}

Return the `RBG{Float32}` color for the (`x`, `y`) pixel in the
given `HDRimage`, obtained through the `pixel_offset` function.
The indexes for a `HDRimage` pixel matrix
with `w` width and `h` height follow this sketch:
```ditaa
|  (0,0)    (0,1)    (0,2)   ...   (0,w-1)  |
|  (1,0)    (1,1)    (1,2)   ...   (1,w-1)  |
|   ...      ...      ...    ...     ...    |
| (h-1,0)  (h-1,1)  (h-1,2)  ...  (h-1,w-1) |
```

See also: [`pixel_offset`](@ref), [`HDRimage`](@ref)
"""
function get_pixel(img::HDRimage, x::Int, y::Int)
    return img.rgb_m[pixel_offset(img, x, y)]
end

function set_pixel(img::HDRimage, x::Int, y::Int, c::RGB{Float32})
    img.rgb_m[pixel_offset(img, x,y)] = c
    return nothing
end

function set_pixel(img::HDRimage, x::Int, y::Int, c::RGB{T}) where {T}
    img.rgb_m[pixel_offset(img,x,y)] = convert(RGB{Float32}, c)
    return nothing
end

"""
    set_pixel(img::HDRimage, x::Int, y::Int, c::RGB{Float32})
    set_pixel(img::HDRimage, x::Int, y::Int, c::RGB{T}) where {T}

Set the new RGB color `c` for the (`x`, `y`) pixel in
the given `HDRimage`, acceded through the `pixel_offset` function.
The indexes for a `HDRimage` pixel
matrix with `w` width and `h` height follow this sketch:
```ditaa
|  (0,0)    (0,1)    (0,2)   ...   (0,w-1)  |
|  (1,0)    (1,1)    (1,2)   ...   (1,w-1)  |
|   ...      ...      ...    ...     ...    |
| (h-1,0)  (h-1,1)  (h-1,2)  ...  (h-1,w-1) |
```

If `c` is of a type `RGB{T}` where `T ≠ Float32`, it's called the
`convert(RGB{Float32}, c)` function, which raises an exception if the
conversion is not possible.

See also: [`pixel_offset`](@ref), [`HDRimage`](@ref)
"""
set_pixel

##########################################################################################92

"""
    write(io::IO, img::HDRimage)

Write the given `HDRimage` image using the given `IO` stream as a PFM file.
The `endianness` used for writing the file is Little Endian (`-1.0`).

See also: [`get_pixel`](@ref), [`HDRimage`](@ref)
"""
function write(io::IO, img::HDRimage)
    endianness=-1.0
    w = img.width
    h = img.height

    # The PFM header, as a Julia string (UTF-8)
    header = "PF\n$w $h\n$endianness\n"

    # Convert the header into a sequence of bytes
    # PS: transcode writes a "raw" sequence of bytes (8bit)
    bytebuf = transcode(UInt8, header)

    # Write on io the header in binary code
    # PS: reinterpret  writes a "raw" sequence of bytes (8bit)
    write(io, reinterpret(UInt8, bytebuf))

    # Write the image (bottom-to-up, left-to-right)
    for y in h-1:-1:0, x in 0:w-1
        color = get_pixel(img, x, y)
        write(io, reinterpret(UInt8,  [color.r]))
        write(io, reinterpret(UInt8,  [color.g]))
        write(io, reinterpret(UInt8,  [color.b]))
    end

end

##########################################################################################92

"""
    InvalidPfmFileFormat <: Exception (var::String)

Self-made exception raised by thne following function
in case a reading operation of a PFM file fails:

- [`parse_img_size`](@ref)

- [`parse_endianness`](@ref)

- [`read_line`](@ref)

- [`read_float`](@ref)

- [`read(::IO, ::Type{HDRimage})`](@ref)

"""
struct InvalidPfmFileFormat <: Exception
    var::String
end

"""
    parse_img_size(line::String) :: (Int64, Int64)

Return the size `(width, height)` parsed from a given `String`, throwing
`InvalidPfmFileFormat` exception if encounters invalid values.
It works inside the `read(::IO, ::Type{HDRimage})` function.

See also: [`read(::IO, ::Type{HDRimage})`](@ref),
[`InvalidPfmFileFormat`](@ref)
"""
function parse_img_size(line::String)
    elements = split(line, " ")
    length(elements) == 2 || throw(InvalidPfmFileFormat(
                                "invalid image size specification: "*
                                "$(length(elements)) instead of 2"))

    try
        width, height = convert.(Int, parse.(Float64, elements))
        (width > 0 && height > 0) || throw(ErrorException)
        return width, height
    catch e
        isa(e, InexactError) || throw(InvalidPfmFileFormat(
                                    "cannot convert width/height "*
                                    "$(elements) to Tuple{Int, Int}")
        )
        isa(e, ErrorException) || throw(InvalidPfmFileFormat(
                                    "width/height cannot be negative,"*
                                    "but in $(elements) at least one of them is <0.")
        )
    end

end

"""
    parse_endianness(ess::String) :: Float64

Return the endianness parsed from a given `String`, throwing
`InvalidPfmFileFormat` exception if encounters an invalid value.
It works inside the `read(::IO, ::Type{HDRimage})` function.

See also: [`read(::IO, ::Type{HDRimage})`](@ref),
[`InvalidPfmFileFormat`](@ref)
"""
function parse_endianness(ess::String)
    try
        val = parse(Float64, ess)
        val in [1.0, -1.0] || throw(InvalidPfmFileFormat(
                                "invalid endianness in PFM file:"*
                                " $(val) instead of +1.0 or -1.0.\n")
        )
        return val
    catch e
        throw(InvalidPfmFileFormat(
            "missing endianness in PFM file: $ess instead of ±1.0"))
    end
end

"""
    read_float(io::IO, ess::Float64) :: Float32

Return a `Float32`, readed from the given IO stream `io` with the
(required) endianness `ess`; `ess` must be `+1.0` (Big Endian) or
`-1.0` (Little Endian).
Controls also if there are enough bits in order to form a `Float32`,
otherwise throw `InvalidPfmFileFormat`.
It works inside the `read(::IO, ::Type{HDRimage})` function.

See also: [`read(::IO, ::Type{HDRimage})`](@ref),
[`InvalidPfmFileFormat`](@ref)
"""
function read_float(io::IO, ess::Float64)
    # controllo che in ingresso abbia una stringa che sia cnovertibile in Float32
    ess == 1.0 || ess == -1.0 || throw(InvalidPfmFileFormat(
                                    "endianness $ess not acceptable."))
    try
        value = read(io, Float32)
        ess == 1.0 ? value=ntoh(value) : value=ltoh(value) # convert machine's endianness
        return value
    catch e
        throw(InvalidPfmFileFormat("color is not Float32, it's a $(typeof(io))"))
    end
end

"""
    read_line(io::IO) :: String

Reads a line from the file whose the given IO object `io` refers to.
Do understand when the file is ended and when a new line begins.
Return the readed line as a `String`.
It works inside the `read(::IO, ::Type{HDRimage})` function.

See also: [`read(::IO, ::Type{HDRimage})`](@ref),
[`InvalidPfmFileFormat`](@ref)
"""
function read_line(io::IO)
    result = b""
    while eof(io) == false
        cur_byte = read(io, UInt8)
        if [cur_byte] in [b"", b"\n"]
            return String(result)
        end
        result = vcat(result, cur_byte)
    end
    return String(result)
end

"""
    read(io::IO, ::Type{HDRimage}) :: HDRimage

Read a PFM image from a stream object `io`, and return a `HDRimage`
object containing the image. If an error occurs, raise a
`InvalidPfmFileFormat` exception.
Calls internally the following functions:

- [`parse_img_size`](@ref)

- [`parse_endianness`](@ref)

- [`read_line`](@ref)

- [`read_float`](@ref)

See also:  [`parse_img_size`](@ref), [`parse_endianness`](@ref),
[`read_line`](@ref), [`read_float`](@ref), [`InvalidPfmFileFormat`](@ref),
[`HDRimage`](@ref)
"""
function read(io::IO, ::Type{HDRimage})
    # magic number
    magic = read_line(io)
    magic == "PF" || throw(InvalidPfmFileFormat("invalid magic number in PFM file:"*
                                                 "$(magic) instead of 'PF'.\n"))

    # image dimension
    img_size = read_line(io)
    typeof(parse_img_size(img_size)) == Tuple{Int,Int} ||
        throw(InvalidPfmFileFormat(
            "invalid img size in PFM file: $(parse_img_size(img_size)) is "*
            "$( typeof(parse_img_size(img_size)) ) instead of 'Tuple{UInt,UInt}'.\n"))
    (width, height) = parse_img_size(img_size)

    # endianness
    endianness = parse_endianness(read_line(io))
    endianness in [1.0, -1.0] || throw(InvalidPfmFileFormat(
                                "invalid endianness in PFM file: "*
                                "$(endianness) instead of +1.0 or -1.0.\n"))

    # color matrix
    result = HDRimage(width, height)
    for y in height-1:-1:0, x in 0:width-1

        (r,g,b) = [read_float(io, endianness) for i in 0:2]
        set_pixel(result, x, y, RGB(r,g,b) )
    end

    return result
end

##########################################################################################92

"""
    parse_command_line(args::Vector{String}) ::
        (String, String, Float64, Float64)

Interpret the command line when the main is executed,
and manage eventual argument errors.

## Input

A `args::Vector{String})` with length 2, 3 or 4.

## Returns

A tuple `(infile, outfile, a, γ)` containing:

- `infile::String` : first string (required), is the input file
  name, which must be a PFM format.

- `outfile::String` : second string (required), is the output filename;
  its format can be PNG or TIFF.

- `a::Float64` : third argument (optional), is the scale factor for
  luminosity correction (default 0.18), passed to `normalize_image!`

- `γ::Float64` : fourth argument (optional), is the gamma factor for
  screen correction (default 1.0), passed to `γ_correction!`

See also : [`normalize_image!`](@ref), [`γ_correction!`](@ref)
"""
function parse_command_line(args::Vector{String})
    (isempty(args) || length(args)==1 || length(args)>4) && throw(Exception)
    infile = nothing; outfile = nothing; a=0.18; γ=1.0
    try
        infile = args[1]
        outfile = args[2]
        open(infile, "r") do io
            read(io, UInt8)
        end
    catch e
        throw(RuntimeError("invalid input file: $(args[1]) does not exist"))
    end

    if length(args)>2
        try
            a = parse(Float64, args[3])
            a > 0. || throw(Exception)
        catch e
            throw(InvalidArgumentError(
                "invalid value for a: $(args[3])  must be a positive number"))
        end

        if length(args) == 4
            try
                γ = parse(Float64, args[4])
                γ > 0. || throw(Exception)
            catch e
                throw(InvalidArgumentError(
                    "invalid value for γ: $(args[4])  must be a positive number"))
            end
        end
    end

    return infile, outfile, a, γ
end

##########################################################################################92

"""
    load_image(path::Union{String, IO}) :: HDRimage

Load an image from the specified `path` in an HDR image format.

See also: [`HDRimage`](@ref)
"""
function load_image(path::Union{String, IO})
    img = load(path)
    img_float32 = float32.(img)
    width, height = size(img_float32)
    vec = reshape(permutedims(img_float32), (width*height,))
    #vec = reshape(reverse(permutedims(img_float32), dims=1), (width*height,))

    return HDRimage(height, width, vec)
end

"""
    ldr2pfm(path::String, outfile::String)

Load an image from the specified `path`, convert it in a pfm file format
and save it as `outfile`.
It works through the `load_image` and the  `write(::IO, ::HDRimage)` function.

See also: [`HDRimage`](@ref), [`load_image`](@ref),
[`write(::IO, ::HDRimage)`](@ref)
"""
function ldr2pfm(path::String, outfile::String)
    img = load_image(path)
    open(outfile, "w") do outf;
        write(outf, img)
    end
    nothing
end


text:
struct LagrangeTensorProductBasis{dim,NF,NF1D,T} <: AbstractBasis{dim,NF}
    basis::LagrangeBasis{NF1D,T}
    points::Matrix{T}
    function LagrangeTensorProductBasis(dim,order;start = -1.0,stop = 1.0)
        @assert 1 <= dim <= 3
        basis1d = LagrangeBasis(order, start = start, stop = stop)
        NF1D = number_of_basis_functions(basis1d)
        NF = NF1D^dim
        T = type_of_interpolation_points(basis1d)
        points = tensor_product_points(dim,basis1d.points')
        new{dim,NF,NF1D,T}(basis1d,points)
    end
end

function Base.show(io::IO, basis::LagrangeTensorProductBasis{dim,NF,NF1D,T}) where {dim,NF,NF1D,T}
    p = order(basis)
    str = "LagrangeTensorProductBasis\n\tDimension  : $dim\n\tOrder      : $p\n\tNum. Funcs.: $NF"
    print(io,str)
end

function order(basis::LagrangeTensorProductBasis)
    return order(basis.basis)
end

function (B::LagrangeTensorProductBasis{1})(x)
    return B.basis(x)
end

function (B::LagrangeTensorProductBasis{1})(x::V) where {V<:AbstractVector}
    @assert length(x) == 1
    return B(x[1])
end

function (B::LagrangeTensorProductBasis{2})(x,y)
    return kron(B.basis(x),B.basis(y))
end

function (B::LagrangeTensorProductBasis{2})(x::V) where {V<:AbstractVector}
    @assert length(x) == 2
    return B(x[1],x[2])
end

function (B::LagrangeTensorProductBasis{3})(x,y,z)
    return kron(B.basis(x),B.basis(y),B.basis(z))
end

function (B::LagrangeTensorProductBasis{3})(x::V) where {V<:AbstractVector}
    @assert length(x) == 3
    return B(x[1],x[2],x[3])
end

function gradient(B::LagrangeTensorProductBasis{1},x)
    return derivative(B.basis,x)
end

function gradient(B::LagrangeTensorProductBasis{1},x::V) where {V<:AbstractVector}
    @assert length(x) == 1
    return gradient(B,x[1])
end

function gradient(B::LagrangeTensorProductBasis{2},dir::Z,x,y) where {Z<:Integer}
    @assert 1 <= dir <= 2
    if dir == 1
        dNx = derivative(B.basis,x)
        Ny = B.basis(y)
        return kron(dNx,Ny)
    else
        Nx = B.basis(x)
        dNy = derivative(B.basis,y)
        return kron(Nx,dNy)
    end
end

function gradient(B::LagrangeTensorProductBasis{2},x,y)
    col1 = gradient(B,1,x,y)
    col2 = gradient(B,2,x,y)
    return hcat(col1,col2)
end

function gradient(B::LagrangeTensorProductBasis{2},dir,x::V) where {V<:AbstractVector}
    @assert length(x) == 2
    @assert 1 <= dir <= 2
    return gradient(B,dir,x[1],x[2])
end

function gradient(B::LagrangeTensorProductBasis{2},x::V) where {V<:AbstractVector}
    @assert length(x) == 2
    return gradient(B,x[1],x[2])
end

function gradient(B::LagrangeTensorProductBasis{3},dir,x,y,z)
    @assert 1 <= dir <= 3
    if dir == 1
        dNx = derivative(B.basis,x)
        Ny = B.basis(y)
        Nz = B.basis(z)
        return kron(dNx,Ny,Nz)
    elseif dir == 2
        Nx = B.basis(x)
        dNy = derivative(B.basis,y)
        Nz = B.basis(z)
        return kron(Nx,dNy,Nz)
    else
        Nx = B.basis(x)
        Ny = B.basis(y)
        dNz = derivative(B.basis,z)
        return kron(Nx,Ny,dNz)
    end
end

function gradient(B::LagrangeTensorProductBasis{3},x,y,z)
    Nx = B.basis(x)
    Ny = B.basis(y)
    Nz = B.basis(z)

    dNx = derivative(B.basis,x)
    dNy = derivative(B.basis,y)
    dNz = derivative(B.basis,z)

    return hcat(kron(dNx,Ny,Nz),kron(Nx,dNy,Nz),kron(Nx,Ny,dNz))
end

function gradient(B::LagrangeTensorProductBasis{3},dir,x::V) where {V<:AbstractVector}
    @assert length(x) == 3
    @assert 1 <= dir <= 3
    return gradient(B,dir,x[1],x[2],x[3])
end

function gradient(B::LagrangeTensorProductBasis{3},x::V) where {V<:AbstractVector}
    @assert length(x) == 3
    return gradient(B,x[1],x[2],x[3])
end

function hessian(B::LagrangeTensorProductBasis{1},x)
    return hessian(B.basis,x)
end

function hessian(B::LagrangeTensorProductBasis{1},x::V) where {V<:AbstractVector}
    @assert length(x) == 1
    return hessian(B,x[1])
end

function hessian(B::LagrangeTensorProductBasis{2},x,y)

    d2Nx = hessian(B.basis,x)
    d1Nx = derivative(B.basis,x)
    d0Nx = B.basis(x)

    d2Ny = hessian(B.basis,y)
    d1Ny = derivative(B.basis,y)
    d0Ny = B.basis(y)

    Nxx = kron(d2Nx,d0Ny)
    Nxy = kron(d1Nx,d1Ny)
    Nyy = kron(d0Nx,d2Ny)

    return hcat(Nxx,Nxy,Nyy)
end

function hessian(B::LagrangeTensorProductBasis{2},x::V) where {V<:AbstractVector}
    @assert length(x) == 2
    return hessian(B,x[1],x[2])
end

function number_of_1d_points(p::M) where {M<:AbstractMatrix}
    m,n = size(p)
    @assert m == 1
    return n
end

function tensor_product_points(N,points)
    @assert 1 <= N <= 3

    np = number_of_1d_points(points)

    if N == 1
        return points
    elseif N == 2
        row1 = repeat(points,inner=(1,np))
        row2 = repeat(points,outer=(1,np))
        return vcat(row1,row2)
    else
        np2 = np^2
        row1 = repeat(points,inner=(1,np2))
        row2 = repeat(repeat(points,inner=(1,np)),outer=(1,np))
        row3 = repeat(points,outer=(1,np2))
        return vcat(row1,row2,row3)
    end
end


text:
"""
	sphdistance(cmd0::String="", arg1=[], kwargs...)

Create Voronoi distance, node, or natural nearest-neighbor grid on a sphere

Full option list at [`sphdistance`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance .html)

Parameters
----------

- **C** : **save_mem** : -- Bool or [] --

    For large data sets you can save some memory (at the expense of more processing).
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#c)
- **E** : **what_quantity** : -- Str --   Flags = d|n|z[dist]

    Specify the quantity that should be assigned to the grid nodes.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#e)
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = sphdistance(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#g)
- **I** : **inc** : -- Str or Number --

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#i)
- **L** : **dist_unit** : -- Str --      Flags = d|e|f|k|M|n|u

    Specify the unit used for distance calculations.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#l)
- **N** : **nodetable** : -- Str --      Flags = nodetable

    Read the information pertaining to each Voronoi polygon (the unique node lon, lat and polygon area)
    from a separate file.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#n)
- **Q** : **voronoi_polyg** : -- Str --     Flags = voronoifile

    Append the name of a file with pre-calculated Voronoi polygons.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance .html#q)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function sphdistance(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0  && return monolitic("sphdistance ", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'C', d, [:C :save_mem])
	cmd = add_opt(cmd, 'E', d, [:E :what_quantity])
    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'L', d, [:L :dist_unit])
	cmd = add_opt(cmd, 'N', d, [:N :nodetable])
	cmd = add_opt(cmd, 'Q', d, [:Q :voronoi_polyg])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "sphdistance", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
sphdistance(arg1=[], cmd0::String=""; kw...) = sphdistance(cmd0, arg1; kw...)


text:
module BenchHeap

using DataStructures
using BenchmarkTools
using Random

function push_heap(h::AbstractHeap, xs::Vector)
    n = length(xs)

    for i = 1 : n
        push!(h, xs[i])
    end
end

function pop_heap(h::AbstractHeap)
    n = length(h)

    for i = 1 : n
        pop!(h)
    end
end

suite = BenchmarkGroup()

heaptypes = [BinaryHeap, MutableBinaryHeap]
aexps = [1,3]
datatypes = [Int, Float64]
baseorderings = Dict(
    "Min" => DataStructures.LessThan,
    #"Max" => DataStructures.GreaterThan,
    )
fastfloatorderings = Dict(
    # These will be enabled upon reordering change
    #"FastMin" => DataStructures.FasterForward(),
    #"FastMax" => DataStructures.FasterReverse(),
    )

for heap in heaptypes
    for aexp in aexps
        for dt in datatypes
            Random.seed!(0)
            a = rand(dt, 10^aexp)

            orderings = baseorderings
            if dt == Float64
                # swap to faster ordering operation
                for (k,v) in orderings
                    if haskey(fastfloatorderings, k)
                        orderings["Slow"*k] = v
                        orderings[k] = fastfloatorderings[k]
                    end
                end
            end

            for (ord_str, ord) in orderings
                prepath = [string(heap)]
                postpath = [string(dt), "10^"*string(aexp), ord_str]
                suite[vcat(prepath, ["make"], postpath)] =
                    @benchmarkable $(heap){$dt,$ord}($a)
                suite[vcat(prepath, ["push"], postpath)] =
                    @benchmarkable push_heap(h, $a) setup=(h=$(heap){$dt,$ord}())
                suite[vcat(prepath, ["pop"], postpath)] =
                    @benchmarkable pop_heap(h) setup=(h=$(heap){$dt,$ord}($a))
            end
        end
    end
end

# Quick check to ensure no Float regressions with Min/Max convenience functions
# These don't fit in well with the above loop, since ordering is hardcoded.
heapalias = Dict(
    "BinaryMinHeap" => BinaryMinHeap,
    "BinaryMaxHeap" => BinaryMaxHeap,
    "BinaryMinMaxHeap" => BinaryMinMaxHeap, # <- no alias issue
)
for (heapname, heap) in heapalias
    for aexp in aexps
        for dt in [Float64]
            Random.seed!(0)
            a = rand(dt, 10^aexp)
            prepath = [heapname]
            postpath = [string(dt), "10^"*string(aexp)]
            suite[vcat(prepath, ["make"], postpath)] =
                @benchmarkable $(heap)($a)
            suite[vcat(prepath, ["push"], postpath)] =
                @benchmarkable push_heap(h, $a) setup=(h=$(heap){$dt}())
            suite[vcat(prepath, ["pop"], postpath)] =
                @benchmarkable pop_heap(h) setup=(h=$(heap)($a))
        end
    end
end

for func in [nlargest, nsmallest]
    for aexp in [4]
        Random.seed!(0);
        a = rand(10^aexp);
        for nexp in [2]
            n = 10^nexp
            suite[[string(func), "a=rand(10^"*string(aexp)*")", "n=10^"*string(nexp)]] =
                @benchmarkable $(func)($n, $a)
        end
    end
end

end  # module

BenchHeap.suite


text:
using Test
using QCMaterial

@testset "util.topylist" begin
    org_array = [1, 2, 3.0]
    pylist = topylist(org_array)
    @test all(org_array .== pylist)
end

@testset "util.doublefunc" begin
    org_array = [1, 2, 3]
    @test all(doublefunc(org_array) .== 2 .* org_array)
end

@testset "util.numerical_grad" begin
    # Some test for numerical_grad
    f(x) = x[1] + 2*x[2]
    deriv = numerical_grad(f, zeros(2))
    @assert deriv ≈ [1.0, 2.0]
end


text:
module MapsMocksTests

using Test
using Gridap
using TensorValues

using ..MapsMocks

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
ao = [  a+pj for pj in p  ]
m = MockMap(a)
test_map(m,p,ao)

ao = [p[i]*j+a  for j in 1:3, i in 1:length(p)]
m = TestMap(a,3)
test_map(m,p,ao)

a = VectorValue(10.1,10.3)
b = VectorValue(15.3,20.5)
ao = [  a+pj for pj in p  ]
m = MockMap(a)
test_map(m,p,ao)

ao = [p[i]*j+a  for j in 1:3, i in 1:length(p)]
m = TestMap(a,3)
test_map(m,p,ao)

end # module


text:
module TGST
# routine tests for the Symmetric Tridiagonal eigensolvers in GenericSchur

using LinearAlgebra
using Test
using GenericSchur
using GenericSchur: geigen!, geigvals!

# some interesting examples from the LAPACK papers

function wilkinson(n,::Type{T}) where T
    if mod(n,2) == 1
        nh = (n-1) >>> 1
        D = T.(vcat(nh:-1:1,[0],1:nh))
    else
        D = abs.(T.(collect(-(n-1):2:(n-1)))/2)
    end
    E = fill(T(1),n-1)
    D,E
end

function clement(n,::Type{T}) where T
    D = zeros(T,n)
    E = [sqrt(T(i*(n+1-i))) for i in 1:n-1]
    D,E
end

function legendre(n,::Type{T}) where T
    D = zeros(T,n)
    E = [T(i)/sqrt(T((2i-1)*(2i)+1)) for i in 2:n]
    D,E
end

function laguerre(n,::Type{T}) where T
    D = T.(collect(3:2:2*n+1))
    E = T.(collect(2:n))
    D,E
end

function hermite(n,::Type{T}) where T
    D = zeros(T,n)
    E = [sqrt(T(i)) for i in 1:n-1]
    D,E
end

# LAPACK-style classes
function st_testmat1(n,itype,::Type{T}=Float64; dtype=:normal) where T
    dtype ∈ [:normal, :unif, :uniform] || throw(ArgumentError("incomprehensible dtype"))
    λ = zeros(T,n)
    ulp = eps(real(T))
    kpar = 1/ulp
    if itype == 1
        λ[1] = 1
        for j=2:n
            λ[j] = 1 / T(j)
        end
    elseif itype == 2
        λ[1:n-1] .= 1
        λ[n] = T(1) / n
    elseif itype == 3
        for j=1:n
            λ[j] = 1 / (kpar^((j-1)/(n-1)))
        end
    elseif itype == 4
        for j=1:n
            λ[j] = T(1) - (T(j-1)/T(n-1)) * (1-1/kpar)
        end
    elseif itype == 5
        λ .= exp.(-log(n) * rand(n))
    elseif itype == 6
        if dtype == :normal
            λ .= randn(n)
        else
            λ .= rand(dtype, n)
        end
    elseif itype == 7
        λ[1:n-1] .= (1:n-1)*ulp
        λ[n] = 1
    elseif itype == 8
        λ[1] = ulp
        for j = 2:n-1
            λ[j] = 1 + sqrt(ulp) * j
        end
        λ[n] = 2
    elseif itype == 9
        λ[1] = 1
        for j=2:n
            λ[j] = λ[j-1] + 100 * ulp
        end
    end
    Q,r = qr(randn(n,n))
    Q = T.(Matrix(Q)) # until LinearAlgebra is fixed
    # FIXME: sign adjustment
    A = Q' * diagm(λ) * Q

    F = hessenberg!(A)
    S = F.H
    return S, λ
end
const itypes1 = 1:9
# downselect so subsetting is reasonable
const itypes1a = [1,3,4,5,6]

function st_testmat2(n,itype,::Type{T}=Float64; dtype=:normal) where T <: Real
    dtype ∈ [:normal, :unif, :uniform] || throw(ArgumentError("incomprehensible dtype"))
    λ = zeros(T,n)
    ulp = eps(real(T))
    if itype == 0
        S = SymTridiagonal(zeros(T,n),zeros(T,n-1))
        λ = zeros(T,n)
    elseif itype == 1
        S = SymTridiagonal(ones(T,n),zeros(T,n-1))
        λ = ones(T,n)
    elseif itype == 2
        S = SymTridiagonal(fill(T(2),n),ones(T,n-1))
        λ = nothing
    elseif itype == 3
        S = SymTridiagonal(wilkinson(n,T)...)
        λ = nothing
    elseif itype == 4
        S = SymTridiagonal(clement(n,T)...)
        λ = T.(collect(-n:2:n))
    elseif itype == 5
        S = SymTridiagonal(legendre(n,T)...)
        λ = nothing
    elseif itype == 6
        S = SymTridiagonal(laguerre(n,T)...)
        λ = nothing
    elseif itype == 7
        S = SymTridiagonal(hermite(n,T)...)
        λ = nothing
    elseif itype == 8
        # exercise the scaling logic, all at once
        small = floatmin(T) / ulp
        large = floatmax(T) * ulp
        S1 = SymTridiagonal(wilkinson(n,T)...)
        S = gluemats(small * S1, S1, large * S1)
    else
        throw(ArgumentError("unimplemented itype"))
    end
    S,λ
end
const itypes2 = 1:8
const broken_cases2 = [(BigFloat, 8)]

function gluemats(mat1::SymTridiagonal{T}, mats...; gval=zero(T)) where T
    D = diag(mat1)
    E = diag(mat1,1)
    for mat in mats
        append!(D, diag(mat))
        push!(E, gval)
        append!(E, diag(mat,1))
    end
    SymTridiagonal(D, E)
end

function batch(n, ::Type{T}=Float64; thresh=50, quiet=true) where {T}
    dmax=0.0; vmax=0.0
    for itype in itypes1
        @testset "class 1 type $itype" begin
            A,_ = st_testmat1(n, itype, T)
            de, ve = runtest(diag(A), diag(A,1))
            dmax = max(dmax, de)
            vmax = max(vmax, ve)
            @test de < thresh
            @test ve < thresh
        end
    end
    quiet || println("peak errors: $dmax $vmax")
    dmax=0.0; vmax=0.0
    for itype in itypes2
        @testset "class 2 type $itype" begin
            A,_ = st_testmat2(n, itype, T)
            de, ve = runtest(diag(A), diag(A,1))
            dmax = max(dmax, de)
            vmax = max(vmax, ve)
            if (T, itype) in broken_cases2
                @test_broken de < thresh
            else
                @test de < thresh
            end
            @test ve < thresh
        end
    end
    quiet || println("peak errors: $dmax $vmax")
    nothing
end

maxnorm(A) = maximum(abs.(vec(A)))

function runtest(D::AbstractVector{T}, E; λ=nothing, V=nothing) where T <: Real
    if λ === nothing
        λ, V = geigen!(SymTridiagonal(copy(D), copy(E)))
    end
    n = length(D)
    m = size(V,2)
    ntol = min(m,n)
    myeps = eps(T)
    if m == n
        resnorm = maxnorm(V*Diagonal(λ)*V' - SymTridiagonal(D,E))
    else
        resnorm = maxnorm(Diagonal(λ) - V' * SymTridiagonal(D,E) * V)
    end
    Tnorm = max(opnorm(SymTridiagonal(D,E),1), floatmin(T))
    if resnorm < Tnorm
        d_err =  (resnorm / Tnorm) / (myeps * ntol)
    elseif Tnorm < 1
        d_err = (min(resnorm, Tnorm * ntol) / Tnorm) / (myeps * ntol)
    else
        d_err = min(resnorm / Tnorm, ntol) / (myeps * ntol)
    end
    if m == n
        v_err = maxnorm(V * V' - I) / (myeps * ntol)
    else
        v_err = maxnorm(V' * V - I) / (myeps * ntol)
    end
    d_err, v_err
end

for T in [Float32, Float64, BigFloat]
    @testset "SymTridiagonal $T" begin
        batch(32, T)
    end
end

end # module


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

using NodeNumbering

@testset "Renumbering the RCM ordered nodes." begin

    neworder = Int[6, 7, 2, 5, 3, 8, 1, 4, 9]

    expected = Dict(
               7 => 2,
               9 => 9,
               4 => 8,
               2 => 3,
               3 => 5,
               5 => 4,
               8 => 6,
               6 => 1,
               1 => 7)
    result = NodeNumbering.renumbering(neworder)
    @test result == expected
end