text:
# MDP problem tools and value iteratation
module mdp
using DelimitedFiles
export State, DetState, Policy, RandomPolicy

# Type definitions ############################################################
struct State
  a::Array{Int, 1}
  a2r::Array{Float64, 1}
  Pa2ns::Array{Float64, 2}
  cPa2ns::Array{Float64, 2}
  ns::Array{Int, 1}
end

struct DetState
  a::Array{Int, 1}
  a2r::Array{Float64, 1}
  ns::Array{Int, 1}
end

mutable struct Policy
  S::Union{Dict{Int, State}, Dict{Int, DetState}}
  U::Dict{Int, Float64}
  Aidx::Dict{Int, Int}
  drate::Float64
end

function Policy(S::Union{Dict{Int, State}, Dict{Int, DetState}},
                drate::Float64=1.0)
  U = Dict{Int, Float64}()
  sizehint!(U, length(S))

  Aidx = Dict{Int, Int}()
  sizehint!(Aidx, length(S))

  for (id, s) in S
    (U[id], Aidx[id]) = findmax(s.a2r)
  end

  return Policy(S, U, Aidx, drate)
end

function RandomPolicy(S::Union{Dict{Int, State}, Dict{Int, DetState}},
                      drate::Float64=1.0)
  U = Dict{Int, Float64}()
  sizehint!(U, length(S))

  Aidx = Dict{Int, Int}()
  sizehint!(Aidx, length(S))

  for (id, s) in S
    Aidx[id] = rand(1:length(s.a))
    U[id] = s.a2r[Aidx[id]]
  end

  return Policy(S, U, Aidx, drate)
end
###############################################################################

# Utility Functions ###########################################################
function optimal_policy(S::Dict{Int, Union{State, DetState}}, drate::Float64)
  P = Policy(S, drate)

  ret = 1.0
  k = 40
  while ret > 0.0 || k > 0
    ret = iterate!(P)
    println(ret)
    k = ret == 0.0 ? k - 1 : k
  end

  return P
end

function read_data(name::AbstractString)::Array{Int, 2}
  # read the data and variable names in
  full_data = readdlm("data/" * name * ".csv", ',')
  D = Int.(full_data[2:end, :])

  return D
end

function write_data(name::AbstractString, D::AbstractArray{Int, 2})
  # read the data and variable names in
  full_data = [["s" "a" "r" "ns"]; D]
  writedlm("data/" * name * ".csv", full_data, ',')
end

function MLE(D::AbstractArray{Int, 2})
  # count the number of states and actions, next states for each state -------#
  us = unique(D[:, 1])
  s2ns = Dict{Int, Set{Int}}()
  s2a = Dict{Int, Set{Int}}()
  for s in us
    s2ns[s] = Set{Int}()
    s2a[s] = Set{Int}()
  end

  for i in 1:size(D, 1)
    sid = D[i, 1]
    a = D[i, 2]
    r = D[i, 3]
    ns = D[i, 4]
    push!(s2ns[sid], ns)
    push!(s2a[sid], a)
  end
  #---------------------------------------------------------------------------#

  # maximum likelihood estimate of state transistions and reward --------------#
  counts = Dict{Int, Array{Int, 2}}()
  rewards_total = Dict{Int, Array{Float64, 1}}()
  rewards_count = Dict{Int, Array{Int, 1}}()
  A = Dict{Int, Array{Int, 1}}()
  NS = Dict{Int, Array{Int, 1}}()
  for id in keys(s2ns)
    A[id] = sort(collect(s2a[id]))
    NS[id] = sort(collect(s2ns[id]))

    counts[id] = fill(0, (length(s2a[id]), length(s2ns[id])))
    rewards_total[id] = fill(0.0, length(s2a[id]))
    rewards_count[id] = fill(0, length(s2a[id]))
  end

  for i in 1:size(D, 1)
    sid = D[i, 1]
    a = D[i, 2]
    r = D[i, 3]
    ns = D[i, 4]

    idxa = 1
    for j in 1:length(A[sid])
      idxa = j
      A[sid][j] == a &&  break
    end

    idxns = 1
    for j in 1:length(NS[sid])
      idxns = j
      NS[sid][j] == ns && break
    end

    counts[sid][idxa, idxns] += 1

    rewards_total[sid][idxa] += r
    rewards_count[sid][idxa] += 1
  end
  #---------------------------------------------------------------------------#

  # allocate the model -------------------------------------------------------#
  chain = Dict{Int, State}()
  for id in keys(s2ns)
    a2r = rewards_total[id] ./ rewards_count[id]

    Pa2ns = convert.(Float64, counts[id])
    for i in 1:size(Pa2ns, 1)
      tot = sum(Pa2ns[i, :]) == 0.0 ? 1.0 : sum(Pa2ns[i, :])
      Pa2ns[i, :] /= tot
    end

    cPa2ns = cumsum(Pa2ns, dims=2)

    chain[id] = State(A[id], a2r, Pa2ns, cPa2ns, NS[id])
  end
  #---------------------------------------------------------------------------#

  return chain
end

function determinateMLE(D::AbstractArray{Int, 2})
  states = MLE(D)

  states_det = Dict{Int, DetState}()
  for (id, s) in states
    a = copy(s.a)
    idx = [cart_idx[2] for cart_idx in findmax(s.Pa2ns, dims=2)[2]][:]
    ns = copy(s.ns[idx])

    states_det[id] = DetState(a, s.a2r, ns)
  end

  return states_det
end

function iterate!(p::Policy)
  nbchange = 0
  change = 0.0

  det = typeof(p.S) == Dict{Int, DetState}

  SID = collect(keys(p.S))
  EUr = fill(0.0, 0)
  for k in 1:length(SID)
    sid = SID[k]
    s = p.S[sid]

    if length(EUr) < length(s.a2r)
      EUr = fill(0.0, length(s.a2r))
    end
    for i in 1:length(EUr)
      EUr[i] = 0.0
    end

    for i in 1:length(s.a2r)
      EUr[i] = 0.0
      if det
        if s.ns[i] > 0 # check if action is not final
          EUr[i] = p.drate * p.U[s.ns[i]]
        end
      else
        for j in 1:length(s.ns)
          if s.ns[j] > 0 # check if action is not final
            EUr[i] += p.drate * s.Pa2ns[i, j] * p.U[s.ns[j]]
          else
            println("Action is final")
          end
        end
      end
    end

    oldaidx = p.Aidx[sid]
    oldu = p.U[sid]

    maxAidx = 1
    maxU = -Inf
    for i in 1:length(s.a2r)
      if s.a2r[i] + EUr[i] >= maxU
        maxU = s.a2r[i] + EUr[i]
        maxAidx = i
      end
    end
    #(p.U[sid], p.Aidx[sid]) = findmax(s.a2r + EUr)
    p.U[sid] = maxU
    p.Aidx[sid] = maxAidx

    change += abs(oldu - p.U[sid])
    if oldaidx != p.Aidx[sid]
      nbchange += 1
    end
  end

  return (change, nbchange)
end

function estimate_score(p::Policy, N::Int=10)
  det = typeof(p.S) == Dict{Int, DetState}

  score_tot = 0.0
  for i in 1:N
    score_part = 0.0
    sid = rand(keys(p.S))

    M = p.drate != 1.0 ? ceil(log(eps()) / log(p.drate)) : (1 << 31)
    #M = 100
    for j in 1:M
      aidx = p.Aidx[sid]
      score_part += p.drate^(j - 1) * p.S[sid].a2r[aidx]
      nsid = 1

      if det
        nsid = p.S[sid].ns[aidx]
      else
        r = rand()
        for k in 1:length(p.S[sid].ns)
          if p.S[sid].cPa2ns[aidx, k] > r
            nsid = p.S[sid].ns[k]
            break
          end
        end
      end

      if nsid < 0 # check if action was final
        #println("Stopping after $(j) actions")
        break
      else
        sid = nsid
      end
    end
    #println("####################################################")

    score_tot += score_part
  end

  return score_tot / N
end

function generate_data(name::AbstractString, S::Dict{Int, State}, N::Int=10^5)
  fp = open("data/" * name * ".csv", "w")
  write(fp, "s,a,r,ns\n")

  k = keys(S)
  for i in 1:N
    s = rand(k)
    aidx = rand(1:length(S[s].a))
    a = S[s].a[aidx]
    r = round(Int, S[s].a2r[aidx])
    ns = -1

    rnd = rand()
    for j in 1:length(S[s].ns)
      if S[s].cPa2ns[aidx, j] > rnd
        ns = S[s].ns[j]
        break
      end
    end

    write(fp, "$(s),$(a),$(r),$(ns)\n")
  end

  close(fp)
end

function write_policy(name::AbstractString, p::Policy)
  fp = open("data/" * name * ".policy", "w")

  for s in 1:maximum(keys(p.S))
    if haskey(p.S, s) && haskey(p.Aidx, s)
      write(fp, string(p.S[s].a[p.Aidx[s]]))
    else
      write(fp, string(1))
    end
    write(fp, "\n")
  end

  close(fp)
end
###############################################################################

dir = @__DIR__
include(dir * "/mdp_vit.jl")

end


text:
module Zeros

export Zero, testzero, One, testone

"A type that stores no data, and holds the value zero."
struct Zero <: Integer
end

"A type that stores no data, and holds the value one."
struct One <: Integer
end

const StaticBool = Union{Zero, One}

Base.promote_rule(::Type{One}, ::Type{T}) where {T<:Number} = T
Base.promote_rule(::Type{Zero}, ::Type{T}) where {T<:Number} = T

# For some reason the promote-rules involving StaticBool and Bool need to be
# defined in both directions, or we get stack overflow.
for S in (Zero, One, Bool), T in (Zero, One, Bool)
    if S ≠ T
        Base.promote_rule(::Type{S}, ::Type{T}) = Bool
    end
end

Base.convert(::Type{T}, ::Zero) where {T<:Number} = T(zero(T))
Base.convert(::Type{T}, ::One) where {T<:Number} = oneunit(T)
Base.convert(::Type{Zero}, ::One) = throw(InexactError(:Zero, Zero, One()))
Base.convert(::Type{One}, ::Zero) = throw(InexactError(:One, One, Zero()))

# Some of the more common constructors that do not default to `convert`
# Note:  We cannot have a (::Type{T})(x::StaticBool) where {T<:Number} constructor
# instead of all of these, because of ambiguities with user-defined types.
for T in (Bool, Integer, AbstractFloat,
          BigInt, Int128, Int16, Int32, Int64, Int8,
          UInt128, UInt16, UInt32, UInt64, UInt8,
          BigFloat, Float16, Float32, Float64)
    (::Type{T})(::Zero) = zero(T)
    (::Type{T})(::One) = one(T)
end

#disambig
for T in (:Number, :Complex, :Rational, :BigFloat, :Zero, :One)
    @eval Zero(x::$T) = iszero(x) ? Zero() : throw(InexactError(:Zero, Zero, x))
    @eval One(x::$T) = isone(x) ? One() : throw(InexactError(:One, One, x))
end

Base.zero(::StaticBool) = Zero()
Base.zero(::Type{<:StaticBool}) = Zero()
Base.one(::StaticBool) = One()
Base.one(::Type{<:StaticBool}) = One()

Base.Complex(x::Real, ::Zero) = x

# Loop over types in order to make methods specific enough to avoid ambiguities.
for T in (Number, Real, Rational, Integer, Complex, Complex{Bool})
    Base.:+(x::T, ::Zero) = x
    Base.:+(::Zero, x::T) = x
    Base.:-(x::T, ::Zero) = x
    Base.:-(::Zero, x::T) = -x
    Base.:*(::T, ::Zero) = Zero()
    Base.:*(::Zero, ::T) = Zero()
    Base.:/(::Zero, ::T) = Zero()
    Base.:/(::T, ::Zero) = throw(DivideError())
    Base.:*(x::T, ::One) = x
    Base.:*(::One, x::T) = x
end

# Division sometimes returns a different type from the arguments, e.g. for Int/Int.
Base.:/(x::AbstractFloat, ::One) = x

Base.inv(::One) = One()
Base.inv(::Zero) = throw(DivideError())

# Loop over rounding modes in order to make methods specific enough to avoid ambiguities.
for R in (RoundingMode, RoundingMode{:Down}, RoundingMode{:Up}, Union{RoundingMode{:Nearest}, RoundingMode{:NearestTiesAway}, RoundingMode{:NearestTiesUp}})
    Base.div(x::Integer, ::One, ::R) = x
end

# These functions are intentionally not type-stable.
"Convert to Zero() if equal to zero. (Use immediately before calling a function.)"
testzero(x::Number) = iszero(x) ? Zero() : x
"Convert to One() if equal to one. (Use immediately before calling a function.)"
testone(x::Number) = isone(x) ? One() : x

# This functions give a strange "of_sametype" error. (See int.jl)
# Hence we overload them all, even when result is not One() or Zero()
for op in (:+, :-, :*, :&, :|, :xor)
    for (T1,x) in ((:Zero, 0), (:One, 1))
        for (T2,y) in ((:Zero, 0), (:One, 1))
            val = @eval $op($x,$y)
            @eval Base.$op(::$T1, ::$T2) = $(testzero(testone(val)))
        end
    end
end

Base.:+(::One, x::Bool) = 1 + x
Base.:+(x::Bool, ::One) = x + 1

Base.:-(::Zero) = Zero()
Base.:-(::One) = -1

Base.:/(::Zero, ::Zero) = throw(DivideError())
Base.:/(::One, ::One) = One()

# Disambig
Base.:/(::One, ::Zero) = throw(DivideError())
Base.:/(::Zero, ::One) = Zero()

Base.:<(::T,::T) where {T<:StaticBool} = false
Base.:<=(::T,::T) where {T<:StaticBool} = true

Base.ldexp(::Zero, ::Integer) = Zero()
Base.copysign(::Zero,::Real) = Zero()
Base.flipsign(::Zero,::Real) = Zero()
Base.modf(::Zero) = (Zero(), Zero())

# Alerady working due to default definitions:
# ==, !=, abs, isinf, isnan, isinteger, isreal, isimag,
# signed, unsigned, float, big, complex, isodd, iseven

for op in (:sign, :round, :floor, :ceil, :trunc, :significand, :sqrt, :cbrt)
    @eval Base.$op(x::StaticBool) = x
end

Base.log(::One) = Zero()
Base.log10(::One) = Zero()
Base.log2(::One) = Zero()
Base.log1p(::Zero) = Zero()
Base.exp(::Zero) = One()
Base.exp2(::Zero) = One()
Base.exp10(::Zero) = One()
Base.expm1(::Zero) = Zero()
Base.sin(::Zero) = Zero()
Base.cos(::Zero) = One()
Base.sincos(::Zero) = (Zero(), One())
Base.tan(::Zero) = Zero()
Base.asin(::Zero) = Zero()
Base.atan(::Zero) = Zero()
Base.sinpi(::Zero) = Zero()
Base.sinpi(::One) = Zero()
Base.cospi(::Zero) = One()
Base.sinh(::Zero) = Zero()
Base.cosh(::Zero) = One()
Base.tanh(::Zero) = Zero()
Base.hypot(::Zero, x::Number) = abs(x)
Base.hypot(x::Number, ::Zero) = abs(x)
Base.hypot(::Zero, ::Zero) = Zero()

# ^ has a lot of very specific methods in Base....
for T in (Float16, Float32, Float64, BigFloat, AbstractFloat, Rational, Complex{<:AbstractFloat}, Complex{<:Integer}, Integer, BigInt)
    Base.:^(::T, ::Zero) = One()
end

Base.to_power_type(::Zero) = false

# # Avoid promotion of triplets
# for op in [:fma :muladd]
#     @eval $op(::Zero, ::Zero, ::Zero) = Zero()
#     for T in (Real, Integer)
#         @eval $op(::Zero, ::$T, ::Zero) = Zero()
#         @eval $op(::$T, ::Zero, ::Zero) = Zero()
#         @eval $op(::Zero, x::$T, y::$T) = convert(promote_type(typeof(x),typeof(y)),y)
#         @eval $op(x::$T, ::Zero, y::$T) = convert(promote_type(typeof(x),typeof(y)),y)
#         @eval $op(x::$T, y::$T, ::Zero) = x*y
#         @eval $op(::One, x::$T, y::$T) = x+y
#         @eval $op(x::$T, ::One, y::$T) = x+y
#     end
# end

for op in (:mod, :rem), T in (:Real, :Rational)
  @eval Base.$op(::Zero, ::$T) = Zero()
end
Base.mod(::Zero, ::Zero) = throw(DivideError()) # disambig
Base.rem(::Zero, ::Zero) = throw(DivideError()) #
Base.mod(::One, ::One) = Zero()
Base.rem(::One, ::One) = Zero()

# Sum up vector of One() to Int
Base.reduce_empty(::typeof(Base.add_sum), ::Type{One}) = zero(Int)
Base.reduce_first(::typeof(Base.add_sum), x::One) = Int(x)

# Sum up arrays of Zero() to Zero()
Base.add_sum(::Zero, ::Zero) = Zero()

# Product of vector of Zero() to Bool
Base.reduce_empty(::typeof(Base.mul_prod), ::Type{Zero}) = true
Base.reduce_first(::typeof(Base.mul_prod), ::Zero) = false

# Product of vector of One() to One()
Base.reduce_empty(::typeof(Base.mul_prod), ::Type{One}) = One()

# Multiply up arrays of One() to One()
Base.mul_prod(::One, ::One) = One()

# Multiply up arrays of Zero() to `false` (when length is nonzero)
Base.mul_prod(::Zero, x::Zero) = false
Base.mul_prod(::Zero, x::Bool) = false
Base.mul_prod(::Bool, x::Zero) = false

Base.:(:)(::One, stop::Integer) = Base.OneTo(stop)

# Display Zero() as `𝟎` ("mathematical bold digit zero")
# so that it looks slightly different from `0` (and same for One()).
Base.show(io::IO, ::Zero) = print(io, "𝟎") # U+1D7CE
Base.show(io::IO, ::One) = print(io, "𝟏") # U+1D7CF

Base.string(z::StaticBool) = Base.print_to_string(z)

Base.Checked.checked_abs(x::StaticBool) = x
Base.Checked.checked_mul(x::StaticBool, y::StaticBool) = x*y
Base.Checked.mul_with_overflow(x::StaticBool, y::StaticBool) = (x*y, false)
Base.Checked.checked_add(x::StaticBool, y::StaticBool) = x+y

if VERSION < v"1.2"
    # Disambiguation needed for older Julia versions
    Base.copysign(::Zero, x::Unsigned) = Zero()
    Base.flipsign(::Zero, x::Unsigned) = Zero()
end

if VERSION ≥ v"1.3"
    # @eval, because this needs to parse in Julia 1.0, even though it will not run
    @eval begin
        export $(Symbol("𝟎")), $(Symbol("𝟏"))
        const $(Symbol("𝟎")) = Zero()  # \bfzero <tab>
        const $(Symbol("𝟏")) = One()   # \bfone <tab>
    end
end

include("pirate.jl")

end # module

ext: jl


text:
using ControlSystems, MonteCarloMeasurements
using Optim

# NOTE: for complex structured perturbations the typical upper bound algorithm appears to be quite tight. For real perturbations maybe the interval method has something to offer?

# function Base.:∈(a::Real, b::Complex{<:AbstractParticles})
#     mi, ma = pextrema(b.re)
#     mi <= a <= ma || return false
#     mi, ma = pextrema(b.im)
#     mi <= a <= ma
# end

# # any0det(D::Matrix{<:Complex{<:Interval}}) = 0 ∈ det(D)

# function any0det(D::Matrix{<:Complex{<:Interval}})
#     # CD = complex.(IntervalArithmetic.sup.(real.(D)), IntervalArithmetic.sup.(imag.(D)))
#     # dmid = det(CD)
#     # real(dmid) < 0 && return true
#     0 ∈ det(D)
# end

# """
#     bisect_a(P, K, w; W = (:), Z = (:), tol = 0.001, N = 32, upper = true, δ = δc, σ=-1)
#     bisect_a(M, D, w; W = (:), Z = (:), tol = 0.001, N = 32, upper = true, δ = δc, σ=-1)

# EXPERIMENTAL AND SUBJECT TO BUGS, BREAKAGE AND REMOVAL

# For each frequency in `w`, find the largest `a` such that the loop with uncertainty elements of norm no greater than `a`, located at the inputs `W` and outputs `Z` of `P`, is stable.

# By default, a conservative lower bound on the disk margin is returned. The conservatism comes from multiple factors
# - Each uncertainty element is represented as an IntervalBox in the complex plane rather than a ball with radius `a`.
# - IntervalArithmetic is inherently conservative due to the "dependency problem", i.e., `x-x ≠ 0`.

# If `upper = true`, a Monte-Carlo approach is used to find an upper bound of the disk margin, this will be optimistic. If you calculate both and they are tight, you have a good indication of the true disk margin.

# # Arguments:
# - `P`: Plant
# - `K`: Controller
# - `w`: Frequency vector
# - `W`: The inputs at which to insert a perturbation, defaults to all.
# - `Z`: The outputs at which to insert a perturbation, defaults to all.
# - `tol`: The tolerance in the bisection and the resolution of the resulting `a`.
# - `N`: The number of samples for the upper bound estimation.
# - `upper`: Calculate upper or lower bound?
# - `δ = δc` for complex perturbations and `δ = δr` for real perturbations.
# """
# function bisect_a(M0::AbstractArray, D::AbstractMatrix; tol=1e-3, kwargs...)
#     @views map(axes(M0, 3)) do i
#         M = @views M0[:,:,i]
#         au = opnorm(M)
#         al = eltype(D) <: Complex ? ρ(M) : 0
#         au-al < tol && return al
#         iters = ceil(Int, log2((au-al)/tol))
#         local a
#         for j = 1:iters
#             a = (au+al)/2
#             if any0det(I - (1/a*M)*D) # 0 ∈ det(I - (a*M0[:,:,i])*D)
#                 al = a
#             else
#                 au = a
#             end
#         end
#         1/al
#     end
# end

# function bisect_a(args...; tol=1e-3, kwargs...)
#     M0, D = get_M(args...; kwargs...)
#     bisect_a(M0, D; tol=1e-3, kwargs...)
# end

# """
#     M,D = get_M(P, K, w; W = (:), Z = (:), N = 32, upper = false, δ = δc, σ=-1)

# Return the frequency response of `M` in the `M-Δ` formulation that arises when individual, complex/real perturbations are introduced on inputs `W` and outputs `Z` (defaults to all).

# # Arguments:
# - `P`: System
# - `K`: Controller
# - `w`: Frequency vector
# - `W`: Input indices that ar perturbed
# - `Z`: Output indices that ar perturbed
# - `N`: Number of samples for the upper bound computation
# - `upper`: Indicate whether an upper or lower bound is to be computed
# - `δ = δc` for complex perturbations and `δ = δr` for real perturbations.
# """
# function get_M(P::LTISystem, K::LTISystem, w::AbstractVector; W = (:), Z = (:), N = 32, upper=false, δ = δc, σ = -1)
#     Z1 = W2 = Z == (:) ? (1:P.ny) : Z
#     W1 = Z2 = W == (:) ? (1:P.nu) : W
#     ny,nu = length(Z1), length(W1)
#     D = Δ(ny+nu, δ)
#     if upper
#         D = rand(D, N)
#     else
#         D = Diagonal([Interval(d) for d in diag(D)])
#     end
#     if isempty(Z1) # No outputs
#         L = (K*P)[W1, W1]
#     elseif isempty(W1) # No inputs
#         L = (P*K)[Z1, Z1]
#     else
#         L = [ss(zeros(P.ny, P.ny), P.timeevol) P;-K ss(zeros(K.ny, K.ny), P.timeevol)]
#         L = L[[Z1; P.ny .+ Z2], [W1; P.nu .+ W2]]
#     end
#     n = L.ny
#     X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
#     M = starprod(X,L)
#     # M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
#     M0 = freqresp(M, w).parent
#     M0, D
# end

# function get_M(L::LTISystem, w::AbstractVector; N = 32, upper=false, δ = δc, σ = -1)
#     ny,nu = size(L)
#     D = Δ(ny+nu, δ)
#     if upper
#         D = rand(D, N)
#     else
#         D = Diagonal([Interval(d) for d in diag(D)])
#     end
#     n = L.ny
#     X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
#     M = starprod(X,L)
#     # M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
#     M0 = freqresp(M, w).parent
#     M0, D
# end

# function muloss(M, d)
#     D = Diagonal(d)
#     opnorm(D\M*D)
# end

#=
If there's a single block only, we have
μ = maximum(abs(e) for e in eigvals(M) if e is real) if the block is real*I
μ = ρ(M) if the block is complex*I
μ = opnorm(M) if the block is full complex

=#

"""
    μ = structured_singular_value(M; tol=1e-4, scalings=false, dynamic=false)

Compute (an upper bound of) the structured singular value μ for diagonal Δ of complex perturbations (other structures of Δ are not yet supported).
`M` is assumed to be an (n × n × N_freq) array or a matrix.

We currently don't have any methods to compute a lower bound, but if all perturbations are complex the spectral radius `ρ(M)` is always a lower bound (usually not a good one).

If `scalings = true`, return also a `n × nf` matrix `Dm` with the diagonal scalings `D` such that
```
D = Diagonal(Dm[:, i])
σ̄(D\\M[:,:,i]*D)
```
is minimized.

If `dynamic = true`, the perturbations are assumed to be time-varying `Δ(t)`. In this case, the same scaling is used for all frequencies and the returned `D` if `scalings=true` is a vector `d` such that `D = Diagonal(d)`.
"""
function structured_singular_value(M::AbstractArray{T}; tol=1e-4, scalings=false, dynamic=false) where T
    Ms1 = similar(M[:,:,1])
    Ms2 = similar(M[:,:,1])
    function muloss(M, d)
        D = Diagonal(d)
        mul!(Ms1, M, D)
        ldiv!(Ms2, D, Ms1)
        opnorm(Ms2)
    end
    n = size(M, 2)
    d0 = ones(real(T), n)
    opts = Optim.Options(
        store_trace       = false,
        show_trace        = false,
        show_every        = 10,
        iterations        = 1000,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 0,
        f_abstol          = 0,
        g_tol             = tol,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
    )
    if dynamic
        res = Optim.optimize(
            d->maximum(i->muloss(@view(M[:,:,i]),d), axes(M, 3)),
            d0,
            # i == 1 ? ParticleSwarm() : BFGS(alphaguess = LineSearches.InitialStatic(alpha=0.9), linesearch = LineSearches.HagerZhang()),
            n == 1 ? BFGS() : ParticleSwarm(),
            opts,
        )
        d0 = res.minimizer
        return scalings ? (res.minimum, d0) : res.minimum
    else
        if scalings
            Dm = Matrix{Float64}(undef, n, size(M,3))
        end
        mu = map(axes(M, 3)) do i
            @views M0 = M[:,:,i]
            res = Optim.optimize(
                d->muloss(M0,d),
                d0,
                # i == 1 ? ParticleSwarm() : BFGS(alphaguess = LineSearches.InitialStatic(alpha=0.9), linesearch = LineSearches.HagerZhang()),
                n == 1 ? BFGS() : i == 1 ? ParticleSwarm() : NelderMead(), # Initialize using Particle Swarm
                opts,
            )
            d0 = res.minimizer # update initial guess
            if scalings
                Dm[:, i] .= d0
            end
            res.minimum
        end
        if scalings
            (M isa AbstractMatrix ? mu[] : mu), Dm
        else
            M isa AbstractMatrix ? mu[] : mu
        end
    end
end

"""
    robstab(M0::UncertainSS, w=exp10.(LinRange(-3, 3, 1500)); kwargs...)

Return the robust stability margin of an uncertain model, defined as the inverse of the structured singular value.
Currently, only diagonal complex perturbations supported.
"""
robstab(M0::UncertainSS, args...; kwargs...) = 1/norm(structured_singular_value(M0, args...; kwargs...), Inf)

"""
    structured_singular_value(M0::UncertainSS, [w::AbstractVector]; kwargs...)

- `w`: Frequency vector, if none is provided, the maximum μ over a brid 1e-3 : 1e3 will be returned.
"""
function structured_singular_value(M0::LTISystem, w::AbstractVector; kwargs...)
    if M0 isa UncertainSS
        if length(M0.Δ) == 1 && M0.Δ[] isa δDyn
            return sigma(M0.M, w)[1][:, 1]
        end
        all(d isa δ{<:Complex} for d in M0.Δ) || error("Structured singular value only supported for diagonal complex perturbations")
        M0 = M0.M # not it's own method due to method ambiguity misery
    end
    M = freqresp(M0, w).parent
    μ = structured_singular_value(M; kwargs...)
end

function structured_singular_value(M0::LTISystem; kwargs...)
    w = exp10.(LinRange(-3, 3, 1500))
    μ = structured_singular_value(M0, w; kwargs...)
    maximum(μ)
end

"""
    sim_diskmargin(L, σ::Real, w::AbstractVector)
    sim_diskmargin(L, σ::Real = 0)

Simultaneuous diskmargin at the outputs of `L`.
Uses should consider using [`diskmargin`](@ref).
"""
function sim_diskmargin(L::LTISystem,σ::Real,w::AbstractVector)
    # S̄ = S+(σ-1)/2*I = lft([(1+σ)/2 -1;1 -1], L)
    n = L.ny
    X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
    M = starprod(X,L)
    M0 = freqresp(M, w).parent
    mu = structured_singular_value(M0)
    imu = inv.(structured_singular_value(M0))
    simultaneous = [Diskmargin(imu, σ; ω0 = w, L) for (imu, w) in zip(imu,w)]
end

"""
    sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real = 0)

Simultaneuous diskmargin at both outputs and inputs of `P`.
Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet
https://arxiv.org/abs/2003.04771
See also [`ncfmargin`](@ref).
"""
function sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real=0)
    L = [ss(zeros(P.ny, P.ny)) P;-C ss(zeros(C.ny, C.ny))]
    sim_diskmargin(L,σ)
end

"""
    sim_diskmargin(L, σ::Real = 0)

Return the smallest simultaneous diskmargin over the grid 1e-3:1e3
See also [`ncfmargin`](@ref).
"""
function sim_diskmargin(L, σ::Real=0)
    m = sim_diskmargin(L, σ, LinRange(-3, 3, 500))
    m = argmin(d->d.α, m)
end

"""
    diskmargin(P::LTISystem, C::LTISystem, σ, w::AbstractVector, args...; kwargs...)

Simultaneuous diskmargin at outputs, inputs and input/output simultaneously of `P`.
Returns a named tuple with the fields `input, output, simultaneous_input, simultaneous_output, simultaneous` where `input` and `output` represent loop-at-a-time margins, `simultaneous_input` is the margin for simultaneous perturbations on all inputs and `simultaneous` is the margin for perturbations on all inputs and outputs simultaneously.
See also [`ncfmargin`](@ref).
"""
function diskmargin(P::LTISystem, C::LTISystem, σ, w::AbstractVector, args...; kwargs...)
    L = C*P
    issiso(L) && return diskmargin(L, σ, w, args...)
    input = loop_diskmargin(L, σ, w)
    simultaneous_input = sim_diskmargin(L,σ,w)
    L = P*C
    output = loop_diskmargin(L, σ, w)
    simultaneous_output = sim_diskmargin(L,σ,w)
    te = P.timeevol
    L = [ss(zeros(P.ny, P.ny), te) P;-C ss(zeros(C.ny, C.ny), te)]
    simultaneous = sim_diskmargin(L,σ,w)
    (; input, output, simultaneous_input, simultaneous_output, simultaneous)
end

##

"""
    loop_diskmargin(L, args...; kwargs...)

Calculate the loop-at-a-time diskmargin for each output of `L`.

See also [`diskmargin`](@ref), [`sim_diskmargin`](@ref).
Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet
https://arxiv.org/abs/2003.04771
"""
function loop_diskmargin(L::LTISystem, args...; kwargs...)
    dms = map(1:L.ny) do i
        open_L = broken_feedback(L, i)
        diskmargin(open_L, args...; kwargs...)
    end
end

"""
    loop_diskmargin(P, C, args...; kwargs...)

Calculate the loop-at-a-time diskmargin for each output and input of `P`.
See also [`diskmargin`](@ref), [`sim_diskmargin`](@ref).
Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet
https://arxiv.org/abs/2003.04771
"""
function loop_diskmargin(P::LTISystem,C::LTISystem,args...; kwargs...)
    input = loop_diskmargin(C*P, args...; kwargs...)
    output = loop_diskmargin(P*C, args...; kwargs...)
    (; input, output)
end

"""
    broken_feedback(L, i)
Closes all loops in square MIMO system `L` except for loops `i`.
Forms L1 in fig 14. of "An Introduction to Disk Margins" https://arxiv.org/abs/2003.04771
"""
function broken_feedback(L::LTISystem, i)
    ny, nu = size(L)
    ny == nu || throw(ArgumentError("Only square loop-transfer functions supported"))
    connection_inds = setdiff(1:ny, i)
    i isa AbstractVector || (i = [i])
    open_L = feedback(
        L,
        ss(I(ny-1), L.timeevol), # open one loop
        U1 = connection_inds,
        Y1 = connection_inds,
        Z1 = i,
        W1 = i,
    )
    @assert issiso(open_L)
    open_L
end

"""
    loop_scaling(M0::Matrix, tol = 0.0001)

Find the optimal diagonal scaling matrix `D` such that `D\\M0*D` has a minimized condition number.
Applicable to square `M0` only. See also [`structured_singular_value`](@ref) with option `dynamic=true`.
Use [`loop_scale`](@ref) to find and apply the scaling to a loop-transfer function.
"""
function loop_scaling(M::Matrix, tol=1e-4)

    Ms1 = similar(M)
    Ms2 = similar(M)
    function scaleloss(d)
        D = Diagonal(d)
        mul!(Ms1, M, D)
        ldiv!(Ms2, D, Ms1)
        cond(Ms2)
    end
    n = size(M, 2)
    d0 = ones(real(eltype(M)), n)
    opts = Optim.Options(
        store_trace       = false,
        show_trace        = false,
        show_every        = 10,
        iterations        = 1000,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 0,
        f_abstol          = 0,
        g_tol             = tol,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
    )

    res = Optim.optimize(
        scaleloss,
        d0,
        # i == 1 ? ParticleSwarm() : BFGS(alphaguess = LineSearches.InitialStatic(alpha=0.9), linesearch = LineSearches.HagerZhang()),
        n == 1 ? BFGS() : NelderMead(),
        opts,
    )
    d = res.minimizer
    abs.(d)
end

"""
    loop_scale(L::LTISystem, w = 0)

Find the optimal diagonal scaling matrix `D` such that `D\\L(iw)*D` has a minimized condition number at frequency `w`. Applicable to square `L` only.
Use [`loop_scaling`](@ref) to obtain `D`.
"""
function loop_scale(L::LTISystem, w = 0)
    M = freqresp(L, w)
    D = loop_scaling(M)
    ny = L.ny
    D = Diagonal(D)
    inv(D)*L*D
end

# function loop_scaling(M0::Matrix, tol=1e-4)
#     M = [
#         0I inv(M0)
#         M0     0I
#     ]
#     Ms1 = similar(M)
#     Ms2 = similar(M)
#     function scaleloss(d)
#         d02[1:n] .= d
#         d02[n+1:end] .= d
#         D = Diagonal(d02)
#         mul!(Ms1, M, D)
#         ldiv!(Ms2, D, Ms1)
#         cond(Ms2)
#     end
#     n = size(M0, 2)
#     d0 = ones(real(eltype(M)), n)
#     d02 = ones(real(eltype(M)), 2n)
#     opts = Optim.Options(
#         store_trace       = false,
#         show_trace        = false,
#         show_every        = 10,
#         iterations        = 1000,
#         allow_f_increases = false,
#         time_limit        = 100,
#         x_tol             = 0,
#         f_abstol          = 0,
#         g_tol             = tol,
#         f_calls_limit     = 0,
#         g_calls_limit     = 0,
#     )

#     res = Optim.optimize(
#         scaleloss,
#         d0,
#         # i == 1 ? ParticleSwarm() : BFGS(alphaguess = LineSearches.InitialStatic(alpha=0.9), linesearch = LineSearches.HagerZhang()),
#         n == 1 ? BFGS() : NelderMead(),
#         opts,
#     )
#     d = res.minimizer
#     abs.(d)
# end


text:
using SolidStateDetectors

#define the charge drift model in your custom json file
cdm = ADLChargeDriftModel("custom_config_file.json")

#define a range of relevant electric field strengths in units of V/m
E = Float32.(0:100:1000000)

#calculate the drift velocities from this in units of m/s
ve100 = map(E -> SolidStateDetectors.Vl(E, cdm.electrons.axis100), E)
ve111 = map(E -> SolidStateDetectors.Vl(E, cdm.electrons.axis111), E)
vh100 = map(E -> SolidStateDetectors.Vl(E, cdm.holes.axis100), E)
vh111 = map(E -> SolidStateDetectors.Vl(E, cdm.holes.axis111), E)

#cross check whether the results look reasonable by plotting (optional)
using Plots
plot(E, [ve100, ve111, vh100, vh111],
    xlabel = "Electric field strength [V/m]", ylabel = "Drift velocity [m/s]")

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
(pwd() != @__DIR__) && cd(@__DIR__) # allow starting app from bin/ dir

using Jyunsim
push!(Base.modules_warned_for, Base.PkgId(Jyunsim))
Jyunsim.main()


text:
# computes cohesion forces

using BenchmarkTools

include("../../src/part/CellLists.jl")   # for Point2D

function LennardJonesForceMagnitude_TEST()
    ϵ  = 1.0
    d  = 1.0
    r  = 1.0
    rc = 2.0

    #@code_warntype LennardJonesForceMagnitude(ϵ,d,r,rc)
    @benchmark LennardJonesForceMagnitude($ϵ,$d,$r,$rc)

    #nothing
end

function ForceCalculation_TEST()
    ϵ  = 1.0
    d  = 1.0
    r  = 1.0
    rc = 1.0

    Δx = 2.0
    Δy = 1.5

    #@code_warntype ForceCalculation(ϵ,d,rc,Δx,Δy)
    @benchmark ForceCalculation($ϵ,$d,$rc,$Δx,$Δy)

    #nothing
end

# N: number of particles
function computeCohesion_TEST(N=100)
    nodeList = [Point2D(2*rand()-1.0,2*rand()-1.0) for i=1:N]
    radiusList = ones(N)

    #@code_warntype computeCohesion(nodeList,radiusList)
    @btime computeCohesion_SMARTER($nodeList,$radiusList)

    nothing
end

function computeCohesion_CL!_TEST(N=100,L=0.1)
    nodeList = [Point2D(2*rand()-1.0,2*rand()-1.0) for i=1:N]
    radiusList = 0.5*L*ones(N)

    # construct cell list
    totalBounds = [-1.0-L,1.0+L,-1.0-L,1.0+L]
    cl = generateCellList(nodeList,totalBounds,L)

    # initialize vector for force assignment
    cfX = zeros(Float64,N)
    cfY = zeros(Float64,N)

    #@code_warntype computeCohesion_CL!(nodeList,radiusList,cl,cfX,cfY)
    @btime computeCohesion_CL!($nodeList,$radiusList,$cl,$cfX,$cfY)

    nothing
end

# computes L (cell list side length) ideal for lennard-jones force calculations
# from particle-particle interactions (cohesion)
function computeLJ_L(rc::T,rList::Vector{T}) where T<:Real
    d = 2*maximum(rList)
    return 1.1*rc*d/sqrt(2)
end
function computeLJ_L(rc::T,radius::T) where T<:Real
    d = 2*radius
    return 1.1*rc*d/sqrt(2)
end

# ---------------
# ACCURACY TESTS
# ---------------

function force_plot_TEST()
    ϵ  = 1.0         # strength
    rc = 1.5         # ratio to cut-off

    N = 30
    L = 0.1
    #nodeList = [Point2D(-L+L/3,0.0),Point2D(L-L/3,0.0)]

    nodeList = [Point2D(2*rand()-1.0,2*rand()-1.0) for i=1:N]
    radiusList = 0.5*L*ones(N)

    # initalize force arrays
    cfX = zeros(N); cfY = zeros(N)

    # direct
    computeCohesion!(nodeList,radiusList,rc,ϵ,cfX,cfY)

    for i=1:N
        println(i,"\t",sqrt(cfX[i]^2 + cfY[i]^2))
    end

    return nodeList,radiusList,cfX,cfY,rc
end

function twoParticleComparison_TEST(N)
    ϵ  = 1.0         # strength
    rc = 1.5         # ratio to cut-off

    nodeList = [Point2D(2*rand()-1.0,2*rand()-1.0) for i=1:N]
    radiusList = 0.01*ones(N)

    L = computeLJ_L(rc,radiusList)
    #println("L: ",L)

    # check that minimum cell list length is satisfied
    if (L < rc*2*radiusList[1]/sqrt(2))
        println("L:   ",L)
        println("min: ",rc*2*radiusList[1]/sqrt(2))
        throw(error("cell list length must be larger!"))
    end

    # construct cell list
    totalBounds = [-1.0-L,1.0+L,-1.0-L,1.0+L]
    cl = generateCellList(nodeList,totalBounds,L)

    # initalize force arrays
    cfXdirect = zeros(N); cfYdirect = zeros(N)
    cfXcl     = zeros(N); cfYcl     = zeros(N)

    # direct
    println("N: ",N)
    print("direct:    ")
    #@code_warntype computeCohesion!(nodeList,radiusList,rc,ϵ,cfXdirect,cfYdirect)
    @btime computeCohesion!($nodeList,$radiusList,$rc,$ϵ,$cfXdirect,$cfYdirect)

    # cell list
    print("cell list: ")
    #@code_warntype computeCohesion_CL!(nodeList,radiusList,rc,ϵ,cl,cfXcl,cfYcl)
    @btime computeCohesion_CL!($nodeList,$radiusList,$rc,$ϵ,$cl,$cfXcl,$cfYcl)
    println()

    if false
        xdiff = zeros(N)
        ydiff = zeros(N)
        for i=1:N
            if cfXdirect[i] != 0.0
                xdiff[i] = abs((cfXdirect[i] - cfXcl[i])/cfXdirect[i])
            else
                xdiff[i] = abs(cfXdirect[i] - cfXcl[i])
            end

            if cfYdirect[i] != 0.0
                ydiff[i] = abs((cfYdirect[i] - cfYcl[i])/cfYdirect[i])
            else
                ydiff[i] = abs(cfYdirect[i] - cfYcl[i])
            end
        end

        if any(xdiff .> 1e-12) || any(ydiff .> 1e-12)
            println("force vectors are different!")
            println("max difference: ",maximum([xdiff;ydiff]))
        else
            println("force vectors match!")
        end

        println()

        println("X: ")
        for i=1:N
            println(i,"\t cfXd:  ",cfXdirect[i])
            println("\t cfXcl: ",cfXcl[i])
        end
        println()
        println("Y: ")
        for i=1:N
            println(i,"\t cfYd:  ",cfYdirect[i])
            println("\t cfYcl: ",cfYcl[i])
        end
    end

    nothing
end

function main()
    for N in [100,200,400,800,1600,3200,6400,12800]
        twoParticleComparison_TEST(N);
    end

    nothing
end

max_stars_repo_name: pseastham/eFEMpart


text:
import Base.Broadcast: Broadcasted, Extruded, BroadcastStyle, ArrayStyle

BroadcastStyle(::Type{<:CuArray}) = ArrayStyle{CuArray}()

function Base.similar(bc::Broadcasted{ArrayStyle{CuArray}}, ::Type{T}) where T
    similar(CuArray{T}, axes(bc))
end

# replace base functions with libdevice alternatives
# TODO: do this with Cassette.jl

cufunc(f) = f
cufunc(::Type{T}) where T = (x...) -> T(x...) # broadcasting type ctors isn't GPU compatible

Broadcast.broadcasted(::ArrayStyle{CuArray}, f, args...) =
  Broadcasted{ArrayStyle{CuArray}}(cufunc(f), args, nothing)

libdevice = :[
  cos, cospi, sin, sinpi, tan, acos, asin, atan,
  cosh, sinh, tanh, acosh, asinh, atanh,
  log, log10, log1p, log2, logb, ilogb,
  exp, exp2, exp10, expm1, ldexp,
  erf, erfinv, erfc, erfcinv, erfcx,
  brev, clz, ffs, byte_perm, popc,
  isfinite, isinf, isnan, nearbyint,
  nextafter, signbit, copysign, abs,
  sqrt, rsqrt, cbrt, rcbrt, pow,
  ceil, floor, saturate,
  lgamma, tgamma,
  j0, j1, jn, y0, y1, yn,
  normcdf, normcdfinv, hypot,
  fma, sad, dim, mul24, mul64hi, hadd, rhadd, scalbn].args

for f in libdevice
  isdefined(Base, f) || continue
  @eval cufunc(::typeof(Base.$f)) = CUDAnative.$f
end

using MacroTools

const _cufuncs = copy(libdevice)
cufuncs() = (global _cufuncs; _cufuncs)

function replace_device(ex)
  global _cufuncs
  MacroTools.postwalk(ex) do x
    x in _cufuncs ? :(CuArrays.cufunc($x)) : x
  end
end

macro cufunc(ex)
  global _cufuncs
  def = MacroTools.splitdef(ex)
  f = def[:name]
  def[:name] = Symbol(:cu, f)
  def[:body] = replace_device(def[:body])
  push!(_cufuncs, f)
  quote
    $(esc(MacroTools.combinedef(def)))
    CuArrays.cufunc(::typeof($(esc(f)))) = $(esc(def[:name]))
  end
end

# ForwardDiff Integration
using ForwardDiff: Dual, value, partials, unary_dual_definition
using DiffRules

for f in libdevice
  if haskey(DiffRules.DEFINED_DIFFRULES, (:Base,f,1))
    f == :tanh && continue
    diffrule = DiffRules.DEFINED_DIFFRULES[(:Base,f,1)]
    DiffRules.DEFINED_DIFFRULES[(:CUDAnative,f,1)] =
      (args...) -> replace_device(diffrule(args...))
    eval(unary_dual_definition(:CUDAnative, f))
  end
end

DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :tanh, 1)] = x ->
  replace_device(:(1-tanh(x)^2))
eval(unary_dual_definition(:CUDAnative, :tanh))


text:
using MinimalPerfectHashes
using Base.Test
using Graphs

months = ["jan", "feb", "mar", "apr", "may", "jun", "jul",
  "aug", "sep", "oct", "nov", "dec"]

mph = MinimalPerfectHash(months)

results = unique(map((x) -> hash(x, mph), months))

@test length(results) == 12
@test maximum(results)    == 12
@test minimum(results)    == 1


text:
using DBaseReader
using Base.Test

data = joinpath(dirname(@__FILE__), "data")

# Basic tests
filename = "$data/cb_2014_us_region_20m.dbf"
d = readdbf(filename)
@test length(d) == 7
@test length(d[:NAME]) == 4
@test typeof(d[:NAME][1]) == String && d[:NAME][1] == "Northeast"
@test typeof(d[:ALAND][1]) == Int64 && d[:ALAND][1] == 419356559348
@test Set(collect(keys(d))) == Set([:REGIONCE,:AFFGEOID,:GEOID,:NAME,:LSAD,:ALAND,:AWATER])

# QGIS 2.10-generated dbf containing both a Date-type field
# and a Real Numeric field
filename = "$data/test.dbf"
d = readdbf(filename)
@test d[:ADate][2] == Date(1927, 10, 15)
@test typeof(d[:AReal][2]) == Float64 && d[:AReal][2] == 6.62607

# Test FoxPro 'F' type
filename = "$data/taz_small.dbf"
d = readdbf(filename)
@test d[:TAZ2K][1] == 501030005

# TODO: test for
# * Logical fields (Have never seen one associated with a shapefile though)


text:
using GNBC
using Documenter

makedocs(;
    modules=[GNBC],
    authors="Daymond Ling",
    repo="https://github.com/DaymondLing/GNBC.jl/blob/{commit}{path}#L{line}",
    sitename="GNBC.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DaymondLing.github.io/GNBC.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DaymondLing/GNBC.jl",
)


text:
# Test transcription_model
@testset "transcription_model" begin
    # initialize model
    m = InfiniteModel()
    # test normal usage
    @test isa(transcription_model(m), Model)
    # test error
    set_optimizer_model(m, Model())
    @test_throws ErrorException transcription_model(m)
end

# Test build_optimizer_model!
@testset "build_optimizer_model!" begin
    # initialize model
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    m = InfiniteModel(mockoptimizer)
    @infinite_parameter(m, par in [0, 1], num_supports = 2)
    @infinite_parameter(m, pars[1:2] in [0, 1], supports = [0, 1])
    @variable(m, 1 >= x >= 0, Infinite(par), Int)
    @variable(m, y == 2, Infinite(par, pars), Bin, start = 0)
    @variable(m, x0, Point(x, 0))
    @variable(m, 0 <= y0 <= 1, Point(y, 0, [0, 0]), Int)
    @variable(m, 0 <= z <= 1, Bin)
    @variable(m, w == 1, Int, start = 1)
    meas1 = support_sum(x - w, par)
    meas2 = integral(y, pars)
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @constraint(m, c4, meas2 - 2y0 + x <= 1, DomainRestrictions(par => [0.5, 1]))
    @constraint(m, c5, meas2 == 0)
    @constraint(m, @deriv(x, par) == 0)
    @constraint(m, sin(w) + integral(x^3, par) == 0)
    @objective(m, Min, x0 + meas1)
    set_silent(m)
    set_time_limit_sec(m, 42.)
    # test normal usage
    @test isa(build_optimizer_model!(m, Val(:TransData)), Nothing)
    @test optimizer_model_ready(m)
    @test num_variables(optimizer_model(m)) == 30
    @test time_limit_sec(optimizer_model(m)) == 42
end


text:
@testset "mapreduce" begin
    println("--------")
    n = 100000000
    a = rand(n)
    println("threaded tmapreduce: $(Threads.nthreads()) threads")
    @static if VERSION < v"0.7-"
        tmapreduce(log, +, 0.0, a)
        @time tmapreduce(log, +, 0.0, a)
    else
        tmapreduce(log, +, a, init=0.0)
        @time tmapreduce(log, +, a, init=0.0)
    end

    println("unthreaded")
    @static if VERSION < v"0.7-"
        mapreduce(log, +, 0.0, a)
        @time mapreduce(log, +, 0.0, a)
    else
        mapreduce(log, +, a, init=0.0)
        @time mapreduce(log, +, a, init=0.0)
    end

    @static if VERSION < v"0.7-"
        @test tmapreduce(log, +, 0.0, ones(1000)) == 0
        srand(1234321)
        a = rand(1000)
        r1 = tmapreduce((x)->(2*x), +, 0.0, a)
        r2 = mapreduce((x)->(2*x), +, 0.0, a)
        @test r1 ≈ r2
    else
        @test tmapreduce(log, +, ones(1000), init=0.0) == 0
        Random.seed!(1234321)
        a = rand(1000)
        r1 = tmapreduce((x)->(2*x), +, a, init=0.0)
        r2 = mapreduce((x)->(2*x), +, a, init=0.)
        @test r1 ≈ r2
    end
end


text:
# A collection of bootstrapped code which should be generated
# dynamically in future.

abstract type TNamed <: ROOTStreamedObject end
struct TNamed_1 <: TNamed end
function readfields!(io, fields, ::Type{TNamed_1})
    parsefields!(io, fields, TObject)
    fields[:fName] = readtype(io, String)
    fields[:fTitle] = readtype(io, String)
end

abstract type TAttLine <: ROOTStreamedObject end
struct TAttLine_1 <: TAttLine end
function readfields!(io, fields, T::Type{TAttLine_1})
    fields[:fLineColor] = readtype(io, Int16)
    fields[:fLineStyle] = readtype(io, Int16)
    fields[:fLineWidth] = readtype(io, Int16)
end
struct TAttLine_2 <: TAttLine end
function readfields!(io, fields, T::Type{TAttLine_2})
    fields[:fLineColor] = readtype(io, Int16)
    fields[:fLineStyle] = readtype(io, Int16)
    fields[:fLineWidth] = readtype(io, Int16)
end

abstract type TAttFill <: ROOTStreamedObject end
struct TAttFill_1 <: TAttFill end
function readfields!(io, fields, T::Type{TAttFill_1})
    fields[:fFillColor] = readtype(io, Int16)
    fields[:fFillStyle] = readtype(io, Int16)
end
struct TAttFill_2 <: TAttFill end
function readfields!(io, fields, T::Type{TAttFill_2})
    fields[:fFillColor] = readtype(io, Int16)
    fields[:fFillStyle] = readtype(io, Int16)
end

abstract type TAttMarker <: ROOTStreamedObject end
struct TAttMarker_2 <: TAttFill end
function readfields!(io, fields, T::Type{TAttMarker_2})
    fields[:fMarkerColor] = readtype(io, Int16)
    fields[:fMarkerStyle] = readtype(io, Int16)
    fields[:fMarkerSize] = readtype(io, Float32)
end

abstract type TAttAxis <: ROOTStreamedObject end
struct TAttAxis_4 <: TAttAxis end
function readfields!(io, fields, T::Type{TAttAxis_4})
    fields[:fNdivisions] = readtype(io, Int32)
    fields[:fAxisColor] = readtype(io, Int16)
    fields[:fLabelColor] = readtype(io, Int16)
    fields[:fLabelFont] = readtype(io, Int16)
    fields[:fLabelOffset] = readtype(io, Float32)
    fields[:fLabelSize] = readtype(io, Float32)
    fields[:fTickLength] = readtype(io, Float32)
    fields[:fTitleOffset] = readtype(io, Float32)
    fields[:fTitleSize] = readtype(io, Float32)
    fields[:fTitleColor] = readtype(io, Int16)
    fields[:fTitleFont] = readtype(io, Int16)
end

abstract type TAxis <: ROOTStreamedObject end
struct TAxis_10 <: TAxis end
function readfields!(io, fields, T::Type{TAxis_10})
    # overrides things like fName,... that were set from the parent TH1 :(
    stream!(io, fields, TNamed)
    stream!(io, fields, TAttAxis)
    fields[:fNbins] = readtype(io, Int32)
    fields[:fXmin] = readtype(io, Float64)
    fields[:fXmax] = readtype(io, Float64)
    fields[:fXbins] = readtype(io, TArrayD)
    fields[:fFirst] = readtype(io, Int16)
    fields[:fLast] = readtype(io, Int16)
    fields[:fBits2] = readtype(io, UInt16)
    fields[:fTimeDisplay] = readtype(io, Bool)
    fields[:fTimeFormat] = readtype(io, String)
end

abstract type TH1 <: ROOTStreamedObject end
struct TH1_8 <: TH1 end
function readfields!(io, fields, T::Type{TH1_8}) end

abstract type TH2 <: ROOTStreamedObject end
struct TH2_4 <: TH2 end
function readfields!(io, fields, T::Type{TH2_4}) end

@with_kw struct ROOT_3a3a_TIOFeatures <: ROOTStreamedObject
    fIOBits
end

function parsefields!(io, fields, T::Type{ROOT_3a3a_TIOFeatures})
    preamble = Preamble(io, T)
    skip(io, 4)
    fields[:fIOBits] = readtype(io, UInt8)
    endcheck(io, preamble)
end

# FIXME maybe this is unpack, rather than readtype?
function readtype(io, T::Type{ROOT_3a3a_TIOFeatures})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

# FIXME this should be generated
@with_kw struct TLeaf
    # FIXME these two come from TNamed
    fName
    fTitle

    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount
end

function parsefields!(io, fields, ::Type{T}) where {T<:TLeaf}
    preamble = Preamble(io, T)
    stream!(io, fields, TNamed)
    fields[:fLen] = readtype(io, Int32)
    fields[:fLenType] = readtype(io, Int32)
    fields[:fOffset] = readtype(io, Int32)
    fields[:fIsRange] = readtype(io, Bool)
    fields[:fIsUnsigned] = readtype(io, Bool)
    fields[:fLeafCount] = readtype(io, UInt32)

    # FIXME this needs to be checked, sometimes the TLeaf is too short
    observed = position(io) - preamble.start
    for _ in 1:(preamble.cnt - observed)
        read(io, 1)
    end
    endcheck(io, preamble)
end

# FIXME this should be generated
@with_kw struct TLeafElement
    # FIXME these two come from TNamed
    fName
    fTitle

    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    fID
    fType
end

function parsefields!(io, fields, ::Type{T}) where {T<:TLeafElement}
    parsefields!(io, fields, TLeaf)
    fields[:fID] = readtype(io, Int32)
    fields[:fType] = readtype(io, Int32)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafElement})
    @initparse
    preamble = Preamble(io, T)
    parsefields!(io, fields, T)
    endcheck(io, preamble)
    T(;fields...)
end

# FIXME this should be generated and inherited from TLeaf
@with_kw struct TLeafI
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end

function parsefields!(io, fields, T::Type{TLeafI})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, Int32)
    fields[:fMaximum] = readtype(io, Int32)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafI})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

primitivetype(l::TLeafI) = l.fIsUnsigned ? UInt32 : Int32

# FIXME this should be generated and inherited from TLeaf
@with_kw struct TLeafL
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end

function parsefields!(io, fields, T::Type{TLeafL})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, Int64)
    fields[:fMaximum] = readtype(io, Int64)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafL})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

primitivetype(l::TLeafL) = l.fIsUnsigned ? UInt64 : Int64

# FIXME this should be generated and inherited from TLeaf
@with_kw struct TLeafO
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end
primitivetype(l::TLeafO) = Bool

function parsefields!(io, fields, T::Type{TLeafO})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, Bool)
    fields[:fMaximum] = readtype(io, Bool)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafO})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

# FIXME this should be generated and inherited from TLeaf
@with_kw struct TLeafF
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end

function parsefields!(io, fields, T::Type{TLeafF})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, Float32)
    fields[:fMaximum] = readtype(io, Float32)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafF})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

primitivetype(l::TLeafF) = Float32

# FIXME this should be generated and inherited from TLeaf
# https://root.cern/doc/master/TLeafB_8h_source.html#l00026
@with_kw struct TLeafB
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end

function parsefields!(io, fields, T::Type{TLeafB})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, UInt8)
    fields[:fMaximum] = readtype(io, UInt8)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafB})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

primitivetype(l::TLeafB) = UInt8
# FIXME this should be generated and inherited from TLeaf
@with_kw struct TLeafD
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end

function parsefields!(io, fields, T::Type{TLeafD})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, Float64)
    fields[:fMaximum] = readtype(io, Float64)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TLeafD})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

primitivetype(l::TLeafD) = Float64

# FIXME this should be generated and inherited from TLeaf
@with_kw struct TLeafC
    # from TNamed
    fName
    fTitle

    # from TLeaf
    fLen
    fLenType
    fOffset
    fIsRange
    fIsUnsigned
    fLeafCount

    # own fields
    fMinimum
    fMaximum
end

function parsefields!(io, fields, ::Type{T}) where {T<:TLeafC}
    preamble = Preamble(io, T)
    parsefields!(io, fields, TLeaf)
    fields[:fMinimum] = readtype(io, Int32)
    fields[:fMaximum] = readtype(io, Int32)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, ::Type{T}) where {T<:TLeafC}
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

abstract type TBranch <: ROOTStreamedObject end
abstract type TBranchElement <: ROOTStreamedObject end
function Base.hash(b::Union{TBranch, TBranchElement}, h::UInt)
    h = hash(b.fFileName, h)
    h = hash(b.fName, h)
    h = hash(b.fEntries, h)
end
Base.length(b::Union{TBranch, TBranchElement}) = b.fEntries
Base.eachindex(b::Union{TBranch, TBranchElement}) = Base.OneTo(b.fEntries)
numbaskets(b::Union{TBranch, TBranchElement}) = findfirst(x->x>(b.fEntries-1),b.fBasketEntry)-1

@with_kw struct TBranch_12 <: TBranch
    cursor::Cursor
    # from TNamed
    fName
    fTitle

    # from TAttFill
    fFillColor
    fFillStyle

    fCompress
    fBasketSize
    fEntryOffsetLen
    fWriteBasket
    fEntryNumber

    fOffset
    fMaxBaskets
    fSplitLevel
    fEntries
    fFirstEntry
    fTotBytes
    fZipBytes

    fBranches
    fLeaves
    fBaskets
    fBasketBytes::Vector{Int64}
    fBasketEntry::Vector{Int64}
    fBasketSeek::Vector{Int64}
    fFileName
end
function readfields!(cursor::Cursor, fields, ::Type{T}) where {T<:TBranch_12}
    io = cursor.io
    tkey = cursor.tkey
    refs = cursor.refs

    stream!(io, fields, TNamed)
    stream!(io, fields, TAttFill)

    fields[:fCompress] = readtype(io, Int32)
    fields[:fBasketSize] = readtype(io, Int32)
    fields[:fEntryOffsetLen] = readtype(io, Int32)
    fields[:fWriteBasket] = readtype(io, Int32)
    fields[:fEntryNumber] = readtype(io, Int64)

    fields[:fOffset] = readtype(io, Int32)
    fields[:fMaxBaskets] = readtype(io, UInt32)
    fields[:fSplitLevel] = readtype(io, Int32)
    fields[:fEntries] = readtype(io, Int64)
    fields[:fFirstEntry] = readtype(io, Int64)
    fields[:fTotBytes] = readtype(io, Int64)
    fields[:fZipBytes] = readtype(io, Int64)

    fields[:fBranches] = unpack(io, tkey, refs, TObjArray)
    fields[:fLeaves] = unpack(io, tkey, refs, TObjArray)
    fields[:fBaskets] = unpack(io, tkey, refs, TObjArray)
    # fields[:fBaskets] = unpack(io, tkey, refs, Undefined)
    speedbump = true  # FIXME speedbump?

    speedbump && skip(io, 1)

    fields[:fBasketBytes] = [readtype(io, Int32) for _ in 1:fields[:fMaxBaskets]]

    speedbump && skip(io, 1)

    # this is also called fBsketEvent, as far as I understood
    fields[:fBasketEntry] = [readtype(io, Int64) for _ in 1:fields[:fMaxBaskets]]

    speedbump && skip(io, 1)

    fields[:fBasketSeek] = [readtype(io, Int64) for _ in 1:fields[:fMaxBaskets]]
    fields[:fFileName] = readtype(io, String)
end
@with_kw struct TBranch_13 <: TBranch
    cursor::Cursor
    # from TNamed
    fName
    fTitle

    # from TAttFill
    fFillColor
    fFillStyle

    fCompress
    fBasketSize
    fEntryOffsetLen
    fWriteBasket
    fEntryNumber

    fIOFeatures

    fOffset
    fMaxBaskets
    fSplitLevel
    fEntries
    fFirstEntry
    fTotBytes
    fZipBytes

    fBranches
    fLeaves
    fBaskets
    fBasketBytes::Vector{Int64}
    fBasketEntry::Vector{Int64}
    fBasketSeek::Vector{Int64}
    fFileName
end
function readfields!(cursor::Cursor, fields, ::Type{T}) where {T<:TBranch_13}
    io = cursor.io
    tkey = cursor.tkey
    refs = cursor.refs

    stream!(io, fields, TNamed)
    stream!(io, fields, TAttFill)

    fields[:fCompress] = readtype(io, Int32)
    fields[:fBasketSize] = readtype(io, Int32)
    fields[:fEntryOffsetLen] = readtype(io, Int32)
    fields[:fWriteBasket] = readtype(io, Int32)
    fields[:fEntryNumber] = readtype(io, Int64)

    fields[:fIOFeatures] = readtype(io, ROOT_3a3a_TIOFeatures)

    fields[:fOffset] = readtype(io, Int32)
    fields[:fMaxBaskets] = readtype(io, UInt32)
    fields[:fSplitLevel] = readtype(io, Int32)
    fields[:fEntries] = readtype(io, Int64)
    fields[:fFirstEntry] = readtype(io, Int64)
    fields[:fTotBytes] = readtype(io, Int64)
    fields[:fZipBytes] = readtype(io, Int64)

    fields[:fBranches] = unpack(io, tkey, refs, TObjArray)
    fields[:fLeaves] = unpack(io, tkey, refs, TObjArray)
    fields[:fBaskets] = unpack(io, tkey, refs, TObjArray)
    # fields[:fBaskets] = unpack(io, tkey, refs, Undefined)
    speedbump = true  # FIXME speedbump?

    speedbump && skip(io, 1)

    fields[:fBasketBytes] = [readtype(io, Int32) for _ in 1:fields[:fMaxBaskets]]

    speedbump && skip(io, 1)

    # this is also called fBsketEvent, as far as I understood
    fields[:fBasketEntry] = [readtype(io, Int64) for _ in 1:fields[:fMaxBaskets]]

    speedbump && skip(io, 1)

    fields[:fBasketSeek] = [readtype(io, Int64) for _ in 1:fields[:fMaxBaskets]]
    fields[:fFileName] = readtype(io, String)
end

@with_kw struct TBranchElement_9 <: TBranchElement
    cursor::Cursor
    # from TNamed
    fName
    fTitle

    # from TAttFill
    fFillColor
    fFillStyle

    fCompress
    fBasketSize
    fEntryOffsetLen
    fWriteBasket
    fEntryNumber

    fOffset
    fMaxBaskets
    fSplitLevel
    fEntries
    fFirstEntry
    fTotBytes
    fZipBytes

    fBranches
    fLeaves
    fBaskets
    fBasketBytes::Vector{Int64}
    fBasketEntry::Vector{Int64}
    fBasketSeek::Vector{Int64}
    fFileName

    # own fields
    fClassName
    fParentName
    fClonesName
    fCheckSum
    fClassVersion
    fID
    fType
    fStreamerType
    fMaximum
    fBranchCount
    fBranchCount2
end

@with_kw struct TBranchElement_10 <: TBranchElement
    cursor::Cursor
    # from TNamed
    fName
    fTitle

    # from TAttFill
    fFillColor
    fFillStyle

    fCompress
    fBasketSize
    fEntryOffsetLen
    fWriteBasket
    fEntryNumber

    fIOFeatures=0x00

    fOffset
    fMaxBaskets
    fSplitLevel
    fEntries
    fFirstEntry
    fTotBytes
    fZipBytes

    fBranches
    fLeaves
    fBaskets
    fBasketBytes::Vector{Int64}
    fBasketEntry::Vector{Int64}
    fBasketSeek::Vector{Int64}
    fFileName

    # own fields
    fClassName
    fParentName
    fClonesName
    fCheckSum
    fClassVersion
    fID
    fType
    fStreamerType
    fMaximum
    fBranchCount
    fBranchCount2
end

function readfields!(cursor::Cursor, fields, ::Type{T}) where {T<:TBranchElement_9}
    io = cursor.io
    tkey = cursor.tkey
    refs = cursor.refs

    stream!(cursor, fields, TBranch)

    fields[:fClassName] = readtype(io, String)
    fields[:fParentName] = readtype(io, String)
    fields[:fClonesName] = readtype(io, String)
    fields[:fCheckSum] = readtype(io, UInt32)
    fields[:fClassVersion] = readtype(io, Int32)
    fields[:fID] = readtype(io, Int32)
    fields[:fType] = readtype(io, Int32)
    fields[:fStreamerType] = readtype(io, Int32)
    fields[:fMaximum] =readtype(io, Int32)
    fields[:fBranchCount] = readobjany!(io, tkey, refs)
    fields[:fBranchCount2] = readobjany!(io, tkey, refs)
end

function readfields!(cursor::Cursor, fields, ::Type{T}) where {T<:TBranchElement_10}
    io = cursor.io
    tkey = cursor.tkey
    refs = cursor.refs

    stream!(cursor, fields, TBranch)

    fields[:fClassName] = readtype(io, String)
    fields[:fParentName] = readtype(io, String)
    fields[:fClonesName] = readtype(io, String)
    fields[:fCheckSum] = readtype(io, UInt32)
    fields[:fClassVersion] = readtype(io, Int16)
    fields[:fID] = readtype(io, Int32)
    fields[:fType] = readtype(io, Int32)
    fields[:fStreamerType] = readtype(io, Int32)
    fields[:fMaximum] =readtype(io, Int32)
    fields[:fBranchCount] = readobjany!(io, tkey, refs)
    fields[:fBranchCount2] = readobjany!(io, tkey, refs)
end

# FIXME preliminary TTree structure
@with_kw struct TTree
    # TNamed
    fName
    fTitle

    # TAttLine
    fLineColor
    fLineStyle
    fLineWidth

    # TAttFill
    fFillColor
    fFillStyle

    # TAttMarker
    fMarkerColor
    fMarkerStyle
    fMarkerSize

    fEntries
    fTotBytes
    fZipBytes
    fSavedBytes
    fFlushedBytes
    fWeight
    fTimerInterval
    fScanField
    fUpdate
    fDefaultEntryOffsetLen
    fNClusterRange
    fMaxEntries
    fMaxEntryLoop
    fMaxVirtualSize
    fAutoSave
    fAutoFlush
    fEstimate

    fClusterRangeEnd
    fClusterSize

    fIOFeatures

    fBranches
    fLeaves

    fAliases
    fIndexValues
    fIndex
    fTreeIndex
    fFriends
end

TH1F(io, tkey::TKey, refs) = TH(io, tkey, refs)
TH2F(io, tkey::TKey, refs) = TH(io, tkey, refs)
TH1D(io, tkey::TKey, refs) = TH(io, tkey, refs)
TH2D(io, tkey::TKey, refs) = TH(io, tkey, refs)

function TH(io, tkey::TKey, refs)
    fields = Dict{Symbol, Any}()

    io = datastream(io, tkey)
    preamble = Preamble(io, Missing)

    is2d = startswith(tkey.fClassName, "TH2")
    if is2d
        stream!(io, fields, TH2, check=false)
    end

    stream!(io, fields, TH1, check=false)
    stream!(io, fields, TNamed)
    stream!(io, fields, TAttLine)
    stream!(io, fields, TAttFill)
    stream!(io, fields, TAttMarker)
    fields[:fNcells] = readtype(io, Int32)

    for axis in ["fXaxis_", "fYaxis_", "fZaxis_"]
        subfields = Dict{Symbol, Any}()
        stream!(io, subfields, TAxis, check=false)
        fields[Symbol(axis, :fLabels)] = readobjany!(io, tkey, refs)
        fields[Symbol(axis, :fModLabs)] = readobjany!(io, tkey, refs)
        for (k,v) in subfields
            fields[Symbol(axis, k)] = v
        end
        # FIXME this line makes non-uniform binned histograms work, but not sure why
        readtype(io, Int32)
    end

    fields[:fBarOffset] = readtype(io, Int16)
    fields[:fBarWidth] = readtype(io, Int16)
    for symb in [:fEntries, :fTsumw, :fTsumw2, :fTsumwx, :fTsumwx2, :fMaximum, :fMinimum, :fNormFactor]
        fields[symb] = readtype(io, Float64)
    end

    fields[:fContour] = readtype(io, TArrayD)
    fields[:fSumw2] = readtype(io, TArrayD)
    fields[:fOption] = readtype(io, String)
    # if user saved after calling h.Fit() with a TF1, then this will error
    fields[:fFunctions] = unpack(io, tkey, refs, TList)
    fields[:fBufferSize] = readtype(io, Int32)
    skip(io, 1) # speedbump
    fields[:fBuffer] = readtype(io, TArrayD)
    fields[:fBinStatErrOpt] = readtype(io, Int16)
    fields[:fStatOverflows] = readtype(io, Int16)

    if is2d
        for symb in [:fScalefactor, :fTsumwy, :fTsumwy2, :fTsumwxy]
            fields[symb] = readtype(io, Float64)
        end
    end

    arraytype = endswith(tkey.fClassName, 'F') ? TArrayF : TArrayD
    fields[:fN] = readtype(io, arraytype)
    fields
end

function TDirectory(io, tkey::TKey, refs)
    fobj = io
    seekstart(fobj, tkey)

    # almost verbatim from L95 to L101 of root.jl
    dir_header = unpack(fobj, ROOTDirectoryHeader)
    seek(fobj, dir_header.fSeekKeys)
    header_key = unpack(fobj, TKey)
    n_keys = readtype(fobj, Int32)
    keys = [unpack(fobj, TKey) for _ in 1:n_keys]

    directory = ROOTDirectory(header_key.fName, dir_header, keys, fobj, refs)
    return directory
end

# FIXME idk what is going on but this just looks like a TTree.....
function TNtuple(io, tkey::TKey, refs)
    io = datastream(io, tkey)
    preamble = Preamble(io, Missing)
    tree = TTree(io, tkey, refs; top=false) #embeded tree
end

# FIXME preliminary TTree implementation
function TTree(io, tkey::TKey, refs; top=true)
    # if embeded in a Ntuple, don't run datastream again
    io = top ? datastream(io, tkey) : io

    @initparse

    preamble = Preamble(io, Missing)
    # @show preamble

    stream!(io, fields, TNamed)

    stream!(io, fields, TAttLine)
    stream!(io, fields, TAttFill)
    stream!(io, fields, TAttMarker)

    fields[:fEntries] = readtype(io, Int64)
    fields[:fTotBytes] = readtype(io, Int64)
    fields[:fZipBytes] = readtype(io, Int64)
    fields[:fSavedBytes] = readtype(io, Int64)
    fields[:fFlushedBytes] = readtype(io, Int64)
    fields[:fWeight] = readtype(io, Float64)
    fields[:fTimerInterval] = readtype(io, Int32)
    fields[:fScanField] = readtype(io, Int32)
    fields[:fUpdate] = readtype(io, Int32)
    # See https://github.com/cbourjau/alice-rs/blob/6af19a78fe5521f5b27466d7d20f7dfacd38a38f/root-io/src/tree_reader/tree.rs#L148
    if preamble.version >= 18
        fields[:fDefaultEntryOffsetLen] = readtype(io, Int32)
    end
    if preamble.version >= 19
        fields[:fNClusterRange] = readtype(io, UInt32)
    end
    fields[:fMaxEntries] = readtype(io, Int64)
    fields[:fMaxEntryLoop] = readtype(io, Int64)
    fields[:fMaxVirtualSize] = readtype(io, Int64)
    fields[:fAutoSave] = readtype(io, Int64)
    fields[:fAutoFlush] = readtype(io, Int64)
    fields[:fEstimate] = readtype(io, Int64)

    # FIXME what about speedbumps??
    speedbump = true

    # See https://github.com/cbourjau/alice-rs/blob/6af19a78fe5521f5b27466d7d20f7dfacd38a38f/root-io/src/tree_reader/tree.rs#L148
    if haskey(fields, :fNClusterRange)
        speedbump && skip(io, 1)
        fields[:fClusterRangeEnd] = [readtype(io, Int64) for _ in 1:fields[:fNClusterRange]]
        speedbump && skip(io, 1)
        fields[:fClusterSize] = [readtype(io, Int64) for _ in 1:fields[:fNClusterRange]]
    end

    # for key in keys(fields)
    #     @show key, fields[key]
    # end

    if preamble.version >= 20
        fields[:fIOFeatures] = readtype(io, ROOT_3a3a_TIOFeatures)
    else
        fields[:fIOFeatures] = missing
    end

    fields[:fBranches] = unpack(io, tkey, refs, TObjArray)
    fields[:fLeaves] = unpack(io, tkey, refs, TObjArray)

    fields[:fAliases] = readobjany!(io, tkey, refs)
    fields[:fIndexValues] = readtype(io, TArrayD)
    fields[:fIndex] = readtype(io, TArrayI)
    fields[:fTreeIndex] = readobjany!(io, tkey, refs)
    fields[:fFriends] = readobjany!(io, tkey, refs)

    # uproot unpacks Undefined instances here, we cannot since
    # the Base.GenericIOBuffer{Array{UInt8,1}} from the compression
    # library throws an EOFError when we read more than available.
    # FIXME this needs to be checked though!
    # while !eof(io)
    #     read(io, 1)
    # end
    # unpack(io, tkey, refs, Undefined)
    # println(fields[:fBranches])

    # endcheck(io, preamble)
    TTree(;fields...)
end

# FIXME what to do with auto.py's massive type translation?
# https://github.com/scikit-hep/uproot3/blob/54f5151fb7c686c3a161fbe44b9f299e482f346b/uproot3/interp/auto.py#L360-L365


text:
using Zygote: pullback
using ProximalCore

function ProximalCore.gradient!(grad, f, x)
    fx, pb = pullback(f, x)
    grad .= pb(one(fx))[1]
    return fx
end


text:
@testset "CCall to SBML works and SBML returns a version" begin
    @test SBMLVersion() isa VersionNumber
end


text:
# TravelTime
let
    model = TravelTime()
    cost = 1.0
    transition_cost = 1.0
    new_cost = accumulate_cost(model,cost,transition_cost)
    @test new_cost == cost + transition_cost
end
# TravelDistance
let
    model = TravelDistance()
    cost = 1.0
    transition_cost = 1.0
    new_cost = accumulate_cost(model,cost,transition_cost)
    @test new_cost == cost + transition_cost
end
# NullCost
let
    model = NullCost()
    cost = 1.0
    transition_cost = 1.0
    new_cost = accumulate_cost(model,cost,transition_cost)
    @test new_cost == cost
end
# DeadlineCost
let
    t_max = 10.0
    model = DeadlineCost(t_max)
    cost = get_initial_cost(model)
    @test cost == 0.0
    transition_cost = 1.0
    new_cost = accumulate_cost(model, cost, transition_cost)
    @test new_cost == 1.0
    h_cost = 8.0
    @test add_heuristic_cost(model, cost, h_cost) == 0.0
    h_cost = 11.0
    @test add_heuristic_cost(model, cost, h_cost) == 1.0
    set_deadline!(model,0.0)
    @test add_heuristic_cost(model, cost, h_cost) == cost + h_cost
    # set_deadline!(model,t_max,t_max-4)
end
# SumOfMakeSpans cost
let
    tF = [1.0,2.0,3.0]
    root_nodes = [2,3]
    weights = [1.0,1.0]
    deadlines = [0.0,0.0]
    let
        model = SumOfMakeSpans(tF,root_nodes,weights,deadlines)
        set_deadline!(model,tF[root_nodes])
        @test model.deadlines == tF[root_nodes]

        cost = get_initial_cost(model)
        @test add_heuristic_cost(model, 0, 2) == 0 + 0
        @test add_heuristic_cost(model, 0, 3) == 1 + 0
        @test add_heuristic_cost(model, 0, 4) == 2 + 1
        for (c,g) in [(0,1),(0,2),(0,3),(0,4)]
            @test add_heuristic_cost(model, c, g) == c + sum(model.weights .* max.(0, g .- model.tF[model.root_nodes]))
        end
    end
    let
        model = MakeSpan(tF,root_nodes,weights,deadlines)
        set_deadline!(model,tF[root_nodes])
        @test model.deadlines == tF[root_nodes]
        cost = get_initial_cost(model)
        for (c,g) in [(0,1),(0,2),(0,3),(0,4)]
            @test add_heuristic_cost(model, c, g) == c + maximum(model.weights .* max.(0, g .- model.tF[model.root_nodes]))
        end
    end
end
# ConflictCostModel
let
    G = Graph(5)
    for v in 1:nv(G)-1
        add_edge!(G,v,v+1)
    end
    T = 10
    n_agents = 2
    h = HardConflictCost(G,T,n_agents)
    set_path!(h,1,[5,4,3],0)
    set_path!(h,2,[1,2,3],0)
    path_id_to_time_idx = get_conflicting_paths(h.model.table)
    @test path_id_to_time_idx[1] == path_id_to_time_idx[2] == 2
    # Agent 1 has a cost of 1 whenever it conflicts with agent 2's path
    @test get_conflict_value(h.model,1,2,1) == 1.0
    @test get_conflict_value(h.model,1,3,2) == 1.0
    @test get_conflict_value(h.model,1,4,0) == 0.0
    @test get_conflict_value(h.model,1,3,1) == 0.0
    # Agent 2 should have zero cost even when it conflicts with its own path
    @test get_conflict_value(h.model,2,2,1) == 0.0
    # Partial path reset - check that later conflicts still remain
    partially_set_path!(h.model,2,[2,1],0)
    @test get_conflict_value(h.model,1,2,1) == 0.0
    @test get_conflict_value(h.model,1,3,2) == 1.0

    set_path!(h.model,2,[3,4,5],0)
    # Now that agent 2's path has been updated, the cost should be zero for a
    # vtx no longer on the path
    @test get_conflict_value(h.model,1,3,2) == 0.0
end
# FatPaths (SoftConflictTable)
let
    for mapf in [
            CRCBS.init_fat_path_mapf(init_mapf_1(),FlatFPCost()),
            CRCBS.init_fat_path_mapf(init_mapf_1(),NormalizedFPCost()),
        ]
        table = mapf.env.cost_model.cost_models[2].model.table
        env = build_env(mapf,initialize_root_node(mapf),2)

        @test get_conflict_value(table,1,env,CBSEnv.State(5,1)) == 0
        @test get_conflict_value(table,2,env,CBSEnv.State(5,1)) >= 0
        @test get_conflict_value(table,2,env,CBSEnv.Action(Edge(1,5),1),1) >= 0
        @test get_conflict_value(table,2,env,CBSEnv.Action(Edge(5,1),1),1) >= 0
        @test get_conflict_value(table,2,env,CBSEnv.Action(Edge(5,1),1),2) == 0
        cost = get_transition_cost(env,CBSEnv.State(5,0),CBSEnv.Action(Edge(5,1),1),CBSEnv.State(1,1))
        @test cost[2] >= 0.0

        CRCBS.clear_fat_path!(table,1)
        @test get_conflict_value(table,1,env,CBSEnv.State(5,1)) == 0
        @test get_conflict_value(table,2,env,CBSEnv.State(5,1)) == 0
        @test get_conflict_value(table,2,env,CBSEnv.Action(Edge(1,5),1),1) == 0
        @test get_conflict_value(table,2,env,CBSEnv.Action(Edge(5,1),1),1) == 0
        @test get_conflict_value(table,2,env,CBSEnv.Action(Edge(5,1),1),2) == 0
        cost = get_transition_cost(env,CBSEnv.State(5,0),CBSEnv.Action(Edge(5,1),1),CBSEnv.State(1,1))
        @test cost[2] == 0.0
    end

end
let
        mapf = CRCBS.init_fat_path_mapf(init_mapf_4(),NormalizedFPCost())
        table = mapf.env.cost_model.cost_models[2].model.table
        env = build_env(mapf,initialize_root_node(mapf),2)
        s = CBSEnv.State(5,1)
        a = CBSEnv.Action(Edge(5,1),1)

        @test get_conflict_value(table,1,env,s) == 0
        @test get_conflict_value(table,2,env,s) == 1.0/2
        @test get_conflict_value(table,2,env,reverse(a),1) == 1.0/2
        @test get_conflict_value(table,2,env,a,1) >= 1.0/2
        @test get_conflict_value(table,2,env,a,2) == 0
        s0 = CBSEnv.State(5,0)
        sp = get_next_state(env,s0,a)
        cost = get_transition_cost(env,s0,a,sp)
        @test cost[2] >= 0.0

end
# test FatPaths
let
    mapf = CRCBS.init_fat_path_mapf(CRCBS.init_mapf_4())
    solver = CBSSolver(AStar{cost_type(mapf)}())
    solution, cost = solve!(solver,mapf)
    paths = convert_to_vertex_lists(solution)
    @test paths[1] == [1,5,9,13,14,15]
    @test paths[2] == [16,12,8,4,3,2]
end
# MetaCost
let
    t0 = 0.0
    n_agents = 3
    cost = MetaCost{Float64}([1.0,3.0,4.0],0.0)
    transition_cost = [1.0,1.0,1.0]
    let
        FullCostModel(MaxCost(),FinalTime())
        FullCostModel(maximum,FinalTime())
        MakeSpan()
        SumOfTravelTime()
        SumOfTravelDistance()
    end
    let
        model = MetaCostModel(MakeSpan(),n_agents)
        new_cost = accumulate_cost(model,cost,transition_cost)
        @test new_cost.total_cost == 5.0
    end
    let
        model = MetaCostModel(SumOfTravelTime(),n_agents)
        new_cost = accumulate_cost(model,cost,transition_cost)
        @test new_cost.total_cost == 11.0
    end
    let
        model = FullCostModel(max,MetaCostModel(SumOfTravelTime(),n_agents))
        get_cost_model(model)
        new_cost = accumulate_cost(model,cost,transition_cost)
        @test new_cost.total_cost == 11.0
    end
    # @test MetaCost{Float64}([1.0],0.0) < MetaCost{Float64}([0.0],1.0)
end
# CompositeCost
let
    model = construct_composite_cost_model(TravelTime(),TravelTime(),TravelDistance())
    cost = (0.0,0.0,0.0)
    transition_cost = (1.0,1.0,1.0)
    @test accumulate_cost(model, cost,transition_cost) == (1.0,1.0,1.0)
    get_initial_cost(model)
    get_infeasible_cost(model)
end
let
    cost_model = construct_composite_cost_model(
        FullDeadlineCost(DeadlineCost(10)),
        FullCostModel(sum,NullCost()),
        SumOfTravelTime()
    )
    get_infeasible_cost(cost_model)
end
let
    p = Path{GraphState,GraphAction,Float64}(s0=GraphState(1,0))
    vtx_grid = initialize_dense_vtx_grid(4,4)
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    env = CBSEnv.LowLevelEnv(graph=graph)
    for e in [
            Edge(0,1),
            Edge(1,2),
            Edge(2,2),
        ]
        s = get_final_state(p)
        a = GraphAction(e,1)
        sp = get_next_state(env,s,a)
        add_to_path!(p,env,s,a,sp)
    end
    for (m,c) in [
        (TravelTime(),3),
        (TravelDistance(),2)]
    end
    @test compute_path_cost(TravelTime(),env,p,1) == 3
    @test compute_path_cost(SumOfTravelTime(),env,p,1) == 3
    @test compute_path_cost(TravelDistance(),env,p,1) == 2
    @test compute_path_cost(SumOfTravelDistance(),env,p,1) == 2
    @test compute_path_cost(construct_composite_cost_model(
        SumOfTravelTime(),
        SumOfTravelDistance(),
        ),env,p,1) == (3,2)
    @test compute_path_cost(DeadlineCost(10),env,p,1) == 3
    @test compute_path_cost(FullDeadlineCost(DeadlineCost(10)),env,p,1) == 3
end


text:
module RAMS
"""
RAMS Julia Library

"""

using HDF5
using Dates
using ProgressMeter
using Statistics

# Get thermo.jl functions
include("thermo.jl")
export temperature, pressure, wsat, rh, mslp

"""
    RAMSDates(flist)

Takes array of file paths and returns datetimes.

# Arguments
- `flist::Array{String,1}`: 1D array of string file paths

# Returns
 - `dtarr::Array{DateTime,1}`: 1D array of datetimes
"""
function RAMSDates(flist::Array{String,1})

    dtregex = r"[0-9]{4}-[0-9]{2}-[0-9]{2}-[0-9]{6}"
    dtarr = DateTime[]
    for f in flist
        dtmatch = match(dtregex, f)
        if dtmatch === nothing
            throw(ErrorException("Filename did not match yyyy-mm-dd-hhmmss"))
        else
            year,month,day,time = split(dtmatch.match,'-')
            year = parse(Int, year)
            month = parse(Int, month)
            day = parse(Int, day)

            hour = parse(Int, time[1:2])
            minute = parse(Int, time[3:4])
            second = parse(Int, time[5:6])

            t = DateTime(Date(year, month, day), Time(hour, minute, second))
            push!(dtarr, t)

        end
    end
    return dtarr
end
export RAMSDates

"""
    RAMSVar(flist, varname)

Function to read a variable from a list of RAMS data files.
If `varname` is a string, the specified variable will be returned.
If `varname` is `Vector{String}`, the listed variables will be added together
and returned.

# Arguments
- `flist::Array{String,1}`: 1D array of string file paths
- `varname::String`: Name of variable.
- Optional: `varname::Vector{String}`: Vector of variable names, the listed variables will be added together

# Returns
- `var::Array`: Output variable
"""
function RAMSVar(flist::Array{String,1}, varname::String; meandims=nothing)
    temp = h5read(flist[1], varname)
    nt = length(flist)

    if meandims !== nothing
        basedims = []
        for (i,d) in enumerate(size(temp))
            if !(i in meandims) append!(basedims, d) end
        end
    else
        basedims = [i for i in size(temp)]
    end

    dims = vcat(basedims, nt)
    t = typeof(temp[1])
    var = zeros(t, dims...)

    @showprogress for (i,f) in enumerate(flist[1:end])
        if meandims !== nothing selectdim(var,length(basedims)+1,i) .+= dropmean(h5read(f, varname), meandims)
        else selectdim(var,length(basedims)+1,i) .= h5read(f, varname) end
    end

    return var
end

# TODO Extend this function too
function RAMSVar(flist::Array{String,1}, varname::Vector{String}; meandims=nothing)
    temp = h5read(flist[1], varname[1])
    nt = length(flist)

    if meandims !== nothing
        basedims = []
        for (i,d) in enumerate(size(temp))
            if !(i in meandims) append!(basedims, d) end
        end
    else
        basedims = [i for i in size(temp)]
    end

    dims = vcat(basedims, nt)
    t = typeof(temp[1])
    var = zeros(t, dims...)

    @showprogress for (i,f) in enumerate(flist[1:end])
        for vname in varname
            if meandims !== nothing selectdim(var,length(basedims)+1,i) .+= dropmean(h5read(f, vname), meandims)
            else selectdim(var,length(basedims)+1,i) .+= h5read(f, vname) end
        end
    end
    return var
end
export RAMSVar

"""
    dropmean(var, meandrop)
Function to drop dimensions used while taking a mean.
Specify `meandrop` as a tuple of dimensions (e.g. `meandrop=(1,2)`)
"""
function dropmean(var, meandrop)
    return dropdims(mean(var, dims=meandrop), dims=meandrop)
end
export dropmean

"""
    vert_int(var, ztn; notime=false)
Vertically integrate `var` along `z` axis with heights specified in `ztn`
using the Trapezoidal rule. Assumes that `var` has dimensions `[x, y, z, t]`,
`[x, y, z]`, or `[z, t]`.

# Returns:
- `int_var::Array`: Integrated array with dimensions `[x,y,t]`
"""
function vert_int(var::Array{Float32, 4}, ztn)
    (nx, ny, nz, nt) = size(var)
    int_var = zeros(typeof(var[1]), (nx, ny, nt))

    @showprogress for t=1:nt
        for x=1:nx, y=1:ny
            s = 0.0
            for z in 2:nz
                s += ((var[x,y,z,t] + var[x,y,z-1,t])/2) * (ztn[z] - ztn[z-1])
            end
            int_var[x,y,t] = s

        end
    end
    return int_var
end

function vert_int(var::Array{Float32, 3}, ztn)
    (nx, ny, nz) = size(var)
    int_var = zeros(typeof(var[1]), (nx, ny))

    @showprogress for x=1:nx, y=1:ny
        s = 0.0
        for z in 2:nz
            s += ((var[x,y,z] + var[x,y,z-1])/2) * (ztn[z] - ztn[z-1])
        end
        int_var[x,y] = s

    end
    return int_var
end

function vert_int(var::Array{Float32, 2}, ztn)
    (nz, nt) = size(var)
    int_var = zeros(typeof(var[1]), (nt))

    @showprogress for t=1:nt
        s = 0.0
        for z in 2:nz
            s += ((var[z,t] + var[z-1, t])/2) * (ztn[z] - ztn[z-1])
        end
        int_var[t] = s

    end
    return int_var
end

export vert_int

"""
    list_files(dir::String)
List all RAMS data files (`*.h5`) in `dir`.
Effectively, this removes the `*.txt` header files from the built-in `readdir()` results.
"""
function list_files(dir::String)
    flist = readdir(dir, join=true)
	outlist = String[]
    for f in flist
    	m = match(r"^.*\.h5$", f)
        if !isnothing(m)
        	push!(outlist, m.match)
		end
	end
    return outlist
end
export list_files

end # module


text:
module Ocean

using ..BalanceLaws
using ..Problems

export AbstractOceanModel,
    AbstractOceanProblem,
    AbstractOceanCoupling,
    Uncoupled,
    Coupled,
    AdvectionTerm,
    NonLinearAdvectionTerm

abstract type AbstractOceanModel <: BalanceLaw end
abstract type AbstractOceanProblem <: AbstractProblem end

abstract type AbstractOceanCoupling end
struct Uncoupled <: AbstractOceanCoupling end
struct Coupled <: AbstractOceanCoupling end

abstract type AdvectionTerm end
struct NonLinearAdvectionTerm <: AdvectionTerm end

function ocean_init_state! end
function ocean_init_aux! end
function ocean_boundary_state! end

function coriolis_parameter end
function kinematic_stress end
function surface_flux end

include("OceanBC.jl")

include("HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl")
include("ShallowWater/ShallowWaterModel.jl")
include("SplitExplicit/SplitExplicitModel.jl")
include("SplitExplicit01/SplitExplicitModel.jl")
include("OceanProblems/SimpleBoxProblem.jl")

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
@testset "662.maximum-width-of-binary-tree.jl" begin
    @test width_of_binary_tree(nothing) == 0
    @test width_of_binary_tree(TreeNode{Int}([1, 3, 2, 5, 3, nothing, 9])) == 4
    @test width_of_binary_tree(TreeNode{Int}([1, 3, nothing, 5, 3])) == 2
    @test width_of_binary_tree(TreeNode{Int}([1, 3, 2, 5])) == 2
    @test width_of_binary_tree(TreeNode{Int}([1, 3, 2, 5, nothing, nothing, 9, 6, nothing,
                                              nothing, 7])) == 8
end


text:
export pairwisefst

# generate a wrapper struct so we can nicely print the results
struct PairwiseFST
    results::DataFrame
    method::String
end

# pretty-printing of FST results
function Base.show(io::IO, data::PairwiseFST)
    show(
        io,
        data.results,
        show_row_number=false,
        rowlabel = Symbol(" "),
        eltypes = false,
        row_names = names(data.results),
        title = "Pairwise FST: " * data.method
    )
end

"""
    pairwisefst(data::PopData; method::String, by::String = "global", iterations::Int64)
Calculate pairwise FST between populations in a `PopData` object. Set `iterations`
to a value greater than `0` to perform a single-tailed permutation test to obtain
P-values of statistical significance. Use `by = "locus"` to perform a locus-by-locus FST for
population pairs (iterations and significance testing ignored).
#### Methods:
- `"Hudson92"`: Hudson et al. (1992) method (only for biallelic data)
- `"Nei87"`: Nei (1987) method
- `"WC84"` : Weir & Cockerham (1984) method (default)

**Examples**
```julia
data = @nancycats
wc = pairwise_fst(data, method = "WC84")
wc_sig = pairwise_fst(data, iterations = 1000)
```
"""
function pairwisefst(data::PopData; method::String = "WC84", by::String = "global", iterations::Int64 = 0)
    if occursin("nei", lowercase(method))
        if by == "locus"
            return _pairwise_Nei_lxl(data)
        else
            iterations > 0 && return _permuted_Nei(data, iterations)
            iterations == 0 && return _pairwise_Nei(data)
        end
    elseif occursin("wc", lowercase(method))
        if by == "locus"
            throw(ArgumentError("locus-by-locus pairwise FST is not yet implemented for the WeirCockerham method." * feature_req()))
            #return _pairwise_WeirCockerham_lxl(data)
        else
            iterations > 0 && return _permuted_WeirCockerham(data, iterations)
            iterations == 0 && return _pairwise_WeirCockerham(data)
        end
    elseif occursin("hudson", lowercase(method))
        if by == "locus"
            return _pairwise_Hudson_lxl(data)
        else
            iterations > 0 && return _permuted_Hudson(data, iterations)
            iterations == 0 && return _pairwise_Hudson(data)
        end
    else
        throw(ArgumentError("please use one of \"WC84\", \"Nei87\", or \"Hudson92\" methods"))
    end
end


text:
using Test
include("cas_theory.jl")
include("cas_simplify.jl")

@test :(4a)         == @simplify 2a + a + a
@test :(a*b*c)      == @simplify a * c * b
@test :(2x)         == @simplify 1 * x * 2
@test :((a*b)^2)    == @simplify (a*b)^2
@test :((a*b)^6)    == @simplify (a^2*b^2)^3
@test :(a+b+d)      == @simplify a + b + (0*c) + d
@test :(a+b)        == @simplify a + b + (c*0) + d - d
@test :(a)          == @simplify (a + d) - d
@test :(a + b + d)  == @simplify a + b * c^0 + d
@test :(a * b * x ^ (d+y))  == @simplify a * x^y * b * x^d
@test :(a * b * x ^ 74103)  == @simplify a * x^(12 + 3) * b * x^(42^3)

@test 1 == @simplify (x+y)^(a*0) / (y+x)^0
@test 2 == @simplify cos(x)^2 + 1 + sin(x)^2
@test 2 == @simplify cos(y)^2 + 1 + sin(y)^2
@test 2 == @simplify sin(y)^2 + cos(y)^2 + 1

@test :(y + sec(x)^2 ) == @simplify 1 + y + tan(x)^2
@test :(y + csc(x)^2 ) == @simplify 1 + y + cot(x)^2

# @simplify ∂(x^2, x)

@time @simplify ∂(x^(cos(x)), x)

@test :(2x^3) == @simplify x * ∂(x^2, x) * x

# @simplify ∂(y^3, y) * ∂(x^2 + 2, x) / y * x

# @simplify (6 * x * x * y)

# @simplify ∂(y^3, y) / y

# # ex = :( ∂(x^(cos(x)), x) )
# ex = :( (6 * x * x * y) )
# g = EGraph(ex)
# saturate!(g, cas; mod=@__MODULE__)
# g.classes
# extract!(g, simplcost; root=g.root)

# params = SaturationParams(
#     scheduler=BackoffScheduler,
#     eclasslimit=5000,
#     timeout=7,
#     schedulerparams=(1000,5),
#     #stopwhen=stopwhen,
# )

# ex = :((x+y)^(a*0) / (y+x)^0)
# g = EGraph(ex)
# @profview println(saturate!(g, cas, params; mod=@__MODULE__))

# ex = extract!(g, simplcost)
# ex = rewrite(ex, canonical_t; clean=false, m=@__MODULE__)


text:
__precompile__()

#
# The MathGL module
#

module MathGL

mgl = MathGL
export mgl

#
# The Capi module
#

module Capi

#providing access to the c api
include("capi/capi.jl")

end #module Capi

# Simulating the mathgl c++ api
# Make the C functions with name mgl_fname accessible by mgl.fname
using MathGL.Capi
cmgl = MathGL.Capi

# Wrap the mathgl classes in julia types
# Up to now only the graph class and the data classes (rudimentary, only
# linked arrays in data) are supported
include("types.jl")

# Add show methods for mathgl graphs, i.e. for display in jupyter
include("multimedia.jl")

# Funktions to set options, defaults, ect.
include("setup.jl")

# Plotting function
# Auxilliary plotting functions, e.g. bounding box, legends, colorbars, ...
include("misc.jl")

# Primitives (Marks, lines, errorboxes, circles...)
include("primitives.jl")

# Text & text-positioning (along directions, curves...)
include("text.jl")

# Plot 1d data
include("1dplot.jl")

# Plot 2d data
include("2dplot.jl")

# Plot 3d data
include("3dplot.jl")

# Plotting functions that are in the "dual plot" section of the mathgl manual
include("dualplot.jl")

# Plot vector fields, flows, ect., "vector fields" section in the mathgl manual
include("vectorplot.jl")

# Special plots, contained in section "other plotting" in the mathgl manual
include("specialplot.jl")

# Macros to simulate the mgl scripting language / maybe more in the future
include("macroplot.jl")

#=# Module to draw in qt windows=#
#=module Qt=#
#=include("qtwindow.jl")=#
#=end=#

# A module for fast and (very, very) dirty usage of MathGL
# Floods the toplevel-namespace...
#=module Shortnames=#
#=include("shortnames.jl")=#
#=end #module Shortnames=#

end #module MathGL


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libigc"
version = v"1.0.10409"

# IGC depends on LLVM, a custom Clang, and a Khronos tool. Instead of building these pieces
# separately, taking care to match versions and apply Intel-specific patches where needed
# (i.e. we can't re-use Julia's LLVM_jll) collect everything here and perform a monolithic,
# in-tree build with known-good versions.

# Collection of sources required to build IGC
# NOTE: these hashes are taken from the release notes in GitHub,
#       https://github.com/intel/intel-graphics-compiler/releases
sources = [
    GitSource("https://github.com/intel/intel-graphics-compiler.git", "d5bef0c991c41e03a567187eb78fe35b6c116847"),
    GitSource("https://github.com/intel/opencl-clang.git", "50bf6d7524bea3e3d595b7c252c9ad2f91340231"),
    GitSource("https://github.com/KhronosGroup/SPIRV-LLVM-Translator.git", "cf681c88ea3f0c40d1a85677232fe6102b7e084f"),
    GitSource("https://github.com/KhronosGroup/SPIRV-Tools.git", "eeb973f5020a5f0e92ad6da879bc4df9f5985a1c"),
    GitSource("https://github.com/KhronosGroup/SPIRV-Headers.git", "ae217c17809fadb232ec94b29304b4afcd417bb4"),
    GitSource("https://github.com/intel/vc-intrinsics.git", "a7270625195501019507017aab17e50f5dac6466" #=v0.1.0=#),
    GitSource("https://github.com/llvm/llvm-project.git", "1fdec59bffc11ae37eb51a1b9869f0696bfd5312" #=llvmorg-11.1.0=#),
    # patches
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
# the build system uses git
export HOME=$(pwd)
git config --global user.name "Binary Builder"
git config --global user.email "your@email.com"

# move everything in places where it will get detected by the IGC build system
mv opencl-clang llvm-project/llvm/projects/opencl-clang
mv SPIRV-LLVM-Translator llvm-project/llvm/projects/llvm-spirv

# Work around compilation failures
# https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86678
atomic_patch -p0 patches/gcc-constexpr_assert_bug.patch
# https://reviews.llvm.org/D64388
sed -i '/add_subdirectory/i add_definitions(-D__STDC_FORMAT_MACROS)' intel-graphics-compiler/external/llvm/llvm.cmake

cd intel-graphics-compiler
install_license LICENSE.md

CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# NOTE: igc currently can't cross compile due to a variety of issues:
# - https://github.com/intel/intel-graphics-compiler/issues/131
# - https://github.com/intel/opencl-clang/issues/91
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING:BOOL=OFF)

# Explicitly use our cmake toolchain file
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})

# Silence developer warnings
CMAKE_FLAGS+=(-Wno-dev)

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux", libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("GenX_IR", :GenX_IR),
    ExecutableProduct(["iga32", "iga64"], :iga),
    LibraryProduct(["libiga32", "libiga64"], :libiga),
    LibraryProduct("libigc", :libigc),
    LibraryProduct("libigdfcl", :libigdfcl),
    # opencl-clang
    LibraryProduct("libopencl-clang", :libopencl_clang),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

# IGC only supports Ubuntu 18.04+, which uses GCC 7.4.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"8", lock_microarchitecture=false)

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil