text:
#=

Predictor

## Introduction

In an ideal world the predictor only uses an order (2,1) Padé approximant to generate
an initial guess for the Newton corrector.
The choice of Padé approximants as predictors is motived by the results in [1].
In particular, the ability to provide a trust region of the prediction method to
significantly decrease the likelihood of path jumping.
additionally, sice Padé approximants are rational functions, these perform very well
for diverging paths and close to singularities.

Now, for an Padé approximant of order (2,1) the first 3 derivatives are needed.
The most robust way to achieve this it to use automatic differentiation.
However, due to significant compile times, we also have to consider the case
that the Padé approximants are computed by numerical differentiation using the scheme
derived in [2]. The numerical differentiation scheme only uses O(h^2) approximations
since we are mostly constrained by the available precision.

However, the smaller the trust region is (i.e. the closer we are to a singularity) the
smaller we need to take the discretization in the numerical differentiation.
Thus, the error in the numerical differentiation is inversely corrolated to the trust
region.
Additionally, the third derivative needs to be sufficiently accurate. Otherwise
the predictor doesn't generate sufficiently accurate intial guesses resulting
in very small step sizes.
Therefore it can necessary to drop down to an Padé approximant of order (1,1).
Usually the arithmetic is sufficient to obtain for the second derivative a sufficiently
accurate approximation.

If also the second derivative is not available we use cubic hermite interpolation.
For diverging paths this works much much better than (explicit) Runge-Kutta methods.
The reason that we do not use expl. RK methods is that, then when numerical differentiation
has to give up, we are in a regime where the ODE of the path is *stiff*. Thus, expl.
RK just don't yield satisfying results.

[1]: A Robust Numerical Path Tracking Algorithm for Polynomial Homotopy Continuation,
 Telen, Van Barel, Verschelde https://arxiv.org/abs/1909.04984

[2]:Mackens, Wolfgang. "Numerical differentiation of implicitly defined space curves."
   Computing 41.3 (1989): 237-260.
=#

module PredictionMethod
@enum methods begin
    Pade21
    Pade11
    Hermite
    Euler
end
end
using .PredictionMethod: PredictionMethod

## Type dispatch on automatic differentiation or numerical differentiation
struct AD{N} end
function AD(N::Int)
    if !(0 ≤ N ≤ 3)
        throw(ArgumentError("`automatic_differentiation` has to be between 0 and 4."))
    end
    AD{N}()
end

"""
    Predict(H::AbstractHomotopy, ::AD{N}) where {N}

Construct a predictor. The predicor uses a Padé approximant of order (2,1) to produce
an initial guess for the Newton corrector.
For this, the predictor tries to compute the local derivatives up to order 3.
This is done using automatic differentiation for the derivatives up to order `3`.
Otherwise numerical differentiation is used.
"""
Base.@kwdef mutable struct Predictor{T<:AD}
    AD::T

    method::PredictionMethod.methods
    order::Int
    trust_region::Float64
    local_error::Float64
    cond_H_ẋ::Float64

    tx⁰::TaylorVector{1,ComplexF64}
    tx¹::TaylorVector{2,ComplexF64}
    tx²::TaylorVector{3,ComplexF64}
    tx³::TaylorVector{4,ComplexF64}
    tx⁴::TaylorVector{5,ComplexF64}
    t::ComplexF64
    tx_norm::Vector{Float64}

    # finite diff
    xtemp::Vector{ComplexF64}
    u::Vector{ComplexF64}
    u₁::Vector{ComplexF64}
    u₂::Vector{ComplexF64}
    # for hermite interpolation
    prev_tx¹::TaylorVector{2,ComplexF64}
    prev_t::ComplexF64
end

function Predictor(H::AbstractHomotopy, autodiff::AD)
    m, n = size(H)
    tx⁴ = TaylorVector{5}(ComplexF64, n)
    tx⁰ = TaylorVector{1}(tx⁴)
    tx¹ = TaylorVector{2}(tx⁴)
    tx² = TaylorVector{3}(tx⁴)
    tx³ = TaylorVector{4}(tx⁴)
    Predictor(
        AD = autodiff,
        method = PredictionMethod.Pade21,
        order = 4,
        trust_region = Inf,
        local_error = Inf,
        cond_H_ẋ = Inf,
        tx⁰ = tx⁰,
        tx¹ = tx¹,
        tx² = tx²,
        tx³ = tx³,
        tx⁴ = tx⁴,
        t = complex(NaN),
        tx_norm = zeros(4),
        xtemp = zeros(ComplexF64, n),
        u = zeros(ComplexF64, m),
        u₁ = zeros(ComplexF64, m),
        u₂ = zeros(ComplexF64, m),
        prev_tx¹ = TaylorVector{2}(ComplexF64, n),
        prev_t = complex(NaN),
    )
end

function init!(predictor::Predictor)
    predictor.cond_H_ẋ = 1.0
    predictor.t = predictor.prev_t = NaN
    predictor.trust_region = predictor.local_error = NaN
    predictor
end

order(predictor::Predictor) = predictor.order
local_error(predictor::Predictor) = predictor.local_error
trust_region(predictor::Predictor) = predictor.trust_region

###################
### FINITE DIFF ###
###################

function g!(u, H, x::Vector, t, h, xtemp)
    evaluate!(u, H, x, t + h)
end
function g!(u, H, tx::TaylorVector{2}, t, h, xtemp)
    @inbounds for i = 1:length(tx)
        x, x¹ = tx[i]
        xtemp[i] = muladd(h, x¹, x)
    end
    evaluate!(u, H, xtemp, t + h)
end
function g!(u, H, tx::TaylorVector{3}, t, h, xtemp)
    @inbounds for i = 1:length(tx)
        x, x¹, x² = tx[i]
        xtemp[i] = muladd(h, muladd(h, x², x¹), x)
    end
    evaluate!(u, H, xtemp, t + h)
end

function g!(u, H, tx::TaylorVector{4}, t, h, xtemp)
    @inbounds for i = 1:length(tx)
        x, x¹, x², x³ = tx[i]
        xtemp[i] = muladd(h, muladd(h, muladd(h, x³, x²), x¹), x)
    end
    evaluate!(u, H, xtemp, t + h)
end

function finite_diff!(u, predictor::Predictor, H::AbstractHomotopy, x, t, h; order::Int)
    @unpack u₁, u₂, xtemp = predictor

    g!(u₁, H, x, t, h, xtemp)
    g!(u₂, H, x, t, -h, xtemp)

    hk = h^order
    if iseven(order)
        u .= 0.5 .* (u₁ .+ u₂) ./ hk
    else
        u .= 0.5 .* (u₁ .- u₂) ./ hk
    end
end

function finite_diff_taylor!(
    u,
    ::Val{N},
    predictor::Predictor,
    H::AbstractHomotopy,
    x,
    t;
    prev_λ::Float64 = NaN,
) where {N}
    # Use not machine precision to account for some error in the evaluation
    ε = 1.4210854715202004e-14
    # λ should be a scaling factor such that the derivatives of x(t/λ) are of of the same
    # order of magnitude. We can either compute this by balancing the derivatives or
    # re-using the trust region size of the previous step
    if N == 1
        λ = isfinite(prev_λ) ? prev_λ : 1.0
    elseif N == 2
        λ = isfinite(prev_λ) ? prev_λ : predictor.tx_norm[1] / predictor.tx_norm[2]
    elseif N == 3
        λ = predictor.tx_norm[2] / predictor.tx_norm[3]
    end
    # truncation err = (1/λ)^2*h^2
    # round-off err   = ε/h^N
    # -->
    # Compute optimal ĥ by
    #      round-off err = trunc err
    # <=>  ελ^-N/ĥ^3 = λ^(-N-2)*ĥ^2
    # <=>  (ε λ^2)^(1/(N+2)) = ĥ
    # To be a little bit pessimistic, reduce λ by a quarter
    λ *= 0.25

    h = (ε * λ^2)^(1 / (N + 2))
    # h should be at most half λ
    h = min(h, 0.5λ)
    # Check that truncation and round-off error are both acceptable
    trunc_err = h^2 / λ^2
    rnd_err = ε / h^N
    err = trunc_err + rnd_err
    if trunc_err + rnd_err > 0.01 || !isfinite(h)
        return false, err
    end

    finite_diff!(u, predictor, H, x, t, h; order = N)

    return true, err^2
end

## Default handling ignores incremental
taylor!(u, v::Val, H::AbstractHomotopy, tx, t, incremental::Bool) = taylor!(u, v, H, tx, t)
@generated function taylor!(u, v::Val{M}, predictor::Predictor{AD{N}}, H, tx, t) where {M,N}
    if M ≤ N
        quote
            taylor!(u, v, H, tx, t, true)
            true, eps()
        end
    else
        quote
            finite_diff_taylor!(u, v, predictor, H, tx, t; prev_λ = predictor.trust_region)
        end
    end
end

"""
    update!(
        P::Predictor,
        H::AbstractHomotopy,
        x,
        t,
        jacobian::Jacobian,
        norm::AbstractNorm,
        x̂ = nothing,
    )

Upadte the predictor with the new solution `(x,t)` of `H(x,t) = 0`.
This computes new local derivatives and chooses an appriopriate prediction method.
"""
function update!(
    predictor::Predictor,
    H::AbstractHomotopy,
    x,
    t,
    J::Jacobian,
    norm::AbstractNorm,
    x̂ = nothing, # the predicted value
)
    @unpack u, tx⁰, tx¹, tx², tx³, tx⁴, xtemp, tx_norm = predictor
    x⁰, x¹, x², x³, x⁴ = vectors(tx⁴)

    @inbounds for i = 1:length(xtemp)
        predictor.prev_tx¹[i, 1] = predictor.tx¹[i, 1]
        predictor.prev_tx¹[i, 2] = predictor.tx¹[i, 2]
    end
    predictor.prev_t, predictor.t = predictor.t, t
    if isnothing(x̂)
        predictor.local_error = NaN
    else
        Δs = fast_abs(t - predictor.prev_t)
        predictor.local_error = norm(x̂, x) / Δs^predictor.order
    end

    x⁰ .= x
    predictor.t = t
    tx_norm[1] = norm(x)

    # Compute ẋ
    trust, err = taylor!(u, Val(1), predictor, H, x, t)
    u .= .-u
    LA.ldiv!(xtemp, J, u)
    # Check error made in the linear algebra
    δ = fixed_precision_iterative_refinement!(xtemp, workspace(J), u, norm)
    predictor.cond_H_ẋ = δ / eps()
    tol_δ₁ = max(1e-10, err)
    δ > tol_δ₁ && iterative_refinement!(xtemp, J, u, InfNorm(); tol = tol_δ₁, max_iters = 5)
    tx_norm[2] = norm(xtemp)
    x¹ .= xtemp

    trust, err = taylor!(u, Val(2), predictor, H, tx¹, t)
    if !trust
        # Use cubic hermite
        @label use_hermite
        if isfinite(predictor.prev_t)
            predictor.method = PredictionMethod.Hermite
            predictor.order = 3
            predictor.trust_region = tx_norm[1] / tx_norm[2]
            if isnan(predictor.local_error)
                predictor.local_error = (tx_norm[2] / tx_norm[1])^3
            end
            return predictor
        else
            predictor.method = PredictionMethod.Euler
            predictor.order = 2
            predictor.trust_region = tx_norm[1] / tx_norm[2]
            if isnan(predictor.local_error)
                predictor.local_error = (tx_norm[2] / tx_norm[1])^2
            end
            return predictor
        end
    end
    u .= .-u
    LA.ldiv!(xtemp, J, u)
    tol_δ₂ = max(1e-8, err)
    δ > tol_δ₂ && iterative_refinement!(xtemp, J, u, norm; tol = tol_δ₂, max_iters = 4)
    tx_norm[3] = norm(xtemp)
    x² .= xtemp

    trust, err = taylor!(u, Val(3), predictor, H, tx², t)
    if !trust
        @label use_pade11
        # only trust tx_norm[2] / tx_norm[3] if tx_norm[1] / tx_norm[2] is of the same order
        if 0.1 ≤ (tx_norm[1] / tx_norm[2]) / (tx_norm[2] / tx_norm[3]) ≤ 10
            # Use Padé (1,1) predictor
            predictor.method = PredictionMethod.Pade11
            predictor.order = 3
            predictor.trust_region = tx_norm[2] / tx_norm[3]
            if isnan(predictor.local_error)
                predictor.local_error = (tx_norm[3] / tx_norm[2])^3
            end
            return predictor
        else
            @goto use_hermite
        end
    end

    u .= .-u
    LA.ldiv!(xtemp, J, u)
    tol_δ₃ = max(1e-6, err)
    δ > tol_δ₃ && iterative_refinement!(xtemp, J, u, norm; tol = tol_δ₃, max_iters = 3)
    tx_norm[4] = norm(xtemp)
    x³ .= xtemp

    if 0.1 ≤ (tx_norm[2] / tx_norm[3]) / (tx_norm[3] / tx_norm[4]) ≤ 10
        # Use Padé (2,1) predictor
        predictor.method = PredictionMethod.Pade21
        predictor.order = 4
        predictor.trust_region = tx_norm[3] / tx_norm[4]
        if isnan(predictor.local_error)
            predictor.local_error = (tx_norm[4] / tx_norm[3])^4
        end
    else
        @goto use_pade11
    end

    predictor
end

predict!(x̂, pred::Predictor, H, x, t, Δt) = predict!(x̂, pred, pred.method, H, x, t, Δt)
function predict!(x̂, pred::Predictor, method::PredictionMethod.methods, H, x, t, Δt)
    if method == PredictionMethod.Pade21
        λ = pred.trust_region
        λ² = λ^2
        λ³ = λ^3
        tol = 1e-8
        for (i, (x, x¹, x², x³)) in enumerate(pred.tx³)
            c, c¹, c², c³ = fast_abs.((x, x¹, x², x³))
            # check if only taylor series is used
            τ = tol * √(c^2 + (c¹ * λ)^2 + (c² * λ²)^2 + (c³ * λ³)^2)
            if c³ * λ³ ≤ τ || c² * λ² ≤ τ
                x̂[i] = x + Δt * (x¹ + Δt * x²)
            else
                δᵢ = 1 - Δt * x³ / x²
                x̂[i] = x + Δt * (x¹ + Δt * x² / δᵢ)
            end
        end
    elseif method == PredictionMethod.Pade11
        λ = pred.trust_region
        λ² = λ^2
        tol = 1e-8
        for (i, (x, x¹, x²)) in enumerate(pred.tx²)
            c, c¹, c² = fast_abs.((x, x¹, x²))
            # check if only taylor series is used
            τ = tol * √(c^2 + (c¹ * λ)^2 + (c² * λ²)^2)
            if c¹ * λ ≤ τ || c² * λ² ≤ τ
                x̂[i] = x + Δt * x¹
            else
                δᵢ = 1 - Δt * x² / x¹
                x̂[i] = x + Δt * x¹ / δᵢ
            end
        end
    elseif method == PredictionMethod.Hermite
        cubic_hermite!(x̂, pred.prev_tx¹, pred.prev_t, pred.tx¹, t, t + Δt)
    elseif method == PredictionMethod.Euler
        for (i, (x, x¹)) in enumerate(pred.tx¹)
            x̂[i] = x + Δt * x¹
        end
    end

    x̂
end

function cubic_hermite!(x̂, tx¹₀, t₀, tx¹₁, t₁, t)
    if isreal(t₀) && isreal(t₁) && isreal(t)
        _cubic_hermite!(x̂, tx¹₀, real(t₀), tx¹₁, real(t₁), real(t))
    else
        _cubic_hermite!(x̂, tx¹₀, t₀, tx¹₁, t₁, t)
    end
end
@inline function _cubic_hermite!(x̂, tx¹₀, t₀, tx¹₁, t₁, t)
    s = (t - t₀) / (t₁ - t₀)
    h₀₀ = (1 + 2s) * (1 - s)^2
    h₁₀ = (t - t₀) * (1 - s)^2
    h₀₁ = s^2 * (3 - 2s)
    h₁₁ = (t - t₀) * s * (s - 1)
    @inbounds for i in eachindex(x̂)
        x̂[i] = h₀₀ * tx¹₀[i, 1] + h₁₀ * tx¹₀[i, 2] + h₀₁ * tx¹₁[i, 1] + h₁₁ * tx¹₁[i, 2]
    end
    x̂
end


text:
"""
    target(target,xlength,ylength,zlength)

Creates a cuboid target of dimensions (xlength, ylength and zlength) at the target position
for shotline intersection computation
"""
function target(t::AbstractTarget,xlength::Float64,ylength::Float64,zlength::Float64)
    target = Cuboid(xlength, ylength, zlength)
    target = OpticSim.transform(target,OpticSim.Geometry.translation(t.position[1], t.position[2], t.position[3]))
    return target
end

"""
    target(size)

Creates a cuboid target of dimensions specified by the vector size
"""
function target(size::Vector{Float64})
    target = Cuboid(size[1], size[2], size[3])
    #target = OpticSim.transform(target,OpticSim.Geometry.translation(t.position[1], t.position[2], t.position[3]))
    return target
end

function targetPosition(target::CSGGenerator{Float64},position::Vector{Float64})
    #target = Cuboid(size[1], size[2], size[3])
    target = OpticSim.transform(target,OpticSim.Geometry.translation(position[1], position[2], position[3]))
    return target
end

ext: jl


text:
"""
Geroch-Held-Penrose (GHP) operators; see

   Geroch, Held, and Penrose,
   Journal of Mathematical Physics 14, 874 (1973)

For the formulas in the coordinates we use, see
(in particular Sec V.C)

   Ripley, Loutrel, Giorgi, and Pretorius,
   Phys.Rev.D 103 (2021) 104018, arXiv:2010.00162,

"""
module GHP

include("Fields.jl")
include("Sphere.jl")
using .Fields: Field
using .Sphere: swal_raising_matrix, swal_lowering_matrix, angular_matrix_mult!

export GHP_ops, Initialize_GHP_ops, set_edth!, set_edth_prime!, set_thorn!, set_thorn_prime!

struct GHP_ops

   pre_edth_DT::Array{ComplexF64,2}
   pre_edth_raised::Array{ComplexF64,2}
   pre_edth::Array{ComplexF64,2}
   pre_edth_prime_DT::Array{ComplexF64,2}
   pre_edth_prime_lowered::Array{ComplexF64,2}
   pre_edth_prime::Array{ComplexF64,2}
   pre_thorn_DT::Array{ComplexF64,2}
   pre_thorn_DR::Array{ComplexF64,2}
   pre_thorn::Array{ComplexF64,2}
   pre_thorn_prime_DT::Array{ComplexF64,2}
   pre_thorn_prime_DR::Array{ComplexF64,2}

   raise::Array{ComplexF64,3}
   lower::Array{ComplexF64,3}

   smap::Dict{Int64,Int64}

   function GHP_ops(;
         Rvals::Vector{Float64},
         Cvals::Vector{Float64},
         Svals::Vector{Float64},
         mval::Int64,
         bhm::Float64,
         bhs::Float64,
         cl::Float64)

      nx = length(Rvals)
      ny = length(Cvals)

      spins = [-2,-1,0,1,2]

      pre_edth_DT            = zeros(ComplexF64,nx,ny)
      pre_edth_raised        = zeros(ComplexF64,nx,ny)
      pre_edth               = zeros(ComplexF64,nx,ny)
      pre_edth_prime_DT      = zeros(ComplexF64,nx,ny)
      pre_edth_prime_lowered = zeros(ComplexF64,nx,ny)
      pre_edth_prime         = zeros(ComplexF64,nx,ny)
      pre_thorn_DT           = zeros(ComplexF64,nx,ny)
      pre_thorn_DR           = zeros(ComplexF64,nx,ny)
      pre_thorn              = zeros(ComplexF64,nx,ny)
      pre_thorn_prime_DT     = zeros(ComplexF64,nx,ny)
      pre_thorn_prime_DR     = zeros(ComplexF64,nx,ny)

      raise = zeros(ComplexF64,ny,ny,length(spins))
      lower = zeros(ComplexF64,ny,ny,length(spins))

      smap = Dict{Int64,Int64}()

      for j=1:ny
         cy = Cvals[j]
         sy = Svals[j]
         for i=1:nx
            R = Rvals[i]
            ##=============
            ## divided by R
            ##=============
            pre_edth_DT[i,j] = (1.0/sqrt(2.0))*(
                                    1.0/((cl^2)
                                    -
                                    im*bhs*R*cy)
                                )*(-im*bhs*sy)

            pre_edth_raised[i,j] = (1.0/sqrt(2.0))*(
                                    1.0/((cl^2)
                                    -
                                    im*bhs*R*cy)
                                   )

            pre_edth[i,j] = (
                  (im*bhs*R*sy/sqrt(2.0))
                  /
                  ((im*(cl^2) + bhs*R*cy)^2)
                 )

            ##=============
            ## divided by R
            ##=============
            pre_edth_prime_DT[i,j] = (1.0/sqrt(2.0))*(
                                       1.0/((cl^2)
                                       +
                                       im*bhs*R*cy)
                                      )*im*bhs*sy

            pre_edth_prime_lowered[i,j] = (1.0/sqrt(2.0))*(
                                             1.0/((cl^2)
                                             +
                                             im*bhs*R*cy)
                                            )

            pre_edth_prime[i,j] = (
                  (im*bhs*R*sy/sqrt(2.0))
                  /
                  ((cl^2 + im*bhs*R*cy)^2)
                 )

            ##=============
            ## divided by R
            ##=============
            pre_thorn_DT[i,j] = (
               (1.0/((cl^4)+((bhs*R*cy)^2)))
               *
               R*2.0*bhm*(2.0*bhm-((bhs/cl)^2)*R)
              )

            ##=============
            ## divided by R
            ##=============
            pre_thorn_DR[i,j] = (
               (1.0/((cl^4)+((bhs*R*cy)^2)))*(
               -  0.5*((cl^2)-(2.0*bhm*R) + ((bhs*R/cl)^2))
               )
              )

            ##=================
            ## NOT divided by R
            ##=================
            pre_thorn[i,j] = (
               (1.0/((cl^4)+((bhs*R*cy)^2)))*R*(im*bhs)
              )

            pre_thorn_prime_DT[i,j] = (
               (2.0 + (4.0*bhm*R/(cl^2)))
              )

            ##=============
            ## divided by R
            ##=============
            pre_thorn_prime_DR[i,j] = ((1.0/cl)^2)*R
         end
      end

      for (i,s) in enumerate(spins)
         tmp_raise = swal_raising_matrix( ny,s,mval)
         tmp_lower = swal_lowering_matrix(ny,s,mval)

         raise[:,:,i] .= tmp_raise
         lower[:,:,i] .= tmp_lower
      end

      for (i,s) in enumerate(spins)
         smap[s]=i
      end

      return new(
         pre_edth_DT,
         pre_edth_raised,
         pre_edth,
         pre_edth_prime_DT,
         pre_edth_prime_lowered,
         pre_edth_prime,
         pre_thorn_DT,
         pre_thorn_DR,
         pre_thorn,
         pre_thorn_prime_DT,
         pre_thorn_prime_DR,
         raise,
         lower,
         smap
      )
   end
end

function Initialize_GHP_ops(;
      Rvals::Vector{Float64},
      Cvals::Vector{Float64},
      Svals::Vector{Float64},
      Mvals::Vector{Int64},
      bhm::Float64,
      bhs::Float64,
      cl::Float64)
   return Dict([
      (mv,GHP_ops(Rvals=Rvals,Cvals=Cvals,Svals=Svals,mval=mv,bhm=bhm,bhs=bhs,cl=cl))
      for mv in Mvals
     ])
end

function set_edth!(;
      edth ::Array{ComplexF64,2},
      spin ::Int64,
      boost::Int64,
      m_ang::Int64,
      f     ::Array{ComplexF64,2},
      DT    ::Array{ComplexF64,2},
      raised::Array{ComplexF64,2},
      Op
     )
   nx, ny = size(edth)
   p = (spin+boost)

   angular_matrix_mult!(raised,f,view(Op.raise,:,:,Op.smap[spin]))

   for j=1:ny
      for i=1:nx
         edth[i,j] = (
            Op.pre_edth_DT[i,j]     *DT[i,j]
            +
            Op.pre_edth_raised[i,j] *raised[i,j]
            +
            p*Op.pre_edth[i,j]*f[i,j]
        )
      end
   end
   return nothing
end

function set_edth_prime!(;
      edth_prime::Array{ComplexF64,2},
      spin ::Int64,
      boost::Int64,
      m_ang::Int64,
      f      ::Array{ComplexF64,2},
      DT     ::Array{ComplexF64,2},
      lowered::Array{ComplexF64,2},
      Op
     )
   nx, ny = size(edth)
   q = (-spin+boost)

   angular_matrix_mult!(lowered,f,view(Op.lower,:,:,Op.smap[spin]))

   for j=1:ny
      for i=1:nx
         edth_prime[i,j] = (
            Op.pre_edth_prime_DT[i,j]     *DT[i,j]
            +
            Op.pre_edth_prime_lowered[i,j]*lowered[i,j]
            +
            q*Op.pre_edth_prime[i,j]*f[i,j]
           )
      end
   end
   return nothing
end

function set_thorn!(;
      thorn::Array{ComplexF64,2},
      spin ::Int64,
      boost::Int64,
      m_ang::Int64,
      falloff::Int64,
      f    ::Array{ComplexF64,2},
      DT   ::Array{ComplexF64,2},
      DR   ::Array{ComplexF64,2},
      ep_0 ::Array{ComplexF64,2},
      R    ::Vector{Float64},
      Op
     )
   nx, ny = size(thorn)
   p = ( spin+boost)
   q = (-spin+boost)

   for j=1:ny
      for i=1:nx
         thorn[i,j] = (
            Op.pre_thorn_DT[i,j]*DT[i,j]
            +
            Op.pre_thorn_DR[i,j]*(R[i]*DR[i,j] + falloff*f[i,j])
            +
            m_ang*Op.pre_thorn[i,j]*f[i,j]
            -
            R[i]*(p*ep_0[i,j] + q*conj(ep_0[i,j]))*f[i,j]
           )
      end
   end
   return nothing
end

##===========================================================================
## no rescaling in R for thorn prime
##===========================================================================
function set_thorn_prime!(;
      thorn_prime::Array{ComplexF64,2},
      falloff::Int64,
      f ::Array{ComplexF64,2},
      DT::Array{ComplexF64,2},
      DR::Array{ComplexF64,2},
      R ::Vector{Float64},
      Op
     )
   for j=1:ny
      for i=1:nx
         thorn_prime_arr[i,j] = (
            Op.pre_thorn_prime_DT[i,j]*DT[i,j]
            +
            Op.pre_thorn_prime_DR[i,j]*(R[i]*DR[i,j] + falloff*f[i,j])
           )
      end
   end
   return nothing
end

end


text:
  mutable struct TestMesh <: AbstractDGMesh{Float64}
    pertNeighborEls::Array{Int, 2}
    dofs::Array{Int, 3}
    neighbor_nums::Array{Int, 2}
    numDof::Int
    numNodesPerElement::Int
    numDofPerNode::Int
    numEl::Int
    coloringDistance::Int
  end

  # type for testing isFieldDefined
  mutable struct TestType
    a::Int
    b::Array{Float64, 1}
    function TestType()
      return new()
    end
  end

function normLinf(A::AbstractArray)

  return maximum(abs.(A))
end

@testset "--- Testing misc.jl ---" begin

  # test mat-vec
  A = rand(3,3)
  x = rand(3)
  b = A*x
  b2 = zeros(3)
  smallmatvec!(A, x, b2)
  b3 = smallmatvec(A, x)

  @test isapprox( b, b2) atol=1e-14
  @test isapprox( b, b3) atol=1e-14

  A = rand(4, 3)
  x = rand(3)
  b = A*x
  b2 = smallmatvec(A, x)

  @test isapprox( b2, b) atol=1e-14

  # test reverse mode
  A = rand(4,3)

  x = rand(3)
  x2 = zeros(Complex128, 3)
  copy!(x2, x)

  b = zeros(4)
  b2 = zeros(Complex128, 4)

  jac = zeros(4, 3)
  jac2 = zeros(4, 3)
  h = 1e-20
  pert = Complex128(0, h)
  for i=1:3
    x2[i] += pert
    fill!(b2, 0)
    smallmatvec!(A, x2, b2)

    for j=1:4
      jac[j, i] = imag(b2[j])/h
    end

    x2[i] -= pert
  end

  for j=1:4
    b[j] = 1  # actuall b_bar
    fill!(x, 0)  # actuall x_bar
    smallmatvec_revv!(A, x, b)

    for i=1:3
      jac2[j, i] = x[i]
    end

    b[j] = 0
  end

  @test isapprox( normLinf(jac - jac2), 0.0) atol=1e-13

  # test reinterpert-BLAS interaction

  println("testing BLAS fallbacks")
  for d=2:20  # size of arrays
    A = rand(Complex128, d, d, 1)
    Av = unsafe_aview(A, :, :, 1)

    x = rand(d, 1)
    xv = unsafe_aview(x, :, 1)
#    xv = x[:, 1]

    b = zeros(Complex128, d, 1)
    bv = unsafe_aview(b, :, 1)
#    bv = b[:, 1]

    # test unsafe_aview
    smallmatvec!(Av, xv, bv)

    b2 = A[:, :, 1]*x[:, 1]

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    # test safe view
    Av2 = aview(A, :, :, 1)
    xv2 = aview(x, :, 1)
    bv2 = aview(b, :, 1)

    smallmatvec!(Av2, xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    # check mixed argument types
    b3 = zeros(b2)
    smallmatvec!(Av2, xv2, b3)

    @test isapprox( normLinf(b3 - b2), 0.0) atol=1e-14

    smallmatvec!(A[:, :, 1], xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    smallmatvec!(ROView(Av), xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatvec!(ROView(Av2), xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatvec!(ROView(Av2), ROView(xv), bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatvec!(ROView(Av2), xv, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    # test transposed mat-vec
    b2 = A[:, :, 1].'*x[:, 1]

    smallmatTvec!(Av, xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatTvec!(Av2, xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

  end

  # test mat-mat
  for d=2:20  # size of arrays
    A = rand(Complex128, d, d, 1)
    Av = unsafe_aview(A, :, :, 1)

    x = rand(d, d, 1)
    xv = unsafe_aview(x, :, :, 1)
#    xv = x[:, 1]

    b = zeros(Complex128, d, d, 1)
    bv = unsafe_aview(b, :, :, 1)
#    bv = b[:, 1]

    b2 = A[:, :, 1]*x[:, :, 1]

    smallmatmat!(Av, xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    # test safe aview
    Av2 = aview(A, :, :, 1)
    xv2 = aview(x, :, :, 1)
    bv2 = aview(b, :, :, 1)

    smallmatmat!(Av2, xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    b3 = zeros(b2)

    smallmatmat!(Av2, xv2, b3)

    @test isapprox( normLinf(b3 - b2), 0.0) atol=1e-14

    smallmatmat!(A[:, :, 1], xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    smallmatmat!(ROView(Av2), xv, bv)

    @test isapprox( normLinf(bv - bv2), 0.0) atol=1e-14

    smallmatmat!(ROView(Av2), xv, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    # test mat-matT
    b2 = A[:, :, 1]*x[:, :, 1].'
    smallmatmatT!(Av, xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatmatT!(Av2, xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    smallmatmatT!(A[:, :, 1], xv2, bv2)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    b3 = zeros(b2)
    smallmatmatT!(Av2, xv2, b3)

    @test isapprox( normLinf(b3 - b2), 0.0) atol=1e-14

    smallmatmatT!(ROView(Av2), xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatmatT!(ROView(Av2), xv, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

    # test matT-mat
    b2 = A[:, :, 1].'*x[:, :, 1]

    smallmatTmat!(Av, xv, bv)

    @test isapprox( normLinf(bv - b2), 0.0) atol=1e-14

    smallmatTmat!(Av2, xv2, bv2)

    @test isapprox( normLinf(bv2 - b2), 0.0) atol=1e-14

  end

  # A is d1 x d2, x = d2 x d3, b = d1 x d3
  for d1=2:20
    for d2=2:20
      for d3=2:20

        # mat-mat
        A = rand(d1, d2)
        x = rand(d2, d3)
        b = A*x
        b2 = rand(d1, d3)
        smallmatmat!(A, x, b2)
        b3 = smallmatmat(A, x)
        b4 = rand(d1, d3)
        b5 = 2*A*x + 3*b4
        smallmatmat_kernel!(A, x, b4, 2, 3)

        @test normLinf(b - b2) < 1e-13
        @test normLinf(b - b3) < 1e-13
        @test normLinf(b4 - b5) < 1e-13

        # mat-matT
        A = rand(d1, d2)
        x = rand(d3, d2)
        b = A*(x.')
        b2 = rand(d1, d3)
        smallmatmatT!(A, x, b2)
        b3 = smallmatmatT(A, x)
        b4 = rand(d1, d3)
        b5 = 2*A*(x.') + 3*b4
        smallmatmatT_kernel!(A, x, b4, 2, 3)

        @test normLinf(b - b2) < 1e-13
        @test normLinf(b - b3) < 1e-13
        @test normLinf(b4 - b5) < 1e-13

        # matT-mat
        A = rand(d2, d1)
        x = rand(d2, d3)
        b = A.'*x
        b2 = rand(d1, d3)
        smallmatTmat!(A, x, b2)
        b3 = smallmatTmat(A, x)
        b4 = rand(d1, d3)
        b5 = 2*A.'*x + 3*b4
        smallmatTmat_kernel!(A, x, b4, 2, 3)

        @test normLinf(b - b2) < 1e-13
        @test normLinf(b - b3) < 1e-13
        @test normLinf(b4 - b5) < 1e-13
      end
    end
  end

  # test mat-vec
  # A = d1 x d2, x = d2, b = d1
  for d1=1:20
    for d2=1:20
      # mat-vec
      A = rand(d1, d2)
      x = rand(d2)
      b = A*x
      b2 = rand(d1)
      smallmatvec!(A, x, b2)
      b3 = smallmatvec(A, x)
      b4 = rand(d1)
      b5 = 2*A*x + 3*b4
      smallmatvec_kernel!(A, x, b4, 2, 3)

      @test normLinf(b - b2) < 1e-13
      @test normLinf(b - b3) < 1e-13
      @test normLinf(b4 - b5) < 1e-13

      # matT-vec
      A = rand(d2, d1)
      x = rand(d2)
      b = A.'*x
      b2 = rand(d1)
      smallmatTvec!(A, x, b2)
      b3 = smallmatTvec(A, x)
      b4 = rand(d1)
      b5 = 2*A.'*x + 3*b4
      smallmatTvec_kernel!(A, x, b4, 2, 3)

      @test normLinf(b - b2) < 1e-13
      @test normLinf(b - b3) < 1e-13
      @test normLinf(b4 - b5) < 1e-13

    end
  end

  branch_name = getBranchName()
  @test  length(branch_name)  > 0

  time_str = getTimeString()
  @test  length(time_str)  > 0

  A = [0.0 0 1; 0 0 0; 1 0 0]
  numz, arr = checkZeroRows(A, eps())
  @test ( numz )== 1
  @test ( arr )== [false, true, false]

  numz, arr = checkZeroColumns(A, eps())
  @test ( numz )== 1
  @test ( arr )== [false, true, false]

  numz, arr = checkIdenticalColumns(A, 1, eps())
  @test ( numz )== 0
  @test ( arr )== [false, false, false]

  q = FIFOQueue{Int}()

  push!(q, 1)
  push!(q, 2)
  push!(q, 3)

  @test ( length(q) )== 3
  @test ( q.head )== 1
  @test ( q.tail )== 3
  val = pop!(q)
  @test ( val )== 1
  @test ( q.head )== 2
  @test ( isempty(q) )== false
  @test ( front(q) )== 2

  resize!(q, 10)
  @test ( length(q.s) )== 10

  empty!(q)
  @test ( isempty(q) )== true

  q2 = FIFOQueue{Int}(size_hint=100)
  for i=1:100
    push!(q2, i)
  end

  @test ( length(q2) )== 100
  for i=1:50
    @test ( pop!(q2) )== i
  end

  push!(q2, 1)
  @test ( length(q2.s) )== 100
  @test ( length(q2) )== 51
  @test ( pop!(q2) )== 51

  # test calcDiffElementArea
  dxidx = [1. 1; 0 1]
  nrm = [1., 0]
  workvec = rand(2)
  calcDiffElementArea(nrm, dxidx, workvec)
  @test ( workvec )== [1.0, 1.0]

  @test ( rmfile("abc.txt") )== nothing

  s = printbacktrace()
#  println("typeof(s) = ", typeof(s))
#  println("fieldnames(s) = ", fieldnames(s))
  @test ( typeof(s) )== String

  A = [1.0 2.0; 3.0 4.0]
  @test ( isSymmetric(A) )== false
  A[2,1] = A[1,2]
  @test ( isSymmetric(A) )== true
  tol = 1e-12
  A[2,1] += tol/2
  @test ( isSymmetric(A, tol) )== true
  A[2,1] += tol
  @test ( isSymmetric(A, tol) )== false

  make_symmetric!(A)
  @test ( isSymmetric(A, eps()) )== true

  # test topology type
  face_verts = [1 1 1 2; 2 2 3 3; 3 4 4 4]
  ElementTopology{3}(face_verts)  # test the assertions didn't fire

  face_verts[3,2] = 3 # duplicate a face
  @test_throws Exception  ElementTopology{3}(face_verts)

  ElementTopology3()
  ElementTopology2()

  face_verts = ElementTopology3().face_verts
  edge_verts = [1 2 3 1 2 3;
                2 3 1 4 4 4]
  topo = ElementTopology{3}(face_verts, edge_verts, topo2=ElementTopology2())

  @test ( topo.face_edges )== [1 1 2 3; 2 5 6 6; 3 4 5 4]
  @test ( topo.face_edges_flipped )== [false false false true;
                                     false false false false;
                                     false true  true  true]

  # test isFieldDefined

  obj = TestType()

  ret = isFieldDefined(obj, :a)
  @test ( ret )== true  # bitstypes fields are always defined

  obj.a = 2
  ret = isFieldDefined(obj, :a)
  @test ( ret )== true

  ret = isFieldDefined(obj, :a, :b)
  @test ( ret )== false

  @test_throws Exception  isFieldDefined(obj, :a, :c)
  @test_throws Exception  isFieldDefined(obj)

  obj.b = rand(3)
  ret = isFieldDefined(obj, :a, :b)
  @test ( ret )== true

  # test file name modification
  fname = "abc.dat"
  fname2 = get_parallel_fname(fname, 1)
  @test ( fname2 )== "abc_1.dat"

  fname = "abc"
  fname2 = get_parallel_fname(fname, 1)
  @test ( fname2 )== "abc_1"

  # test append_path, prepend_path
  pth = ""
  new_entry = "/dir1/dir2"
  new_pth = append_path(pth, new_entry)
  @test ( new_pth )== new_entry

  pth = "/dir0"
  new_pth = append_path(pth, new_entry)
  @test ( new_pth )== "/dir0:/dir1/dir2"

  pth = ""
  new_pth = prepend_path(pth, new_entry)
  @test ( new_pth )== new_entry

  pth = "/dir0"
  new_pth = prepend_path(pth, new_entry)
  @test ( new_pth )== "/dir1/dir2:/dir0"

  str = joinpath_ascii(pwd(), "hello")
  @test ( str )== joinpath(pwd(), "hello")

  # test splitting fnames
  pth = "test.dat"
  fstem, fext = split_fname(pth)
  @test ( fstem )== "test"
  @test ( fext )== ".dat"

  pth = "test"
  fstem, fext = split_fname(pth)
  @test ( fstem )== "test"
  @test ( fext )== ""

end


text:
msg_notfound(dir, filename) = "The MNIST file \"$filename\" was not found in \"$dir\". You can download the dataset at http://yann.lecun.com/exdb/mnist/, or alternatively use MNIST.download_helper(directory) to do it for you."

msg_prompt(dir, files) = """
Interactive session detected. MNIST.download_helper initiated.

Dataset: THE MNIST DATABASE of handwritten digits
Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
Website: http://yann.lecun.com/exdb/mnist/

[LeCun et al., 1998a]
    Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. "Gradient-based learning applied to document recognition." Proceedings of the IEEE, 86(11):2278-2324, November 1998

The specified directory \"$dir\" is missing the files $(join(map(f->"\"$f\"", files), ", ", " and ")) of the full data set.

The files are available for download at the offical website linked above.
We can download these files for you if you wish, but that doesn't free
you from the burden of using the data responsibly and respect copyright.
The authors of MNIST aren't really explicit about any terms of use,
so please read the website to make sure you want to download the dataset.

    http://yann.lecun.com/exdb/mnist/

Did you visit the website and want to download the dataset to \"$dir\"? [y/n] """

function downloaded_file(dir, filename)
    path = joinpath(dir, filename)
    if !isfile(path)
        if isinteractive()
            warn(msg_notfound(dir, filename))
            download_helper(dir)
        else
            error(msg_notfound(dir, filename))
        end
    end
    path
end

"""
    download_helper([dir]; [i_accept_the_terms_of_use = false])

Check if the MNIST dataset is contained in the specified `dir`,
or if any of the four files are missing. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`.

In the case that any of the four files are missing and
`i_accept_the_terms_of_use=false` the function will raise a
warning or an error depending on if julia is run in an
interactive session. If an interactive session is detected the
user will be presented with information and the option to
download the dataset to the specified `dir`.

If the download should happen automatically, please first visit
the website at http://yann.lecun.com/exdb/mnist, before setting
`i_accept_the_terms_of_use=true`.
"""
function download_helper(dir; i_accept_the_terms_of_use = false)
    files = filter(file->!isfile(joinpath(dir, file)),
                   [TRAINIMAGES, TRAINLABELS, TESTIMAGES, TESTLABELS])
    if !isempty(files)
        if !i_accept_the_terms_of_use && isinteractive()
            print(msg_prompt(dir, files))
            answer = first(readline())
            if answer == 'y'
                i_accept_the_terms_of_use = true
            end
        end
        if i_accept_the_terms_of_use
            mkpath(dir)
            for file in files
                url = "http://yann.lecun.com/exdb/mnist/$file"
                path = joinpath(dir, file)
                info("downloading $file from $url to $dir")
                run(download_cmd(url, path))
            end
        else
            error("Unable to download the dataset. Please visit the website at http://yann.lecun.com/exdb/mnist and download the files manually.")
        end
    else
        info("Nothing to do.")
    end
    nothing
end


text:
module PathLoss

# Exports
export
# models
FreeSpaceModel,
OkumuraHataModel,
Ecc33Model,
Cost231Model,
Cost231HataModel,
EricssonModel,
SuiModel,
LeeModel,
FlatEarthModel,
# functions
distanceInKm,
pathloss,
azimuthLineAtoB,
#enums
CityKind,
AreaKind,
TerrainKind,
LeeArea

# package code goes here
include(joinpath(Pkg.dir(),"PathLoss","src","Geo","GeoFunctions.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","Model.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Enums.jl"))
include(joinpath(Pkg.dir(),"PathLoss", "src","Models","Sui.jl"))
include(joinpath(Pkg.dir(),"PathLoss", "src","Models","FreeSpace.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","OkumuraHata.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","Cost231.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","Ecc33.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","Cost231Hata.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","Ericsson.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","LeeModel.jl"))
include(joinpath(Pkg.dir(),"PathLoss","src","Models","FlatEarthModel.jl"))
end # module


text:
#=
system.jl

    State space system abstract type and constructors

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/02/17
=#
"""
    StateSpaceSystem

Abstract type for state space systems.
"""
abstract type StateSpaceSystem end

"""
    LinearTimeInvariant <: StateSpaceSystem

Constructor for a linear time invariant instance of the state space system type
with system matrices ``Z``, ``T``, ``d``, ``c``, ``H``, and ``Q`` and
initialization of the Kalman filter ``a₁`` and ``P₁``.
"""
struct LinearTimeInvariant{Ty, TZ, TT, Td, Tc, TH, TQ, Ta, TP} <: StateSpaceSystem
	y::Ty	# data
	Z::TZ	# system matrix Z
	T::TT	# system matrix T
	d::Td	# mean adjustment observation
	c::Tc	# mean adjustment state
	H::TH	# system matrix H
	Q::TQ	# system matrix Q
	a1::Ta	# initial state
	P1::TP	# initial state variance
end
# Constructor
function LinearTimeInvariant{Ty, TZ, TT, Td, Tc, TH, TQ, Ta, TP}(n::Integer, p::Integer, T_len::Integer) where {Ty, TZ, TT, Td, Tc, TH, TQ, Ta, TP}
	# data
	y= Ty(undef, n, T_len)

    # Initialize system components
    Z= TZ(undef, n, p)
    T= TT <: Diagonal ? TT(undef, p) : TT(undef, p, p)
	d= Td(undef, n)
	c= Tc(undef, p)
    H= TH <: Diagonal ? TH(undef, n) : Symmetric(Matrix{eltype(TH)}(undef, n, n))
    Q= TQ <: Diagonal ? TQ(undef, p) : Symmetric(Matrix{eltype(TQ)}(undef, p, p))

	# Initial conditions
	a1= similar(Ta, p)
	P1= TP <: Diagonal ? TP(undef, p) : Symmetric(Matrix{eltype(TP)}(undef, p, p))

    return LinearTimeInvariant(y, Z, T, d, c, H, Q, a1, P1)
end

"""
    LinearTimeVariant <: StateSpaceSystem

Constructor for a linear time variant instance of the state space system type
with system matrices ``Zₜ``, ``Tₜ``, ``dₜ``, ``cₜ``, ``Hₜ``, and ``Qₜ`` and
initialization of the Kalman filter ``a₁`` and ``P₁``.
"""
struct LinearTimeVariant{Ty, TZ, TT, Td, Tc, TH, TQ, Ta, TP} <: StateSpaceSystem
	y::Ty	# data
	Z::TZ	# system matrices Zₜ
	T::TT	# system matrices Tₜ
	d::Td	# mean adjustments observation
	c::Tc	# mean adjustments state
	H::TH	# system matrices Hₜ
	Q::TQ	# system matrices Qₜ
	a1::Ta	# initial state
	P1::TP	# initial state variance
end
# Constructor
function LinearTimeVariant{Ty, TZ, TT, Td, Tc, TH, TQ, Ta, TP}(n::Integer, p::Integer, T_len::Integer) where {Ty, TZ, TT, Td, Tc, TH, TQ, Ta, TP}
	# data
	y= Ty(undef, n, T_len)

    # Initialize system components
    Z= [TZ(undef, n, p) for _ in 1:T_len]
    T= TT <: Diagonal ? [TT(undef, p) for _ in 1:T_len] :
                        [TT(undef, p, p) for _ in 1:T_len]
	d= Td(undef, n, T_len)
	c= Tc(undef, p, T_len)
    H= TH <: Diagonal ? [TH(undef, n) for _ in 1:T_len] :
                        [Symmetric(Matrix{eltype(TH)}(undef, n, n)) for _ in 1:T_len]
    Q= TQ <: Diagonal ? [TQ(undef, p) for _ in 1:T_len] :
                        [Symmetric(Matrix{eltype(TQ)}(undef, p, p)) for _ in 1:T_len]

	# Initial conditions
	a1= similar(Ta, p)
	P1= TP <: Diagonal ? TP(undef, p) : Symmetric(Matrix{eltype(TP)}(undef, p, p))

    return LinearTimeVariant(y, Z, T, d, c, H, Q, a1, P1)
end


text:
module Cyclotomics

import Primes
using SparseArrays

export coeffs, conductor, E, Cyclotomic

include("zumbroich.jl")
include("cycl.jl")
include("predicates.jl")
include("arithmetic.jl")
include("normalform.jl")
include("io.jl")

end # of module Cyclotomics


text:
#  Calling a function that requires no arguments:
f() = print("Hello world!")
f()

#  Calling a function with a fixed number of arguments:
function f(x, y, z)
    x*y - z^2
end

f(3, 4, 2)

#  Calling a function with optional arguments:
#  Note Julia uses multiple dispatch based on argument number and type, so
# f() is always different from f(x) unless default arguments are used, as in:

pimultiple(mult=1.0) = pi * mult # so pimultiple() defaults to pi * (1.0) or pi

#  Calling a function with a variable number of arguments:

f(a,b,x...) = reduce(+, 0, x) - a - b

# here a and b are single arguments, but x is a tuple of x plus whatever follows x, so:
a = b = c = d = e = 3
f(a,b,c)           # x within the function is (c) so == 0 + c - a - b
f(a,b,c,d,e)      # x is a tuple == (c,d,e) so == (0 + c + d + e) - a - b
f(a,b)             # x is () so == 0 - a - b

#  Calling a function with named arguments:
# Functions with keyword arguments are defined using a semicolon in the function signature,
#  as in
#             function plot(x, y; style="solid", width=1, color="black")
#
# When the function is called, the semicolon is optional, so plot here can be
# either called with plot(x, y, width=2) or less commonly as plot(x, y; width=2).

#  Using a function in statement context:
#  Any function can be used as a variable by its name.

circlearea(x) = x^2 * pi
map(circlearea, [r1, r2, r3, r4])

#  Using a function in first-class context within an expression:
cylindervolume = circlearea(r) * h

#  Obtaining the return value of a function:
radius = 2.5
area = circlearea(2.5)

#  Distinguishing built-in functions and user-defined functions:
#  Julia does not attempt to distinguish these in any special way,
#  but at the REPL command line there is ? help available for builtin
#  functions that would not generally be available for the user-defined ones.

#  Distinguishing subroutines and functions:
#  All subroutines are called functions in Julia, regardless of whether they return values.

#  Stating whether arguments are passed by value or by reference:
#  As in Python, all arguments are passed by pointer reference, but assignment to a passed argument
#  only changes the variable within the function. Assignment to the values referenced by the argument
## DOES however change those values. For instance:

a = 3
b = [3]
c = [3]

function f(x, y)
    a = 0
    b[1] = 0
    c = [0]
end # a and c are now unchanged but b = [0]

#  Is partial application possible and how:
#  In Julia, there are many different ways to compose functions. In particular,
#  Julia has an "arrow" operator -> that may be used to curry other functions.

f(a, b) = a^2 + a + b
v = [4, 6, 8]
map(x -> f(x, 10), v)  # v = [30, 52, 82]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using MPIReco

# multi-gradient is a special case of multi-patch
@testset "multi-gradient in-memory reconstruction" begin
  # low gradient with patch no 3
  b = MultiMPIFile(["./data/MG_G1", "./data/MG_G2_03"])
  bSFs = MultiMPIFile(["./data/SF_MG_G1", "./data/SF_MG_G2"])
  names = names = (:color, :x, :y, :z, :time)

  c1 = reconstruction(bSFs, b;
			    SNRThresh=2, frames=1, lambd=0.003, minFreq=80e3,
			    recChannels=1:2,iterations=3, roundPatches=false)
  @test axisnames(c1) == names
  @test axisvalues(c1) == (1:1, -26.0u"mm":1.0u"mm":26.0u"mm", -26.75u"mm":1.0u"mm":26.25u"mm", 0.0u"mm":1.0u"mm":0.0u"mm", 0.0u"ms":0.6528u"ms":0.0u"ms")
  im1 = reverse(c1[1,:,:,1,1]',dims=1)
  exportImage("./img/MultiGradient1.png", im1)

  #Iabs = abs.(im1)
  #Iabs[Iabs .< 0.2*maximum(Iabs)] .= 0
  #Icolored = colorview(Gray, Iabs./maximum(Iabs))
  #FileIO.save("./img/MultiGradient1.png", Icolored )

  # low gradient with all 3 high gradient patches
  b = MultiMPIFile(["./data/MG_G1", "./data/MG_G2_01", "./data/MG_G2_02", "./data/MG_G2_03"])

  c2 = reconstruction(bSFs, b;
			    SNRThresh=2, frames=1, lambd=0.003, minFreq=80e3,
			    recChannels=1:2,iterations=3, roundPatches=false)
  @test axisnames(c2) == names
  # TODO wo kommt dieser Floatingpoint Fehler rein? Lässt der sich vermeiden?
  @test axisvalues(c2) == (1:1, -26.0u"mm":1.0u"mm":26.0u"mm", -26.500000000000004u"mm":1.0u"mm":30.499999999999996u"mm", 0.0u"mm":1.0u"mm":0.0u"mm", 0.0u"ms":0.6528u"ms":0.0u"ms")
  im2 = reverse(c2[1,:,:,1,1]',dims=1)
  exportImage("./img/MultiGradient2.png", im2)

  #Iabs = abs.(im2)
  #Iabs[Iabs .< 0.2*maximum(Iabs)] .= 0
  #Icolored = colorview(Gray, Iabs./maximum(Iabs))
  #FileIO.save("./img/MultiGradient2.png", Icolored )

  # low gradient only
  b = MultiMPIFile(["./data/MG_G1"])
  bSFs = MultiMPIFile(["./data/SF_MG_G1"])

  c3 = reconstruction(bSFs, b;
			    SNRThresh=2, frames=1, lambd=0.003, minFreq=80e3,
			    recChannels=1:2,iterations=3, roundPatches=false)
  im3 = reverse(c3[1,:,:,1,1]',dims=1)
  exportImage("./img/MultiGradient3.png", im3)

  #Iabs = abs.(im3)
  #Iabs[Iabs .< 0.2*maximum(Iabs)] .= 0
  #Icolored = colorview(Gray, Iabs./maximum(Iabs))
  #FileIO.save("./img/MultiGradient3.png", Icolored )
end


text:
# Testing of samplers

using Distributions
using Base.Test

import Distributions:
    CategoricalDirectSampler,
    AliasTable,
    BinomialGeomSampler,
    BinomialTPESampler,
    BinomialPolySampler,
    BinomialAliasSampler,
    PoissonADSampler,
    PoissonCountSampler,
    ExponentialSampler,
    GammaMTSampler

n_tsamples = 10^6

## Categorical samplers

@test_throws ErrorException CategoricalDirectSampler(Float64[])
@test_throws ErrorException AliasTable(Float64[])

for S in [CategoricalDirectSampler, AliasTable]
    println("    testing $S")
    for p in Any[[1.0], [0.3, 0.7], [0.2, 0.3, 0.4, 0.1]]
        test_samples(S(p), Categorical(p), n_tsamples)
    end
end

## Binomial samplers

binomparams = [(0, 0.4), (0, 0.6), (5, 0.0), (5, 1.0),
               (1, 0.2), (1, 0.8), (3, 0.4), (4, 0.6),
               (40, 0.5), (100, 0.4), (300, 0.6)]

for (S, paramlst) in [
    (BinomialGeomSampler, [(0, 0.4), (0, 0.6), (5, 0.0), (5, 1.0), (1, 0.2), (1, 0.8), (3, 0.4), (4, 0.6)]),
    (BinomialTPESampler, [(40, 0.5), (100, 0.4), (300, 0.6)]),
    (BinomialPolySampler, binomparams),
    (BinomialAliasSampler, binomparams) ]

    println("    testing $S")
    for pa in paramlst
        n, p = pa
        test_samples(S(n, p), Binomial(n, p), n_tsamples)
    end
end

## Poisson samplers

for (S, paramlst) in [
    (PoissonCountSampler, [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0]),
    (PoissonADSampler, [5.0, 10.0, 15.0, 20.0, 30.0])]

    println("    testing $S")
    for μ in paramlst
        test_samples(S(μ), Poisson(μ), n_tsamples)
    end
end

## Exponential samplers

for S in [ExponentialSampler]
    println("    testing $S")
    for scale in [1.0, 2.0, 3.0]
        test_samples(S(scale), Exponential(scale), n_tsamples)
    end
end

## Gamma samplers

for S in [GammaMTSampler]
    for pa in [(1.0, 1.0), (2.0, 1.0), (3.0, 1.0), (0.5, 1.0),
               (1.0, 2.0), (3.0, 2.0), (0.5, 2.0)]
        test_samples(S(pa...), Gamma(pa...), n_tsamples)
    end
end


text:
using Plots
using DifferentialEquations

labs = ["Prey", "Predator"];
function f(du,u,p,t)
    du[1] = (p[1]*u[1]) - (p[2]*u[1]*u[2])
    du[2] = (p[2]*u[1]*u[2]) - (p[3]*u[2])
end
u0 = [10.0;4.0]
p = [0.5;0.05;0.05]
tspan = (0.0, 500.0)
prob = ODEProblem(f,u0,tspan,p)

# You need to solve the problem.  i.e. you need to run the simulation:
res = solve(prob)

# Then you need to get the simulated values in a nice format for plotting
t = res.t
prey = [values[1] for values in res.u]
pred = [values[2] for values in res.u]

# Plot the dynamics as usual
p1 = plot(t,pred)
plot!(t,prey)

# Plot the phase-space diagram
p2 = plot(pred,prey)

# Combine the two plots side by side
plot(p1,p2,layout=(1,2),legend=false,linewidth=3 )


text:
begin
    fwL = zeros(3)
    fwR = zeros(3)
    γ = 5 / 3
    prim = [1.0, 0.0, 1.0]
    w = KitBase.prim_conserve(prim, γ)
    dx = 1e-2
    dt = 1e-3
    res = zeros(3)
    avg = zeros(3)
end

KitBase.step!(fwL, w, prim, fwR, γ, dx, res, avg)

begin
    ffL = zeros(16)
    ffR = zeros(16)
    f = rand(16)
    u = randn(16)
    ω = ones(16)
end

KitBase.step!(fwL, ffL, w, prim, f, fwR, ffR, u, ω, γ, 1e-3, 0.72, 1.0, dx, dt, res, avg)

KitBase.step!(
    zeros(5),
    zeros(16, 16, 16),
    KitBase.prim_conserve([1.0, 0.0, 0.0, 0.0, 1.0], γ),
    [1.0, 0.0, 0.0, 0.0, 1.0],
    rand(16, 16, 16),
    zeros(5),
    zeros(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    ones(16, 16, 16),
    γ,
    1e-3,
    0.72,
    1.0,
    dx,
    dt,
    zeros(5),
    zeros(5),
)

KitBase.step!(
    fwL,
    ffL,
    ffL,
    w,
    prim,
    f,
    f,
    fwR,
    ffR,
    ffR,
    u,
    ω,
    2,
    γ,
    1e-3,
    0.72,
    1.0,
    dx,
    dt,
    res,
    avg,
)
KitBase.step!(
    zeros(3, 2),
    zeros(16, 2),
    zeros(16, 2),
    hcat(w, w),
    hcat(prim, prim),
    hcat(f, f),
    hcat(f, f),
    zeros(3, 2),
    zeros(16, 2),
    zeros(16, 2),
    hcat(u, u),
    hcat(ω, ω),
    2,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    1e-2,
    1.0,
    dx,
    dt,
    hcat(res, res),
    hcat(avg, avg),
)

# fsm
KitBase.step!(
    zeros(5),
    zeros(16, 16, 16),
    [1.0, 0.0, 0.0, 0.0, 1.0],
    KitBase.conserve_prim([1.0, 0.0, 0.0, 0.0, 1.0], 5 / 3),
    rand(16, 16, 16),
    zeros(5),
    zeros(16, 16, 16),
    5 / 3,
    1.0,
    5,
    ones(16, 16, 16, 20),
    ones(16, 16, 16, 20),
    ones(16, 16, 16),
    0.1,
    0.01,
    zeros(5),
    zeros(5),
    :fsm,
)
KitBase.step!(
    [1.0, 0.0, 0.0, 0.0, 1.0],
    KitBase.conserve_prim([1.0, 0.0, 0.0, 0.0, 1.0], 5 / 3),
    rand(16, 16, 16),
    zeros(5),
    zeros(16, 16, 16),
    zeros(5),
    zeros(16, 16, 16),
    zeros(5),
    zeros(16, 16, 16),
    zeros(5),
    zeros(16, 16, 16),
    5 / 3,
    1.0,
    5,
    ones(16, 16, 16, 20),
    ones(16, 16, 16, 20),
    ones(16, 16, 16),
    0.1,
    0.01,
    zeros(5),
    zeros(5),
    :fsm,
)

# Rykov
KitBase.step!(
    zeros(4),
    ffL,
    ffL,
    ffL,
    [1.0, 0.0, 1.0, 0.1],
    KitBase.conserve_prim([1.0, 0.0, 1.0, 0.1], 2, 2),
    f,
    f,
    f,
    zeros(4),
    ffR,
    ffR,
    ffR,
    u,
    ω,
    2,
    2,
    1e-3,
    0.81,
    0.72,
    91.5 / 273,
    18.1,
    1 / 1.55,
    0.2354,
    0.3049,
    dx,
    dt,
    zeros(4),
    zeros(4),
    :rykov,
)


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
pkg.jl defines package management functions.
"""

# --- Constants

const cli_tools = ["jltest", "jlcoverage", "jlcodestyle"]

# --- CLI installer functions

const default_julia_flags = ["--startup-file=no", "-q", "--compile=min", "-O0"]

"""
    TestTools.install(; <keyword arguments>)

Install all of the CLI utilities.

# Keyword arguments

* `julia::AbstractString`: path to julia executable. Default: path of the current running
    julia

* `bin_dir::AbstractString`: directory to install CLI utilities into.
    Default: `~/.julia/bin`

* `julia_flags::Vector{<:AbstractString}`: command line flags for CLI executables.
    Default: `["--startup-file=no", "-q", "--compile=min", "-O0"]`

* `force::Bool`: flag used to indicate that existing CLI executables should be
    overwritten. Default: `false`
"""
function install(;
    julia::AbstractString=joinpath(Sys.BINDIR, Base.julia_exename()),
    bin_dir::AbstractString=joinpath(DEPOT_PATH[1], "bin"),
    julia_flags::Vector{<:AbstractString}=default_julia_flags,
    force::Bool=false,
)
    # --- Install CLI utilities

    for cli in cli_tools
        install_cli(cli; julia=julia, bin_dir=bin_dir, julia_flags=julia_flags, force=force)
    end

    # --- Emit informational message

    @info """
          Make sure that `$(bin_dir)` is in PATH, or manually add a
          symlink from a directory in PATH to the installed program file.
          """
end

"""
    TestTools.install_cli(name::AbstractString; <keyword arguments>)

Install executable for CLI named `cli`.

Valid values for `name`: "jltest", "jlcoverage", "jlcodestyle".

# Keyword arguments

* `julia::AbstractString`: path to julia executable. Default: path of the current running
    julia

* `bin_dir::AbstractString`: directory to install CLI executable into.
    Default: `~/.julia/bin`

* `julia_flags::Vector{<:AbstractString}`: command line flags for CLI executable.
    Default: `["--startup-file=no", "-q", "--compile=min", "-O0"]`

* `force::Bool`: flag used to indicate that existing CLI executable should be
    overwritten. Default: `false`
"""
function install_cli(
    cli::AbstractString;
    julia::AbstractString=joinpath(Sys.BINDIR, Base.julia_exename()),
    bin_dir::AbstractString=joinpath(DEPOT_PATH[1], "bin"),
    julia_flags::Vector{<:AbstractString}=default_julia_flags,
    force::Bool=false,
)
    # --- Check arguments

    if !(cli in cli_tools)
        throw(ArgumentError("Invalid `cli`: $(cli)"))
    end

    # --- Preparations

    # Get OS
    os = Sys.iswindows() ? :windows : :unix

    # Set name of CLI
    if os == :windows
        cli *= ".cmd"
    end

    # Get absolute path to installation directory
    bin_dir = abspath(expanduser(bin_dir))

    # Get absolute path to executable to be installed
    exec_path = joinpath(bin_dir, cli)

    # Check if the executable already exists
    if ispath(exec_path) && !force
        error(
            "File `$(Base.contractuser(exec_path))` already exists. " *
            "Use `TestTools.install(force=true)` to overwrite.",
        )
    end

    # Create installation directory
    mkpath(bin_dir)

    # --- Install executable

    open(exec_path, "w") do io
        if os == :windows
            # TODO: Find a way to embed the script in the file
            print(
                io,
                """
:: -----------------------------------------------------------------------------------------
::  Copyright 2022 Velexi Corporation
::
::  Licensed under the Apache License, Version 2.0 (the "License");
::  you may not use this file except in compliance with the License.
::  You may obtain a copy of the License at
::
::      http://www.apache.org/licenses/LICENSE-2.0
::
::  Unless required by applicable law or agreed to in writing, software
::  distributed under the License is distributed on an "AS IS" BASIS,
::  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
::  See the License for the specific language governing permissions and
::  limitations under the License.
:: -----------------------------------------------------------------------------------------
@ECHO OFF
$(julia) $(join(julia_flags, ' ')) $(abspath(@__DIR__, cli, "cli", "main.jl")) %*
""",
            )

        else # unix
            print(
                io,
                """
#!/usr/bin/env bash
#=
exec $(julia) $(join(julia_flags, ' ')) "\${BASH_SOURCE[0]}" "\$@"
=#

#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

""",
            )
            open(abspath(@__DIR__, cli, "cli", "main.jl"), "r") do main
                line_count = 1
                for line in eachline(main)
                    if line_count > 10
                        println(io, line)
                    end
                    line_count += 1
                end
            end
        end
    end

    # Set permissions on executable
    chmod(exec_path, 0o0100755) # equivalent to -rwxrwxr-x (chmod +x exec_path)

    # --- Emit informational message

    @info "Installed $(cli) to `$(Base.contractuser(exec_path))`."

    return nothing
end

# --- CLI uninstaller functions

"""
    TestTools.uninstall(; <keyword arguments>)

Unnstall all of the CLI utilities.

# Keyword arguments

* `bin_dir::AbstractString`: directory containing CLI executables to uninstall.
    Default: `~/.julia/bin`
"""
function uninstall(; bin_dir::AbstractString=joinpath(DEPOT_PATH[1], "bin"))
    for cli in cli_tools
        uninstall_cli(cli; bin_dir=bin_dir)
    end
end

"""
    TestTools.uninstall_cli(name::AbstractString; <keyword arguments>)

Uninstall executable for CLI named `cli`.

Valid values for `name`: "jltest", "jlcoverage", "jlcodestyle".

# Keyword arguments

* `bin_dir::AbstractString`: directory containing CLI executable to uninstall.
    Default: `~/.julia/bin`
"""
function uninstall_cli(
    cli::AbstractString; bin_dir::AbstractString=joinpath(DEPOT_PATH[1], "bin")
)
    # --- Check arguments

    if !(cli in cli_tools)
        throw(ArgumentError("Invalid `cli`: $(cli)"))
    end

    # --- Preparations

    # Get OS
    os = Sys.iswindows() ? :windows : :unix

    # Set name of CLI
    if os == :windows
        cli *= ".cmd"
    end

    # Get absolute path to installation directory
    bin_dir = abspath(expanduser(bin_dir))

    # Get absolute path to executable to be installed
    exec_path = joinpath(bin_dir, cli)

    # --- Uninstall CLI

    rm(exec_path; force=true)

    # --- Emit informational message

    @info "Uninstalled `$(Base.contractuser(exec_path))`."

    return nothing
end


text:
function randuint() :: UInt32
    Base.llvmcall((
"""
define i32 @randuint() #0 {
  %1 = tail call { i32, i32 } @llvm.x86.rdrand.32() #1
  %2 = extractvalue { i32, i32 } %1, 0
  ret i32 %2
}

; Function Attrs: nounwind
declare { i32, i32 } @llvm.x86.rdrand.32() #1

attributes #0 = { nounwind ssp uwtable "less-precise-fpmad"="false" "no-frame-pointer-elim"="true" "no-frame-pointer-elim-non-leaf" "no-infs-fp-math"="true" "no-nans-fp-math"="true" "stack-protector-buffer-size"="8" "target-cpu"="core2" "target-features"="+ssse3,+cx16,+rdrnd,+sse,+sse2,+sse3" "unsafe-fp-math"="true" "use-soft-float"="false" }
attributes #1 = { nounwind }
"""
,
"""
%1 = call i32 @randuint()
ret i32 %1
"""),
        UInt32, Tuple{}
        )
end


text:
import Statistics
using Random

include("typedefs.jl")

function get_centroid(neighbourhood, id)::Node
    xdev = 40
    ydev = 40
    mean_x = Statistics.mean([node.cx for node in neighbourhood]) + rand(-xdev:xdev)
    mean_y = Statistics.mean([node.cx for node in neighbourhood]) + rand(-ydev:ydev)
    distances = [hypot(mean_x-node.cx, mean_y-node.cy) for node in neighbourhood]
    return Node(id, "provider", mean_x, mean_y)
end

function clusterise(nodelist::Vector{Node}, k::Int)::Vector{Node}
    return random_radius_clustering(nodelist, k)
end

function random_radius_clustering(nodelist::Vector{Node}, k::Int)::Vector{Node}
    n = length(nodelist)
    m = n - k
    min_r = 5
    max_r = 30
    centroids = []
    last_id = maximum(x -> x.id, nodelist)
    closed_set = copy(nodelist)
    for i in range(1, stop=m)
        open_set = copy(nodelist)
        random_point = open_set[rand(1:m)]
        distances = [hypot(random_point.cx-node.cx, random_point.cy-node.cy) for node in open_set]
        neighbours = []
        eligible_nodes = [node for (i, node) in enumerate(open_set) if distances[i] >= min_r && distances[i] <= max_r]
        eligible_distances = [distances[i] for (i, node) in enumerate(open_set) if distances[i] >= min_r && distances[i] <= max_r]
        # for j in range(1, stop=k)
        #     min_dist_ind = argmin(eligible_distances)
        #     min_dist_node = popat!(eligible_nodes, min_dist_ind)
        #     push!(neighbours, min_dist_node)
        # end
        push!(centroids, get_centroid(eligible_nodes, last_id))
        last_id += 1
    end

    append!(closed_set, centroids)
    return closed_set
end

max_stars_repo_name: dwirya/nodegen

max_issues_repo_name: dwirya/nodegen

max_forks_repo_name: dwirya/nodegen


text:
type GLFunc
    p::Ptr{Void}
end
# based on getCFun macro
macro glfunc(opengl_func)
    arguments = map(opengl_func.args[1].args[2:end]) do arg
        isa(arg, Symbol) && return Expr(:(::), arg)
        arg
    end
    # Get info out of arguments of `opengl_func`
    arg_names       = map(arg->arg.args[1], arguments)
    return_type     = opengl_func.args[2]
    input_types     = map(arg->arg.args[2], arguments)
    func_name       = opengl_func.args[1].args[1]
    func_name_sym   = Expr(:quote, func_name)
    func_name_str   = string(func_name)
    ptr_expr        = :(getprocaddress_e($func_name_str))

    if is_windows() # windows has some function pointers statically available and some not, this is how we deal with it:
        ptr = Libdl.dlsym_e(gl_lib, func_name)
        if (ptr != C_NULL)
            ptr_expr = :(($func_name_sym, "opengl32"))
            ret = quote
                function $func_name($(arg_names...))
                    ccall($ptr_expr, $return_type, ($(input_types...),), $(arg_names...))
                end
                $(Expr(:export, func_name))
            end
            return esc(ret)
        end
    end
    ptr_sym = gensym("$(func_name)_func_pointer")
    ret = quote
        const $ptr_sym = GLFunc(C_NULL)
        function $func_name($(arg_names...))
            if $ptr_sym.p::Ptr{Void} == C_NULL
                $ptr_sym.p::Ptr{Void} = $ptr_expr
            end
            ccall($ptr_sym.p::Ptr{Void}, $return_type, ($(input_types...),), $(arg_names...))
        end
        $(Expr(:export, func_name))
        end
    return esc(ret)
end
if is_windows()
    const gl_lib = Libdl.dlopen("opengl32")
end
include("glFunctions.jl")
if is_windows()
    Libdl.dlclose(gl_lib)
end


text:
export channel_match

channel_match( S::GphysData, j::Integer, fs::Float64) = (j == 0 ? 0 : S.fs[j] == fs ? j : 0)
function channel_match( S::GphysData,
                        j::Integer,
                        fs::Float64,
                        gain::AbstractFloat,
                        loc::InstrumentPosition,
                        resp::InstrumentResponse,
                        units::String
                      )
  (j == 0) && (return 0)
  return (
    try
      @assert(S.fs[j] == fs)
      @assert(S.gain[j] == gain)
      @assert(S.units[j] == units)
      @assert(loc == S.loc[j])
      @assert(resp == S.resp[j])
      j
    catch err
      0
    end
  )
end

function channel_match( S::GphysData, C::GphysChannel )
  j = findid(S, C.id)
  j = channel_match( S, j, C.fs, C.gain, C.loc, C.resp, C.units)
  return j
end

function channel_match( C::GphysChannel, D::GphysChannel ; use_gain::Bool=true )
  ff = use_gain ? (:id, :fs, :gain, :loc, :resp, :units) : (:id, :fs, :loc, :resp, :units)
  return all([isequal(getfield(C, f), getfield(D, f)) for f in ff ])
end

# This will seek out a match and correct unset values if a partial match is found
function cmatch_p!( C::GphysChannel, D::GphysChannel )
  fs_match    = max(    C.fs == D.fs   ,    C.fs == default_fs   ,    D.fs == default_fs   )
  gain_match  = max(  C.gain == D.gain ,  C.gain == default_gain ,  D.gain == default_gain )
  loc_match   = max(   C.loc == D.loc  ,   C.loc == default_loc  ,   D.loc == default_loc  )
  resp_match  = max(  C.resp == D.resp ,  C.resp == default_resp ,  D.resp == default_resp )
  units_match = max( C.units == D.units, C.units == ""           , D.units == ""           )
  m = min(C.id == D.id, fs_match, gain_match, loc_match, resp_match, units_match)
  if m
    # Fill any "unset" values in D from C
    (    D.fs == default_fs   ) && (    D.fs = C.fs               )
    (  D.gain == default_gain ) && (  D.gain = C.gain             )
    (   D.loc == default_loc  ) && (   D.loc = deepcopy(C.loc)    )
    (  D.resp == default_resp ) && (  D.resp = deepcopy(C.resp)   )
    ( D.units == ""           ) && ( D.units = identity(C.units)  )

    # Fill any "unset" values in C from D
    (    C.fs == default_fs   ) && (    C.fs = D.fs               )
    (  C.gain == default_gain ) && (  C.gain = D.gain             )
    (   C.loc == default_loc  ) && (   C.loc = deepcopy(D.loc)    )
    (  C.resp == default_resp ) && (  C.resp = deepcopy(D.resp)   )
    ( C.units == ""           ) && ( C.units = identity(D.units)  )
  end
  return m
end


text:
#=
    Rule based controller
=#
mutable struct RBCOptions
    policy_selection::Int64

    RBCOptions(; policy_selection = 1) = new(policy_selection)
end

mutable struct RBC <: AbstractController
    options::RBCOptions
    decisions::NamedTuple
    history::AbstractScenarios

    RBC(; options = RBCOptions()) = new(options)
end

### Policies
function π_1(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::RBC)
    # Utils to simplify the writting
    Δh = mg.parameters.Δh
    liion, tes, h2tank = mg.storages[1], mg.storages[2], mg.storages[3]
    heater, elyz, fc = mg.converters[1], mg.converters[2], mg.converters[3]

    # Net power elec
    p_net_E = mg.demands[1].carrier.power[h,y,s] - mg.generations[1].carrier.power[h,y,s]

    # Liion
    _, _, u_liion = compute_operation_dynamics(liion, (Erated = liion.Erated[y,s], soc = liion.soc[h,y,s], soh = liion.soh[h,y,s]), p_net_E, Δh)

    if p_net_E < 0.
        # Elyz
        _, u_elyz_E, elyz_H, elyz_H2 = compute_operation_dynamics(elyz, (powerMax = elyz.powerMax[y,s], soh = elyz.soh[h,y,s]), p_net_E - u_liion, Δh)
        # H2 tank
        _, u_h2tank = compute_operation_dynamics(h2tank, (Erated = h2tank.Erated[y,s], soc = h2tank.soc[h,y,s]), -elyz_H2, Δh)
        # Test H2
        elyz_H2 == - u_h2tank ? nothing : u_elyz_E = elyz_H = elyz_H2 = u_h2tank = 0.
        # FC
        u_fc_E, fc_H, fc_H2 = 0., 0., 0.
        # Heater
        u_heater_E, heater_H = compute_operation_dynamics(heater, (powerMax = heater.powerMax[y,s],), p_net_E - u_liion - u_elyz_E, Δh)
    else
        # FC
        _, u_fc_E, fc_H, fc_H2 = compute_operation_dynamics(fc, (powerMax = fc.powerMax[y,s], soh = fc.soh[h,y,s]), p_net_E - u_liion, Δh)
        # H2 tank
        _, u_h2tank = compute_operation_dynamics(h2tank, (Erated = h2tank.Erated[y,s], soc = h2tank.soc[h,y,s]), -fc_H2, Δh)
        # Test H2
        fc_H2 == - u_h2tank ? nothing : u_fc_E = fc_H = fc_H2 = u_h2tank = 0.
        # Elyz
        u_elyz_E, elyz_H, elyz_H2 = 0., 0., 0.
        # Heater
        u_heater_E, heater_H = 0., 0.
    end

    # Net heating power post H2
    p_net_H = mg.demands[2].carrier.power[h,y,s] - fc_H - elyz_H - heater_H

    # TES
    _, u_tes = compute_operation_dynamics(tes, (Erated = tes.Erated[y,s], soc = tes.soc[h,y,s]), p_net_H, Δh)

    # Heater
    if p_net_H < 0.
        _u_heater_E = 0.
    else
        _u_heater_E, _ = compute_operation_dynamics(heater, (powerMax = heater.powerMax[y,s],), - (p_net_H - u_tes) / heater.η_E_H, Δh)
    end

    # Store values
    controller.decisions.storages[1][h,y,s] = u_liion
    controller.decisions.storages[2][h,y,s] = u_tes
    controller.decisions.storages[3][h,y,s] = u_h2tank
    controller.decisions.converters[1][h,y,s]  = u_heater_E + _u_heater_E
    controller.decisions.converters[2][h,y,s] = u_elyz_E
    controller.decisions.converters[3][h,y,s] = u_fc_E
end
function π_2(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::RBC)
    controller.decisions.storages[1][h,y,s] = mg.demands[1].carrier.power[h,y,s] - mg.generations[1].carrier.power[h,y,s]
end
function π_3(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::RBC)
    # Compute the heater electrical power based on the simple model
    controller.decisions.converters[1][h,y,s] = - max(min(mg.demands[2].carrier.power[h,y,s] / mg.converters[1].η_E_H, mg.converters[1].powerMax[y,s]), 0.)
    # Compute the liion decision from the power balance
    controller.decisions.storages[1][h,y,s] = mg.demands[1].carrier.power[h,y,s] - mg.generations[1].carrier.power[h,y,s] - controller.decisions.converters[1][h,y,s]
end

### Offline
function initialize_controller!(mg::Microgrid, controller::RBC, ω::AbstractScenarios)
    # Preallocation
    preallocate!(mg, controller)

    return controller
end

### Online
function compute_operation_decisions!(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::RBC)
    # Chose policy
    if controller.options.policy_selection == 1
        return π_1(h, y, s, mg, controller)
    elseif controller.options.policy_selection == 2
        return π_2(h, y, s, mg, controller)
    elseif controller.options.policy_selection == 3
        return π_3(h, y, s, mg, controller)
    else
        println("Policy not defined !")
    end
end


text:
doc"""
Various parameters for the solve! call.
```
  maxiter := maximum total number of iterations
  erroriter := number of iterations between updating the SolverStats object and checking stopping conditions
  breakiter := number of iterations between the `onbreak` callback
  stopepochs := number of epochs (1 epoch == 1 update to stats) without improvement for early stopping
  minerror := when validation error drops below this, stop
  onbreak := callback function, run once every breakiter iterations. function args: (net::NeuralNet, solverParams::SolverParams, stats::SolverStats)
```
"""
type SolverParams
  maxiter::Int   # maximum total number of iterations
  erroriter::Int  # number of iterations be
  breakiter::Int
  plotiter::Int
  plotfields::Vector{Symbol}
  stopepochs::Int
  minerror::Float64
  onbreak::Function
end

function SolverParams(; maxiter = 1000,
                        erroriter = 1000,
                        breakiter = 10000,
                        plotiter = -1,
                        plotfields = Symbol[:x, :xhat, :y, :Σ, :a],
                        stopepochs = 100,
                        minerror = 1e-5,
                        onbreak = donothing)
  SolverParams(maxiter, erroriter, breakiter, plotiter, plotfields, stopepochs, minerror, onbreak)
end

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

@enum SolverStatus RUNNING CONVERGED STOPPEDEARLY MAXITER

type SolverStats
  numiter::Int
  trainError::Float64
  validationError::Float64
  bestValidationError::Float64
  epochSinceImprovement::Int
  status::SolverStatus
  bestModel
  plotter
end
SolverStats() = SolverStats(0, Inf, Inf, Inf, 0, RUNNING, nothing, nothing)

Base.string(stats::SolverStats) = "SolverStats{$(stats.status) n=$(stats.numiter), trainerr=$(stats.trainError), valerr=$(stats.validationError), besterr=$(stats.bestValidationError), epochSinceImprovement=$(stats.epochSinceImprovement)}"
Base.print(io::IO, stats::SolverStats) = print(io, string(stats))
Base.show(io::IO, stats::SolverStats) = print(io, string(stats))

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

doc"""
Batch stochastic gradient descent solver.  Sample from training data to update net, stop if converged
or if validation error is not improving.  Returns SolverStats summary object.
"""
function solve!(net::NetStat, traindata::DataSampler, validationdata::DataSampler, transformY::Bool = false)

  stats = SolverStats()

  # optionally set up the plotter (only if plotiter ≥ 0)
  if net.solverParams.plotiter >= 0
    stats.plotter = NetProgressPlotter(net, stats, net.solverParams.plotfields)
  end

  println("\nsolve: $net\n")

  # loop through maxiter times
  for i in 1:net.solverParams.maxiter

    stats.numiter += 1

    # randomly sample one data item and update the network
    data = sample(traindata)
    fit!(net, data, transformY)

    # update the plot?
    piter = net.solverParams.plotiter
    if (piter > 0 && i % net.solverParams.plotiter == 0) || piter == 0
      fit!(stats.plotter)
    end

    # # check for convergence
    if i % net.solverParams.erroriter == 0
      stats.trainError = totalCost(net, traindata)
      stats.validationError = totalCost(net, validationdata)
      println("Status: $stats")

      # check for improvement in validation error
      if stats.validationError < stats.bestValidationError
        stats.bestValidationError = stats.validationError
        stats.bestModel = copy(net)
        stats.epochSinceImprovement = 0
      else
        stats.epochSinceImprovement += 1

        # early stopping... no improvement
        if stats.epochSinceImprovement >= net.solverParams.stopepochs
          println("Early stopping: $stats")
          stats.status = STOPPEDEARLY
          return stats
        end
      end

      # check if our error is low enough
      if stats.validationError <= net.solverParams.minerror
        println("Converged, breaking: $stats")
        stats.status = CONVERGED
        return stats
      end
    end

    # take a break?
    if i % net.solverParams.breakiter == 0
      net.solverParams.onbreak(net, stats)
    end
  end

  println("Maxiter reached: $stats")
  stats.status = MAXITER
  stats
end


text:
using Ipopt
using PowerModelsPrivacyPreserving
using PowerModels
using Distributions
using Random

const PM = PowerModels
const PMPP = PowerModelsPrivacyPreserving

ipopt = Ipopt.Optimizer
file =  "test/data/matpower/case5.m"
data_unpert = PM.parse_file(file)

for (i,gen) in data_unpert["gen"]
    gen["eta"] = 0.95
end
# data_min_loss = deepcopy(data_unpert)
# data_min_cost = deepcopy(data_unpert)

result_unpert_cost = PMPP.run_opf_bf_dvorkin(data_unpert, BFAPowerModel, ipopt)

##
# https://github.com/mlubin/RobustCCOPFSupplement/
# https://github.com/mlubin/RobustCCOPFSupplement/blob/master/codejl/ccopfmodel_simulation.jl
# Lubin, M., Dvorkin, Y., & Backhaus, S. (2015). A robust approach to chance constrained optimal power flow with renewable generation. IEEE Trans. Power Syst., 11. Cryptography and Security. Retrieved from http://arxiv.org/abs/1504.0601
# https://jumpchance.readthedocs.io/en/latest/
# https://github.com/mlubin/JuMPChance.jl
result_unpert_cost_cc = PMPP.run_opf_bf_dvorkin_cc(data_unpert, BFAPowerModel, ipopt)


text:
# Autogenerated wrapper script for FFTW_jll for x86_64-apple-darwin
export libfftw3, libfftw3f

JLLWrappers.@generate_wrapper_header("FFTW")
JLLWrappers.@declare_library_product(libfftw3, "@rpath/libfftw3.3.dylib")
JLLWrappers.@declare_library_product(libfftw3f, "@rpath/libfftw3f.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfftw3,
        "lib/libfftw3.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3f,
        "lib/libfftw3f.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Ferrite
using Tensors
using Test
using ForwardDiff
import SHA
using Random
using LinearAlgebra
using SparseArrays

include("test_utils.jl")
include("test_interpolations.jl")
include("test_cellvalues.jl")
include("test_facevalues.jl")
include("test_quadrules.jl")
include("test_assemble.jl")
include("test_dofs.jl")
include("test_constraints.jl")
include("test_grid_dofhandler_vtk.jl")
include("test_abstractgrid.jl")
include("test_mixeddofhandler.jl")
include("test_l2_projection.jl")
include("test_pointevaluation.jl")
# include("test_notebooks.jl")
include("test_apply_rhs.jl")
include("test_examples.jl")
@test all(x -> isdefined(Ferrite, x), names(Ferrite))  # Test that all exported symbols are defined


text:
using Dates, Test

include("clock.jl")

@testset "Create a new clock with an initial time" begin
    # on the hour
    @test Clock(8, 0) == Clock(8, 0)

    # past the hour
    @test Clock(11, 9) == Clock(11, 9)

    # midnight is zero hours
    @test Clock(24, 0) == Clock(0, 0)

    # hour rolls over
    @test Clock(25, 0) == Clock(1, 0)

    # hour rolls over continuously
    @test Clock(100, 0) == Clock(4, 0)

    # sixty minutes is next hour
    @test Clock(1, 60) == Clock(2, 0)

    # minutes roll over
    @test Clock(0, 160) == Clock(2, 40)

    # minutes roll over continuously
    @test Clock(0, 1723) == Clock(4, 43)

    # hour and minutes roll over
    @test Clock(25, 160) == Clock(3, 40)

    # hour and minutes roll over continuously
    @test Clock(201, 3001) == Clock(11, 1)

    # hour and minutes roll over to exactly midnight
    @test Clock(72, 8640) == Clock(0, 0)

    # negative hour
    @test Clock(-1, 15) == Clock(23, 15)

    # negative hour rolls over
    @test Clock(-25, 0) == Clock(23, 0)

    # negative hour rolls over continuously
    @test Clock(-91, 0) == Clock(5, 0)

    # negative minutes
    @test Clock(1, -40) == Clock(0, 20)

    # negative minutes roll over
    @test Clock(1, -160) == Clock(22, 20)

    # negative minutes roll over continuously
    @test Clock(1, -4820) == Clock(16, 40)

    # negative sixty minutes is previous hour
    @test Clock(2, -60) == Clock(1, 0)

    # negative hour and minutes both roll over
    @test Clock(-25, -160) == Clock(20, 20)

    # negative hour and minutes both roll over continuously
    @test Clock(-121, -5810) == Clock(22, 10)
end

@testset "Add minutes" begin
    # add minutes
    @test Clock(10, 0) + Dates.Minute(3) == Clock(10, 3)

    # add no minutes
    @test Clock(6, 41) + Dates.Minute(0) == Clock(6, 41)

    # add to next hour
    @test Clock(0, 45) + Dates.Minute(40) == Clock(1, 25)

    # add more than one hour
    @test Clock(10, 0) + Dates.Minute(61) == Clock(11, 1)

    # add more than two hours with carry
    @test Clock(0, 45) + Dates.Minute(160) == Clock(3, 25)

    # add across midnight
    @test Clock(23, 59) + Dates.Minute(2) == Clock(0, 1)

    # add more than one day (1500 min = 25 hrs)
    @test Clock(5, 32) + Dates.Minute(1500) == Clock(6, 32)

    # add more than two days
    @test Clock(1, 1) + Dates.Minute(3500) == Clock(11, 21)
end

@testset "Subtract minutes" begin
    # subtract minutes
    @test Clock(10, 3) - Dates.Minute(3) == Clock(10, 0)

    # subtract to previous hour
    @test Clock(10, 3) - Dates.Minute(30) == Clock(9, 33)

    # subtract more than an hour
    @test Clock(10, 3) - Dates.Minute(70) == Clock(8, 53)

    # subtract across midnight
    @test Clock(0, 3) - Dates.Minute(4) == Clock(23, 59)

    # subtract more than two hours
    @test Clock(0, 0) - Dates.Minute(160) == Clock(21, 20)

    # subtract more than two hours with borrow
    @test Clock(6, 15) - Dates.Minute(160) == Clock(3, 35)

    # subtract more than one day (1500 min = 25 hrs)
    @test Clock(5, 32) - Dates.Minute(1500) == Clock(4, 32)

    # subtract more than two days
    @test Clock(2, 20) - Dates.Minute(3000) == Clock(0, 20)
end

@testset "Compare two clocks for equality" begin
    # clocks with same time
    @test Clock(15, 37) == Clock(15, 37)

    # clocks a minute apart
    @test Clock(15, 36) != Clock(15, 37)

    # clocks an hour apart
    @test Clock(14, 37) != Clock(15, 37)

    # clocks with hour overflow
    @test Clock(10, 37) == Clock(34, 37)

    # clocks with hour overflow by several days
    @test Clock(3, 11) == Clock(99, 11)

    # clocks with negative hour
    @test Clock(22, 40) == Clock(-2, 40)

    # clocks with negative hour that wraps
    @test Clock(17, 3) == Clock(-31, 3)

    # clocks with negative hour that wraps multiple times
    @test Clock(13, 49) == Clock(-83, 49)

    # clocks with minute overflow
    @test Clock(0, 1) == Clock(0, 1441)

    # clocks with minute overflow by several days
    @test Clock(2, 2) == Clock(2, 4322)

    # clocks with negative minute
    @test Clock(2, 40) == Clock(3, -20)

    # clocks with negative minute that wraps
    @test Clock(4, 10) == Clock(5, -1490)

    # clocks with negative minute that wraps multiple times
    @test Clock(6, 15) == Clock(6, -4305)

    # clocks with negative hours and minutes
    @test Clock(7, 32) == Clock(-12, -268)

    # clocks with negative hours and minutes that wrap
    @test Clock(18, 7) == Clock(-54, -11513)

    # full clock and zeroed clock
    @test Clock(24, 0) == Clock(0, 0)
end

@testset "displaying clocks" begin
    @test sprint(show, Clock(8, 0)) == "\"08:00\""
    @test sprint(show, Clock(11, 9)) == "\"11:09\""
    @test sprint(show, Clock(24, 0)) == "\"00:00\""
    @test sprint(show, Clock(25, 0)) == "\"01:00\""
    @test sprint(show, Clock(100, 0)) == "\"04:00\""
    @test sprint(show, Clock(1, 60)) == "\"02:00\""
    @test sprint(show, Clock(0, 160)) == "\"02:40\""
    @test sprint(show, Clock(0, 1723)) == "\"04:43\""
    @test sprint(show, Clock(25, 160)) == "\"03:40\""
    @test sprint(show, Clock(201, 3001)) == "\"11:01\""
    @test sprint(show, Clock(72, 8640)) == "\"00:00\""
    @test sprint(show, Clock(-1, 15)) == "\"23:15\""
    @test sprint(show, Clock(-25, 0)) == "\"23:00\""
    @test sprint(show, Clock(-91, 0)) == "\"05:00\""
    @test sprint(show, Clock(1, -40)) == "\"00:20\""
    @test sprint(show, Clock(1, -160)) == "\"22:20\""
    @test sprint(show, Clock(1, -4820)) == "\"16:40\""
    @test sprint(show, Clock(2, -60)) == "\"01:00\""
    @test sprint(show, Clock(-25, -160)) == "\"20:20\""
    @test sprint(show, Clock(-121, -5810)) == "\"22:10\""
end


text:
using Images
using JSON
using ArgParse
using HDF5

SPLITS = ["train", "restval", "val", "test"]

include("lib/data.jl")
include("lib/imgproc.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "Convert common image captioning datasets to JLD format."

    @add_arg_table s begin
        ("--images"; required=true; help="images dir")
        ("--captions"; required=true; help="captions archive file (karpathy)")
        ("--savefile"; required=true; help="output file in HDF5 format")
        ("--imsize"; arg_type=Int; nargs=2; default=[224,224])
        ("--feedback"; arg_type=Int; default=0)
        ("--debug"; action=:store_true)
        ("--seed"; arg_type=Int; default=0; help="random seed")
        ("--nocrop"; action=:store_true)
        ("--randomcrop"; action=:store_true)
        ("--extradata"; action=:store_true)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); display(o)
    println(); flush(STDOUT)
    o[:seed] > 0 && srand(o[:seed])

    imgpath = abspath(o[:images])
    capfile = abspath(o[:captions])
    newsize = tuple(o[:imsize]...)
    crop = !o[:nocrop]
    savefile = abspath(o[:savefile])
    f = h5open(savefile, "w")

    # process images
    for splitname in SPLITS
        @printf("Processing %s split... [%s]\n", splitname, now())
        flush(STDOUT)

        randomcrop = false
        if splitname == "train" || (splitname == "restval" && o[:extradata])
            randomcrop = o[:randomcrop]
        end

        entries = first(get_entries(capfile, [splitname]))
        data = Any[]
        for i = 1:length(entries)
            entry = entries[i]
            if o[:debug]
                @printf("Image: %s\n", entry["filename"]); flush(STDOUT)
            end

            # processing
            img = read_image(entry["filename"], imgpath)
            img = process_image(img, newsize; crop=crop, randomcrop=randomcrop)

            write(f, entry["filename"], img)
            entries[i] = 0

            # feedback
            if o[:feedback] > 0 && i % o[:feedback] == 0
                @printf("Processed %d images by so far...\n", i)
                flush(STDOUT)
            end
        end
    end
    close(f)
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
#! /usr/bin/julia

# Rosetta Code, Hostname

println(gethostname())

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# DynamicBounds.jl
# A package for compute bounds and relaxations of the solutions of
# parametric differential equations.
# See https://github.com/PSORLab/DynamicBounds.jl
#############################################################################
# src/library.jl
# A module containing the access functions and library utilities for the
# library of relaxation problem benchmarks.
#############################################################################

export STANDARD_LIBRARY, reinitialize_standard_lib!, fetch_instance

const STANDARD_LIB_PATH = joinpath(@__DIR__(), "library")
const PROTECTED_LIBS = ["pODEs"]

mutable struct LibraryProblem
    prob
    id::Symbol
    desc::String
    source::String
    url::String
end

const STANDARD_LIBRARY = LibraryProblem[]

function reinitialize_standard_lib!()
    for folder in PROTECTED_LIBS
        for f in readdir(joinpath(STANDARD_LIB_PATH, folder))
            include(joinpath(STANDARD_LIB_PATH, folder, f))
            lib_prob = LibraryProblem(prob, id, desc, source, url)
            push!(STANDARD_LIBRARY, lib_prob)
        end
    end
end

pdelib_dir = joinpath(@__DIR__(), "library")
reinitialize_standard_lib!()

"""
    fetch_instance
"""
function fetch_instance(lib::String, instance::String)
    simple_instance = replace(instance, Pair(".jl", ""))
    if isfile(joinpath(pdelib_dir, lib, "$(simple_instance).jl"))
        include(joinpath(pdelib_dir, lib, "$(simple_instance).jl"))
        return LibraryProblem(prob, id, desc, source, url)
    else
        @warn "No instance detected..."
        return nothing
    end
end


text:
module TestSparseMatrix

import SparseArrays
using Test

import MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

"""
    test_empty_size()

Julia 1.7 is going to outlaw constructing empty sparse matrices. Test that we
can still do so.
"""
function test_empty_size()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.ZeroBasedIndexing,
    }()
    MOI.empty!(A)
    @test A.rowval == Int[]
    @test A.nzval == Float64[]
    @test A.colptr == Int[0]
end

function test_empty_large()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.ZeroBasedIndexing,
    }()
    MOI.empty!(A)
    MOI.Utilities.add_column(A)
    MOI.Utilities.add_column(A)
    MOI.Utilities.set_number_of_rows(A, 5)
    B = convert(SparseArrays.SparseMatrixCSC{Float64,Int}, A)
    @test size(B) == (5, 2)
end

function test_VectorAffine_ZeroBased()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.ZeroBasedIndexing,
    }()
    MOI.empty!(A)
    x = MOI.VariableIndex.(1:3)
    f = MOI.VectorAffineFunction(
        vcat(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(1.0, x)),
            MOI.VectorAffineTerm.(2, MOI.ScalarAffineTerm.([2.0, 3.0], x[2:3])),
        ),
        [0.5, 1.2],
    )

    index_map = MOI.Utilities.IndexMap()
    for i in 1:3
        MOI.Utilities.add_column(A)
        index_map[x[i]] = x[i]
    end
    MOI.Utilities.allocate_terms(A, index_map, f)
    MOI.Utilities.set_number_of_rows(A, 2)
    MOI.Utilities.load_terms(A, index_map, f, 0)
    MOI.Utilities.final_touch(A)
    B = convert(SparseArrays.SparseMatrixCSC{Float64,Int}, A)
    @test B == [1.0 1.0 1.0; 0.0 2.0 3.0]
    @test A.rowval == [0, 0, 1, 0, 1]
    @test A.nzval == [1.0, 1.0, 2.0, 1.0, 3.0]
    @test A.colptr == [0, 1, 3, 5]
end

function test_extract_function()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.ZeroBasedIndexing,
    }()
    MOI.empty!(A)
    x = MOI.VariableIndex.(1:3)
    f = MOI.VectorAffineFunction(
        vcat(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(1.0, x)),
            MOI.VectorAffineTerm.(2, MOI.ScalarAffineTerm.([2.0, 3.0], x[2:3])),
        ),
        [0.5, 1.2],
    )

    index_map = MOI.Utilities.IndexMap()
    MOI.Utilities.add_columns(A, 3)
    for i in 1:3
        index_map[x[i]] = x[i]
    end
    MOI.Utilities.allocate_terms(A, index_map, f)
    MOI.Utilities.set_number_of_rows(A, 2)
    MOI.Utilities.load_terms(A, index_map, f, 0)
    MOI.Utilities.final_touch(A)
    row_1 = MOI.Utilities.extract_function(A, 1, 0.5)
    @test row_1 ≈ MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.5)
    @test MOI.Utilities.extract_function(A, 1:2, [0.5, 1.2]) ≈ f
    empty_f = MOI.Utilities.extract_function(A, 1:0, Float64[])
    @test empty_f ≈
          MOI.VectorAffineFunction(MOI.VectorAffineTerm{Float64}[], Float64[])
    return
end

function test_VectorAffine_OneBased()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.OneBasedIndexing,
    }()
    MOI.empty!(A)
    x = MOI.VariableIndex.(1:3)
    f = MOI.VectorAffineFunction(
        vcat(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(1.0, x)),
            MOI.VectorAffineTerm.(2, MOI.ScalarAffineTerm.([2.0, 3.0], x[2:3])),
        ),
        [0.5, 1.2],
    )

    index_map = MOI.Utilities.IndexMap()
    for i in 1:3
        MOI.Utilities.add_column(A)
        index_map[x[i]] = x[i]
    end
    MOI.Utilities.allocate_terms(A, index_map, f)
    MOI.Utilities.set_number_of_rows(A, 2)
    MOI.Utilities.load_terms(A, index_map, f, 0)
    MOI.Utilities.final_touch(A)
    B = convert(SparseArrays.SparseMatrixCSC{Float64,Int}, A)
    @test B == [1.0 1.0 1.0; 0.0 2.0 3.0]
    @test A.rowval == [1, 1, 2, 1, 2]
    @test A.nzval == [1.0, 1.0, 2.0, 1.0, 3.0]
    @test A.colptr == [1, 2, 4, 6]
end

function test_ScalarAffine_ZeroBased()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.ZeroBasedIndexing,
    }()
    MOI.empty!(A)
    x = MOI.VariableIndex.(1:3)
    f1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.5)
    f2 =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0], x[2:3]), 1.2)
    index_map = MOI.Utilities.IndexMap()
    for i in 1:3
        MOI.Utilities.add_column(A)
        index_map[x[i]] = x[i]
    end
    MOI.Utilities.allocate_terms(A, index_map, f1)
    MOI.Utilities.allocate_terms(A, index_map, f2)
    MOI.Utilities.set_number_of_rows(A, 2)
    MOI.Utilities.load_terms(A, index_map, f1, 0)
    MOI.Utilities.load_terms(A, index_map, f2, 1)
    MOI.Utilities.final_touch(A)
    B = convert(SparseArrays.SparseMatrixCSC{Float64,Int}, A)
    @test B == [1.0 1.0 1.0; 0.0 2.0 3.0]
    @test A.rowval == [0, 0, 1, 0, 1]
    @test A.nzval == [1.0, 1.0, 2.0, 1.0, 3.0]
    @test A.colptr == [0, 1, 3, 5]
end

function test_ScalarAffine_OneBased()
    A = MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        Int,
        MOI.Utilities.OneBasedIndexing,
    }()
    MOI.empty!(A)
    x = MOI.VariableIndex.(1:3)
    f1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.5)
    f2 =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0], x[2:3]), 1.2)
    index_map = MOI.Utilities.IndexMap()
    for i in 1:3
        MOI.Utilities.add_column(A)
        index_map[x[i]] = x[i]
    end
    MOI.Utilities.allocate_terms(A, index_map, f1)
    MOI.Utilities.allocate_terms(A, index_map, f2)
    MOI.Utilities.set_number_of_rows(A, 2)
    MOI.Utilities.load_terms(A, index_map, f1, 0)
    MOI.Utilities.load_terms(A, index_map, f2, 1)
    MOI.Utilities.final_touch(A)
    B = convert(SparseArrays.SparseMatrixCSC{Float64,Int}, A)
    @test B == [1.0 1.0 1.0; 0.0 2.0 3.0]
    @test A.rowval == [1, 1, 2, 1, 2]
    @test A.nzval == [1.0, 1.0, 2.0, 1.0, 3.0]
    @test A.colptr == [1, 2, 4, 6]
end

end

TestSparseMatrix.runtests()