text:
#=

This file contains helper functions that are used throughout the module

=#

""" Convenience function to perform (-1)^x using x's parity """
exp_m1(x) = iseven(x) ? 1 : -1

"""
    $(FUNCTIONNAME)(size_parameter)
Computes the number of required Legendre functions  for a given size parameter.
See eq 6 in Sanghavi 2014
- `size_parameter` size parameter of the aerosol (2πr/λ)
The function returns a rounded integer, following conventions by BH, Rooj/Stap, Siewert
"""
get_n_max(size_parameter) = round(Int, size_parameter + 4.05 * size_parameter^(1/3) + 10)

"""
    $(FUNCTIONNAME)(size_parameter,refractive_idx::Number,an,bn,Dn)
Computes Mie coefficients `an` and `bn` as a function of size parameter and complex
refractive index. See eq 4.88 in Bohren and Huffman
- `size_parameter` size parameter of the aerosol (2πr/λ)
- `refractive_idx` refractive index of the aerosol (complex number)
- `an` and `bn` pre-allocated arrays, need to match at least n_max for the given size parameter
- `Dn` pre-allocated array for the logarithmic derivative (see BH, eq 4.88)
(need to check whether Dn can be created internally without causing too many allocations)

The function returns a rounded integer, following conventions by BH, Rooj/Stap, Siewert
"""
function compute_mie_ab!(size_param, refractive_idx::Number, an, bn, Dn)
    FT = eltype(refractive_idx)
    FT2 = eltype(an)

    # Compute y
    y = size_param * refractive_idx

    # Maximum expansion (see eq. A17 from de Rooij and Stap, 1984)
    n_max = get_n_max(size_param)

    # Make sure downward recurrence starts higher up
    # (at least 15, check eq. A9 in de Rooij and Stap, 1984, may need to check what is needed)
    nmx = length(Dn)
    @assert size(an)[1] >= n_max
    @assert size(an) == size(bn)
    fill!(Dn, 0);

    # Dn as in eq 4.88, Bohren and Huffman, to calculate an and bn
    # Downward Recursion, eq. 4.89, Bohren and Huffman
    [Dn[n] = ((n+1) / y) - (1 / (Dn[n+1] + (n+1) / y)) for n = (nmx - 1):-1:1]

    # Get recursion for bessel functions ψ and ξ
    ψ₀, ψ₁, χ₀, χ₁ =  (cos(size_param), sin(size_param), -sin(size_param), cos(size_param))
    ξ₁ = FT(ψ₁, -χ₁)

    # This solves Bohren and Huffman eq. 4.88 for an and bn, computing updated ψ and ξ on the fly
    for n = 1:n_max
        fn = (2n + 1) / (n * (n + 1))
        ψ  = (2n - 1) * ψ₁ / size_param - ψ₀
        χ  = (2n - 1) * χ₁ / size_param - χ₀
        ξ   = FT(ψ, -χ)
        t_a = Dn[n] / refractive_idx + n / size_param
        t_b = Dn[n] * refractive_idx + n / size_param

        an[n] = (t_a * ψ - ψ₁) / (t_a * ξ - ξ₁)
        bn[n] = (t_b * ψ - ψ₁) / (t_b * ξ - ξ₁)

        ψ₀ = ψ₁
        ψ₁ = ψ
        χ₀ = χ₁
        χ₁ = χ
        ξ₁ = FT(ψ₁, -χ₁)
    end
end

"""
    $(FUNCTIONNAME)(model::MieModel, λ, radius)
Compute all an, bn using compute_mie_ab!
Input: MieModel, wavelength (λ), radius
Output: an, bn. Both of shape (aerosol.nquad_radius, N_max) (N_max from aerosol.r_max)
"""
function compute_anbn(model::MieModel, λ, radius)

    @unpack computation_type, aerosol, r_max, nquad_radius, λ, polarization_type, truncation_type, wigner_A, wigner_B = model
    @unpack size_distribution, nᵣ, nᵢ = aerosol

    FT = eltype(λ)
    FT2 = eltype(nᵣ)

    # Find overall N_max from the maximum radius
    N_max = Scattering.get_n_max(2 * π * r_max / λ)

    # Where to store an, bn, computed over size distribution
    an = zeros(Complex{FT2}, nquad_radius, N_max)
    bn = zeros(Complex{FT2}, nquad_radius, N_max)

    # Loop over the size distribution, and compute an, bn, for each size
    for i in 1:nquad_radius

        # Get current radius and size parameter
        r = radius[i]
        size_param = 2 * π * r / λ

        # Pre-allocate Dn:
        y = size_param * (nᵣ - nᵢ);
        nmx = round(Int, max(N_max, abs(y)) + 51)
        Dn = zeros(Complex{FT2}, nmx)

        # Compute an, bn
        Scattering.compute_mie_ab!(size_param, nᵣ + nᵢ * im,
                                      view(an, i, :),
                                      view(bn, i, :), Dn)
    end

    return an, bn;
end

"""
    $(FUNCTIONNAME)(an, bn, ab_pairs, w, Nmax, N_max_)
From the an, bn matrices, precompute all (an✶)am, (an✶)bm, (bn✶)am, (bn✶)bm
This allows quick computation of (an✶ + bn✶) × (am + bm)
"""
function compute_avg_anbns!(an, bn, ab_pairs, w, Nmax, N_max_)
    FT2 = eltype(an)

    # Unpack ab_pairs
    mat_anam, mat_anbm, mat_bnam, mat_bnbm = ab_pairs

    # Fill all matrices with 0
    [fill!(mat, 0) for mat in [mat_anam, mat_bnbm, mat_anbm, mat_bnam]]

    @inbounds for n in 1:Nmax, m in n:Nmax
        anam = bnbm = anbm = bnam = FT2(0);
        @inbounds for i = 1:size(an, 1)
            if m < N_max_[i] && n < N_max_[i]
                anam += w[i] * an[i,n]' * an[i,m]
                bnbm += w[i] * bn[i,n]' * bn[i,m]
                anbm += w[i] * an[i,n]' * bn[i,m]
                bnam += w[i] * bn[i,n]' * an[i,m]
            end
        end
        @inbounds mat_anam[m,n] = anam;
        @inbounds mat_bnbm[m,n] = bnbm;
        @inbounds mat_anbm[m,n] = anbm;
        @inbounds mat_bnam[m,n] = bnam;
    end
end

"""
    $(FUNCTIONNAME)(an, bn, π_, τ_, S₁, S₂)
Determines the amplitude functions `S₁`,`S₂` in Mie theory
- `an` and `bn` pre-calculated Mie coefficients `an` and `bn`, see [`compute_mie_ab!`](@ref) function
- `π` and `τ` pre-calculated associated Legendre functions `π` and `τ`, see [`compute_mie_π_τ!`](@ref) function
The function returns `S₁`,`S₂` as a function of the cosine of the scattering angle `ξ`.
Users need to make sure `an` and `bn`, `π` and `τ` are pre-computed.
"""
function compute_mie_S₁S₂!(an, bn, π_, τ_, S₁, S₂)
    FT = eltype(an)
    nmax = size(an, 1);
    nμ   = size(π_, 1);

    # Verify sizes
    @assert size(S₁) == size(S₂)
    @assert length(S₁) == nμ

    for l in 1:nmax, iμ in 1:nμ
            S₁[iμ] += (2l + 1) / (l * (l + 1)) * (an[l] * τ_[iμ,l] + bn[l] * π_[iμ,l])
            S₂[iμ] += (2l + 1) / (l * (l + 1)) * (an[l] * π_[iμ,l] + bn[l] * τ_[iμ,l])
    end
end

"""
    $(FUNCTIONNAME)(n,xmin,xmax; norm=false)
Returns the `n` Gauss-Legendre quadrature points and weights with a change of interval between xmin and xmax
- `n` number of quadrature points
- `xmin`,`xmax` lower and upper bound of integral
- `norm`: if `true`, normalizes the weights so that a mean can be computed instead of full integration
The function returns `n` quadrature points ξ within [xmin,xmax] with associated weightes `w`
"""
function gauleg(n, xmin, xmax; norm=false)
    ξ, w = gausslegendre(n)
    ξ = (xmax - xmin) / 2 * ξ .+ (xmin + xmax) / 2
    norm ? w /= sum(w) : w *= (xmax - xmin) / 2
    return ξ, w
end

@doc raw"""
    $(FUNCTIONNAME)(greek_coefs, μ; returnLeg = false)
Returns the reconstructed elements of the 4x4 scattering matrix at positions
f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄ from the greek coefficients

f₁₁ represents the phase function p for the Intensity (first Stokes Vector element) and is normalized as follows:
```math
\frac{1}{4\pi}\int_0^{2\pi}d\phi \int_{-1}^1 p(\mu) d\mu  = 1
```

- `greek_coefs` greek coefficients (Domke Type)
- `returnLeg` if `false` (default), just return `f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄`, if `true`,
- return `f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄, P, P²` (i.e. also the two legendre polynomials as matrices)
"""
function reconstruct_phase(greek_coefs, μ; returnLeg=false)

    FT = eltype(greek_coefs.α)
    l_max = length(greek_coefs.α);
    nμ = length(μ)

    # Compute legendre polynomials
    P, P², R², T² = compute_legendre_poly(μ, l_max)

    # To stay general, we also don't assume f₂₂=f₁₁ or f₄₄=f₃₃
    # which only holds for spherical
    f₁₁, f₃₃, f₁₂, f₃₄, f₂₂, f₄₄ = (zeros(FT, nμ), zeros(FT, nμ), zeros(FT, nμ),
                                    zeros(FT, nμ), zeros(FT, nμ), zeros(FT, nμ))

    # Compute prefactor
    fac = zeros(l_max);
    [fac[l + 1] = sqrt(1 / ((l-1) * l * (l+1) * (l+2))) for l = 2:(l_max-1)]

    # In matrix form:
    f₁₁[:] = P * greek_coefs.β                                           # a₁ in Rooij notation
    f₄₄[:] = P * greek_coefs.δ                                           # a₄ in Rooij notation
    f₁₂[:] = P² * (fac .* greek_coefs.γ)                                 # b₁ in Rooij notation
    f₃₄[:] = P² * (fac .* greek_coefs.ϵ)                                 # b₂ in Rooij notation
    f₂₂[:] = R² * (fac .* greek_coefs.α) .+ T² * (fac .* greek_coefs.ζ)  # a₂ in Rooij notation
    f₃₃[:] = R² * (fac .* greek_coefs.ζ) .+ T² * (fac .* greek_coefs.α)  # a₃ in Rooij notation

    # Put elements into a struct
    scattering_matrix = ScatteringMatrix(f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄)

    # For truncation in δ-BGE, we need P and P² as well, convenient to return here:
    return returnLeg ? (scattering_matrix, P, P²) : scattering_matrix
end

"""
    $(FUNCTIONNAME)(depol)
Returns the greek coefficients (as [`GreekCoefs`](@ref)) of the Rayleigh phase function given
depolarization value.
- `depol` Depolarization (best use 0 as default )
"""
function get_greek_rayleigh(depol::Number)
    FT = eltype(depol)
    # Rayleigh Greek Parameters
    dpl_p = (1 - depol)  / (1 + depol / 2)
    dpl_q = (1 + depol)  / (1 - depol)
    dpl_r = (1 - 2depol) / (1 - depol)

    α  =  FT[0.0, 0.0,             3dpl_p]
    β  =  FT[1.0, 0.0,             0.5 * dpl_p]
    γ  =  FT[0.0, 0.0,             dpl_p * sqrt(1.5)]
    δ  =  FT[0.0, dpl_p * dpl_r * 1.5, 0.0]
    ϵ  =  FT[0.0, 0.0,             0.0]
    ζ  =  FT[0.0, 0.0,             0.0]
    return GreekCoefs(α, β, γ, δ, ϵ, ζ)
end

"""
    $(FUNCTIONNAME)(k, an, bn, w)
Calculate the average Scattering and Extinction Cross Section
Eqn. 1, averaged over size distribution
"""
function compute_avg_C_scatt_ext(k, an, bn, w)
    n_ = collect(1:size(an)[2]);
    n_ = 2n_ .+ 1
    coef = 2π / k^2 * n_'
    return (coef * (w' * (abs2.(an') + abs2.(bn'))')', coef * (w' * real(an + bn))')
end

""" Compute probability weights of radii """
function compute_wₓ(size_distribution, wᵣ, r, r_max)

    wₓ = pdf.(size_distribution,r)      # Weights from distribution
    wₓ .*= wᵣ                           # pre multiply with wᵣ to get proper means eventually:

    # normalize (could apply a check whether cdf.(size_distribution,r_max) is larger than 0.99:
    #println("Test")
    @info "Fraction of size distribution cut by max radius: $((1-cdf.(size_distribution,r_max))*100) %"
    wₓ /= sum(wₓ)
    return wₓ
end

#####
##### Π-matrix construction methods (Sanghavi 2014, eq. 15)
#####

"""
    $(FUNCTIONNAME)(mo::Stokes_IQUV, P, R, T, l::Int, m::Int; sign_change=false)
Compute Π matrix for all stokes vector elements used in computations of the phase matrix
See Sanghavi 2014, eq. 15
"""
function construct_Π_matrix(mo::Stokes_IQUV, P, R, T, l::Int, m::Int; sign_change=false)
    if sign_change # (basically gets it for -μ due to symmetries on P,R,T)
        if isodd(l-m)
            Π = [SMatrix{4,4}([-P[i,l,m] 0 0 0 ; 0 -R[i,l,m] -T[i,l,m] 0; 0 -T[i,l,m] -R[i,l,m] 0; 0 0 0 -P[i,l,m]]) for i in 1:size(P, 1)]
        else
            Π = [SMatrix{4,4}([P[i,l,m] 0 0 0 ; 0 R[i,l,m] T[i,l,m] 0; 0 T[i,l,m] R[i,l,m] 0; 0 0 0 P[i,l,m]]) for i in 1:size(P, 1)]
        end
    else
        Π = [SMatrix{4,4}([P[i,l,m] 0 0 0 ; 0 R[i,l,m] -T[i,l,m] 0; 0 -T[i,l,m] R[i,l,m] 0; 0 0 0 P[i,l,m]]) for i in 1:size(P, 1)]
    end
    return Π
end

"""
    $(FUNCTIONNAME)(mo::Stokes_IQU, P, R, T, l::Int, m::Int; sign_change=false)
Compute Π matrix for  stokes vector elements I,Q,U used in computations of the phase matrix
See Sanghavi 2014, eq. 15
"""
function construct_Π_matrix(mo::Stokes_IQU, P, R, T, l::Int, m::Int; sign_change=false)
    if sign_change # (basically gets it for -μ due to symmetries on P,R,T)
        if isodd(l-m)
            Π = [SMatrix{3,3}([-P[i,l,m] 0 0  ; 0 -R[i,l,m] -T[i,l,m] ; 0 -T[i,l,m] -R[i,l,m] ]) for i in 1:size(P, 1)]
        else
            Π = [SMatrix{3,3}([P[i,l,m] 0 0  ; 0 R[i,l,m] T[i,l,m] ; 0 T[i,l,m] R[i,l,m] ]) for i in 1:size(P, 1)]
        end
    else
        Π = [SMatrix{3,3}([P[i,l,m] 0 0  ; 0 R[i,l,m] -T[i,l,m] ; 0 -T[i,l,m] R[i,l,m] ]) for i in 1:size(P, 1)]
    end
    return Π
end

"""
    $(FUNCTIONNAME)(mo::Stokes_I, P, R, T, l::Int, m::Int; sign_change=false)
Compute Π matrix for  stokes vector elements I used in computations of the phase matrix

"""
construct_Π_matrix(mod::Stokes_I, P, R, T, l::Int, m::Int; sign_change=false) = sign_change ? -P[:,l,m] : P[:,l,m]

#####
##### B-matrix construction methods (Sanghavi 2014, eq. 16)
#####

"""
    $(FUNCTIONNAME)(mod::Stokes_IQUV, α, β, γ, δ, ϵ, ζ, l::Int)
Compute B matrix for all stokes vector elements used in computations of the phase matrix
See Sanghavi 2014, eq. 16
"""
construct_B_matrix(mod::Stokes_IQUV, α, β, γ, δ, ϵ, ζ, l::Int) = SMatrix{4,4}([β[l] γ[l] 0 0 ; γ[l] α[l] 0 0; 0 0 ζ[l] ϵ[l]; 0 0 -ϵ[l] δ[l]])

"""
    $(FUNCTIONNAME)(mod::Stokes_IQU, α, β, γ, δ, ϵ, ζ, l::Int)
Compute B matrix for stokes vector elements I,Q,U used in computations of the phase matrix
    See Sanghavi 2014, eq. 16
"""
construct_B_matrix(mod::Stokes_IQU, α, β, γ, δ, ϵ, ζ, l::Int) = SMatrix{3,3}([β[l] γ[l] 0 ; γ[l] α[l] 0 ; 0 0 ζ[l]])

"""
$(FUNCTIONNAME)(mod::Stokes_I, α, β, γ, δ, ϵ, ζ, l::Int)
Compute Π matrix for stokes vector elements I used in computations of the phase matrix
See Sanghavi 2014, eq. 16
"""
construct_B_matrix(mod::Stokes_I, α, β, γ, δ, ϵ, ζ, l::Int) = β[l]

"""
    $(FUNCTIONNAME)(mod::AbstractPolarizationType, μ, α, β, γ, δ, ϵ, ζ, m::Int)
Compute moments of the phase matrix
"""
function compute_Z_moments(mod::AbstractPolarizationType, μ, greek_coefs::GreekCoefs, m::Int ; arr_type = Array)
    @unpack α, β, γ, δ, ϵ, ζ = greek_coefs
    FT = eltype(β)
    n = length(μ)

    # Change from 0-index to 1-index (i.e. the lowest m is 0 ),
    # make more logical later to avoid confusion later (m=0 has a meaning!!)
    m = m+1

    # Set prefactor for moments (note 1-notation for `m` here):
    fact = (m == 1) ? 0.5 : 1.0

    # get l_max just from length of array:
    l_max = length(β)

    # Check that all μ are positive here ([0,1])
    # @show μ
    @assert all(0 .< μ .≤ 1) "all μ's within compute_Z_moments have to be ∈ ]0,1]"

    # Compute legendre Polynomials at μ and up to lmax
    P, R, T    = Scattering.compute_associated_legendre_PRT(μ, l_max)
    P⁻, R⁻, T⁻ = Scattering.compute_associated_legendre_PRT(-μ, l_max)

    # Pre-compute all required B matrices
    𝐁_all = [construct_B_matrix(mod, α, β, γ, δ, ϵ, ζ, i) for i in 1:l_max]

    # Get dimension of square matrix (easier for Scalar/Stokes dimensions)
    B_dim = Int(sqrt(length(𝐁_all[1])))

    # Create matrices:
    nb = B_dim * n
    𝐙⁺⁺, 𝐙⁻⁺ = (zeros(FT, nb, nb), zeros(FT, nb, nb))
    A⁺⁺, A⁻⁺ = (zeros(FT, B_dim, B_dim, n, n), zeros(FT, B_dim, B_dim, n, n))

    # Iterate over l
    for l = m:l_max

        # B matrix for l
        𝐁 = 𝐁_all[l];

        # Construct Π matrix for l,m pair (change to in place later!)
        # See eq. 15 in Sanghavi 2014, note that P,R,T are already normalized
        Π  = construct_Π_matrix(mod, P, R, T, l, m)
        Π⁻ = construct_Π_matrix(mod, P⁻, R⁻, T⁻, l, m)

        # Iterate over angles
        for i in eachindex(μ), j in eachindex(μ)
            if B_dim == 1
                A⁺⁺[B_dim,B_dim,i,j] += Π[i] * 𝐁 * Π[j]
                A⁻⁺[B_dim,B_dim,i,j] += Π[i] * 𝐁 * Π⁻[j]
            else
                A⁺⁺[:,:,i,j] += Π[i] * 𝐁 * Π[j]
                A⁻⁺[:,:,i,j] += Π[i] * 𝐁 * Π⁻[j]
            end
        end
    end

    # Now get to the Z part:
    for imu in eachindex(μ), jmu in eachindex(μ)

        # Indices adjusted for size of A
        ii, jj = ((imu - 1) * B_dim, (jmu - 1) * B_dim)

        # This is equivalent to Z̄ = 1/(1+δ) * C̄m+S̄m = 1/(1+δ) * (A+DAD+AD-DA)
        # (see eq 11 in Sanghavi et al, 2013)
        for i in 1:B_dim, j in 1:B_dim
            𝐙⁺⁺[ii + i,jj + j] = 2fact * A⁺⁺[i,j,imu,jmu]
            if i <= 2 && j >= 3
                𝐙⁻⁺[ii + i,jj + j] = -2fact * A⁻⁺[i,j,imu,jmu]
            elseif i >= 3 && j <= 2
                𝐙⁻⁺[ii + i,jj + j] = -2fact * A⁻⁺[i,j,imu,jmu]
            else
                𝐙⁻⁺[ii + i,jj + j] = 2fact * A⁻⁺[i,j,imu,jmu]
            end
        end
    end

    # Return Z-moments
    return arr_type(𝐙⁺⁺), arr_type(𝐙⁻⁺)
end


text:
"""
    abstract type AEdge end

An abstract edge type.
"""
abstract type AEdge end

"""
    src(e)

Returns the source of an edge.
"""
src(e::AEdge) = error("Method not defined")

"""
    dst(e)

Returns the destination of an edge.
"""
dst(e::AEdge) = error("Method not defined")

"""
    reverse(e)

Returns an edge with swapped `src(e)` and `dst(e)`.
"""
reverse(e::AEdge) = error("Method not defined")

"""
    abstract type AGraph end

Abstract undirected graph type
"""
abstract type AGraph end

"""
    abstract ADiGraph

Abstract directed graph type
"""
abstract type ADiGraph end

"""
    const AGraphOrDiGraph = Union{AGraph, ADiGraph}

Union of [`AGraph`](@ref) and [`ADiGraph`](@ref).
"""
const AGraphOrDiGraph = Union{AGraph, ADiGraph}

####### Required interface for concrete types ########################

"""
    nv(g)

The number of vertices in `g`.
"""
nv(g::AGraphOrDiGraph) = error("Method not defined")

"""
    ne(g)

The number of edges in `g`.

Time Complexity: O(1)
"""
ne(g::AGraphOrDiGraph) = error("Method not defined")

"""
    add_edge!(g, e) -> (ok, new_edge)

Add to `g` the edge `e`.

    add_edge!(g, u, v) -> (ok, new_edge)

Add to `g` an edge from `u` to `v`.

`ok=false` if add fails (e.g. if vertices are not in the graph or the edge
is already present) and `true` otherwise.
`new_edge` is the descriptor of the new edge.
"""
add_edge!(g::AGraphOrDiGraph, u, v) = error("Method not defined")

"""
    rem_edge!(g, e)

Remove the edge `e`.

    rem_edge!(g, u, v)

Remove the edge from `u` to `v`.

Returns false if edge removal fails (e.g., if the edge does not exist) and true otherwise.
"""
rem_edge!(g::AGraphOrDiGraph, u, v) = error("Method not defined")

"""
    add_vertex!(g)

Add a new vertex to the graph `g`.
"""
add_vertex!(g::AGraphOrDiGraph) = error("Method not defined")

# length() has to be appliable to the result
"""
    in_neighbors(g, v)

Returns an iterable to all neighbors connected to vertex `v` by an incoming edge.

NOTE: it may return a reference, not a copy. Do not modify result.
"""
in_neighbors(g::ADiGraph, v) = error("Method not defined")

# length() has to be appliable to the result
"""
    out_neighbors(g::AGraphOrDiGraph, v)

Returns an iterable to all neighbors connected to vertex `v` by an outgoing edge.

NOTE: it may return a reference, not a copy. Do not modify result.
"""
out_neighbors(g::AGraphOrDiGraph, v) = error("Method not defined")

"""
    edge(g, u, v)

Returns an edge from 'u' to 'v'. The edge doesn't necessarily exists
in `g`.
"""
edge(g::AGraphOrDiGraph, u, v) = Edge{Int}(u, v)

#TODO check consistency
"""
    edgetype(g)
    edgetype(G)

Returns the type of edges of graph `g` (or graph type `G`), i. e.
the element type returned of the iterator `edges(g)`.
"""
edgetype(::Type{G}) where {G<:AGraphOrDiGraph} = Edge{Int}

"""
    vertextype(g)
    vertextype(G)

Returns the integer type of vertices of graph `g` (or graph type `G`).
"""
vertextype(::Type{G}) where {G<:AGraphOrDiGraph} = Int

"""
    graphtype{G<:AGraphOrDiGraph}(::Type{G})

The graph type corresponding to `G`. If `G<:AGraph` returns `G`,
if `G<:ADiGraph` returns a type `H<:AGraph`.
"""
graphtype(::Type{G}) where {G<:AGraphOrDiGraph} = error("Method not defined")

"""
    digraphtype{G<:AGraphOrDiGraph}(::Type{G})

The digraph type corresponding to `G`. If `G<:ADiGraph` returns `G`,
if `G<:AGraph` returns a type `H<:ADiGraph`.
"""
digraphtype(::Type{G}) where {G<:AGraphOrDiGraph} = error("Method not defined")

"""
    pop_vertex!(g)

Remove the last vertex of `g`. Equivalent to rem_vertex!(g, nv(g)).
"""
pop_vertex!(g::AGraphOrDiGraph) = error("Method not defined")

"""
    swap_vertices!(g, u, v)

Swap the labels of vertices `u` and `v`.
In the new graph all old neighbors of vertex `u` will be neighbors of `v` and
viceversa.
"""
swap_vertices!(g::AGraphOrDiGraph, u::Integer, v::Integer) = error("Method not defined")

ext: jl


text:
"""
  generateCellList(nodeList,square,L)

  Generates list of cells, where L is the width of each cell block

INPUT
nodeList
totalBounds
L

OUTPUT
cl
"""
function generate_cell_list(particleList,totalBounds::Vector{T},L::T) where T<:Real
  X0=totalBounds[1]; X1=totalBounds[2]; Y0=totalBounds[3]; Y1=totalBounds[4]
  Nx = Int(ceil((X1 - X0)/L))
  Ny = Int(ceil((Y1 - Y0)/L))
  n_nodes = length(particleList)

  neighborMat = zeros(Int,8)
  cl = CellList(Array{Cell}(undef,Nx*Ny),totalBounds,L)

  tempZeroArr = Array{Int}(undef,0)

  # compute neighbors of cells
  k=1
  for i=0:(Nx-1), j=1:Ny
    matInd = 1

    # zero-out neighborMat
    fill!(neighborMat,zero(Int))

    # compute bounding square
    x0 = X0 + i*L; x1 = X0 + (i+1)*L; y0 = Y1-j*L; y1 = Y1-(j-1)*L

    # compute neighbors
    if i>0 neighborMat[matInd] = (i-1)*Nx+j; matInd+=1 end                  # left
    if j>1 neighborMat[matInd] = i*Nx+j-1; matInd+=1 end                    # top
    if i<(Nx-1) neighborMat[matInd] = (i+1)*Nx+j; matInd+=1 end             # right
    if j<Ny neighborMat[matInd] = i*Nx+j+1; matInd+=1 end                   # bottom
    if (j>1 && i<(Nx-1)) neighborMat[matInd] = (i+1)*Nx+j-1; matInd+=1 end  # top-right
    if (j>1 && i>0) neighborMat[matInd] =(i-1)*Nx+j-1; matInd+=1 end        # top-left
    if (j<Ny && i<(Nx-1)) neighborMat[matInd] = (i+1)*Nx+j+1; matInd+=1 end # bottom-right
    if (j<Ny && i>0) neighborMat[matInd] = (i-1)*Nx+j+1; matInd+=1 end      # bottom-left

    # generate cell
    cl.cells[k] = Cell(copy(tempZeroArr),[x0,x1,y0,y1],neighborMat[1:(matInd-1)])
    k+=1
  end

  # assign nodes to respective cells
  for nodeInd=1:n_nodes
    for cellInd=1:Nx*Ny
      if is_inside_rect(cl.cells[cellInd].square,particleList[nodeInd].pos)
        push!(cl.cells[cellInd].particleIDList,nodeInd)
        break
      end
    end
  end

  return cl
end

"""

"""
function update_cell_list!(cl::CellList,particleList)
  X0=cl.bounds[1]; X1=cl.bounds[2]; Y0=cl.bounds[3]; Y1=cl.bounds[4]
  Nx = Int(ceil((X1 - X0)/cl.sideLength))
  Ny = Int(ceil((Y1 - Y0)/cl.sideLength))
  nNodes = length(particleList)

  # 'zero-out' cell node lists
  for cellInd=1:Nx*Ny
    cl.cells[cellInd].particleIDList = Array{Int}(undef,0)
  end

  # assign nodes to respective cells
  for nodeInd=1:nNodes
    for cellInd=1:Nx*Ny
      if is_inside_rect(cl.cells[cellInd].square,particleList[nodeInd].pos)
        push!(cl.cells[cellInd].particleIDList,nodeInd)
        break
      end
    end
  end

  nothing
end

"""

updates cell list with assumption that most nodes are in same cell as before, or have gone to neighbors

should be much faster, but isn't......

DO NOT USE

function updateCellList_SMART!(cl::CellList,nodeList::Vector{Point2D{T}}) where T<:Real
  X0=cl.bounds[1]; X1=cl.bounds[2]; Y0=cl.bounds[3]; Y1=cl.bounds[4]
  Nx = Int(ceil((X1 - X0)/cl.sideLength))
  Ny = Int(ceil((Y1 - Y0)/cl.sideLength))
  nNodes = length(nodeList)

  nodesToFind = Array{Int}(undef,0)

  # assign nodes to respective cells
  for cellInd=1:Nx*Ny
    numDeleted = 0
    for ti=1:length(cl.cells[cellInd].nodeList)
      found = false
      if !(is_inside_rect(cl.cells[cellInd].square,nodeList[cl.cells[cellInd].nodeList[ti-numDeleted]]))
        # search neighbors of cell for node
        for nbrs in cl.cells[cellInd].neighborList
          if is_inside_rect(cl.cells[nbrs].square,nodeList[cl.cells[cellInd].nodeList[ti-numDeleted]])
            push!(cl.cells[nbrs].nodeList,cl.cells[cellInd].nodeList[ti-numDeleted])
            found = true
            break
          end
        end

        if !(found)
          push!(nodesToFind,cl.cells[cellInd].nodeList[ti-numDeleted])
        end
        # delete node from cell nodelist
        #println("node deleted")
        deleteat!(cl.cells[cellInd].nodeList,ti-numDeleted)
        numDeleted += 1
      end
    end
  end

  # for nodes not found, search through all elements
  #println(length(nodesToFind))
  for nodeInd in nodesToFind
    for cellInd=1:Nx*Ny
      if is_inside_rect(cl.cells[cellInd].square,nodeList[nodeInd])
        push!(cl.cells[cellInd].nodeList,nodeInd)
        break
      end
    end
  end

  nothing
end
"""


text:
#!/usr/bin/env julia
# Tim Sterne-Weiler 2016

const dir = abspath( splitdir(@__FILE__)[1] )
const ver = chomp(readline(open(dir * "/VERSION")))

tic()
println( STDERR, "Whippet $ver loading and compiling... " )

using ArgParse
using Glob

push!( LOAD_PATH, dir * "/../src" )
using Whippet

function parse_cmd()
  s = ArgParseSettings()
  # TODO finish options...
  @add_arg_table s begin
    "--a", "-a"
      help     = "Replicates for Set A -- Could be: pattern to glob.psi (common-filename-segment [*.psi*]), or comma delimited list of filenames. ie. (-a sample_a) would work for sample_a-rep1.psi.gz,sample_a-rep2.psi.gz,..."
      arg_type = String
#      required = true
    "--b", "-b"
      help     = "Replicates for Set B -- Same rules as for (-a)"
      arg_type = String
#      required = true
    "--out", "-o"
      help     = "Core file name to send .diff.gz output to!"
      arg_type = String
      default  = fixpath( "./output" )
    "--directory", "-d"
      help     = "Directory to search for file patterns or list in -a and -b"
      arg_type = String
      default  = "."
    "--min-reads", "-r"
      help     = "Minimum number of reads for a single event to be included!"
      arg_type = Int64
      default  = 5
    "--min-samples", "-s"
      help     = "Minimum number of samples in a or in b for each event to be considered!"
      arg_type = Int64
      default  = 1
#=    "--min-delta-psi", "-m"
      help     = "Calculate max probability of |deltaPsi| greater than this value (default is 0.0, it is not advisable to change this)."
      arg_type = Float64
      default  = 0.0 =#
    "--emperical-size", "-e"
      help     = "Emperical distribution size to sample from."
      arg_type = Int64
      default  = 1000
    "--seed", "-g"
      help     = "Seed the RNG (Int) for reproducible results on successive runs"
      arg_type = Int64
      default  = 123456
  end
  return parse_args(s)
end

function retrievefilelist( pattern::String, dir::String )
   list = Vector{String}()
   if search(pattern, ',') > 0
      tmp = split( pattern, ',', keep=false )
   else
      tmp = glob( "*" * pattern * "*.psi*", dir )
   end
   # now clean the return
   for file in tmp
      push!( list, string(file) )
   end
   list
end

function main()
   args  = parse_cmd()
   println(STDERR, " $( round( toq(), 6 ) ) seconds" )
   srand( args["seed"] )
   dir   = fixpath( args["directory"] )
   lista = retrievefilelist( args["a"], dir )
   listb = retrievefilelist( args["b"], dir )
   println(STDERR, "Sample A: $(join(map(basename, lista), ','))")
   println(STDERR, "Sample B: $(join(map(basename, listb), ','))")
   if length(lista) <= 0 || length(listb) <= 0
      error("Unable to match files! length(a) == $(length(lista)), length(b) == $(length(listb))!")
   end

   astreams = open_streams( lista )
   bstreams = open_streams( listb )

   println(STDERR, "Now processing files and calculating posterior distributions...")
   @timer process_psi_files( args["out"] * ".diff.gz", astreams, bstreams,
                             min_samp=args["min-samples"],
                             min_reads=args["min-reads"],
                             amt=0.0,
                             size=args["emperical-size"] )
   println(STDERR, "Whippet $ver done." )
end

@timer main()


text:
module GMT

using Printf
using Dates

# Need to know what GMT version is available or if none at all to warn users on how to install GMT.
try
	#global const GMTver = Meta.parse(readlines(`gmt --version`)[1][1:3])
	v =split(readlines(`gmt --version`)[1][1:5],'.')
	global const GMTver = Meta.parse(v[1] * '.' * v[2]) + (Meta.parse(v[3]) * 0.01)
catch
	global const GMTver = 0.0
end

global legend_type  = nothing
global current_cpt  = nothing		# To store the current palette
const global img_mem_layout = [""]			# "TCP"	 For Images.jl. The default is "TRBa"
const global grd_mem_layout = [""]			# "BRP" is the default for GMT PS images.
const global current_view   = [""]			# To store the current viewpoint (-p)
const global multi_col   = Array{Bool,1}(undef,1)		# To allow plottig multiple columns at once (init to false)
const global IamModern   = Array{Bool,1}(undef,1)		# To know if we are in modern mode
const global FirstModern = Array{Bool,1}(undef,1)		# To know
const global IamSubplot  = Array{Bool,1}(undef,1)		# To know if we are in subplot mode
const global usedConfPar = Array{Bool,1}(undef,1)		# Hacky solution for the session's memory trouble
const global convert_syntax = Array{Bool,1}(undef,1);convert_syntax[1] = false	# To only convert to hard core GMT syntax (like Vd=2)
const global show_kwargs = Array{Bool,1}(undef,1);show_kwargs[1] = false	# To just print the kwargs of a option call)
const global FMT = ["ps"]
const def_fig_size  = "12c/8c"              # Default fig size for plot like programs
const def_fig_axes  = " -Baf -BWSen"        # Default fig axes for plot like programs
const def_fig_axes3 = " -Baf -Bza"  		#		"" but for 3D views

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
	@eval Base.Experimental.@optlevel 1
end

export
	GMTver, FMT, gmt,
	arrows, arrows!, bar, bar!, bar3, bar3!, hlines, hlines!, lines, lines!, legend, legend!, vlines, vlines!,
	basemap, basemap!, blockmean, blockmedian, blockmode, clip, clip!, coast, coast!, colorbar, colorbar!,
	colorscale, colorscale!, contour, contour!, contourf, contourf!, events, filter1d, fitcircle, gmt2kml,
	gmtconnect, gmtconvert, gmtinfo, gmtmath, gmtregress, gmtread, gmtselect, gmtset, gmtsimplify, gmtspatial,
	gmtvector, gmtwrite, gmtwhich, grd2cpt, grd2kml, grd2xyz, grdblend, grdclip, grdcontour, grdcontour!, grdcut,
	grdedit, grdfft, grdfilter, grdgradient, grdhisteq, grdimage, grdimage!, grdinfo, grdfill, grdlandmask, grdmath,
	grdmask, grdpaste, grdproject, grdsample, grdtrack, grdtrend, grdvector, grdvector!, grdview, grdview!, grdvolume,
	greenspline, mat2grid, mat2img, histogram, histogram!, image, image!, image_alpha!, imshow, kml2gmt, logo, logo!,
	makecpt, mask, mask!, mapproject, movie, nearneighbor, plot, plot!, plot3, plot3!, plot3d, plot3d!, project, pscontour,
	pscontour!, psconvert, psbasemap, psbasemap!, psclip, psclip!, pscoast, pscoast!, psevents, pshistogram, pshistogram!,
	psimage, psimage!, pslegend, pslegend!, psmask, psmask!, psrose, psrose!, psscale, psscale!, pssolar, pssolar!,
	psternary, psternary!, pstext, pstext!, pswiggle, pswiggle!, psxy, psxy!, psxyz, psxyz!, regress, resetGMT, rose,
	rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, spectrum1d, sphdistance, sphinterpolate,
	sphtriangulate, surface, ternary, ternary!, text, text!, text_record, trend1d, trend2d, triangulate, splitxyz,
	decorated, vector_attrib, wiggle, wiggle!, xyz2grd, gmtbegin, gmtend, gmthelp, subplot, gmtfig, inset, showfig,
	image_alpha!, mat2ds, mat2grid, mat2img, linspace, logspace, contains, fields, tic, toc

include("common_docs.jl")
include("libgmt_h.jl")
if (GMTver >= 6)
	include("libgmt.jl")
end
include("gmt_main.jl")
include("grd_operations.jl")
include("common_options.jl")
include("gmtbegin.jl")
include("blocks.jl")
include("contourf.jl")
include("filter1d.jl")
include("fitcircle.jl")
include("gmt2kml.jl")
include("gmtconnect.jl")
include("gmtconvert.jl")
include("gmtinfo.jl")
include("gmtlogo.jl")
include("gmtreadwrite.jl")
include("gmtset.jl")
include("gmtselect.jl")
include("gmtsimplify.jl")
include("gmtspatial.jl")
include("gmtregress.jl")
include("gmtvector.jl")
include("gmtwich.jl")
include("grd2cpt.jl")
include("grd2kml.jl")
include("grd2xyz.jl")
include("grdblend.jl")
include("grdclip.jl")
include("grdcut.jl")
include("grdedit.jl")
include("grdcontour.jl")
include("grdfft.jl")
include("grdfill.jl")
include("grdfilter.jl")
include("grdgmtmath.jl")
include("grdhisteq.jl")
include("grdinfo.jl")
include("grdimage.jl")
include("grdgradient.jl")
include("grdlandmask.jl")
include("grdmask.jl")
include("grdpaste.jl")
include("grdproject.jl")
include("grdsample.jl")
include("grdtrack.jl")
include("grdtrend.jl")
include("grdvector.jl")
include("grdview.jl")
include("grdvolume.jl")
include("greenspline.jl")
include("imshow.jl")
include("kml2gmt.jl")
include("makecpt.jl")
include("mapproject.jl")
include("movie.jl")
include("nearneighbor.jl")
include("plot.jl")
include("project.jl")
include("psbasemap.jl")
include("psclip.jl")
include("pscoast.jl")
include("pscontour.jl")
include("psconvert.jl")
include("pshistogram.jl")
include("psimage.jl")
include("pslegend.jl")
include("psmask.jl")
include("psscale.jl")
include("psrose.jl")
include("pssolar.jl")
include("pstext.jl")
include("psxy.jl")
include("pswiggle.jl")
include("sample1d.jl")
include("spectrum1d.jl")
include("sphdistance.jl")
include("sphinterpolate.jl")
include("sphtriangulate.jl")
include("splitxyz.jl")
include("surface.jl")
include("subplot.jl")
include("triangulate.jl")
include("trend1d.jl")
include("trend2d.jl")
include("xyz2grd.jl")

function __init__()
	if (5 <= GMTver < 6.0)  println("\n\tGMT version 5 is no longer supported (support ended at 0.23). Must uptdate."); return  end

	if (GMTver == 0.0)
		println("\n\nYou don't seem to have GMT installed and I don't install it automatically,\nso you will have to do it yourself.")
		t = "\n\t\t https://github.com/GenericMappingTools/gmt/releases"
		if (Sys.iswindows() && Sys.WORD_SIZE == 64)
			println("Download and install the official version at (the '..._win64.exe':" * t)
		elseif (Sys.iswindows() && Sys.WORD_SIZE == 32)
			println("Download and install the official version at (the '..._win32.exe':" * t)
		elseif (Sys.isapple())
			println("Install GMT with Homebrew: brew install gmt ghostscript ffmpeg")
		else
			println("https://github.com/GenericMappingTools/gmt/blob/master/INSTALL.md#linux")
		end
		return
	end
	clear_sessions(3600)		# Delete stray sessions dirs older than 1 hour
	global API = GMT_Create_Session("GMT", 2, GMT_SESSION_NOEXIT + GMT_SESSION_EXTERNAL + GMT_SESSION_COLMAJOR)
	if (API == C_NULL)  error("Failure to create a GMT Session")  end
	if haskey(ENV,"JULIA_GMT_IMGFORMAT")  FMT[1] = ENV["JULIA_GMT_IMGFORMAT"]  end
end

if (GMTver >= 6)  include("get_enums.jl")  end	# Needed to cheat the autoregister autobot

include("precompile_GMT_i.jl")
_precompile_()

end # module


text:
# Interface function: nnmf

function nnmf(X::AbstractMatrix{T}, k::Integer;
              init::Symbol=:nndsvdar,
              alg::Symbol=:greedycd,
              maxiter::Integer=100,
              tol::Real=cbrt(eps(T)/100),
              replicates::Integer=1,
              W0::Union{AbstractMatrix{T}, Nothing}=nothing,
              H0::Union{AbstractMatrix{T}, Nothing}=nothing,
              update_H::Bool=true,
              verbose::Bool=false) where T

    eltype(X) <: Number && all(t -> t >= zero(T), X) || throw(ArgumentError("The elements of X must be non-negative."))

    p, n = size(X)
    k <= min(p, n) || throw(ArgumentError("The value of k should not exceed min(size(X))."))

    replicates >= 1 || throw(ArgumentError("The value of replicates must be positive."))

    if !update_H && init != :custom
        @warn "Only W will be updated."
    end

    if init == :custom
        W0 !== nothing && H0 !== nothing || throw(ArgumentError("To use :custom initialization, set W0 and H0."))
        eltype(W0) <: Number && all(t -> t >= zero(T), W0) || throw(ArgumentError("The elements of W0 must be non-negative."))
        p0, k0 = size(W0)
        p == p0 && k == k0 || throw(ArgumentError("Invalid size for W0."))
        eltype(H0) <: Number && all(t -> t >= zero(T), H0) || throw(ArgumentError("The elements of H0 must be non-negative."))
        k0, n0 = size(H0)
        k == k0 && n == n0 || throw(ArgumentError("Invalid size for H0."))
    else
        W0 === nothing && H0 === nothing || @warn "Ignore W0 and H0 except for :custom initialization."
    end

    # determine whether H needs to be initialized
    initH = alg != :projals

    # perform initialization
    if init == :random
        W, H = randinit(X, k; zeroh=!initH, normalize=true)
    elseif init == :nndsvd
        W, H = nndsvd(X, k; zeroh=!initH)
    elseif init == :nndsvda
        W, H = nndsvd(X, k; variant=:a, zeroh=!initH)
    elseif init == :nndsvdar
        W, H = nndsvd(X, k; variant=:ar, zeroh=!initH)
    elseif init == :spa
        W, H = spa(X, k)
    elseif init == :custom
        W, H = W0, H0
    else
        throw(ArgumentError("Invalid value for init."))
    end

    # choose algorithm
    if alg == :projals
        alginst = ProjectedALS{T}(maxiter=maxiter, tol=tol, verbose=verbose, update_H=update_H)
    elseif alg == :alspgrad
        alginst = ALSPGrad{T}(maxiter=maxiter, tol=tol, verbose=verbose, update_H=update_H)
    elseif alg == :multmse
        alginst = MultUpdate{T}(obj=:mse, maxiter=maxiter, tol=tol, verbose=verbose, update_H=update_H)
    elseif alg == :multdiv
        alginst = MultUpdate{T}(obj=:div, maxiter=maxiter, tol=tol, verbose=verbose, update_H=update_H)
    elseif alg == :cd
        alginst = CoordinateDescent{T}(maxiter=maxiter, tol=tol, verbose=verbose, update_H=update_H)
    elseif alg == :greedycd
        alginst = GreedyCD{T}(maxiter=maxiter, tol=tol, verbose=verbose, update_H=update_H)
    elseif alg == :spa
        if init != :spa
            throw(ArgumentError("Invalid value for init, use :spa instead."))
        end
        alginst = SPA{T}(obj=:mse)
    else
        throw(ArgumentError("Invalid algorithm."))
    end

    # run optimization
    ret = solve!(alginst, X, W, H)

    # replicates
    minobjv = ret.objvalue
    for _ in 2:replicates
        W, H = randinit(X, k; zeroh=!initH, normalize=true)
        tmp = solve!(alginst, X, W, H)
        if minobjv > tmp.objvalue
            ret = tmp
            minobjv = tmp.objvalue
        end
    end

    return ret
end


text:
using ZigZagBoomerang
using Random
using LinearAlgebra
using SparseArrays
using ForwardDiff
using Trajectories

Random.seed!(0)

n = 10 # 10-dimensional problem
# Create sparse 10-dimensional precision matrix
GS = [1=>2:10, 2=>3:5, 3=>5:10, 4=>[5], 5=>[], 6=>7:8, 7=>[], 8=>[], 9=>[10], 10=>[]] # start from conditional dependence tree
S = Matrix(1.0I, n, n)
for (i, nb) in GS, j in nb
    S[i, j] = 0.3randn()
end
Γ = sparse(S * S')

ϕ(x) = 0.5*x'*Γ*x # potential or -log(target density)
∇ϕ!(y, x) = ForwardDiff.gradient!(y, ϕ, x)

λref_bps = 1.0
x0, θ0 = randn(n), randn(n)
c_bps = 0.001
T = 1000.0
B = BouncyParticle(Γ, x0*0, λref_bps)
out1, uT, acc = pdmp(∇ϕ!, 0.0, x0, θ0, T, c_bps, B; adapt=false, factor=2.0)
using Plots
dt = 0.1
xx = trajectory(discretize(out1, dt))
p1 = plot(getindex.(xx.x, 1), getindex.(xx.x, 2), linewidth=0.4)

λref_boom = 1.0
c_boom = 3.0
B = Boomerang(Γ, x0*0, λref_boom)
out2, uT, acc = pdmp(∇ϕ!, 0.0, x0, θ0, T, c_boom, B; adapt=false, factor=2.0)
dt = 0.1
xx = trajectory(discretize(out1, dt))
p2 = plot(getindex.(xx.x, 1), getindex.(xx.x, 2), linewidth=0.4)


text:
const _prune_manifest_kwargs_docstring = """
## Required Keyword Arguments

You must specify one (and exactly one) of `project` and `project_filename`.
Similarly, you must specify one (and exactly one) of `manifest` and
`manifest_filename`.

- `project::Union{AbstractString, IO}`: the contents of the input `Project.toml` file
- `project_filename::AbstractString`: the filename of the input `Project.toml` file
- `manifest::Union{AbstractString, IO}`: the contents of the input `Manifest.toml` file
- `manifest_filename::AbstractString`: the filename of the input `Manifest.toml` file
"""

"""
    prune_manifest(; kwargs...) --> new_manifest::AbstractString

Parse the given project and manifest, and generate a new manifest that only
includes packages that are direct or indirect (recursive) dependencies of the
given project. The new manifest is returned as an `AbstractString`.

$(_prune_manifest_kwargs_docstring)
"""
function prune_manifest(; kwargs...)
    return sprint(io -> prune_manifest(io; kwargs...))
end

"""
    prune_manifest(io::IO; kwargs...)

Parse the given project and manifest, and generate a new manifest that only
includes packages that are direct or indirect (recursive) dependencies of the
given project. The new manifest is printed to the given `IO`.

$(_prune_manifest_kwargs_docstring)
"""
function prune_manifest(io::IO;
                        project::Union{AbstractString, IO, Nothing}       = nothing,
                        project_filename::Union{AbstractString, Nothing}  = nothing,
                        manifest::Union{AbstractString, IO, Nothing}      = nothing,
                        manifest_filename::Union{AbstractString, Nothing} = nothing)
    if (project !== nothing) && (project_filename !== nothing)
        throw(ArgumentError("You may not specify both `project` or `project_filename`; you must only specify one"))
    elseif (project !== nothing) && (project_filename === nothing)
        project_dict = TOML.parse(project)
    elseif (project === nothing) && (project_filename !== nothing)
        project_dict = TOML.parsefile(project_filename)
    else
        throw(ArgumentError("You must specify either `project` or `project_filename`"))
    end
    if (manifest !== nothing) && (manifest_filename !== nothing)
        throw(ArgumentError("You may not specify both `manifest` or `manifest_filename`; you must only specify one"))
    elseif (manifest !== nothing) && (manifest_filename === nothing)
        manifest_dict = TOML.parse(manifest)
    elseif (manifest === nothing) && (manifest_filename !== nothing)
        manifest_dict = TOML.parsefile(manifest_filename)
    else
        throw(ArgumentError("You must specify either `manifest` or `manifest_filename`"))
    end

    project_struct = ProjectTOMLDict(project_dict)
    manifest_struct = ManifestTOMLDict(manifest_dict)
    return prune_manifest(io, project_struct, manifest_struct)
end

function prune_manifest(io::IO,
                        project::ProjectTOMLDict,
                        manifest::ManifestTOMLDict)
    name_to_uuid = Dict{String, Base.UUID}()

    for (name, infos_vec) in pairs(manifest.manifest)
        info = only(infos_vec)
        uuid_string = info["uuid"]
        uuid = Base.UUID(uuid_string)
        condition = !haskey(name_to_uuid, name)
        msg = "duplicate package with name `$(name)`"
        condition || throw(ErrorException(msg))
        name_to_uuid[name] = uuid
    end

    for (name, uuid_string_from_project) in pairs(project.project["deps"])
        uuid_from_project = Base.UUID(uuid_string_from_project)
        condition_1 = haskey(name_to_uuid, name)
        msg_1 = "Manifest does not have a dep with name $(name)"
        condition_1 || throw(ErrorException(msg_1))
        uuid_from_manifest = name_to_uuid[name]
        condition_2 = uuid_from_project == uuid_from_manifest
        msg_2 = "For dep $(name), project UUID $(uuid_from_project) does not match manifest UUID $(uuid_from_manifest)"
        condition_2 || throw(ErrorException(msg_2))
    end

    graph = MetaGraphs.MetaDiGraph()
    MetaGraphs.set_indexing_prop!(graph, :uuid)
    for uuid in values(name_to_uuid)
        MetaGraphs.add_vertex!(graph, :uuid, uuid)
    end

    for (name, infos_vec) in pairs(manifest.manifest)
        info = only(infos_vec)
        uuid_string = info["uuid"]
        uuid = Base.UUID(uuid_string)
        deps_names = get(info, "deps", String[])
        for dep_name in deps_names
            dep_uuid = name_to_uuid[dep_name]
            a = graph[uuid, :uuid]
            b = graph[dep_uuid, :uuid]
            MetaGraphs.add_edge!(
                graph,
                graph[uuid, :uuid],
                graph[dep_uuid, :uuid],
            )
        end
    end

    project_recursive_dependencies_uuids = Set{Base.UUID}()

    for (name, uuid_string_from_project) in pairs(project.project["deps"])
        uuid_from_project = Base.UUID(uuid_string_from_project)
        out_pars_indices = MetaGraphs.bfs_parents(graph, graph[uuid_from_project, :uuid]; dir = :out)
        nonzero_out_pars_indices  = findall(x -> x != 0, out_pars_indices)
        nonzero_out_pars_uuids = getindex.(Ref(graph), nonzero_out_pars_indices, :uuid)
        for nonzero_out_par_uuid in nonzero_out_pars_uuids
            push!(project_recursive_dependencies_uuids, nonzero_out_par_uuid)
        end
    end

    output_manifest_dict = Dict{String, Any}()
    for (name, infos_vec) in pairs(manifest.manifest)
        info = only(infos_vec)
        uuid_string = info["uuid"]
        uuid = Base.UUID(uuid_string)
        if uuid in project_recursive_dependencies_uuids
            output_manifest_dict[name] = infos_vec
        end
    end

    write_manifest(io, output_manifest_dict)

    return nothing
end


text:
# Standalone test file that tests spectra visually
#using Plots # uncomment when using the plotting code below

using ClimateMachine.ConfigTypes

using ClimateMachine.Spectra:
    compute_gaussian!,
    compute_legendre!,
    SpectralSphericalMesh,
    trans_grid_to_spherical!,
    power_spectrum_1d,
    power_spectrum_2d,
    compute_wave_numbers
using FFTW

include("spherical_helper_test.jl")

FT = Float64
# -- TEST 1: power_spectrum_1d(AtmosGCMConfigType(), var_grid, z, lat, lon, weight)
nlats = 32

# Setup grid
sinθ, wts = compute_gaussian!(nlats)
yarray = asin.(sinθ) .* 180 / π
xarray = 180.0 ./ nlats * collect(FT, 1:1:(2nlats))[:] .- 180.0
z = 1

# Setup variable
mass_weight = ones(Float64, length(z));
var_grid =
    1.0 * reshape(
        sin.(xarray / xarray[end] * 5.0 * 2π) .* (yarray .* 0.0 .+ 1.0)',
        length(xarray),
        length(yarray),
        1,
    ) +
    1.0 * reshape(
        sin.(xarray / xarray[end] * 10.0 * 2π) .* (yarray .* 0.0 .+ 1.0)',
        length(xarray),
        length(yarray),
        1,
    )
nm_spectrum, wave_numbers = power_spectrum_1d(
    AtmosGCMConfigType(),
    var_grid,
    z,
    yarray,
    xarray,
    mass_weight,
);

# Check visually
plot(wave_numbers[:, 16, 1], nm_spectrum[:, 16, 1], xlims = (0, 20))
contourf(var_grid[:, :, 1])
contourf(nm_spectrum[2:20, :, 1])

# -- TEST 2: power_spectrum_gcm_2d
# Setup grid
sinθ, wts = compute_gaussian!(nlats)
yarray = asin.(sinθ) .* 180 / π
xarray = 180.0 ./ nlats * collect(FT, 1:1:(2nlats))[:] .- 180.0
z = 1

# Setup variable: use an example analytical P_nm function
P_32 = sqrt(105 / 8) * (sinθ .- sinθ .^ 3)
var_grid =
    1.0 * reshape(
        sin.(xarray / xarray[end] * 3.0 * π) .* P_32',
        length(xarray),
        length(yarray),
        1,
    )

mass_weight = ones(Float64, z);
spectrum, wave_numbers, spherical, mesh =
    power_spectrum_2d(AtmosGCMConfigType(), var_grid, mass_weight)

# Grid to spherical to grid reconstruction
reconstruction = trans_spherical_to_grid!(mesh, spherical)

# Check visually
contourf(var_grid[:, :, 1])
contourf(reconstruction[:, :, 1])
contourf(var_grid[:, :, 1] .- reconstruction[:, :, 1])

# Spectrum
contourf(
    collect(0:1:(mesh.num_fourier - 1))[:],
    collect(0:1:(mesh.num_spherical - 1))[:],
    (spectrum[:, :, 1])',
    xlabel = "m",
    ylabel = "n",
)

# Check magnitude
println(0.5 .* sum(spectrum))

dθ = π / length(wts)
cosθ = sqrt.(1 .- sinθ .^ 2)
area_factor = reshape(cosθ .* dθ .^ 2 / 4π, (1, length(cosθ)))

println(sum(0.5 .* var_grid[:, :, 1] .^ 2 .* area_factor))

# NB: can verify against published packages, e.g., https://github.com/jswhit/pyspharm


text:
#############################
##  Numerical experiments associated with section 5.1.2
#############################
## While the paper does not explicitly treat third kind Volterra equations, an extension to these
## cases is relatively straightforward. Note that while in general there won't be exponential convergence
## unless the involved powers are favorably smooth, our method still achieves higher accuracy than many competitor methods.
##

using ApproxFun, MultivariateOrthogonalPolynomials, BandedMatrices, BlockBandedMatrices, SpecialFunctions, Plots
using SparseVolterraExamples

#############################
## Problem in Equation (17)
#############################
####
## The following block computes the coefficient vector of the approximation for a given polynomial order of approximation n.
## The accuracy obtained depends on both the polynomial order for the multiplication as well as the solution.
## Choosing the same order is fine in most cases but better accuracy can be obtained by adjusting for a given problem.
function solveSec512A(n,multin)
    gf(x) = x^(2/3)*(10/3*x^(7/3)-3/16*x^(14/3));
    Kfun(x,y) = y;
    V = triVolterraFullKernelOpP01(Kfun,n,true);
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V;
        V = Conversion(Jacobi(0,0,0..1),Jacobi(1,2,0..1))[1:n,1:n]*V;
        V = Multiplication(Fun(x->x^(2/3),Jacobi(1,2,0..1),multin),Jacobi(1,2,0..1))[1:n,1:n]*Derivative(Jacobi(0,1,0..1),1)[1:n,1:n]-V
    coeff = [DirectEvalLHSP10at0(n);V[1:n-1,1:n]] \ [0;pad(Fun(x->gf(x),Jacobi(1,2, 0..1)).coefficients,n-1)]
    return coeff
end
####
## Plot the analytic solution against the computed solution
coeff = solveSec512A(50,40)
plot(Fun(Jacobi(0,1,0..1),coeff),grid=false,xlabel="x",ylabel="u(x)",label="sparse method")
plot!(x->(x)^(10/3),0,1,grid=false,xlabel="x",ylabel="u(x)",label="analytic solution")
####
## Check the numerical error.
plot(x->(Fun(Jacobi(0,1,0..1),coeff)(x)-(x)^(10/3)),0,1,grid=false,xlabel="x",ylabel="error",label=false)

#############################
## Problem in Equation (18)
#############################
####
## The following block computes the coefficient vector of the approximation for a given polynomial order of approximation n.
## The accuracy obtained depends on both the polynomial order for the multiplication as well as the solution.
## Choosing the same order is fine in most cases but better accuracy can be obtained by adjusting for a given problem.
function solveSec512B(n,multin)
    gf(x) = 9/2*x^4-1/20*x^(11/2)-1/6*x^6;
    Kfun(x,y) = sqrt(y);
    V = triVolterraFullKernelOpP01(Kfun,n,true,n);
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V;
        V = Conversion(Jacobi(0,0,0..1),Jacobi(1,2,0..1))[1:n,1:n]*V;
        V = Multiplication(Fun(x->sqrt(x),Jacobi(1,2,0..1),multin),Jacobi(1,2,0..1))[1:n,1:n]*Derivative(Jacobi(0,1,0..1),1)[1:n,1:n]-(1/20)*Conversion(Jacobi(0,1,0..1),Jacobi(1,2,0..1))[1:n,1:n]*Multiplication(Fun(x->x,Jacobi(0,1,0..1),n),Jacobi(0,1,0..1))[1:n,1:n]-V
    coeff = [DirectEvalLHSP10at0(n);V[1:n-1,1:n]] \ [0;pad(Fun(x->gf(x),Jacobi(1,2, 0..1)).coefficients,n-1)]
    return coeff
end
####
## Plot the analytic solution against the computed solution
coeff = solveSec512B(80,30)
plot(Fun(Jacobi(0,1,0..1),coeff),grid=false,xlabel="x",ylabel="u(x)",label="sparse method")
plot!(x->(x)^(9/2),0,1,grid=false,xlabel="x",ylabel="u(x)",label="analytic solution")
####
## Check the numerical error.
plot(x->(Fun(Jacobi(0,1,0..1),coeff)(x)-(x)^(9/2)),0,1,grid=false,xlabel="x",ylabel="error",label=false)

max_stars_repo_path: examples/Section 5-1-2 - Methods for third kind VIDEs.jl


text:
#--------------------------------------------------------------
# MultilevelTT package
# example code for https://arxiv.org/abs/1802.09062
#--------------------------------------------------------------
module MultilevelTT
push!(LOAD_PATH, @__DIR__) # expose all other modules in this directory
end # module


text:
using Base.LinAlg.BLAS

N = 4000
srand(0)
A = rand(N,N)
B = rand(N,N)
result = zeros(N,N)

function substep{T}(a::Matrix{T}, beta::T, b::Matrix{T}, result::Matrix{T})
    for j=1:size(a,2)
        for i=1:size(b,1)
            @inbounds result[i,j] = a[i,j] + beta*b[i,j]
        end
    end
end

function f0{T}(A::Matrix{T}, alpha::T, B::Matrix{T}, result::Matrix{T}, repeat=1)
    for i=1:repeat
        substep(A, alpha, B, result)
    end
    return nothing
end

@time f0(A, 0.1, B, result)
@time f0(A, 0.1, B, result)

A_vec = reshape(A, length(A))
result_vec = reshape(result, length(result))
N = length(A_vec)
@time BLAS.axpy!(N, 0.1, A, 1, B, 1)
@time BLAS.axpy!(N, 0.1, A, 1, B, 1)

# @time BLAS.blascopy!(N, A_vec, 1, result_vec, 1)
# @time BLAS.blascopy!(N, A_vec, 1, result_vec, 1)
# @time BLAS.scale!(0.1, A)
# @time BLAS.scale!(0.1, A)


text:
module HydrophoneCalibrations

using Interpolations
using Unitful

export volt_to_pressure, volt_to_pressure_and_phase, calibration
export parse_onda

const CALIBRATION_LOADED = Ref{Bool}() # initializes to false

include("configuration_loading.jl")
include("volt_to_pressure.jl")
include("onda_file_parsing.jl")

_load_config()

end


text:
abstract type AbstractAffectFeedforward end

get_device_type(x::AbstractAffectFeedforward) = x.device_type


text:
function DataDynamicDrill(
        u, v, _zchars::ExogTimeVars, _ichars::Vector{<:Tuple},
        ddm::DynamicDrillModel, theta;
        minmaxleases=0:3,
        nper_initial=50:50,
        tstart=5:15
    )

    num_i = length(u)
    @assert length(v) == num_i
    num_zt = length(_zchars)

    data = DataDrill(u, v, _zchars, _ichars, ddm, theta;
        minmaxleases=minmaxleases,
        nper_initial=nper_initial,
        nper_development=0:0,
        tstart=tstart,
        xdomain=0:0
    )

    @assert num_i == length(data)

    # for balanced panel of states + decisions
    lease_states    = zeros(Int, num_zt, j1length(data))
    lease_decisions = zeros(Int, num_zt, j1length(data))

    # for start/end times of each lease
    tdrill  = zeros(Int,      j1length(data)) # when lease j is 1st drilled
    texpire = fill( num_zt+1, j1length(data)) # when lease j did expire

    # for start/end times of constructed data
    tstop1 = zeros(Int, j1length(data)) # end of exploratory period
    tstop2 = zeros(Int, j1length(data)) # end of infill period

    # state space
    wp = statespace(ddm)
    last_state = length(wp)
    expired_state = exploratory_terminal(wp)

    # choose one lease per unit
    selected_initial_leases = [j1_sample(unit) for unit in data]

    # create balanced "panel"
    for unit in data
        for lease in InitialDrilling(unit)

            ztrng = ztrange(lease)
            j = _i(lease)
            lease_states[   ztrng, j] .= _x(lease)
            lease_decisions[ztrng, j] .= _y(lease)

            # if ever drilled
            if lease_ever_drilled(lease)
                obs_drill_first = findfirst(d -> d > 0, _y(lease))
                tdrill[j] = ztrng[obs_drill_first]
            end

            # if ever expired
            if lease_expired(lease)
                obs_expire = findfirst(s -> s == exploratory_terminal(wp), _x(lease))
                texpire[j] = ztrng[obs_expire]
            end

        end
    end

    # now replace decisions with those from 1 particular lease in each unit
    for unit in data

        jselect = selected_initial_leases[_i(unit)] # pick this lease

        y_jselect      = lease_decisions[:,jselect]  # pick these decisions
        x_jselect      = lease_states[:,jselect]     # pick this set of states
        tdrill_jselect = tdrill[jselect]             # 1st drilling

        for lease in InitialDrilling(unit)
            j = _i(lease)

            decisions = view(lease_decisions, :, j)
            states = view(lease_states, :, j)

            # Boolean vars to test for events
            drilled  = 0 < tdrill_jselect
            drilled_b4_j_leased = tdrill_jselect < jtstart(data,j)
            expired_b4_j_drilled = texpire[j] <= tdrill_jselect

            #------------------------------
            # update lease_j with selected lease
            #------------------------------

            # if lease_j starts AFTER first drill, we drop
            if drilled && drilled_b4_j_leased
                decisions .= states .= -1

            # if lease_j expires BEFORE we first drill
            elseif expired_b4_j_drilled
                decisions .= states .= -1

            # update lease_j
            else
                decisions .= y_jselect   # update other decisions
                tdrill[j] = tdrill_jselect
                for t in jtstart(data,j):num_zt-1
                    states[t+1] = ssprime(wp, states[t], decisions[t])
                end
            end

            # unit drilled but lease_j has expired at 1st drill. then drop
            if drilled
                if lease_states[tdrill_jselect, j] == expired_state
                    decisions .= states .= -1
                end
            end

            #------------------------------
            # new start/end times for lease_j
            #------------------------------

            # if ever drilled
            if sum(decisions) > 0
                tstop1[j] = findfirst(d -> d > 0, decisions)

                # totally drilled
                if last_state in states
                    tstop2[j] = findfirst(s -> s == last_state, states)-1

                # not totally drilled
                else
                    tstop2[j] = length(states)
                end

            # else never drilled
            else
                # lease expires
                if expired_state in states
                    tstop1[j] = findfirst(s -> s == expired_state, states)
                # or not
                else
                    tstop1[j] = length(states)
                end
                tstop2[j] = tstop1[j]
            end
        end
    end

    # exploratory & infill drilling periods for each lease_j
    jt1 = view(jtstart(data), OneTo(j1length(data)))
    ttimes1 = [a:b for (a,b) in zip(jt1, tstop1)]
    ttimes2 = [a+1:b for (a,b) in zip(tstop1, tstop2)]

    # number of wells drilled for each lease_j
    num_wells  = vec(mapslices(sum, lease_decisions; dims=1))

    # do we keep this lease?
    keep_lease = vec(mapslices( x -> all(x.>=0), lease_decisions; dims=1))
    all(keep_lease[selected_initial_leases]) || throw(error("not keeping selected leases. boo!"))

    #------------------------------
    # checks that data are correct
    #------------------------------
    for unit in data
        i = _i(unit)
        jrng = j1_range(data, i)
        tstop1s = tstop1[jrng]
        keeps = keep_lease[jrng]
        drilled = num_wells[jrng] .> 0

        # for leases we have drilled & keep, exploratory periods end at same time
        if any(drilled[keeps])
            @assert all_same_value(tstop1s[keeps])

        # otherwise we don't have infill periods
        else
            @assert all(length.( ttimes2[jrng][keeps] ) .== 0)
        end
    end

    #------------------------------
    # create new dataset
    #------------------------------

    # @show ttimes2[selected_initial_leases]
    _tptr  = 1 .+ cumsum(vcat(0, length.(ttimes1[keep_lease]), length.(ttimes2[selected_initial_leases])))
    _j1ptr = 1 .+ cumsum(vcat(0, [sum(keep_lease[j1_range(data,i)]) for i in 1:num_i]))
    _j2ptr = (last(_j1ptr)-1) .+ (1:num_i)
    _jtstart = vcat(first.(ttimes1[keep_lease]), first.(ttimes2[selected_initial_leases]))
    _j1chars = j1chars(data)[keep_lease]

    _jtend = vcat(last.(ttimes1[keep_lease]), last.(ttimes2[selected_initial_leases]))
    @assert all(_jtend .<= num_zt)

    nobs = last(_tptr)-1
    xnew = zeros(Int, nobs)
    ynew = zeros(Int, nobs)

    for i = 1:num_i
        jnewrng = range_i_to_ip1(_j1ptr, i)

        # ensure these sum to 1
        _j1chars[jnewrng] ./= sum(_j1chars[jnewrng])

        keeps = keep_lease[j1_range(data,i)]
        joldrng = collect(j1_range(data,i))[keeps]
        ttimes1_jold = ttimes1[joldrng]

        selected_j = selected_initial_leases[i]

        for (jidx, jnew) in enumerate(jnewrng)
            tnewrng = range_i_to_ip1(_tptr, jnew)
            toldrng = ttimes1_jold[jidx]
            jolddidx = joldrng[jidx]
            length(tnewrng) == length(toldrng) || throw(error("i=$i, jnew=$jnew, jidx=$jidx, tnew=$tnewrng vs $toldrng"))
            xnew[tnewrng] .= lease_states[   toldrng, jolddidx]
            ynew[tnewrng] .= lease_decisions[toldrng, jolddidx]
        end

        let tnewrng = range_i_to_ip1(_tptr, _j2ptr[i])
            tinfill_i = ttimes2[selected_j]
            @assert length(tinfill_i) == length(tnewrng)
            xnew[tnewrng] .= lease_states[   tinfill_i, selected_j]
            ynew[tnewrng] .= lease_decisions[tinfill_i, selected_j]
        end
    end

    @assert num_i ≈ sum(_j1chars)

    new_data = DataDrill(
        ddm, _j1ptr, _j2ptr, _tptr, _jtstart,
        _j1chars, ichars(data), ynew, xnew, zchars(data)
    )

    return new_data
end


text:
# Quantum Key Distribution (BB84)
# You can run this from a terminal as `julia --color=yes bb84.jl`
# Try enabling Eve (line 28) to see the key distribution fail.

using Probably

N = 100

a2b = Channel{Any}(32)
b2a = Channel{Any}(32)

eve_bits = []
eve(x) = push!(eve_bits, measure(x))

# Alice
@async begin
  bits = []
  bases = []
  for i = 1:N
    # Create a qubit in |0⟩, |1⟩, |+⟩ or |-⟩
    value, basis = rand(Bool), rand([:standard, :hadamard])
    push!(bits, value)
    push!(bases, basis)
    x = QBool()
    value && not(x)
    basis == :hadamard && hadamard(x)
    # Optionally show the qubit to eve (if you think you can trust her)
    # eve(x)
    # Send it to bob
    put!(a2b, x)
  end
  # Share measurement bases used (reveals which bits are correct without
  # revealing the values)
  bob_bases = take!(b2a)
  put!(a2b, bases)
  bits = bits[bases .== bob_bases]
  # Share a subset of the bits to verify they were not tampered with
  subset = [rand(Bool) for i = 1:length(bits)]
  put!(a2b, (subset, bits[subset]))
  bits = bits[.!subset]
end

# Bob
let
  bits = []
  bases = []
  for i = 1:N
    # Take each bit and measure it in a random basis
    # (will be correct ~50% of the time)
    basis = rand([:standard, :hadamard])
    push!(bases, basis)
    x = take!(a2b)
    basis == :hadamard && hadamard(x)
    push!(bits, measure(x))
  end
  # Share measurement bases used
  put!(b2a, bases)
  alice_bases = take!(a2b)
  bits = bits[bases .== alice_bases]
  # Check a subset of the bits
  subset, alice_bits = take!(a2b)
  if bits[subset] != alice_bits
    @warn("QKD failed, eavesdropping detected")
  else
    bits = bits[.!subset]
    @info("QKD success. Key is:")
    println(Int.(bits))
  end
end


text:
using AbstractBijections

@show f = Bijection(["X1", "X2"], ["Y1", "Y2"])
@show f["X1"]
@show inv(f)["Y1"]
@show domain(f)
@show image(f)

@show f = bijection(["X1" => "Y1", "X2" => "Y2"])
@show f["X1"] == "Y1"
@show inv(f)["Y1"] == "X1"

@show f = Bijection(["Y1", "Y2"])
@show f[1] == "Y1"
@show inv(f)["Y1"] == 1


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_ellipse(fname)
    pagewidth, pageheight = 1200, 1400
    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setopacity(0.5)
    setline(0.3)

    pagetiles = Partition(pagewidth, pageheight, 200, 200)
    for (pos, n) in pagetiles
        randomhue()
        ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :strokepreserve)
        clip()
        for i in 1:10
            ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :fill)
            randomhue()
            ellipse(pos, Point(pos.x + rand(50:pagetiles.tilewidth), pos.y), rand(50:pagetiles.tileheight), :stroke)
        end
        setline(1)
        sethue("black")
        fillstroke()
        ellipse(pos, 5, 5, :fill)
        clipreset()
        # test three point version
        circle.((pos - (50, 0), pos + (50, 0), pos + (0, 2n)), 4, :fill)
        ellipse(pos - (50, 0), pos + (50, 0), pos + (0, 2n), :stroke)
    end
    @test finish() == true
    end

fname = "ellipse-test1.pdf"
test_ellipse(fname)
println("...finished test: output in $(fname)")


text:
using Test

const testdir = dirname(@__FILE__)

println("Type '1' if you want to test GraphBLASJulia, type any other character if you
    want to test GraphBLASJuliaCvdlab")

choice = readline()

if choice == "1"
    tests = [
        "GraphBLASJulia",
        "PrimitiveJulia"
    ]
else
    tests = [
    "GraphBLASJuliaCvdlab",
    "PrimitiveJulia"
    ]
end

println("testing the following functionalities:", tests)

@testset "BFS" begin
    for t in tests
        testpath = joinpath(testdir, "$(t)/$(t).jl")
        println(testpath)
        include(testpath)
    end
end


text:
using JuMP
using GLPK
include("Data.jl")

# Int8: -128 -> 127
# Int16 : -32.768 -> 32.767
# Int32 : -2.147.483.648 -> 2.147.483.647
# Int64 : - 2^63 -> 2^63 - 1
# Float32, Float64, UINT8, UInt16, etccdddd

file = map(file->string(file), ARGS)

if(length(file) != 1)
  println("Wrong Pattern")
  println("julia main.jl instances/questao2.txt")
end

filepath = file[1]

function solver(filepath)
  capacity = Array{Int64,2}
  d1 = Main.DataReader.Data(0,0,0,0,[],[])

  capacity, d1 = Main.DataReader.readData(filepath)

  costs = zeros(Int64,(d1.N_vertex, d1.N_vertex))

  for i in 1:d1.N_vertex
    if i in d1.start_nodes
      costs[i,d1.end_node] = 1
    end
  end

  mcfp = Model(GLPK.Optimizer)

  n = d1.N_vertex

  @variable(mcfp, x[1:n,1:n] >= 0) #ok

  # capacity  --------- taok
  for i in 1:n
    for j in 1:n
      cap = capacity[i,j]
      @constraint(mcfp, x[i,j] - capacity[i,j] <= 0)
    end
  end

  #sum 0 // demand / offer
  for i in 1:n
    if i in d1.start_nodes
      @constraint(mcfp, sum(x[i,j] - x[j,i] for j in 1:n) == d1.max_possible_flow[i])
    elseif i == d1.end_node
      global max_end_flow = 0
      for node in d1.start_nodes
        max_end_flow += d1.max_possible_flow[node]
      end
      @constraint(mcfp, sum(x[i,j] - x[j,i] for j in 1:n) == max_end_flow * -1)
    else
      @constraint(mcfp, sum(x[i,j] - x[j,i] for j in 1:n) == 0)
    end
  end

  @objective(mcfp, Min, sum(costs[i,j] * x[i,j] for i in 1:n, j in 1:n))

  #println(mcfp)
  function solve_mcfp(mcfp)
    JuMP.optimize!(mcfp)
    @show JuMP.termination_status(mcfp)
    println("Max Flow = ", max_end_flow - objective_value(mcfp))

  end
  @time solve_mcfp(mcfp)
end

# run 3 times to compile the func in the first one
solver(filepath)
solver(filepath)
solver(filepath)

max_stars_repo_name: ItamarRocha/Operations-Research

max_issues_repo_name: ItamarRocha/Operations-Research

max_forks_repo_name: ItamarRocha/Operations-Research


text:
using Flux, Gym

# Load game environment
env = make("Pendulum-v0")
ctx = Ctx(env, :human_window)
reset!(env)

actions = [sample(env._env.action_space) for i=1:1000]
i = 1
done = false
reset!(env)
while i <= length(actions) && !done
    global i, done
    print("Iteration: $i, ")
    a, b, done, d = step!(env, actions[i])
    println("θ = $(env._env.state[1])")

    render!(env)
    i += 1
    sleep(0.01) # to see an animation
end


text:
#=
Utility functions used in the QuantEcon library

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-14

=#
meshgrid(x::AbstractVector, y::AbstractVector) = (repeat(x, 1, length(y))',
                                                  repeat(y, 1, length(x)))

fix(x::Real) = x >= 0 ? floor(Int, x) : ceil(Int, x)
fix!(x::AbstractArray{T}, out::Array{Int}) where {T<:Real} = map!(fix, out, x)
fix(x::AbstractArray{T}) where {T<:Real} = fix!(x, similar(x, Int))

"""
`fix(x)`

Round `x` towards zero. For arrays there is a mutating version `fix!`
"""
fix

ckron(A::AbstractArray, B::AbstractArray) = kron(A, B)
ckron(arrays::AbstractArray...) = reduce(kron, arrays)

"""
`ckron(arrays::AbstractArray...)`

Repeatedly apply kronecker products to the arrays. Equilvalent to
`reduce(kron, arrays)`
"""
ckron

"""
`gridmake!(out::AbstractMatrix, arrays::AbstractVector...)`

Like `gridmake`, but fills a pre-populated array. `out` must have size
`prod(map(length, arrays), dims = length(arrays))`
"""
function gridmake!(out, arrays::Union{AbstractVector,AbstractMatrix}...)
    lens = Int[size(e, 1) for e in arrays]

    n = sum(_i -> size(_i, 2), arrays)
    l = prod(lens)
    @assert size(out) == (l, n)

    reverse!(lens)
    repititions = cumprod(vcat(1, lens[1:end-1]))
    reverse!(repititions)
    reverse!(lens)  # put lens back in correct order

    col_base = 0

    for i in 1:length(arrays)
        arr = arrays[i]
        ncol = size(arr, 2)
        outer = repititions[i]
        inner = floor(Int, l / (outer * lens[i]))
        for col_plus in 1:ncol
            row = 0
            for _1 in 1:outer, ix in 1:lens[i], _2 in 1:inner
                out[row+=1, col_base+col_plus] = arr[ix, col_plus]
            end
        end
        col_base += ncol
    end
    return out
end

@generated function gridmake(arrays::AbstractArray...)
    T = reduce(promote_type, eltype(a) for a in arrays)
    quote
        l = 1
        n = 0
        for arr in arrays
            l *= size(arr, 1)
            n += size(arr, 2)
        end
        out = Matrix{$T}(undef, l, n)
        gridmake!(out, arrays...)
        out
    end
end

function gridmake(t::Tuple)
    all(map(x -> isa(x, Integer), t)) ||
        error("gridmake(::Tuple) only valid when all elements are integers")
    gridmake(map(x->1:x, t)...)::Matrix{Int}
end

"""
`gridmake(arrays::Union{AbstractVector,AbstractMatrix}...)`

Expand one or more vectors (or matrices) into a matrix where rows span the
cartesian product of combinations of the input arrays. Each column of the input
arrays will correspond to one column of the output matrix. The first array
varies the fastest (see example)

##### Example

```jlcon
julia> x = [1, 2, 3]; y = [10, 20]; z = [100, 200];

julia> gridmake(x, y, z)
12x3 Array{Int64,2}:
 1  10  100
 2  10  100
 3  10  100
 1  20  100
 2  20  100
 3  20  100
 1  10  200
 2  10  200
 3  10  200
 1  20  200
 2  20  200
 3  20  200
```
"""
gridmake

@doc doc"""
General function for testing for stability of matrix ``A``. Just
checks that eigenvalues are less than 1 in absolute value.

#### Arguments

- `A::Matrix` The matrix we want to check

#### Returns

- `stable::Bool` Whether or not the matrix is stable

"""
function is_stable(A::AbstractMatrix)

    # Check for stability by testing that eigenvalues are less than 1
    stable = true
    d = eigvals(A)
    if maximum(abs, d) > 1.0
        stable = false
    end
    return stable

end

"""
The total number of m-part compositions of n, which is equal to

(n + m - 1) choose (m - 1)

##### Arguments

- `m`::Int : Number of parts of composition
- `n`::Int : Integer to decompose

##### Returns

- ::Int - Total number of m-part compositions of n
"""
function num_compositions(m, n)
    return binomial(n+m-1, m-1)
end

@doc doc"""
Construct an array consisting of the integer points in the
(m-1)-dimensional simplex $\{x \mid x_0 + \cdots + x_{m-1} = n
\}$, or equivalently, the m-part compositions of n, which are listed
in lexicographic order. The total number of the points (hence the
length of the output array) is L = (n+m-1)!/(n!*(m-1)!) (i.e.,
(n+m-1) choose (m-1)).

##### Arguments

- `m`::Int : Dimension of each point. Must be a positive integer.
- `n`::Int : Number which the coordinates of each point sum to. Must
             be a nonnegative integer.

##### Returns
- `out`::Matrix{Int} : Array of shape (m, L) containing the integer
                       points in the simplex, aligned in lexicographic
                       order.

##### Notes

A grid of the (m-1)-dimensional *unit* simplex with n subdivisions
along each dimension can be obtained by `simplex_grid(m, n) / n`.

##### Examples

>>> simplex_grid(3, 4)

3×15 Array{Int64,2}:
 0  0  0  0  0  1  1  1  1  2  2  2  3  3  4
 0  1  2  3  4  0  1  2  3  0  1  2  0  1  0
 4  3  2  1  0  3  2  1  0  2  1  0  1  0  0

##### References

A. Nijenhuis and H. S. Wilf, Combinatorial Algorithms, Chapter 5,
   Academic Press, 1978.
"""
function simplex_grid(m, n)
    # Get number of elements in array and allocate
    L = num_compositions(m, n)
    out = Matrix{Int}(undef, m, L)

    x = zeros(Int, m)
    x[m] = n

    # Fill in first column
    copyto!(out, 1, x, 1, m)

    h = m
    for i in 2:L
        h -= 1

        val = x[h+1]
        x[h+1] = 0
        x[m] = val - 1
        x[h] += 1

        copyto!(out, m*(i-1) + 1, x, 1, m)

        if val != 1
            h = m
        end
    end

    return out
end

@doc doc"""
Return the index of the point x in the lexicographic order of the
integer points of the (m-1)-dimensional simplex $\{x \mid x_0 +
\cdots + x_{m-1} = n\}$.

##### Arguments

- `x`::Array{Int,1} : Integer point in the simplex, i.e., an array of
                      m nonnegative integers that sum to n.
- `m`::Int : Dimension of each point. Must be a positive integer.
- `n`::Int : Number which the coordinates of each point sum to. Must be a
             nonnegative integer.

##### Returns
- `idx`::Int : Index of x.

"""
function simplex_index(x, m, n)
    # If only one element then only one point in simplex
    if m==1
        return 1
    end

    decumsum = reverse(cumsum(reverse(x[2:end])))
    idx = binomial(n+m-1, m-1)
    for i in 1:m-1
        if decumsum[i] == 0
            break
        end

        idx -= num_compositions(m - (i-1), decumsum[i]-1)
    end

    return idx
end

"""
    next_k_array!(a)

Given an array `a` of k distinct positive integers, sorted in
ascending order, return the next k-array in the lexicographic
ordering of the descending sequences of the elements, following
[Combinatorial number system]
(https://en.wikipedia.org/wiki/Combinatorial_number_system). `a` is
modified in place.

# Arguments

- `a::Vector{<:Integer}`: Array of length k.

# Returns

- `a::Vector{<:Integer}`: View of `a`.

# Examples

```julia
julia> n, k = 4, 2;

julia> a = collect(1:2);

julia> while a[end] <= n
           @show a
           next_k_array!(a)
       end
a = [1, 2]
a = [1, 3]
a = [2, 3]
a = [1, 4]
a = [2, 4]
a = [3, 4]
```
"""
function next_k_array!(a::Vector{<:Integer})

    k = length(a)
    if k == 1 || a[1] + 1 < a[2]
        a[1] += 1
        return a
    end

    a[1] = 1
    i = 2
    x = a[i] + 1

    while i < k && x == a[i+1]
        i += 1
        a[i-1] = i - 1
        x = a[i] + 1
    end
    a[i] = x

    return a
end

"""
    k_array_rank([T=Int], a)

Given an array `a` of k distinct positive integers, sorted in
ascending order, return its ranking in the lexicographic ordering of
the descending sequences of the elements, following
[Combinatorial number system]
(https://en.wikipedia.org/wiki/Combinatorial_number_system).

# Notes

`InexactError` exception will be thrown, or an incorrect value will be
returned without warning if overflow occurs during the computation.
It is the user's responsibility to ensure that the rank of the input
array fits within the range of `T`; a sufficient condition for it is
`binomial(BigInt(a[end]), BigInt(length(a))) <= typemax(T)`.

# Arguments

- `T::Type{<:Integer}`: The numeric type of ranking to be returned.
- `a::Vector{<:Integer}`: Array of length k.

# Returns

- `idx::T`: Ranking of `a`.
"""
function k_array_rank(T::Type{<:Integer}, a::Vector{<:Integer})
    if T != BigInt
        binomial(BigInt(a[end]), BigInt(length(a))) ≤ typemax(T) ||
        throw(InexactError(:Binomial, T, a[end]))
    end
    k = length(a)
    idx = one(T)
    for i = 1:k
        idx += binomial(T(a[i])-one(T), T(i))
    end

    return idx
end

k_array_rank(a::Vector{<:Integer}) = k_array_rank(Int, a)


text:
using Flux.Tracker, Base.Test, NNlib
using Flux.Tracker: gradcheck

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(f(xs...)), xs...)
gradtest(f, dims...) = gradtest(f, rand.(dims)...)

@testset "Tracker" begin

@test gradtest((x, W, b) -> σ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> σ.(W*x .+ b), (5,3), (2,5), 2)

@test gradtest(x -> sin.(sum(x, (2, 3))), (3,4,5))

@test gradtest(x -> softmax(x).*(1:3), 3)
@test gradtest(x -> softmax(x).*(1:3), (3,5))

@test gradtest(Flux.mse, rand(5,5), rand(5, 5))
@test gradtest(Flux.crossentropy, rand(5,5), rand(5, 5))

@test gradtest(x -> x', rand(5))

@test gradtest(vcat, rand(5), rand(3))
@test gradtest(vcat, rand(2,3), rand(3,3))

@testset "mean" begin
  @test gradtest(mean, rand(2, 3))

  @test gradtest(x -> mean(x, 1), rand(2, 3))
  @test gradtest(x -> mean(x, 2), rand(2, 3))
  @test gradtest(x -> mean(x, 3), rand(2, 3, 4))

  @test gradtest(x -> mean(x, [1, 2]), rand(2, 3, 4))
end

@test gradtest(rand(5)) do x
  y = x.^2
  2y + x
end

for T in [Float32, Float64]
    @test isa(param(T(1)), TrackedArray{T, 0})
    @test isa(param(rand(T, 2)), TrackedArray{T, 1})
    @test isa(param(rand(T, 2,2)), TrackedArray{T, 2})
end

# TODO: do we wand this behaviour ??
F = typeof(AbstractFloat(1))
for T in [Int32, Int64]
    @test isa(param(T(1)), TrackedArray{F, 0})
    @test isa(param(rand(T, 2)), TrackedArray{F, 1})
    @test isa(param(rand(T, 2,2)), TrackedArray{F, 2})
end

end #testset


text:
function testsatisfiability(formula, issatisfiable; maxsearchdepth=5)
    if issatisfiable == false
        @customtest is_satisfiable(formula, maxsearchdepth=maxsearchdepth) == false
    else
        @customtest_throws OverMaxSearchDepthError is_satisfiable(formula, maxsearchdepth=maxsearchdepth)
    end
end

@testset "is_satisfiable" begin
    # Basic tests, easy to verify by hand.
    testsatisfiability("P(x) ∧ ¬P(x)", false)
    testsatisfiability("P(x) & ~P(y)", true)
    testsatisfiability("(P(x) & ~P(y)) & (~P(x) & P(y))", false)
    testsatisfiability("(P(x) & ~P(y)) ∨ (C(y) & ~C(y))", true)
    testsatisfiability("A() & (B() | ~C()) & ~B() & C()", false)
    testsatisfiability("*{x}(P(x) & ~P(f(x)))", false)
    testsatisfiability("(S(x,p) | G(p)) & (~L(x) | G(p)) & (~S(x,p)) & (~G(p))", false)
    testsatisfiability("(S(x,p) | G(p)) & (~L(x) | G(p))", true)

    # Test based on "Logik für Informatiker", exercise 76.
    testsatisfiability(CNF([
        Clause([Literal(false, Formula("S(x,p)")), Literal(false, Formula("G(p)"))]),
        Clause([Literal(true, Formula("L(x)")), Literal(false, Formula("G(p)"))]),
        Clause([Literal(true, Formula("S(x,p)"))]),
        Clause([Literal(true, Formula("G(p)"))]),
        ]), false)

    # Test based on "Logik für Informatiker", exercise 76.
    testsatisfiability("(*{x}(S(x, p) > L(x)) > G(p)) & ((~?{x}S(x,p)) & (~G(p)))", false)

    # Test based on "Logik für Informatiker", exercise 76.
    testsatisfiability("(*{x}(S(x, p) > L(x)) > G(p)) & ((~?{x}S(x,p)) & (G(p)))", true)

    # Test based on "Logik für Informatiker", exercise 80.
    testsatisfiability(CNF([
        Clause([Literal(false, Formula("P(u,u)"))]),
        Clause([Literal(true, Formula("P(f(x), x)"))])]), true)

    # Test based on "Logik für Informatiker", exercise 85.
    testsatisfiability(CNF([
        Clause([Literal(true, Formula("B(x)")), Literal(false, Formula("R(y,y)")), Literal(false, Formula("R(x,y)"))]),
        Clause([Literal(true, Formula("B(z)")), Literal(true, Formula("R(u,u)")), Literal(true, Formula("R(z,u)"))]),
        Clause([Literal(false, Formula("B(a)"))]),
        ]), false)

    # Test based on "Logik für Informatiker", exercise 85.
    testsatisfiability(CNF([
        Clause([Literal(true, Formula("B(x)")), Literal(false, Formula("R(y,y)")), Literal(false, Formula("R(x,y)"))]),
        Clause([Literal(true, Formula("B(z)")), Literal(true, Formula("R(u,u)")), Literal(true, Formula("R(z,u)"))]),
        Clause([Literal(true, Formula("B(a)"))]),
    ]), true, maxsearchdepth=3)
end


text:
## --------------  GradRule Logistic
# Grad prob is logistic in h/hGrad[college].

Base.@kwdef mutable struct GradRuleSwitchesLogistic <: GradRuleSwitches
    "Students can graduate at the END of this year or thereafter"
    tFirstGrad :: TimeInt = TimeInt(4)
    "Number of courses required for graduation (in DATA courses)"
    minNcForGrad :: ncInt = TimeInt(40)
    "Number of 2 year colleges (cannot graduate from those)"
    n2Year :: ncInt = ncInt(1)
    "Number of colleges"
    nColleges :: ncInt = ncInt(4)
    # Min and max grad prob when graduation is feasible
    gradProbMin :: Double = Double(0.05)
    gradProbMax :: Double = Double(0.95)
    # Grad prob depends on levels or gains
    useHcLevel :: Bool = true
    byCollege :: Bool = true

    # Logistic parameters. Slope is redundant.
    slope :: Double = Double(1.0)
    calSlope :: Bool = false
    lrShifter :: Double = Double(1.0)
    calLrShifter :: Bool = false
    twister :: Double = Double(1.0)
    calTwister :: Bool = false
end

Base.@kwdef mutable struct GradRuleSetLogistic <: GradRuleSet
    @grad_rule_set_common
    # objId :: ObjectId
    # The h0 term in `exp(-k * (h / h0))`. Higher values reduce grad prob.
    hGradV :: IncreasingVector{Double}
    # The slope is redundant b/c of h0
    slope :: Double = 1.0
    lrShifter :: Double = 1.0
    twister :: Double = 1.0
    switches :: GradRuleSwitchesLogistic
    # pvec :: ParamVector
end

## For one college
struct GradRuleLogistic <: GradRule
    # Can graduate after at least this many years
    tFirstGrad :: TimeInt
    # Grad prob depends on h/hGrad
    hGrad :: Double
    # Min no of courses attempted (cumulative)
    nMin :: ncInt
    # Generalized logistic
    glf :: GeneralizedLogistic
    # Grad prob depends on levels or gains
    useHcLevel :: Bool
end

## ----------- Logistic GradRule

# function make_test_gradrule(gs :: GradRuleSwitchesLogistic)
#     glf = GeneralizedLogistic();
#     return GradRuleLogistic(CollInt(4), Double(2.0), ncInt(20), glf,
#         use_hc_level(gs))
# end

function settings_table(gs :: GradRuleSwitchesLogistic)
    nc = n_colleges(gs);
    tFirstGrad = Int(t_first_grad(gs, nc));
    minNcToGrad = Int(min_courses_to_grad(gs, nc; modelUnits = false));
    if gs.calLrShifter
        lrShifterStr = "calibrated";
    else
        lrShifterStr = "fixed at $(gs.lrShifter)";
    end
    boundsStr = "$(gs.gradProbMin) to $(gs.gradProbMax)";
    return [
        "Graduation rule"  "Logistic";
        "Time to graduation"  "$tFirstGrad";
        "Min courses for graduation"  "$minNcToGrad";
        "Grad prob bounds"  boundsStr;
        "Shifter"  lrShifterStr
    ]
end

calibrate_lr_shifter!(gs :: GradRuleSwitchesLogistic) =
    gs.calLrShifter = true;

function fix_lr_shifter!(gs  :: GradRuleSwitchesLogistic,
    lrShifter :: Double)
    gs.calLrShifter = false;
    gs.lrShifter = lrShifter;
end

# There is no hMin property here. Ignore.
function set_hmin!(switches :: GradRuleSwitchesLogistic, hMin :: Double) end

"""
	$(SIGNATURES)

For all colleges.
Assumes college 1 is a 2 year college. Others are 4 year colleges.
"""
function make_grad_rule_set(objId :: ObjectId,
    grSwitches :: GradRuleSwitchesLogistic)

    nc = grSwitches.nColleges;
    hGradV = init_hgrad(objId, grSwitches);

    slope = grSwitches.slope;
    pSlope = Param(:slope, "GradRule slope", "slope",
        slope, slope, Double(0.1) * slope, Double(4.0) * slope, grSwitches.calSlope);

    # Lower bounds below 0.5 give very flat graduation rules (and grad rates by quality)
    twister = grSwitches.twister;
    pTwister = Param(:twister, "GradRule twister", "twister",
        twister, twister, Double(0.5) * twister, Double(2.0) * twister,
        grSwitches.calTwister);

    # LrShifter is redundant b/c we use h/hGrad in grad_prob
    lrShifter = grSwitches.lrShifter;
    lb = max(0.2, min(0.5 * lrShifter, 0.5));
    pLrShifter = Param(:lrShifter, "GradRule shifter", "lrShifter",
        lrShifter, lrShifter, Double(lb), Double(20.0) * lrShifter,
        grSwitches.calLrShifter);

    pvec = ParamVector(objId,  [pSlope, pTwister, pLrShifter]);
    return GradRuleSetLogistic(objId = objId,
        hGradV = hGradV, switches = grSwitches, pvec = pvec,
        slope = pSlope.value,  twister = pTwister.value,  lrShifter = pLrShifter.value)
end

"""
	$(SIGNATURES)

Make grad rule for one college
"""
function make_grad_rule(gs :: GradRuleSetLogistic, iCollege :: Integer)

    glf = GeneralizedLogistic(lrShifter = gs.lrShifter,  slope = gs.slope,
        twister = gs.twister,  lb = grad_prob_min(gs),  ub = grad_prob_max(gs));
    tFirstGrad = t_first_grad(gs, iCollege);
    return GradRuleLogistic(tFirstGrad, h_grad(gs, iCollege),
        min_courses_to_grad(gs, iCollege),  glf, use_hc_level(gs))
end

## Graduation probability at end of t
# One state
function grad_prob(g :: GradRuleLogistic, t :: Integer, h :: Double,
    n :: T1, h0 :: Double) where T1 <: Integer

    if can_graduate(g, t)
        if (n >= min_courses_to_grad(g, t))
            use_hc_level(g)  ?  (dh = h)  :  (dh = h - h0);
            gProb = logistic(g.glf,  dh / g.hGrad);
        else
            gProb = zero(Double);
        end
    else
        gProb = zero(Double);
    end
    return gProb
end

# --------------------


text:
# NOTE: Not using @enum because it craps out when displaying unknown values
const LINKTYPE_NULL = UInt32(0)
const LINKTYPE_ETHERNET = UInt32(1)

struct PcapHeader
    magic::UInt32
    version_major::UInt16
    version_minor::UInt16
    thiszone::Int32
    sigfigs::UInt32
    snaplen::UInt32
    linktype::UInt32
end

function Base.bswap(x::PcapHeader)
    PcapHeader(
        bswap(x.magic),
        bswap(x.version_major),
        bswap(x.version_minor),
        bswap(x.thiszone),
        bswap(x.sigfigs),
        bswap(x.snaplen),
        bswap(x.linktype))
end

function process_header(x::PcapHeader)
    if x.magic == 0xa1b2c3d4
        bswapped = false
        nanotime = false
    elseif x.magic == 0xd4c3b2a1
        bswapped = true
        nanotime = false
    elseif x.magic == 0xa1b23c4d
        bswapped = false
        nanotime = true
    elseif x.magic == 0x4d3cb2a1
        bswapped = true
        nanotime = true
    else
        throw(ArgumentError("Invalid pcap header"))
    end
    if bswapped
        x = bswap(x)
    end
    x, bswapped, nanotime
end


text:
"""
    module TiledFactorization

Tiled algorithms for factoring dense matrices.
"""
module TiledFactorization

using LinearAlgebra
using LinearAlgebra: BlasInt
using LoopVectorization
using TriangularSolve
using RecursiveFactorization
using Octavian

using DataFlowTasks
using DataFlowTasks: R,W,RW

function schur_complement!(C,A,B,tturbo::Val{T}) where {T}
    # RecursiveFactorization.schur_complement!(C,A,B,tturbo) // usually slower than Octavian
    if T
        Octavian.matmul!(C,A,B,-1,1)
    else
        Octavian.matmul_serial!(C,A,B,-1,1)
    end
end

const TILESIZE = Ref(256)

include("tiledmatrix.jl")
include("cholesky.jl")
include("lu.jl")

export
    PseudoTiledMatrix

end


text:
using Thebes

include(dirname(pathof(Thebes)) * "/../data/moreobjects.jl")

moreobjects = [boxcube, boxtorus, concave, cone, crossshape, cube, cuboctahedron, dodecahedron , geodesic, helix2,
icosahedron, icosidodecahedron, octahedron, octtorus, rhombicosidodecahedron,
rhombicuboctahedron, rhombitruncated_cubeoctahedron, rhombitruncated_icosidodecahedron,
snub_cube, snub_dodecahedron, sphere2, tet3d, tetrahedron, triangle, truncated_cube,
truncated_dodecahedron, truncated_icosahedron, truncated_octahedron, truncated_tetrahedron]

function main()
    Drawing(800, 800, "moreobjects.svg")
    origin()
    t = Table(6, 6, 100, 80)

    background("azure")
    setopacity(0.95)
    setline(0.1)

    for (pos, n) in t
        @layer begin
            translate(pos)
            carpet(10)
            object = make(moreobjects[mod1(n, length(moreobjects))])
            scaleby!(object, 10, 10, 10)
            rotateby!(object, 0.1, 0.1,  0.1)
            sortfaces!(object)
            pin(object)
        end
    end
    finish()
end

main()


text:
using Test, CUDA, SpinGlassPEPS

my_tests = []
if CUDA.functional() && CUDA.has_cutensor() && false
    CUDA.allowscalar(false)
    include("cuda/test_helpers.jl")
    push!(
        my_tests,
        "cuda/base.jl",
        "cuda/contractions.jl",
        "cuda/compressions.jl",
        "cuda/spectrum.jl"
    )
end

for my_test in my_tests
    include(my_test)
end


text:
#################################################################
### Just an example, there is not much to wrap around VideoIO ###
#################################################################

#import VideoIO, Images, ImageInTerminal
#
# fh = VideoIO.open("videofile")
# f = VideoIO.open(fh)
#
## or
#
# h=VideoIO.openvideo("/home/zgornel/projects/DATA/video/pontsync.avi")
#
# img = read(h) # reads stuff
#

# Examples from the High-Level interface:
#=
using Images
import ImageView
import VideoIO

io = VideoIO.open(video_file)
f = VideoIO.openvideo(io)

# As a shortcut for just video, you can upen the file directly
# with openvideo
#f = VideoIO.openvideo(video_file)

# Alternatively, you can open the camera with opencamera().
# The default device is "0" on Windows, "/dev/video0" on Linux,
# and "Integrated Camera" on OSX.  If using something other
# than the default, pass it in as the first parameter (as a string).
#f = VideoIO.opencamera()

# One can seek to an arbitrary position in the video
seek(f,2.5)  ## The second parameter is the time in seconds
img = read(f, Image)
canvas, _ = ImageView.view(img)

while !eof(f)
	read!(f, img)
	ImageView.view(canvas, img)
	#sleep(1/30)
end

=#


text:
"""
    Op

An MPI reduction operator, for use with [Reduce/Scan collective operations](@ref Reduce/Scan) to wrap binary operators. MPI.jl will perform this conversion automatically.

## Usage

    Op(op, T=Any; iscommutative=false)

Wrap the Julia reduction function `op` for arguments of type `T`. `op` is assumed to be
associative, and if `iscommutative` is true, assumed to be commutative as well.

## See also

- [`Reduce!`](@ref)/[`Reduce`](@ref)
- [`Allreduce!`](@ref)/[`Allreduce`](@ref)
- [`Scan!`](@ref)/[`Scan`](@ref)
- [`Exscan!`](@ref)/[`Exscan`](@ref)
"""
mutable struct Op
    val::MPI_Op
    fptr
    Op(val::MPI_Op, fptr) = new(val, fptr)
end
Base.:(==)(a::Op, b::Op) = a.val == b.val
Base.cconvert(::Type{MPI_Op}, op::Op) = op
Base.unsafe_convert(::Type{MPI_Op}, op::Op) = op.val
Base.unsafe_convert(::Type{Ptr{MPI_Op}}, op::Op) = convert(Ptr{MPI_Op}, pointer_from_objref(op))

const OP_NULL = Op(Consts.MPI_OP_NULL[], nothing)
const BAND    = Op(Consts.MPI_BAND[], nothing)
const BOR     = Op(Consts.MPI_BOR[], nothing)
const BXOR    = Op(Consts.MPI_BXOR[], nothing)
const LAND    = Op(Consts.MPI_LAND[], nothing)
const LOR     = Op(Consts.MPI_LOR[], nothing)
const LXOR    = Op(Consts.MPI_LXOR[], nothing)
const MAX     = Op(Consts.MPI_MAX[], nothing)
const MIN     = Op(Consts.MPI_MIN[], nothing)
const PROD    = Op(Consts.MPI_PROD[], nothing)
const REPLACE = Op(Consts.MPI_REPLACE[], nothing)
const SUM     = Op(Consts.MPI_SUM[], nothing)
const NO_OP   = Op(Consts.MPI_NO_OP[], nothing)
add_load_time_hook!(() -> OP_NULL.val = Consts.MPI_OP_NULL[])
add_load_time_hook!(() -> BAND.val    = Consts.MPI_BAND[]   )
add_load_time_hook!(() -> BOR.val     = Consts.MPI_BOR[]    )
add_load_time_hook!(() -> BXOR.val    = Consts.MPI_BXOR[]   )
add_load_time_hook!(() -> LAND.val    = Consts.MPI_LAND[]   )
add_load_time_hook!(() -> LOR.val     = Consts.MPI_LOR[]    )
add_load_time_hook!(() -> LXOR.val    = Consts.MPI_LXOR[]   )
add_load_time_hook!(() -> MAX.val     = Consts.MPI_MAX[]    )
add_load_time_hook!(() -> MIN.val     = Consts.MPI_MIN[]    )
add_load_time_hook!(() -> PROD.val    = Consts.MPI_PROD[]   )
add_load_time_hook!(() -> REPLACE.val = Consts.MPI_REPLACE[])
add_load_time_hook!(() -> SUM.val     = Consts.MPI_SUM[]    )
add_load_time_hook!(() -> NO_OP.val   = Consts.MPI_NO_OP[]  )

Op(::typeof(min), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint}} = MIN
Op(::typeof(max), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint}} = MAX
Op(::typeof(+), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint,MPIComplex}} = SUM
Op(::typeof(*), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint,MPIComplex}} = PROD
Op(::typeof(&), ::Type{T}; iscommutative=true) where {T<:MPIInteger} = BAND
Op(::typeof(|), ::Type{T}; iscommutative=true) where {T<:MPIInteger} = BOR
Op(::typeof(⊻), ::Type{T}; iscommutative=true) where {T<:MPIInteger} = BXOR

function free(op::Op)
    if op != OP_NULL && !Finalized()
        # int MPI_Op_free(MPI_Op *op)
        @mpichk ccall((:MPI_Op_free, libmpi), Cint, (Ptr{MPI_Op},), op)
    end
    op.fptr = nothing
    return nothing
end

struct OpWrapper{F,T}
    f::F
end

function (w::OpWrapper{F,T})(_a::Ptr{Cvoid}, _b::Ptr{Cvoid}, _len::Ptr{Cint}, t::Ptr{MPI_Datatype}) where {F,T}
    len = unsafe_load(_len)
    @assert isconcretetype(T)
    a = Ptr{T}(_a)
    b = Ptr{T}(_b)
    for i = 1:len
        unsafe_store!(b, w.f(unsafe_load(a,i), unsafe_load(b,i)), i)
    end
    return nothing
end

function Op(f, T=Any; iscommutative=false)
    @static if MPI_LIBRARY == MicrosoftMPI && Sys.WORD_SIZE == 32
        error("User-defined reduction operators are not supported on 32-bit Windows.\nSee https://github.com/JuliaParallel/MPI.jl/issues/246 for more details.")
    elseif Sys.ARCH ∈ (:aarch64, :ppc64le, :powerpc64le) || startswith(lowercase(String(Sys.ARCH)), "arm")
        error("User-defined reduction operators are currently not supported on non-Intel architectures.\nSee https://github.com/JuliaParallel/MPI.jl/issues/404 for more details.")
    end
    w = OpWrapper{typeof(f),T}(f)
    fptr = @cfunction($w, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint}, Ptr{MPI_Datatype}))

    op = Op(OP_NULL.val, fptr)
    # int MPI_Op_create(MPI_User_function* user_fn, int commute, MPI_Op* op)
    @mpichk ccall((:MPI_Op_create, libmpi), Cint,
                  (Ptr{Cvoid}, Cint, Ptr{MPI_Op}),
                  fptr, iscommutative, op)

    finalizer(free, op)
    return op
end