text:
"""
    pzeros(M, N; fast = false, atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the finite and infinite zeros of the matrix pencil `M-λN` in `val`, the multiplicities of infinite zeros in `iz` and
information on the  Kronecker-structure in the `KRInfo` object.

The information on the multiplicities of infinite zeros is provided in the vector `iz`,
where each `i`-th element `iz[i]` is equal to `k-1`, where `k` is the order of an infinite elementary divisor with `k > 0`.
The number of infinite zeros contained in `val` is the sum of the components of `iz`.

The information on the  Kronecker-structure consists of the right Kronecker indices `rki`, left Kronecker indices `lki`,
infinite elementary divisors `id`, the number of finite eigenvalues `nf`, and the normal rank `nrank` and can be obtained
from `KRInfo` as `KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`pkstruct`](@ref).

The computation of the zeros is performed by reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF)
exhibiting the spliting of the infinite and finite eigenvalue structures of the pencil `M-λN`.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguements `atol1`, `atol2`  and `rtol` specify the absolute tolerance for the nonzero
elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
"""
function pzeros(M::AbstractMatrix, N::Union{AbstractMatrix,Nothing,UniformScaling{Bool}}; fast::Bool = true,
   atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
   rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2)))

   mM, nM = size(M)
   if N === nothing
      r = rank(M,atol = atol1, rtol = rtol)
      return eltype(M)[], Int[], KRInfo(zeros(Int,nM-r), zeros(Int,mM-r), zeros(Int,0), 0, r)
   end

   if isequal(N,I)
      mM == nM || error("M must be a square matrix")
      return eigvalsnosort(M), Int[], KRInfo(Int[], Int[], Int[], mM, mM)
   end

   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast, withQ = false, withZ = false)

   maxmn = max(mM,nM)
   rki = Vector{Int}(undef,maxmn)
   id = Vector{Int}(undef,maxmn)
   lki = Vector{Int}(undef,maxmn)

   tol1 = max(atol1, rtol*opnorm(M,1))
   mrinf = 0
   nrinf = 0
   niz = 0
   i = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      ired = mrinf+1:mM
      jred = nrinf+1:nM
      τ, ρ = _preduce1!(n, m, p, view(M1,ired,jred), view(N1,ired,jred), Q, Z, tol1;
                        fast = fast, withQ = false, withZ = false)
      i += 1
      mui = ρ+τ
      rki[i] = m - mui
      id[i]  = τ
      mrinf += mui
      nrinf += m
      i > 1 && (niz += τ*(i-1))
      n -= ρ
      m = ρ
      p -= τ
   end
   r = mrinf + n

   rtrail = 0
   ctrail = 0
   j = 0
   while p > 0
      # Step 3: Particular case of the dual PREDUCE algorithm
      ired = mrinf+1:mM-rtrail
      jred = nrinf+1:nM-ctrail
      ρ = _preduce4!(n, 0, p, view(M1,ired,jred),view(N1,ired,jred), Q, Z, tol1, fast = fast, withQ = false, withZ = false)
      j += 1
      rtrail += p
      ctrail += ρ
      lki[j] = p - ρ
      n -= ρ
      p = ρ
   end
   if1 = mrinf+1:mrinf+n
   jf1 = nrinf+1:nrinf+n
   return [eigvalsnosort!(view(M1,if1,jf1),view(N1,if1,jf1)); Inf*ones(real(eltype(M1)),niz) ], minf(id[2:i]),
           KRInfo(kroni(rki[1:i]), kroni(lki[1:j]), minf(id[1:i]), n, r)
end
"""
    peigvals(M, N; fast = false, atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> (val, KRInfo)

Return the (finite and infinite) eigenvalues of the matrix pencil `M-λN` in `val` and
information on the  Kronecker-structure in the `KRInfo` object.

The information on the  Kronecker-structure consists of the right Kronecker indices `rki`, left Kronecker indices `lki`,
infinite elementary divisors `id`, the number of finite eigenvalues `nf`, and the normal rank `nrank` and can be obtained
from `KRInfo` as `KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`pkstruct`](@ref).

The computation of the eigenvalues
is performed by reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF) exhibiting the spliting
of the infinite and finite eigenvalue structures of the pencil `M-λN`.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguements `atol1`, `atol2`  and `rtol` specify the absolute tolerance for the nonzero
elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
"""
function peigvals(M::AbstractMatrix, N::Union{AbstractMatrix,Nothing}; fast::Bool = true,
                  atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
                  rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2)))

   mM, nM = size(M)
   if N === nothing
      r = rank(M,atol = atol1, rtol = rtol)
      return eltype(M)[], KRInfo(zeros(Int,nM-r), zeros(Int,mM-r), zeros(Int,0), 0, r)
   end

   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast, withQ = false, withZ = false)

   maxmn = max(mM,nM)
   rki = Vector{Int}(undef,maxmn)
   id = Vector{Int}(undef,maxmn)
   lki = Vector{Int}(undef,maxmn)

   tol1 = max(atol1, rtol*opnorm(M,1))
   mrinf = 0
   nrinf = 0
   ni = 0
   i = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      ired = mrinf+1:mM
      jred = nrinf+1:nM
      τ, ρ = _preduce1!(n, m, p, view(M1,ired,jred), view(N1,ired,jred), Q, Z, tol1;
                        fast = fast, withQ = false, withZ = false)
      i += 1
      mui = ρ+τ
      rki[i] = m - mui
      id[i]  = τ
      mrinf += mui
      nrinf += m
      ni += τ*i
      n -= ρ
      m = ρ
      p -= τ
   end

   r = mrinf + n

   rtrail = 0
   ctrail = 0
   j = 0
   while p > 0
      # Step 3: Particular case of the dual PREDUCE algorithm
      ired = mrinf+1:mM-rtrail
      jred = nrinf+1:nM-ctrail
      ρ = _preduce4!(n, 0, p, view(M1,ired,jred),view(N1,ired,jred), Q, Z, tol1, fast = fast, withQ = false, withZ = false)
      rtrail += p
      ctrail += ρ
      j += 1
      lki[j] = p - ρ
      n -= ρ
      p = ρ
   end
   if1 = mrinf+1:mrinf+n
   jf1 = nrinf+1:nrinf+n
   return [eigvalsnosort!(view(M1,if1,jf1),view(N1,if1,jf1)); Inf*ones(real(eltype(M1)),ni) ],
          KRInfo(kroni(rki[1:i]), kroni(lki[1:j]), minf(id[1:i]), n, r)
end
"""
    KRInfo

Kronecker-structure object definition.

If `info::KRInfo` is the Kronecker-structure object, then:

`info.rki` is a vector whose components contains the right Kronecker indices;

`info.lki` is a vector whose components contains the left Kronecker indices;

`info.id` is a vector whose components contains the orders of the infinite elementary divisors (i.e., the
multiplicities of infinite eigenvalues).

`info.nf` is the number of finite eigenvalues.

`info.nrank` is the normal rank.

Destructuring via iteration produces the components `info.rki`, `info.lki`, `info.id`, `info.nf` and `info.nrank`.
"""
mutable struct KRInfo
   rki::Vector{Int}
   lki::Vector{Int}
   id::Vector{Int}
   nf::Int
   nrank::Int
   function KRInfo(rki,lki,id,nf,nrank)
      (any(rki .< 0) || any(lki .< 0) || any(id .< 0) || nf < 0 || nrank < 0) && error("no negative components allowed")
      new(rki,lki,id,nf,nrank)
   end
end
# iteration for destructuring into components
Base.iterate(info::KRInfo) = (info.rki, Val(:lki))
Base.iterate(info::KRInfo, ::Val{:lki}) = (info.lki, Val(:id))
Base.iterate(info::KRInfo, ::Val{:id}) = (info.id, Val(:nf))
Base.iterate(info::KRInfo, ::Val{:nf}) = (info.nf, Val(:nrank))
Base.iterate(info::KRInfo, ::Val{:nrank}) = (info.nrank, Val(:done))
#Base.iterate(info::KRInfo, ::Val{:done}) = nothing
"""
    pkstruct(M, N; fast = false, atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> KRInfo

Determine the Kronecker-structure information of the matrix pencil `M-λN` and return an `KRInfo` object.

The right Kronecker indices `rki`, left Kronecker indices `lki`, infinite elementary divisors `id`, the
number of finite eigenvalues `nf` and the normal rank `nrank` can be obtained from
`KRInfo` as `KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
The determination of the Kronecker-structure information is performed by reducing the pencil `M-λN` to an
appropriate Kronecker-like form (KLF) exhibiting all structural elements of the pencil `M-λN`.
The reduction is performed using orthogonal similarity transformations and involves rank decisions based
on rank revealing QR-decompositions with column pivoting, if `fast = true`, or, the more reliable,
SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The right Kronecker indices are provided in the integer vector `rki`, where each `i`-th element `rki[i]` is
the row dimension `k` of an elementary Kronecker block of size `k x (k+1)`.
The number of elements of `rki` is the dimension of
the right nullspace of the pencil `M-λN` and their sum is the degree of a right minimal polynomial nullspace basis.

The left Kronecker indices are provided in the integer vector `lki`, where each `i`-th element `lki[i]` is
the column dimension `k` of an elementary Kronecker block of size `(k+1) x k`.
The number of elements of `lki` is the dimension of
the left nullspace of the pencil `M-λN` and their sum is the degree of a left minimal polynomial nullspace basis.

The multiplicities of infinite eigenvalues are provided in the integer vector `id`, where each `i`-th element `id[i]` is
the order of an infinite elementary divisor (i.e., the multiplicity of an infinite eigenvalue).
The number of elements of `id` is the number of infinite elementary divisors of the pencil `M-λN`.

The normal rank `nrank` is the maximum rank of `M-λN` for all values of `λ`. If `N = nothing`, then `nrank` is the rank of `M`.

The keyword arguements `atol1`, `atol2`  and `rtol` specify the absolute tolerance for the nonzero
elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
"""
function pkstruct(M::AbstractMatrix, N::Union{AbstractMatrix,Nothing}; fast = false, atol1::Real = zero(eltype(M)), atol2::Real = zero(eltype(M)),
   rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2)))

   mM, nM = size(M)
   if N === nothing
      r = rank(M,atol = atol1, rtol = rtol)
      return KRInfo(zeros(Int,nM-r), zeros(Int,mM-r), zeros(Int,0), 0, r)
   end
   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast, withQ = false, withZ = false)

   maxmn = max(mM,nM)
   rki = Vector{Int}(undef,maxmn)
   id = Vector{Int}(undef,maxmn)
   lki = Vector{Int}(undef,maxmn)

   tol1 = max(atol1, rtol*opnorm(M,1))
   mrinf = 0
   nrinf = 0
   ni = 0
   i = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      ired = mrinf+1:mM
      jred = nrinf+1:nM
      τ, ρ = _preduce1!(n, m, p, view(M1,ired,jred), view(N1,ired,jred), Q, Z, tol1;
                        fast = fast, withQ = false, withZ = false)
      i += 1
      mui = ρ+τ
      rki[i] = m - mui
      id[i]  = τ
      mrinf += mui
      nrinf += m
      n -= ρ
      m = ρ
      p -= τ
   end

   r = mrinf + n

   rtrail = 0
   ctrail = 0
   j = 0
   while p > 0
      # Step 3: Particular case of the dual PREDUCE algorithm
      ired = mrinf+1:mM-rtrail
      jred = nrinf+1:nM-ctrail
      ρ = _preduce4!(n, 0, p, view(M1,ired,jred),view(N1,ired,jred), Q, Z, tol1, fast = fast, withQ = false, withZ = false)
      j += 1
      rtrail += p
      ctrail += ρ
      lki[j] = p - ρ
      n -= ρ
      p = ρ
   end
   return KRInfo(kroni(rki[1:i]), kroni(lki[1:j]), minf(id[1:i]), n, r)
end
function kroni(ind)
   # Kronecker indices evaluated from the ranks of blocks of the staircase form
   k = findlast(!iszero,ind)
   k === nothing && (return ind[1:0])
   ni = sum(ind[1:k])
   ki = similar(ind,ni)
   ii = 0
   for i = 1:k
       iip = ii+ind[i]
       for j = ii+1:iip
         ki[j] = i-1
       end
       ii = iip
   end
   return ki
end
function minf(ind)
   # multiplicities of infinite eigenvalues evaluated from the blocksizes of the staircase form
   k = findlast(!iszero,ind)
   k === nothing && (return ind[1:0])
   ni = sum(ind[1:k])
   mi = similar(ind,ni)
   ii = 0
   for i = 1:k
       iip = ii+ind[i]
       for j = ii+1:iip
         mi[j] = i
       end
       ii = iip
   end
   return mi
end
"""
    prank(M::AbstractMatrix, N::AbstractMatrix; fastrank = true, atol1::Real=0, atol2::Real=0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ)

Compute the normal rank of a matrix pencil `M - λN`. If `fastrank = true`, the rank is evaluated
by counting how many singular values of `M - γ N` have magnitude greater than `max(max(atol1,atol2), rtol*σ₁)`,
where `σ₁` is the largest singular value of `M - γ N` and `γ` is a randomly generated value. If `fastrank = false`,
the rank is evaluated as `nr + ni + nf + nl`, where `nr` and `nl` are the sums of right and left Kronecker indices,
respectively, while `ni` and `nf` are the number of infinite and finite eigenvalues, respectively. The sums `nr+ni` and
`nf+nl`, are determined from an appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures
of the pencil `M - λN`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`. The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
For efficiency purpose, the reduction to the relevant KLF is only partially performed
using rank decisions based on rank revealing SVD-decompositions.

!!! compat "Julia 1.1"
    The use of `atol` and `rtol` keyword arguments in rank determinations requires at least Julia 1.1.
    To enforce compatibility with Julia 1.0, the newer function rank in Julia 1.1 has been explicitly included.
"""
function prank(M::AbstractMatrix, N::Union{AbstractMatrix,Nothing}; fastrank::Bool = true,
   atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
   rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2)))

   mM, nM = size(M)
   if N === nothing
      return rank(M, atol = atol1, rtol = rtol)
   end

   (mM,nM) == size(N) || throw(DimensionMismatch("M and N must have the same dimensions"))

   if fastrank
      T = promote_type(eltype(M), eltype(N))
      T <: BlasFloat || (T = promote_type(Float64,T))
      M1 = copy_oftype(M,T)
      N1 = copy_oftype(N,T)
      nrmM = opnorm(M1,1)
      nrmM == zero(nrmM) && (return rank(N1, atol = atol2, rtol = rtol))
      nrmN = opnorm(N1,1)
      nrmN == zero(nrmN) && (return rank(M1, atol = atol1, rtol = rtol))
      scale = nrmM/nrmN*rand()
      return rank(M1+scale*N1, atol = max(atol1,atol2), rtol = rtol)
   else
      # Step 0: Reduce to the standard form
      M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = false, withQ = false, withZ = false)

      n == min(mM,nM) && (return n)
      tol1 = max(atol1, rtol*opnorm(M,1))
      mrinf = 0
      nrinf = 0
      while m > 0
         # Steps 1 & 2: Standard algorithm PREDUCE
         ired = mrinf+1:mM
         jred = nrinf+1:nM
         τ, ρ = _preduce1!(n,m,p,view(M1,ired,jred),view(N1,ired,jred),Q,Z,tol1; fast = false, withQ = false, withZ = false)
         mrinf += ρ+τ
         nrinf += m
         n -= ρ
         m = ρ
         p -= τ
      end
      return mrinf+n
   end
end


text:
"""
Attach an MOI.ModelLike to an optimizer,
solve it and retrieve the termination status and objective value
"""
function solve_abstract_model(
    model::MOI.ModelLike,
    optimizer_constructor,
) where {T}
    JuMP_model = JuMP.Model()
    MOI.copy_to(JuMP.backend(JuMP_model), model)
    set_optimizer(JuMP_model, optimizer_constructor)
    optimize!(JuMP_model)
    termination_status = JuMP.termination_status(JuMP_model)
    obj_val = JuMP.objective_value(JuMP_model)
    return termination_status, obj_val
end

"""
Test if strong duality holds for a problem
"""
function test_strong_duality(
    primal_model::MOI.ModelLike,
    dual_model::MOI.ModelLike,
    optimizer_constructor,
    atol::Float64,
    rtol::Float64,
)
    primal_term_status, primal_obj_val =
        solve_abstract_model(primal_model, optimizer_constructor)
    dual_term_status, dual_obj_val =
        solve_abstract_model(dual_model, optimizer_constructor)

    if primal_term_status == dual_term_status == MOI.OPTIMAL
        return isapprox(primal_obj_val, dual_obj_val; atol = atol, rtol = rtol)
    elseif (primal_term_status == MOI.INFEASIBLE) &&
           (dual_term_status == MOI.DUAL_INFEASIBLE)
        return true
    elseif (primal_term_status == MOI.DUAL_INFEASIBLE) &&
           (dual_term_status == MOI.INFEASIBLE)
        return true
    elseif (primal_term_status == MOI.DUAL_INFEASIBLE) &&
           (dual_term_status == MOI.INFEASIBLE_OR_UNBOUNDED)
        return true
    elseif (primal_term_status == MOI.INFEASIBLE_OR_UNBOUNDED) &&
           (dual_term_status == MOI.DUAL_INFEASIBLE)
        return true
    end
    return false # In case strong duality doesn't hold
end

function test_strong_duality(
    primal_problems::Array{Function},
    optimizer_constructor;
    atol = 1e-6,
    rtol = 1e-4,
)
    for primal_problem in primal_problems
        dual_problem = dualize(primal_problem())
        @testset "$primal_problem" begin
            @test test_strong_duality(
                primal_problem(),
                dual_problem.dual_model,
                optimizer_constructor,
                atol,
                rtol,
            )
        end
    end
end

ext: jl


text:
# run benchmarks on integration scenarios

BenchmarkTools.DEFAULT_PARAMETERS.seconds = timelimit
Base.CoreLogging.disable_logging(Base.CoreLogging.Info) #remove info-level outputs to avoid console spam
@warn "Setting up integration benchmarks..."

#%% general parameters
dp_est = 10. #psi
error_tolerance = 0.1 #psi
outlet_pressure = 220 - 14.65 #WHP in psig
oil_API = 35
sg_gas = 0.65
CO2 = 0.005
sg_water = 1.07
roughness = 0.0006500
BHT = 165
id = 2.441

#%% scenario parameters
scenarios = (:A, :B, :C, :D, :E)

parameters =    (rate = (A = 500, B = 250, C = 1000, D = 3000, E = 50),
                 WC = (A = 0.5, B = 0.25, C = 0.75, D = 0.85, E = 0.25),
                 GLR = (A = 4500, B = 6000, C = 3000, D = 1200, E = 10000),
                 WHT = (A = 100, B = 90, C = 105, D = 115, E = 80))

#%% load test Wellbore
testpath = joinpath(dirname(dirname(pathof(PressureDrop))), "test/testdata/Sawgrass_9_32")
surveypath = joinpath(testpath, "Test_survey_Sawgrass_9.csv")
testwell = read_survey(path = surveypath, id = id)

#%% generate linear temperature profiles
function create_temps(scenario)
        WHT = parameters[:WHT][scenario]

        return linear_wellboretemp(WHT = WHT, BHT = BHT, wellbore = testwell)
end

temps = [create_temps(s) for s in scenarios]
temp_profiles = NamedTuple{scenarios}(temps) #temp profiles labelled by scenario

@warn "Benchmarking Beggs & Brill..."
corr = BeggsAndBrill
timings = Array{Float64,1}(undef, length(scenarios))
for (index, scenario) in enumerate(scenarios)
    WC = parameters[:WC][scenario]
    rate = parameters[:rate][scenario]
    q_o = rate * (1 - WC)
    q_w = rate * WC
    GLR = parameters[:GLR][scenario]
    temps = temp_profiles[scenario]

    timings[index] = @belapsed begin traverse_topdown(wellbore = testwell, roughness = roughness, temperatureprofile = $temps,
                                pressurecorrelation = corr, dp_est = dp_est, error_tolerance = error_tolerance,
                                q_o = $q_o, q_w = $q_w, GLR = $GLR, APIoil = oil_API, sg_water = sg_water, sg_gas = sg_gas,
                                WHP = outlet_pressure, molFracCO2 = CO2)
                        end

end

@warn "$corr timing | min: $(minimum(timings)) s | max: $(maximum(timings)) s | mean: $(sum(timings)/length(timings)) s"

@warn "Benchmarking Hagedorn & Brown..."
corr = HagedornAndBrown
timings = Array{Float64,1}(undef, length(scenarios))
for (index, scenario) in enumerate(scenarios)
    WC = parameters[:WC][scenario]
    rate = parameters[:rate][scenario]
    q_o = rate * (1 - WC)
    q_w = rate * WC
    GLR = parameters[:GLR][scenario]
    temps = temp_profiles[scenario]

    timings[index] = @belapsed begin traverse_topdown(wellbore = testwell, roughness = roughness, temperatureprofile = $temps,
                                pressurecorrelation = corr, dp_est = dp_est, error_tolerance = error_tolerance,
                                q_o = $q_o, q_w = $q_w, GLR = $GLR, APIoil = oil_API, sg_water = sg_water, sg_gas = sg_gas,
                                WHP = outlet_pressure, molFracCO2 = CO2)
                        end
end

@warn "$corr timing | min: $(minimum(timings)) s | max: $(maximum(timings)) s | mean: $(sum(timings)/length(timings)) s"


text:
using Test
using CartesianMesh
using PolynomialBasis
# using Revise
using CutCell
include("useful_routines.jl")

basis = TensorProductBasis(2,1)
mesh = UniformMesh([0.,0.],[2.,1.],[2,2])
cellmaps = CutCell.cell_maps(mesh)

nodalcoordinates = CutCell.dg_nodal_coordinates(cellmaps,basis.points)
testnodalcoordinates = [0.  0.  1.  1.  0.  0.  1.  1.  1.  1.  2.  2.  1.  1.  2.  2.
                        0.  .5  0.  .5  .5  1.  .5  1.  0.  .5  0.  .5  .5  1.  .5  1.]
@test allapprox(nodalcoordinates,testnodalcoordinates)

nodalconnectivity = CutCell.dg_nodal_connectivity(4,4)
testnodalconnectivity = [1  5  9   13
                         2  6  10  14
                         3  7  11  15
                         4  8  12  16]
@test allequal(nodalconnectivity,testnodalconnectivity)

dgmesh = CutCell.DGMesh([0.,0.],[2.,1.],[2,2],basis)


text:
using Distributed
@time begin
    include("required.jl")
    using DecisionTree
    include(joinpath("lib", "model-evaluation.jl"))
end

## Load data
spe = CSV.read(joinpath("data", "proc", "distributions_spe.csv"); header=true, delim="\t")
spa = CSV.read(joinpath("data", "proc", "distributions_spa.csv"); header=true, delim="\t")
env = CSV.read(joinpath("data", "proc", "distributions_env.csv"); header=true, delim="\t")
var = hcat(env, spa)

## Create training and validation sets, with partitioning propotional to class frequency
import MLJ.partition
function partition_by_class(values, fraction; shuffle=true, rng=42)
    # Get indices for each class
    idx0 = findall(iszero, values)
    idx1 = findall(isone, values)

    # Partition each class into trn, vld
    idx0_trn, idx0_vld = partition(idx0, fraction; shuffle=shuffle, rng=rng)
    idx1_trn, idx1_vld = partition(idx1, fraction; shuffle=shuffle, rng=rng)

    # Combine trn, vld partitions
    idx_trn = vcat(idx0_trn, idx1_trn)
    idx_vld = vcat(idx0_vld, idx1_vld)
    return (trn=idx_trn, vld=idx_vld)
end

# Create empty Arrays
trn_idx = Array{Array{Int64,1},1}()
vld_idx = Array{Array{Int64,1},1}()
# Collect training & validation indices
spe = spe[:, Not(60)] # sp60 has no observations...
@time for sp in eachcol(spe)
    part = partition_by_class(sp, 0.7)
    push!(trn_idx, part.trn)
    push!(vld_idx, part.vld)
end

# Get training and validation sets
function partition_labels_features(labels, features, indices)
    labels_set = Array{Array{Int64,1},1}()
    features_set = Array{Array{Float64,2},1}()
    for i in 1:ncol(labels)
        push!(labels_set, Int64.(labels[indices[i], i]))
        push!(features_set, Array(features[indices[i], :]))
    end
    return labels_set, features_set
end
trn_labels, trn_features = partition_labels_features(spe, var, trn_idx);
vld_labels, vld_features = partition_labels_features(spe, var, vld_idx);

## Train Random Forests for all species
# Set parameters
ntrees = 1000
nsubfeatures = -1
nfolds = 3

# Pair corresponding labels and features
trn_pairs = Pair.(trn_labels, trn_features)
vld_pairs = Pair.(vld_labels, vld_features)

# Train models
@time models = [build_forest(lab, feat, nsubfeatures, ntrees) for (lab, feat) in trn_pairs];
# @save joinpath("..", "julia-rf.jld2") models
# @load joinpath("..", "julia-rf.jld2") models

# Apply models
@time predictions = apply_forest.(models, vld_features)
@time predictions_proba = [
    apply_forest_proba(m, f, [0, 1]) for (m, f) in Pair.(models, vld_features)
]

## Evaluate model performance
# Get accuracy measures
@time acc_mes = accuracy_measures.(predictions, vld_labels)
# Get AUC measures
@time auc_stats = auc.(models, vld_features, vld_labels)
# Combine evaluation measures
res = DataFrame(;
    spe=string.("sp", eachindex(acc_mes)),
    freq_abs=Int.(map(sum, eachcol(spe))),
    accuracy=map(x -> x.accuracy, acc_mes),
    sensitivity=map(x -> x.sensitivity, acc_mes),
    specificity=map(x -> x.specificity, acc_mes),
    kappa=map(x -> x.kappa, acc_mes),
    auc=map(x -> x.score, auc_stats),
    plot=map(x -> x.plot, auc_stats),
)

# Explore measures
auc_plot(1 .- res.specificity, res.sensitivity, NaN)
quantile(res.accuracy)
quantile(filter(!isnan, res.sensitivity))
quantile(filter(!isnan, res.specificity))

# Explore measures for models with > 100 obs only
res100 = filter(x -> x.freq_abs > 100, res)
quantile(res100.accuracy)
quantile(filter(!isnan, res100.sensitivity))
quantile(filter(!isnan, res100.specificity))

# Combine all ROC curves in single plot
p = auc_plot([0], [0], NaN; dpi=200)
for auc in auc_stats
    plot!(p, auc.FP_rates, auc.sensitivities)
end
p
savefig(p, "fig/random-forests/08_rf_auc.png")

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using RRSTExperiments
using Test

@testset "RRSTExperiments.jl" begin
    # Write your tests here.
end


text:
"""
```julia
C_R2V(tau::Float64, delta::Float64)
```
Computes the transformation matrix from the R-frame to the V-frame.
"""
function C_R2V(tau::Float64, delta::Float64)
  return [[-sin(delta)*cos(tau) -sin(delta)*sin(tau) cos(delta)];
         [-sin(tau) cos(tau) 0];
         [-cos(delta)*cos(tau) -cos(delta)*sin(tau) -sin(delta)]]
end

"""
```julia
C_V2R(tau::Float64, delta::Float64)
```
Computes the transformation matrix from the V-frame to the R-frame.
"""
function C_V2R(tau::Float64, delta::Float64)
  return [[-sin(delta)*cos(tau) -sin(tau) -cos(delta)*cos(tau)];
          [-sin(delta)*sin(tau) cos(tau) -cos(delta)*sin(tau)];
          [cos(delta) 0 -sin(delta)]]
end

"""
```julia
C_I2R(theta::Float64)
```
Computes the transformation matrix from the I-frame to the R-frame using the
Greenwich Hour Angle in radians.
"""
function C_I2R(theta::Float64)
  return [[cos(theta) sin(theta) 0];
         [-sin(theta) cos(theta) 0];
         [0 0 1]]
end

"""
```julia
C_R2I(theta::Float64)
```
Computes the transformation matrix from the R-frame to the I-frame using the
Greenwich Hour Angle in radians.
"""
function C_R2I(theta::Float64)
  return C_I2R(-theta)
end

"""
```julia
C_V2I(GMST::T, tau::B, delta::B) where {T <: Real, B <: Real}
```
Computes the transformation matrix from the V-frame to the I-frame using the
Greenwich Hour Angle in radians, and the longitude and latitude in radians.
"""
function C_V2I(GMST::T, tau::B, delta::B) where {T <: Real, B <: Real}
  return C_V2R(tau, delta)*C_R2I(GMST)
end

"""
```julia
VelR = VelocityI2R(VelI::Vector{Float64}, PosI::Vector{Float64}, ω::Float64,
  CI2R::Array{Float64}(3,3))
```
Computes the velocity vector in the R-frame from the I-frame.
"""
function VelocityI2R(VelI::Vector{T}, PosI::Vector{P}, ω::S, CI2R::Array{R,2}) where {T,P,S,R <: Real}
  return CI2R*(VelI + [ω*PosI[2], -ω*PosI[1], 0])
end

"""
```julia
PosR = PositionI2R(PosI::Vector{Float64}, CI2R::Array{Float64}(3,3))
```
Computes the position vector in the R-frame from the I-frame.
"""
function PositionI2R(PosI::Vector{Float64}, CI2R::Array{Float64,2})
  return CI2R*PosI
end

"""
```julia
C_V2B(VehState::Vector{Float64}, GMST::Float64)
```
Computes the transformation matrix from the V-frame to the B-frame, using the
vehicle state and the GMST as inputs.
"""
function C_V2B(VehState::Vector{Float64}, GMST::Float64)
  CI2R = C_I2R(GMST)
  PosSph = PositionCart2Sph(CI2R*VehState[4:6])
  return transpose(C_B2I(VehState[10:13])*(CI2R*C_R2V(PosSph[2], PosSph[3])))
end

"""
```julia
function C_B2I(quatB2I::Vector{Float64})
```
Computes the transformation matrix from the B-frame to the I-frame, using the
the quaternion describing the same orientation.
"""
function C_B2I(quatB2I::Vector{Float64})
  return quat2DCM(quatB2I)
end

"""
```julia
C_B2V(VehState::Vector{Float64}, GMST::Float64)
```
Computes the transformation matrix from the B-frame to the V-frame, using the
vehicle state and the GMST as inputs.
"""
function C_B2V(VehState::Vector{Float64}, GMST::Float64)
  return transpose(C_V2B(VehState, GMST))
end

"""
```julia
C_V2TA(γ_a::Float64, χ_a::Float64)
```
Computes the transformation matrix from the V-frame to the airspeed based
trajectory frame.
"""
function C_V2TA(γ_a::T, χ_a::S) where {T,S <: Real}
  return [[cos(χ_a)*cos(γ_a) sin(χ_a)*cos(γ_a) -sin(γ_a)];
            [-sin(χ_a) cos(χ_a) 0.0];
            [cos(χ_a)*sin(γ_a) sin(χ_a)*sin(γ_a) cos(γ_a)]]
end

"""
```julia
C_B2TA(C_B2V::Array{Float64,2}, C_V2TA::Array{Float64,2})
```
Computes the transformation matrix from the B-frame to the airspeed based
trajectory frame.
"""
function C_B2TA(C_B2V::Array{Float64,2}, C_V2TA::Array{Float64,2})
  return C_B2V*C_V2TA
end

"""
```julia
[α, β, σ] =
  AirAttitude(C_B2V::Array{Float64,2}, C_V2TA::Array{Float64,2})
```
Computes the angle of attack, sideslip angle, bank angle, and the body to
airspeed based trajectory reference frame transformation matrix.
"""
function AirAttitude(C_B2V::Array{T,2}, C_V2TA::Array{S,2}) where {T,S <: Real}
  CB2TA = C_B2TA(C_B2V, C_V2TA)
  return [atan2(CB2TA[1,3],CB2TA[1,1]), asin(CB2TA[1,2]),
    atan2(-CB2TA[3,2], CB2TA[2,2])]
end

"""
```julia
C_B2AA(α::Float64, β::Float64)
```
Computes the transformation matrix from the B-frame to the airspeed based
aerodynamic frame
"""
function C_B2AA(α::Float64, β::Float64)
  return [[cos(α)*cos(β) sin(β) sin(α)*cos(β)];
          [-cos(α)*sin(β) cos(β) -sin(α)*sin(β)];
          [-sin(α) 0.0 cos(α)]]
end

"""
```julia
C_AA2B(α::Float64, β::Float64)
```
Computes the transformation matrix from the airspeed based
aerodynamic frame to the B-frame
"""
function C_AA2B(α::Float64, β::Float64)
  return transpose(C_B2AA(α,β))
end

#=
"""
```julia
PosRsph, C_AA2I, C_V2I, V_VA =
    TransformUtils(state::Vector{Float64}, simtime::Float64, ω::Float64 = 0.)
```
Computes the transformation matrix from the airspeed based
aerodynamic frame to the I-frame
"""
function TransformUtils(state::Vector{Float64}, simt::Float64, ω::Float64 = 0.)
  GMST = calc_GMST(MJD(simt))
  CI2R = C_I2R(GMST)
  VR = VelocityI2R(state[1:3], state[4:6], ω, CI2R)
  PosR = CI2R*state[4:6]
  PosRSph = PositionCart2Sph(PosR)
  CR2V = C_R2V(PosRSph[2], PosRSph[3])
  VV = CR2V*VR
  Vsph = VelocityCart2Sph(VV)
  Angles = AirAttitude(C_B2V(state, GMST), C_V2TA(Vsph[2], Vsph[3]))
  return PosRSph, C_B2I(state[10:13])*C_AA2B(Angles[1], Angles[2]), transpose(CI2R*CR2V),
            Vsph[1]
end
=#


text:
export UnaryOperation, BinaryOperation, NaryOperation
export AbstractOperation, Gradient
export to_expr

import Base: +, *, -, √, tanh, sin, cos, tan, ^, exp, convert, promote_rule

# Unary Operators, (name, symbol)
unary_operators = []
push!(unary_operators, ["Negative", "-"])
push!(unary_operators, ["SquareRoot", "√"])
push!(unary_operators, ["Tanh", "tanh"])
push!(unary_operators, ["Sin", "sin"])
push!(unary_operators, ["Cos", "cos"])
push!(unary_operators, ["Tan", "tan"])
push!(unary_operators, ["Exp", "exp"])

# Binary Operators, (name, symbol)
binary_operators = []
push!(binary_operators, ["Add", "+"])
push!(binary_operators, ["Multiply", "*"])
push!(binary_operators, ["Exponentiation", "^"])

nary_operators = []
push!(nary_operators, ["Sum", "+"])
push!(nary_operators, ["Product", "*"])

abstract type AbstractOperation <: AbstractExpression end
abstract type UnaryOperation  <: AbstractOperation end
abstract type BinaryOperation <: AbstractOperation end
abstract type NaryOperation <: AbstractOperation end

# Define to_expression function
to_expr(x) = x
# Define Struct and Symbol Overload for Unary Operators
for unary_operator in unary_operators
    b_name, b_symbol = Meta.parse.(unary_operator)
    @eval struct $b_name{𝒯} <: UnaryOperation
        term::𝒯
    end
    @eval $b_symbol(a::AbstractExpression) = $b_name(a)
    @eval to_expr(a::$b_name) = Expr(:call, $b_symbol, to_expr(a.term))
end

# Define Struct and Symbol Overload for Binary Operators
for binary_operator in binary_operators
    b_name, b_symbol = Meta.parse.(binary_operator)
    @eval struct $b_name{𝒯, 𝒮} <: BinaryOperation
        term1::𝒯
        term2::𝒮
    end
    @eval $b_symbol(a::AbstractExpression, b::AbstractExpression) = $b_name(a, b)
    @eval to_expr(a::$b_name) = Expr(:call, $b_symbol, to_expr(a.term1), to_expr(a.term2))
end

# Define Struct and Symbol Overload for n-ary Operators
for nary_operator in nary_operators
    # Defining the nary name and symbol
    n_name, n_symbol = Meta.parse.(nary_operator)

    # Defining the struct, along with outer constructor
    @eval struct $n_name{a} <: NaryOperation
        terms::Vector{a}
    end

    # Linking the defined symbol to a cuntion that creates the struct
    @eval $n_symbol(a::Vector{AbstractExpression}) = $n_name(a...)
end

# Special Structs
struct Gradient{𝒯, 𝒰} <: AbstractExpression
    operand::𝒯
    metadata::𝒰
end

# Special Rules
-(a::AbstractExpression, b::AbstractExpression) = a + -b


text:
#=
check_permutation_count:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-26
=#

module CountPermuChecker
    # parameters
    const max_size = 256

    # methods
    function check_permu_count(str1, str2)
        println("")
        println("Check $(str1) and $(str2) are Permutaton or not")

        # initialize count array
        count1 = [0 for i in 1:max_size]
        count2 = [0 for i in 1:max_size]

        # increment count of each character
        for s in str1
            count1[Int(s)] += 1
        end
        for s in str2
            count2[Int(s)] += 1
        end

        # string length is not same
        # can not be permutation
        if length(str1) != length(str2)
            println("Not Permutation")
            return false
        end

        # compare count array
        for i in 1:max_size
            if count1[i] != count2[i]
                println("Not Permutation")
                return false
            end
        end

        println("Permutation")
        return true
    end

    function main()
        @time println(check_permu_count("test", "sett"))
        @time println(check_permu_count("geeksforgeeks", "forgeeksgeeks"))
        @time println(check_permu_count("abcdefg", "aabbccd"))
        @time println(check_permu_count("hogehoge", "hogeehoge"))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .CountPermuChecker
    CountPermuChecker.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
"""
http://www.trackvis.org/docs/?subsect=fileformat

`id_string::NTuple{6,UInt8}`: ID string for track file. The first 5 characters must be \"TRACK\".
`dim::NTuple{3,Int16}`: Dimension of the image volume
`voxel_size::NTuple{3,Float32}`: Voxel size of the image volume
`origin::Ntuple{3,Float32}`: Origin of the image volume. This field is not yet
                             being used by TrackVis. That means the origin is
                             always (0, 0, 0).
`n_scalars::Int16`: Number of scalars saved at each track point (besides x, y and z coordinates).
`scalar_name::NTuple{10,NTuple{20,UInt8}}`: Name of each scalar. Can not be longer than 20 characters each. Can only store up to 10 names.
`n_properties::UInt16`: Number of properties saved at each track.
`property_name::NTuple{10,NTuple{20,UInt8}}`: Name of each property. Can not be longer than 20 characters each. Can only store up to 10 names.
`vox_to_ras::Matrix{Float32}`: 4x4 matrix for voxel to RAS (crs to xyz) transformation. If vox_to_ras[3][3] is 0, it means the matrix is not recorded. This field is added from version 2.
`reserved[444]::NTuple{444,UInt8}`: Reserved space for future version.
`voxel_order::NTuple{4,UInt8}`: Storing order of the original image data. Explained here.
`pad2::NTuple{2,UInt8}`: Paddings
`image_orientation_patient::NTuple{6,Float32}`: Image orientation of the original image. As defined in the DICOM header.
`pad1::NTuple{2,UInt8}`: Paddings
`invert_x::UInt8`,`invert_y::UInt8`,`invert_x::UInt8`,`swap_xy::UInt8`,`swap_yz::UInt8`,`swap_zx::UInt8`: Inversion/rotation flags used to generate this track file. For internal use only.
`n_count::Int32`: Number of tracks stored in this track file. 0 means the number was NOT stored.
`version::Int32`: Version number. Current version is 2.
`hdr_size::Int32`: Size of the header. Used to determine byte swap. Should be 1000.

141531076

f = "../tmp/ABCDStudy/test/data/DTI1_PA1_AP1/streamlines/full_interface_prob_pft.trk"
io = open(f)
s = load(io, TrkStream);
mytracks = load(s, VecPath);

trk = load(s, PathMeta);

"""
struct TrkHeader
    dim::NTuple{3,Int16}
    voxel_size::NTuple{3,Float32}
    origin::NTuple{3,Float32}
    nscalars::Int16
    scalar_name::NTuple{10,NTuple{20,UInt8}}
    nproperties::UInt16
    property_name::NTuple{10,NTuple{20,UInt8}}
    vox_to_ras::NTuple{4,NTuple{4,Float32}}
    reserved::NTuple{444,UInt8}
    voxel_order::NTuple{4,UInt8}
    pad2::NTuple{4,UInt8}
    image_orientation_patient::NTuple{6,Float32}
    pad1::NTuple{2,UInt8}
    invert_x::UInt8
    invert_y::UInt8
    invert_z::UInt8
    swap_xy::UInt8
    swap_yz::UInt8
    swap_zx::UInt8
    ncount::Int32
    version::Int32
    hdr_size::Int32
end

@generated function read_header(io)
    out = Expr(:new, TrkHeader)
    for p in fieldtypes(TrkHeader)
        if p <: Tuple
            t = Expr(:tuple)
            for p_i in p.parameters
                if p_i <: Tuple
                    t_i = Expr(:tuple)
                    for p_j in p_i.parameters
                        push!(t_i.args, :(read(io, $p_j)))
                    end
                    push!(t.args, t_i)
                else
                    push!(t.args, :(read(io, $p_i)))
                end
            end
            push!(out.args, t)
        else
            push!(out.args, :(read(io, $p)))
        end
    end
    return out
end

@generated function read_header_swap(io)
    out = Expr(:new, TrkHeader)
    for p in fieldtypes(TrkHeader)
        if p <: Tuple
            t = Expr(:tuple)
            for p_i in p.parameters
                if p_i <: Tuple
                    t_i = Expr(:tuple)
                    for p_j in p_i.parameters
                        push!(t_i.args, :(bswap(read(io, $p_j))))
                    end
                    push!(t.args, t_i)
                else
                    push!(t.args, :(bswap(read(io, $p_i))))
                end
            end
            push!(out.args, t)
        else
            push!(out.args, :(bswap(read(io, $p))))
        end
    end
    return out
end

read_trk(file::AbstractString; mode="r") = read_trk(open(file, mode))
function read_trk(io)
    b1, b2, b3, b4, b5, b6 = read(io, 6)
    if (b1, b2, b3, b4, b5) === (0x54, 0x52, 0x41, 0x43, 0x4b)
        return read_header(io)
    elseif (b6, b5, b4, b3, b2) === (0x4b, 0x43, 0x41, 0x52, 0x54)
        return read_header_swap(io)
    else
        error("This is not a Track file.")
    end
end

function load(io::IO, ::Type{TrkStream})
    needswap = checkfile(io)
    hdr = read(io, TrkHeader, needswap)
    TrkStream(hdr, io, needswap)
end

function load(s::TrkStream, ::Type{VecPoint})
    nscalars = convert(Int64, s.hdr.nscalars)
    if s.needswap
        VecPoint((bswap.(read!(s.io, Vector{Float32}(undef, 3)))...,),
                 (bswap.(read!(s.io, Vector{Float32}(undef, nscalars)))...,))
    else
        VecPoint((read!(s.io, Vector{Float32}(undef, 3))...,),
                 (read!(s.io, Vector{Float32}(undef, nscalars))...,))
    end
end

"""
#C: Path count
S: Size of a path
N: Dimensions
Tp: Type (e.g., Int, Float32) representing position of each VecPoint
Ns: Number of scalars in a VecPoint
Ts: Type (e.g., Int, Float32) of scalars in a VecPoint

#  Path{3,Float32,nscalars,Float32}
"""

function load(s::TrkStream, ::Type{VecPath})
    nprops = s.hdr.nproperties
    ncount = s.hdr.ncount
    nscalars = convert(Int64, s.hdr.nscalars)
    #vecpointtype = VecPath{3,Float32,nscalars,Float32}

    vecdata = Vector{Path}(undef, ncount);
    props = Dict{String,Any}();

    propnames = []
    for i in 1:nprops
        append!(propnames, s.hdr.property_name[i])
        if propnames[i] == 0
            props[String("property$i")] = []
        else
            props[String([propnames[i]...])] = []
        end
    end

    # load each path
    # i = each pathway
    # p = each point
    # j = each property
    for pathway_i in 1:ncount
        npoints = read!(s.io, Vector{Int32}(undef, 1))[1]

        if s.needswap
            npoints = bswap(npoints)
        end

        if s.needswap

            data[pathway_i] = [bswap.([load(s, VecPoint) for p in OneTo(npoints)])...]
        else
            vecdata[pathway_i] = [(load(s, VecPoint) for p in OneTo(npoints))...]
        end

        for j in 1:nprops
            props[propnames[j]][pathway_i] = read!(s.io, Vector{Float32}(undef, 1))[1]
        end
    end

    if s.needswap
        props = bswap.(props)
    end
    return (vecdata, props)  # this should be ImageMeta
end

#### FINAL

@inline function read_trk_point_swap(io, nscalars::Integer)
    return PointNode(
        map(bswap, read_type_swap(io, NTuple{3,Float32})),
        if nscalars == 0
            Float32[]
        else
            for i in OneTo(nscalars)
            [bswap(read(io, Float32)) for i in ]
        end
    )
end

@inline function read_trk_point_noswap(io, nscalars::Integer)
    return PointNode(
        read_type_noswap(io, NTuple{3,Float32}),
        if nscalars == 0
            Float32[]
        else
            [read(io, Float32) for i in OneTo(nscalars)]
        end
    )
end

function read_tract_noswap(io, nscalars::Integer, nprops::Integer)
    npoints = bswap(read(io, Int32))

    out = NodeFiber(Vector{TrkPoint}(undef, npoints), Vector{Float32}(undef, nprops))
    for i in OneTo(npoints)
        @inbounds out[i] = read_trk_point_noswap(io, nscalars)
    end

    for i in OneTo(nprops)
        @inbounds out.metadata[i] = read(io, Float32)
    end
    return out
end

function read_tract_noswap(io, nscalars::Integer, nprops::Integer)
    npoints = read(io, Int32)
    out = Tract(Vector{TractPoint}(undef, npoints), Vector{Float32}(undef, nprops))
    for i in OneTo(npoints)
        out.points[i] = read_trk_point_noswap(io, nscalars)
    end

    for i in OneTo(nprops)
        out.points[i] = read(io, Float32)
    end
    return out
end

function readtrk_swap(io)
end

function readtrk_noswap(io)
    hdr = NamedTuple{TRK_FIELDS,TRK_TUPLE}(read_type_noswap(io, TRK_TUPLE))
    propnames = map(OneTo(hdr.nproperties)) do i
        n_i = hdr.property_name[i]
        if n_i == 0
            Symbol("property_$i")
        else
            Symbol(n_i)
        end
    end

    return BraintTracts(
        hdr,
        [read_tract_noswap(io, hdr.nscalars, hdr.nproperties)],
        propnames
    )
end

function readtrk(io::IO)
    magic = read_type_noswap(io, NTuple{5,UInt8})
    if magic === MAGIC_BYTES
        hdr = NamedTuple{TRK_FIELDS,TRK_TUPLE}(read_tract_noswap(io, TRK_TUPLE))
        propnames = map(OneTo(hdr.nproperties)) do i
            n_i = hdr.property_name[i]
            if n_i == 0
                Symbol("property_$i")
            else
                Symbol(n_i)
            end
        end

        fibers = Vector{TrkFiber}(undef, npoints)

        for fiber_i in OneTo(hdr.ncount)
            @inbounds fibers[fibers_i] = NodeFiber(Vector{TrkPoint}(undef, npoints), Vector{Float32}(undef, nprops))
            npoints = read(io, Int32)
            for point_i in OneTo(npoints)
                @inbounds fibers[fibers_i][point_i] = PointNode(read_type_noswap(io, NTuple{3,Float32}), read!(io, Vector{Float32}(undef, hdr.nscalars)))
            end
            props = properties(@inbounds(fibers[fibers_i]))
            for property_i in OneTo(hdr.nproperties)
                @inbounds props[property_i] = read(io, Float32)
            end
        end
        return BraintTracts(hdr, fibers, propnames)
    elseif map(bswap, magic) === MAGIC_BYTES
        return readtrk_swap(io)
    else
        error("Not a track file.")
    end
end

@inline function read_trk_point_swap(io, nscalars::Integer)
    return PointNode(
        map(bswap, read_type_swap(io, NTuple{3,Float32})),
        if nscalars == 0
            Float32[]
        else
            for i in OneTo(nscalars)
            [bswap(read(io, Float32)) for i in ]
        end
    )
readtrk(trkfile) = readtrk(open(trkfile))


text:
import KnetLayers: arrtype, Activation, Filtering

abstract type Model end
struct ResNet <: Model; w; end
function (M::ResNet)(m,imgurl::String,avgimg;stage=3)
    img = imgdata(imgurl, avgimg)
    return M(w,m,arrtype(img);stage=stage);
end
function ResNet(atype::Type;stage=3)
    w,m,meta = ResNetLib.resnet101init(;trained=true,stage=stage)
    global avgimg = meta["normalization"]["averageImage"]
    global descriptions = meta["classes"]["description"]
    return w,m,meta,avgimg
end

ResNet() = ResNet(nothing);

struct CNN <: Model
    layer1::Filtering{typeof(conv4)}
    layer2::Filtering{typeof(conv4)}
    drop::Dropout
end
function (m::CNN)(x)
    x1 = m.layer1(m.drop(x))
    x2 = m.layer2(m.drop(x1))
    h,w,c,b = size(x2)
    permutedims(reshape(x2,h*w,c,b),(2,3,1))
end
CNN(h::Int,w::Int,c::Int,d::Int) = CNN(Conv(height=h,width=w,inout=c=>d,padding=1,activation=ELU()),
                                       Conv(height=h,width=w,inout=d=>d,padding=1,activation=ELU()),
                                       Dropout(0.18))

struct mRNN  <: Model
    rnn::LSTM
end

function (m::mRNN)(x;batchSizes=[1])
    B = first(batchSizes)
    if last(batchSizes)!=B
        out = m.rnn(x;batchSizes=batchSizes,hy=true,cy=false)
    else
        x   = reshape(x,size(x,1),B,div(size(x,2),B))
        out = m.rnn(x;hy=true,cy=false)
    end
    return out.y, out.hidden
end
mRNN(input::Int,hidden::Int;o...) = mRNN(LSTM(input=input, hidden=hidden;o...))

struct QUnit  <: Model
    embed::Embed
    rnn::mRNN
    linear::Linear
    drop1::Dropout
    drop2::Dropout
end
function (m::QUnit)(x;batchSizes=[1],train=false)
    xe = m.drop1(m.embed(x))
    y,hyout = m.rnn(xe;batchSizes=batchSizes)
    q  = m.drop2(vcat(hyout[:,:,1],hyout[:,:,2]))
    B = batchSizes[1]

    if ndims(y) == 2
        indices      = bs2ind(batchSizes)
        lngths       = length.(indices)
        Tmax         = maximum(lngths)
        td,B         = size(q)
        d            = div(td,2)
        cw           = Any[];
        for i=1:length(indices)
            y1 = y[:,indices[i]]
            df = Tmax-lngths[i]
            if df > 0
                cpad = zeros(Float32,2d*df) # zeros(Float32,2d,df)
                kpad = arrtype(cpad)
                ypad = reshape(cat1d(y1,kpad),2d,Tmax) # hcat(y1,kpad)
                push!(cw,ypad)
            else
                push!(cw,y1)
            end
        end
        cws_2d =  reshape(vcat(cw...),2d,B*Tmax)
    else
        d      = div(size(y,1),2)
        Tmax   = size(y,3)
        cws_2d = reshape(y,2d,B*Tmax)
    end
    cws_3d =  reshape(m.linear(cws_2d),(d,B,Tmax))
    return q,cws_3d;
end
QUnit(vocab::Int,embed::Int,hidden::Int;bidir=true) = QUnit(Embed(input=vocab, output=embed; winit=rand),
                                                            mRNN(embed,hidden;bidirectional=bidir),
                                                            Linear(input=2hidden,output=hidden),
                                                            Dropout(0.15), Dropout(0.08))

function bs2ind(batchSizes)
    B = batchSizes[1]
    indices = Any[]
    for i=1:B
        ind = i.+cumsum(filter(x->(x>=i),batchSizes)[1:end-1])
        push!(indices,append!(Int[i],ind))
    end
    return indices
end

struct Control  <: Model
    cq::Linear
    att::Linear
end
function (m::Control)(c,q,cws,pad;train=false,tap=nothing)
      d,B,T = size(cws)
      cqi   = reshape(m.cq(vcat(c,q)),(d,B,1))
      cvis  = reshape(cqi .* cws,(d,B*T))
      cvis_2d = reshape(m.att(cvis),(B,T)) #eq c2.1.2
      if pad != nothing
          cvi = reshape(softmax(cvis_2d .- pad,dims=2),(1,B,T)) #eq c2.2
      else
          cvi = reshape(softmax(cvis_2d,dims=2),(1,B,T)) #eq c2.2
      end
      tap!=nothing && get!(tap,"w_attn_$(tap["cnt"])",Array(reshape(cvi,B,T)))
      cnew = reshape(sum(cvi.*cws;dims=3),(d,B))
end
Control(d::Int) = Control(Linear(input=2d,output=d),Linear(input=d,output=1))

struct Read  <: Model
    me::Linear
    Kbe::Linear
    Kbe2::Linear
    Ime
    att::Linear
    drop::Dropout
end

function (m::Read)(mp,ci,cws,KBhw′,KBhw′′;train=false,tap=nothing)
    d,B,N = size(KBhw′); BN = B*N
    mi_3d = reshape(m.me(mp),(d,B,1))
    ImKB  = reshape(mi_3d .* KBhw′,(d,BN)) # eq r1.2
    ImKB′ = reshape(elu.(m.Ime*ImKB .+ KBhw′′),(d,B,N)) #eq r2
    ci_3d = reshape(ci,(d,B,1))
    IcmKB_pre = elu.(reshape(ci_3d .* ImKB′,(d,BN))) #eq r3.1.1
    IcmKB_pre = m.drop(IcmKB_pre)
    IcmKB = reshape(m.att(IcmKB_pre),(B,N)) #eq r3.1.2
    mvi = reshape(softmax(IcmKB,dims=2),(1,B,N)) #eq r3.2
    tap!=nothing && get!(tap,"KB_attn_$(tap["cnt"])",Array(reshape(mvi,B,N)))
    mnew = reshape(sum(mvi.*KBhw′;dims=3),(d,B)) #eq r3.3
end
Read(d::Int) = Read(Linear(input=d,output=d),Linear(input=d,output=d),
                    Linear(input=d,output=d),param(d,d; atype=arrtype, init=xavier),
                    Linear(input=d,output=1), Dropout(0.15))

struct Write  <: Model
    me::Linear
    cproj::Union{Linear,Nothing}
    att::Union{Linear,Nothing}
    mpp
    gating::Union{Linear,Nothing}
end

function (m::Write)(m_new,mi₋1,mj,ci,cj;train=false,selfattn=true,gating=true,tap=nothing)
    d,B        = size(m_new)
    T          = length(mj)
    mi         = m.me(vcat(m_new,mi₋1))
    !selfattn && return mi
    ciproj     = m.cproj(ci)
    ci_3d      = reshape(ciproj,d,B,1)
    cj_3d      = reshape(cat1d(cj...),(d,B,T)) #reshape(hcat(cj...),(d,B,T)) #
    sap        = reshape(ci_3d.*cj_3d,(d,B*T)) #eq w2.1.1
    sa         = reshape(m.att(sap),(B,T)) #eq w2.1.2
    sa′        = reshape(softmax(sa,dims=2),(1,B,T)) #eq w2.1.3
    mj_3d      = reshape(cat1d(mj...),(d,B,T)) #reshape(hcat(mj...),(d,B,T)) #
    mi_sa      = reshape(sum(sa′ .* mj_3d;dims=3),(d,B))
    mi′′       = m.mpp*mi_sa .+ mi #eq w2.3
    !gating && return mi′′
    σci′       = sigm.(m.gating(ci))  #eq w3.1
    mi′′′      = (σci′ .* mi₋1) .+  ((1 .- σci′) .* mi′′) #eq w3.2
end

function Write(d::Int;selfattn=true,gating=true)
    if selfattn
        if gating
            Write(Linear(input=2d,output=d),Linear(input=d,output=d),
                  Linear(input=d,output=1),param(d,d;atype=arrtype, init=xavier),
                  Linear(input=d,output=1))
        else
            Write(Linear(input=2d,output=d),Linear(input=d,output=d),
                  Linear(input=d,output=1),param(d,d;atype=arrtype, init=xavier),
                  nothing)
        end
    else
        Write(Linear(input=2d,output=d),nothing,nothing,nothing,nothing)
    end
end

struct MAC <: Model
    control::Control
    read::Read
    write::Write
end
function (m::MAC)(qi,cws,mi,mj,ci,cj,KBhw′,KBhw′′,pad;train=false,selfattn=true,gating=true,tap=nothing)
    cnew = m.control(ci,qi,cws,pad;train=train,tap=tap)
    ri   = m.read(mi,ci,cws,KBhw′,KBhw′′;train=train,tap=tap)
    mnew = m.write(ri,mi,mj,ci,cj;train=train,selfattn=selfattn,gating=gating)
    return cnew,mnew
end
MAC(d::Int;selfattn=false,gating=false) = MAC(Control(d),Read(d),Write(d))

struct Output <: Model
    qe::Dense
    l1::Dense
    l2::Linear
end

(m::Output)(q,mp) = m.l2(m.l1(cat(m.qe(q),mp;dims=1)))

Output(d::Int) = Output(Dense(input=2d,output=d,activation=ELU()),
                        Dense(input=2d,output=d,activation=ELU()),
                        Linear(input=d,output=28))

struct MACNetwork <: Model
    resnet::ResNet
    cnn::CNN
    qunit::QUnit
    qindex::Linear
    mac::MAC
    output::Output
    c0
    m0
    drop::Dropout
end

function (M::MACNetwork)(qs,batchSizes,xS,xB,xP;answers=nothing,p=12,selfattn=false,gating=false,tap=nothing,allsteps=false)
    train         = answers!=nothing
    #STEM Processing
    KBhw          = M.cnn(xS)
    #Read Unit Precalculations
    d,B,N         = size(KBhw)
    KBhw_2d       = M.drop(reshape(KBhw,(d,B*N)))
    KBhw′_pre     = M.mac.read.Kbe(KBhw_2d) # look if it is necessary
    KBhw′′        = M.mac.read.Kbe2(KBhw′_pre)
    KBhw′         = reshape(KBhw′_pre,(d,B,N))

    #Question Unit
    q,cws         = M.qunit(qs;batchSizes=batchSizes,train=train)
    qi_c          = M.qindex(q)
    #Memory Initialization
    ci            = M.c0*xB
    mi            = M.m0*xB

    if selfattn
        cj=[ci]; mj=[mi]
    else
        cj=nothing; mj=nothing
    end

    for i=1:p
        qi        = qi_c[(i-1)*d+1:i*d,:]
        ci = M.drop(ci);
        mi = M.drop(mi);
        ci,mi = M.mac(qi,cws,mi,mj,ci,cj,KBhw′,KBhw′′,xP;train=train,selfattn=selfattn,gating=gating,tap=tap)
        if selfattn; push!(cj,ci); push!(mj,mi); end
        tap!=nothing && (tap["cnt"]+=1)
    end

    y = M.output(q,mi)

    if answers==nothing
        predmat = convert(Array{Float32},y)
        tap!=nothing && get!(tap,"y",predmat)
        predictions = mapslices(argmax,predmat,dims=1)[1,:]
        if allsteps
            outputs = []
            for i=1:p-1
                yi = M.output(q,mj[i])
                yi = convert(Array{Float32},yi)
                push!(outputs,mapslices(argmax,yi,dims=1)[1,:])
            end
            push!(outputs,predictions)
            return outputs
        end
        return predictions
    else
        return nll(y,answers)
    end
end

function MACNetwork(o::Dict)
           MACNetwork(ResNet(),
                      CNN(3,3,1024,o[:d]),
                      QUnit(o[:vocab_size],o[:embed_size],o[:d]),
                      Linear(input=2*o[:d],output=o[:p]*o[:d]),
                      MAC(o[:d];selfattn=o[:selfattn],gating=o[:gating]),
                      Output(o[:d]),
                      param(o[:d],1;atype=arrtype, init=xavier),
                      param(o[:d],1;atype=arrtype, init=randn),
                      Dropout(0.15))
end

setoptim!(m::MACNetwork,o) =
    for param in params(m); param.opt = Adam(;lr=o[:lr]); end

function benchmark(M::MACNetwork,feats,o;N=10)
    getter(id) = view(feats,:,:,:,id)
    B=32;L=25
    @time for i=1:N
        ids  = randperm(128)[1:B]
        xB   = arrtype(ones(Float32,1,B))
        xS   = arrtype(batcher(map(getter,ids)))
        xQ   = [rand(1:84) for i=1:B*L]
        answers = [rand(1:28) for i=1:B]
        batchSizes = [B for i=1:L]
        xP   = nothing
        y    = @diff M(xQ,batchSizes,xS,xB,xP;answers=answers,p=o[:p],selfattn=o[:selfattn],gating=o[:gating])
    end
end
function benchmark(feats,o;N=30)
    M     = MACNetwork(o);
    benchmark(M,feats,o;N=N)
end

const feats_L = 200704;
const feats_H = 14;
const feats_C = 1024;
function batcher1(feats,args)
    B = length(args)
    totlen = feats_L*B
    result = similar(Array{Float32}, totlen)
    starts = (0:B-1) .*feats_L .+ 1; ends = starts .+ feats_L .- 1;
    for i=1:B
        result[starts[i]:ends[i]] = view(feats,:,:,:,args[i])
    end
    return reshape(result,feats_H,feats_H,feats_C,B)
end

function inplace_batcher(result,data,args)
     B = length(args)
     totlen = feats_L*B
     starts = (0:B-1) .* feats_L .+ 1; ends = starts .+ feats_L .- 1;
     for i=1:B
          result[starts[i]:ends[i]] = view(data,:,:,:,args[i])
     end
     return reshape(result,feats_H,feats_H,feats_C,B)
end


text:
module CooperativeGames

    using IterTools
    using Base.Iterators
    using Combinatorics
    using LightGraphs

    include("types.jl")
    include("utils.jl")
    include("games.jl")
    include("graphs.jl")
    include("solutions/harsanyi.jl")
    include("solutions/shapley.jl")
    include("solutions/banzhaf.jl")
    include("solutions/core.jl")
    include("solutions/fairness.jl")

    export Players, SimpleGame, GraphGame,

        # Methods
        isconvex, isfair, graphtoMyerson,

        # Solutions
        Δₕ, fₛ, fₛⁱ, fᵦⁱ, degreeHarsanyi

end


text:
using EarthSphere
myhome = Earth_LatLong_Coord(40.753,73.9725,"NW")
mywork = Earth_LatLong_Coord(40.690,73.965,"NW")
dist = get_distance(myhome,mywork)
angle = EarthSphere.haversine(myhome,mywork)
println(dist)
println(angle)


text:
function is_convex_mesh_oriented_outward(eM::eMesh{Tri,T2}) where {T2}
	for k = 1:n_tri(eM)
		p3 = vertex_pos_for_tri_ind(eM, k)
		n̂ = triangleNormal(p3)
		c = centroid(p3)
		c = normalize(c)
	end
	return true
end

@testset "swept_mesh" begin
	# constant radius
	eM_box = create_swept_mesh(f_swept_triv, [-1.0, 1.0], 1.0, 4, true)
	@test is_convex_mesh_oriented_outward(eM_box)
	@test 8.0 ≈ volume(eM_box)
	@test verify_mesh(eM_box) == nothing

	# variable radius
	eM_box_2 = create_swept_mesh(f_swept_triv, [-1.0, 1.0], [1.0, 1.0], 4, true)
	@test is_convex_mesh_oriented_outward(eM_box_2)
	@test 8.0 ≈ volume(eM_box_2)
	@test verify_mesh(eM_box_2) == nothing

	for field_k = fieldnames(typeof(eM_box))
		@test getfield(eM_box, field_k) == getfield(eM_box_2, field_k)
	end

	eM_box_check = eMesh_box()
	@test is_convex_mesh_oriented_outward(eM_box_check)
	@test 8.0 ≈ volume(eM_box_check)
end