text:
module BunchKaufmanUtils

using LinearAlgebra, StaticArrays

export pseudosolve, ispossemidef

# Piracy, move to StaticArrays
struct SBunchKaufman{N,T,DT,UT} <: Factorization{T}
    p::SVector{N,Int}
    D::DT
    U::UT
end

Base.size(F::SBunchKaufman{N}) where N = (N, N)
Base.size(F::SBunchKaufman{N}, dim::Integer) where N = dim < 3 ? N : 1

function Base.invperm(a::StaticVector{N,Int}) where N
    n = length(a)
    b = MVector{N,Int}(zeros(Int, n))
    @inbounds for (i, j) in enumerate(a)
        ((1 <= j <= n) && b[j] == 0) ||
            throw(ArgumentError("argument is not a permutation"))
        b[j] = i
    end
    return b
end

function LinearAlgebra.bunchkaufman(A::SMatrix{N,N,T}, rook::Bool=false; check=true) where {N,T}
    F = bunchkaufman(Matrix(A), rook; check=check)
    D, U = F.D, F.U
    return SBunchKaufman{N,eltype(F),typeof(D),typeof(U)}(SVector{N,Int}(F.p), D, U)
end

function Base.:*(F::Union{BunchKaufman,SBunchKaufman}, B::AbstractVecOrMat)
    D, U, p = F.D, F.U, F.p
    return permrows(U*(D*(U'*permrows(B, F.p))), invperm(F.p))
end

function LinearAlgebra.:\(F::SBunchKaufman{N}, B::Union{StaticVector{N},StaticMatrix{N}}) where N
    D, U, p = F.D, F.U, F.p
    X = U' \ (D \ (U \ permrows(B, invperm(p))))
    return similar_type(B, eltype(X))(permrows(X, p))
end

function pseudosolve(F::Union{BunchKaufman,SBunchKaufman}, B::AbstractVecOrMat; tol=eps(real(eltype(F)))*10*size(B, 1))
    D, U, p = F.D, F.U, F.p
    n = size(D, 1)
    Y = U \ permrows(B, p)
    dthresh = tol*max(maximum(abs.(D.d)), isempty(D.du) ? zero(eltype(D.du)) : maximum(abs.(D.du)))
    i = 1
    while i <= n
        if i == n || iszero(D.du[i])
            solve1!(Y, i, D.d[i], dthresh)
            i += 1
        else
            solve2!(Y, i, D.d[i], D.du[i], D.d[i+1], dthresh)
            i += 2
        end
    end
    X = U' \ Y
    return simtype(B, eltype(X))(permrows(X, invperm(p)))
end

function ispossemidef(F::Union{BunchKaufman,SBunchKaufman}; tol=eps(real(eltype(F)))*10*size(F, 1))
    D, U, p = F.D, F.U, F.p
    dthresh = tol*max(maximum(abs.(D.d)), maximum(abs.(D.du)))
    i, n = 1, size(F, 1)
    while i <= n
        if i == n || iszero(D.du[i])
            D.d[i] < -dthresh && return false
            i += 1
        else
            λ1, λ2, _ = symeig(D.d[i], D.du[i], D.d[i+1])
            (λ1 < -dthresh || λ2 < -dthresh) && return false
            i += 2
        end
    end
    return true
end

## Utilities
permrows(v::AbstractVector, p) = v[p]
permrows(M::AbstractMatrix, p) = M[p,:]

simtype(B, ::Type{T}) where T = identity
simtype(B::StaticArray, ::Type{T}) where T = similar_type(B, T)

function solve1!(y::AbstractVector, i, d, dthresh)
    if abs(d) > dthresh
        y[i] /= d
    else
        y[i] = 0
    end
    return y
end

function solve1!(Y::AbstractMatrix, i, d, dthresh)
    if abs(d) > dthresh
        for j in axes(Y, 2)
            Y[i,j] /= d
        end
    else
        for j in axes(Y, 2)
            Y[i,j] = 0
        end
    end
    return Y
end

function solve2!(y::AbstractVector, i, di, dui, di1, dthresh)
    λ1, λ2, V = symeig(di, dui, di1)
    y2 = SVector(y[i], y[i+1])
    vy = V*y2
    dinvvy = SVector(abs(λ1) >= dthresh ? vy[1]/λ1 : zero(eltype(vy)),
                     abs(λ2) >= dthresh ? vy[2]/λ2 : zero(eltype(vy)))
    x = V'*dinvvy
    y[i], y[i+1] = x
    return y
end

function solve2!(Y::AbstractMatrix, i, di, dui, di1, dthresh)
    λ1, λ2, V = symeig(di, dui, di1)
    for j in axes(Y, 2)
        y2 = SVector(Y[i,j], Y[i+1,j])
        vy = V*y2
        dinvvy = SVector(abs(λ1) >= dthresh ? vy[1]/λ1 : zero(eltype(vy)),
                         abs(λ2) >= dthresh ? vy[2]/λ2 : zero(eltype(vy)))
        x = V'*dinvvy
        Y[i,j], Y[i+1,j] = x
    end
    return Y
end

function symeig(a::Real, b::Real, d::Real)
    @assert !iszero(b)
    T, D = a + d, a*d-b^2
    if T >= 0
        λ1 = T/2 + sqrt((a-d)^2/4 + b^2)
        λ2 = D/λ1
    else
        λ2 = T/2 - sqrt((a-d)^2/4 + b^2)
        λ1 = D/λ2
    end
    n1, n2 = sqrt((λ1-d)^2 + b^2), sqrt(b^2 + (λ2-a)^2)
    V = @SMatrix [(λ1-d)/n1 b/n2; b/n1 (λ2-a)/n2]
    return λ1, λ2, V
end

end # module


text:
# fbp/z-test.jl

using Test: @test

@test cuboid_im(:test)
@test disk_phantom_params(:test)
@test ellipse_im(:test)
@test ellipsoid_im(:test)
@test ellipse_sino(:test)
@test image_geom(:test)
@test rect_im(:test)
@test rect_sino(:test)
@test rotate2d(:test)
@test sino_geom_test()

ext: jl


text:
using Documenter, Flax

push!(LOAD_PATH,  "../../src")

makedocs()


text:
@assert "this file should not be included" begin
  @assert 1 == 2
end


text:
"""
    save_div(x, y)

Savely divde `x` by `y`. If `y` is zero, return `x` directly.
"""
save_div(x, y) = ifelse(iszero(y), x, x/y)


text:
module Compiler

using Dates
import LinearAlgebra

import ..OptimizingIR
import ..BusinessDays
import ..BS
import ..CAL_BRL
import ..InterestRates
import ..Observable
import ..Konst
import ..RiskFactor
import ..ContractAttributes
import ..PricingModel
import ..DiscountFactor
import ..DiscountFactorForward
import ..ObservableAt
import ..BufferedObservable
import ..LiftObs
import ..LiftObs2
import ..ReduceObs
import ..PricingDate
import ..AdjustedDate
import ..HistoricalValue
import ..Scale
import ..Worthless
import ..Give
import ..WhenAt
import ..Both
import ..Either
import ..Unit
import ..Contract
import ..get_carry_type
import ..CARRY_TYPE_NONE
import ..CARRY_TYPE_CURVE
import ..has_carry_curve
import ..get_carry_curve_symbol
import ..Currencies
import ..StaticHedgingModel
import ..FutureValueModel
import ..SpotCurrency
import ..AbstractPricer
import ..AbstractCashflowPricer
import ..Scenario
import ..ExposureResult
import ..ForwardDiff
import ..riskfactors
import ..FixedIncomeEvent
import ..FixedIncomeContract
import ..AbstractFixedIncomeContract
import ..SwapContract
import ..price
import ..exposures
import ..MarketData
import ..Stock
import ..BlackScholesModel
import ..BinomialModelDaily
import ..FixedCashRiskFactor
import ..FixedNonCashRiskFactor
import ..exch
import ..risk_factor_symbol
import ..market_data_symbol
import ..get_functional_currency
import ..foreachcashflow
import ..event_symbol
import ..get_horizon
import ..Cashflow
import ..CashRiskFactor
import ..NonCashRiskFactor
import ..Volatility
import ..Anytime
import ..SecondOrderRiskFactor
import ..DeltaNormalExposuresMethod
import ..DeltaGammaApproxExposuresMethod
import ..EuropeanOption
import ..is_riskfree_discountfactor

include("types.jl")
include("opcodes.jl")
include("compilation.jl")
include("pricer_impl.jl")
include("model_general.jl")
include("model_statichedging.jl")
include("model_futurevalue.jl")
include("model_binomial_daily.jl")
include("model_black_scholes.jl")

end # module Compiler


text:
srand(1)
@time for i = 1:8
    n = prod(1:i)
    a = sprand(n, n, 0.1 / n)
    for j in 1:i
        @test sumabs2(SparseMatrixCSC(SparseMatrixBSC(a, j, i)) - a) == 0
        v = rand(size(a, 2))
        @test SparseMatrixBSC(a, j, i) * v == a * v
    end
end


text:
@parameters t
D = Differential(t)

@named plf = plant_face()
equations(plf)

function rep_plf(pl;name, simplify=true)
   sts = @variables x(t)
   tmp = compose(ODESystem([
       D(x) ~ pl.i_L
       ],t,sts,[];name), pl)
   simplify ? structural_simplify(tmp) : tmp
end
@named plf_rep = rep_plf(plf)

p_plf = Dict(
    plf.i_L0 => 400.0,         # g/m2 input per year (half NPP)
    plf.β_Ni0 => 25,
    plf.i_IN0 => 0,   ##<< input of mineral N,
    #
    # P from plant model parameters not used in CN-Sesam soil model
    plf.β_Pi0 => Inf, #25*20, ## leaf litter N:P ~20(massratio Kang10)
    plf.i_IP0 => Inf, #0.65,   ##<< input of mineral P, weathering: Table3 mixed sedimentary rocks 0.65g/m2/yr Hartmann14 10.1016/j.chemgeo.2013.10.025
    plf.s_EP0 => Inf, # 0.5, # plant 1/20 of typical total microbial enzyme synthesis flux
    plf.u_PlantPmax0 => Inf,
    plf.k_PlantP0 => Inf,
)
tspan = (0, 120)
prob = ODEProblem(plf_rep, [0.0], tspan, p_plf)
sol = solve(prob, saveat = 0:1:tspan[2])

@testset "increased input" begin
    is_face_period = 20 .<= sol.t .< 20+50
    @test all(sol[plf.i_L][.!is_face_period] .== 400.0)
    @test all(sol[plf.β_Ni][.!is_face_period] .== 25.0)
    @test all(sol[plf.i_L][is_face_period] .== 1.2*400.0)
    @test all(sol[plf.β_Ni][is_face_period] .== 1.2*25.0)
end;

# using Plots
# plot(sol)
# plot(sol, vars=[plf.i_L])
# plot(sol, vars=[plf.β_Ni])


text:
# Этот файл --- первый рабочий скрипт, который получил начальное приближение.
# Оставлен в оригинальном виде.
# Сейчас этот функционал предоставляется одной функцией: `initial_guess`,
# но там внутри все разбито на более мелкие функции.
# Здесь же, всё в одном месте!

using NonLinearReactionAdvectionDiffusionWithFrontData
using Dierckx
using LinearAlgebra;
using Plots; pyplot();

# -----------------------------------------------------------------------------
u_l(t) = -8;
u_r(t) =  4;
qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления,
ε = 0.01;                       # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.45;                # Область по T
N, M = 1000, 800;               # Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =     qf.(Xₙ);               # Сеточные значения коэффициента лин. усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
u₀ = u_init.(Xₙ, ε=ε, x_tp = 0.01);               # Начальные условия
nothing #hide
#
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
nothing # hide
# -----------------------------------------------------------------------------

# Набросаем эскиз прямой задачи.
# В дальнейшем, нарисуем на нем скорость фронта, которой он обладал проходя
# данный узел пространственной сетки.
psol = plot(title="Эскиз решения");
plot!(Xₙ, u[:, 1], label="u(x, 0)");
plot!(Xₙ, u[:, div(end,2)], label="u(x, 0.5)");
plot!(Xₙ, u[:, end], label="u(x, 1)");
plot!(Xₙ, ϕl[:,1], label="");
plot!(Xₙ, ϕr[:,1], label="");
plot!(Xₙ, ϕ[:,1], label="")

# Вычислим ``\frac{d f_1}{dt}`` Скорость движения фронта, как функцию времени.
# -----------------------------------------------------------------------------
df1dt = zero(f1_data);
df1dt[1] = (-3 * f1_data[1] + 2 * f1_data[2] - f1_data[3]) / (2 * τ);
for m in 2:M-2
    df1dt[m] = (f1_data[m+1] - f1_data[m-1]) / (2 * τ);
end
df1dt[M+1] = (f1_data[M+1 - 2] - 2 * f1_data[M+1 - 1] + 3 * f1_data[M+1]) / (2 * τ);
pvel = plot(title="Скорость движения переходного слоя");
plot!(Tₘ, f1_data, label="f_1");
plot!(Tₘ, df1dt, label="v_m")
#
# Последние точки некрасивые, обрежем их
x = f1_data[2:end-1];
y = df1dt[2:end-1];
tm = Tₘ[2:end-1];
#
pvel = plot(title="Скорость движения переходного слоя");
plot!(tm, x, label="f_1");
plot!(tm, y, label="v_m")
# -----------------------------------------------------------------------------

# Найдем минимальный и максимальный номер узлов сетки Xₙ, которые пересекает
# переходный слой
# -----------------------------------------------------------------------------
# .     Узлы сетки
# x     f1(t) положение переходного слоя
#                     k                                   K
# .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
#                   x x x x x x x x x x x x x x x x x x x x x
k = findfirst(z -> z > minimum(x), Xₙ)
K = findfirst(z -> z > maximum(x), Xₙ) - 1
@info "Выбор точек, которые прошел фронт" minimum(x) Xₙ[k] maximum(x) Xₙ[K]
# -----------------------------------------------------------------------------

# Найдем скорость, которой обладал фронт, проходя узлы сетки,
# т.е. найдем ``v_n`` скорость, как функцию пространственной координаты.
# -----------------------------------------------------------------------------
vspl = Spline1D(x, y);
Vn = zero(Xₙ);          # Скорость в узла, которые фронт не проходил --- 0
for n in k:K            # В остальных, найдем интерполяцией.
        Vn[n] = vspl(Xₙ[n])
end
scatter!(psol, Xₙ, Vn, label="v_n")
# -----------------------------------------------------------------------------

# Матрица системы
# -----------------------------------------------------------------------------
Al = LowerTriangular(h * ones(K-k+1, K-k+1));
Ar = UpperTriangular(-h * ones(K-k+1, K-k+1));
Al[:, 1]    .= (k-1)*h;
Ar[:, end]  .= (N+1-K)*h;
# Повторим чертеж для объяснения двух прошлых строк
# N = 17, тогда нумерация массива будет следующей
# k = 6, K = 16
# Интервалов лежащих слева от k -- 5, (k-1)
# Интервалов лежащих справа от K -- 2, (N+1 - K)
# 1   2   3   4   5   6   7   8   9   11  12  13  14  15  16  17  18
#                     k                                   K
# .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
#                   x x x x x x x x x x x x x x x x x x x x x
A = -(Al + Ar) / 2.0;
# -----------------------------------------------------------------------------

α = 0.001;      # Эмпирический параметр регуляризации
B = Vn[k:K];    # Вектор правой части
B = B .- ( -(ulₘ[1] + ulₘ[1])  / 2.0);          # Константы от интегралов
qx = ( transpose(A) * A + α*I) \ transpose(A)*B
l = [qx[1] for i in 1:k-1];         # Дополняем слева
r = [qx[end-1] for i in K+1:N+1];   # и справа.
qx[end] = qx[end-1];                # Эти точки всегда плохо восстанавливаются,
qx[2] = qx[1];                      # обрежем их.
q_guess = [l; qx; r];
plot(Xₙ, qₙ, label="Истинное")
plot!(Xₙ, q_guess, label="Начальное приближение")

savefig("guess.png")


text:
__precompile__()

module SimpleDiffEq

using Reexport
@reexport using DiffEqBase
using StaticArrays

include("functionmap.jl")
include("euler_maruyama.jl")
include("tsit5/atsit5_cache.jl")
include("tsit5/tsit5.jl")
include("tsit5/atsit5.jl")
include("tsit5/gpuatsit5.jl")

end # module


text:
"""
    num_chunks(::Val{B})

Determine the number and type of chunks needed to store `B` bits.
"""
function num_chunks(::Val{B}) where {B}
    if B ≤ 0
        throw(ArgumentError("`B` must be positive!"))
    elseif B ≤ 8
        return 1, UInt8
    elseif B ≤ 16
        return 1, UInt16
    elseif B ≤ 32
        return 1, UInt32
    else
        return (B - 1) ÷ 64 + 1, UInt64
    end
end

"""
    check_bitstring_typeparams(::Val{B}, ::Val{N})

Check if number of bits `B` is consistent with number of chunks `N`. Throw an error if not.
"""
function check_bitstring_typeparams(::Val{B}, ::Val{N}, ::Type{UInt64}) where {B,N}
    if B > N * 64
        error("$B bits do not fit into $N 64-bit chunks")
    elseif B ≤ (N - 1) * 64
        error("$B bits fit into $(N - 1) 64-bit chunks, but $N chunks were provided")
    end
end
function check_bitstring_typeparams(::Val{B}, ::Val{1}, ::Type{T}) where {B,T}
    if B > sizeof(T) * 8
        error("$B bits do not fit into a $(sizeof(T) * 8)-bit chunk")
    end
end
function check_bitstring_typeparams(::Val{B}, ::Val{1}, ::Type{UInt64}) where {B}
    if B > 64
        error("$B bits do not fit into a 64-bit chunk")
    end
end
function check_bitstring_typeparams(::Val{B}, ::Val{N}, ::Type{T}) where {B,N,T}
    error("Only `UInt64` is supported for multi-bit chunks")
end

"""
    BitString{B,N,T<:Unsigned}

Type for storing bitstrings of static size. Holds `B` bits in `N` chunks, where each chunk is
of type `T`.

`N` is chosen automatically to accommodate `B` bits as efficiently as possible.

# Constructors

* `BitString{B,N,T}(::SVector{N,T})`: unsafe constructor. Does not check for ghost bits.

* `BitString{B,N,T}(i::T)`: as above, but sets `i` as the rightmost chunk.

* `BitString{B}(::Integer)`: Convert integer to `BitString`. Integer is truncated to the
  correct number of bits.

"""
struct BitString{B,N,T<:Unsigned}
    chunks::SVector{N,T}

    # This constructor is only to be used internally. It doesn't check for ghost bits.
    function BitString{B,N,T}(s::SVector{N,T}) where {B,N,T}
        check_bitstring_typeparams(Val(B), Val(N), T)
        return new{B,N,T}(s)
    end
    function BitString{B,N,T}(i::T) where {B,N,T<:Unsigned}
        check_bitstring_typeparams(Val(B), Val(N), T)
        return new{B,N,T}(setindex(zero(SVector{N,UInt64}), i, N))
    end
end

###
### Basic properties.
###
"""
    num_chunks(::Type{<:BitString})
    num_chunks(s::BitString)

Number of chunks in bitstring. Equivalent to `length(chunks(s))`.
"""
num_chunks(::Type{<:BitString{<:Any,N}}) where {N} = N

"""
    chunk_type(::Type{<:BitString})
    chunk_type(s::BitString)

Type of unsigned integer used to store the chunks.
"""
chunk_type(::Type{<:BitString{<:Any,<:Any,T}}) where {T} = T

"""
    num_bits(::Type{<:BitString})
    num_bits(s::BitString)

Total number of bits stored in bitstring.
"""
num_bits(::Type{<:BitString{B}}) where {B} = B

"""
    top_chunk_bits(::Type{<:BitString})
    top_chunk_bits(s::BitString)

Number of bits stored in top chunk. Equivalent to `chunk_bits(s, 1)`.
"""
function top_chunk_bits(::Type{<:BitString{B}}) where B
    return B % 64 == 0 ? 64 : B % 64
end

for f in (:num_chunks, :chunk_type, :num_bits, :top_chunk_bits)
    @eval $f(s::BitString) = $f(typeof(s))
end

"""
    chunks(s::BitString)

`SVector` that stores the chunks of `s`.
"""
chunks(s::BitString) = s.chunks

"""
    chunks_bits(::Type{<:BitString}, i)
    chunks_bits(s, i)

Number of bits in the `i`-th chunk of `s`.
"""
chunk_bits(s, i) = chunk_bits(typeof(s), i)
chunk_bits(::Type{<:BitString{B,1}}, _) where {B} = B
function chunk_bits(::Type{S}, i) where {S<:BitString}
    return ifelse(i == 1, top_chunk_bits(S), 64)
end

function ghost_bit_mask(::Type{S}) where S<:BitString
    T = chunk_type(S)
    unused_bits = sizeof(T) * 8 - top_chunk_bits(S)
    return ~zero(T) >>> unused_bits
end

"""
    remove_ghost_bits(s::BitString)

Remove set bits outside data field if any are present.

See also: [`has_ghost_bits`](@ref).
"""
function remove_ghost_bits(s::S) where {S<:BitString}
    mask = ghost_bit_mask(S)
    return S(setindex(s.chunks, s.chunks[1] & mask, 1))
end

@inline function remove_ghost_bits(s::S) where {S<:BitString{<:Any,1}}
    mask = ghost_bit_mask(S)
    return S(chunks(s) .& mask)
end

"""
    has_ghost_bits(s::BitString)

Check for bits outside data field.

See also: [`remove_ghost_bits`](@ref).
"""
function has_ghost_bits(s::S) where {S<:BitString}
    top = first(chunks(s))
    mask = ~zero(UInt64) << top_chunk_bits(S)
    return top & mask > 0
end

###
### Alternative/useful constructors. These are not super efficient, but they are safe.
###
function BitString{B}(i::Union{Int128,Int64,Int32,Int16,Int8}) where {B}
    return remove_ghost_bits(BitString{B}(unsigned(i)))
end
function BitString{B}(i::Union{UInt64,UInt32,UInt16,UInt8}) where {B}
    N, T = num_chunks(Val(B))
    s = setindex(zero(SVector{N,T}), T(i), N)
    return remove_ghost_bits(BitString{B,N,T}(s))
end
function BitString{B}(i::UInt128) where {B}
    N, T = num_chunks(Val(B))
    left = i >>> 0x40 % T # left will only be used if T == UInt64 and N > 1
    right = i  % T
    s = ntuple(Val(N)) do i
        i == N ? right : i == N - 1 ? left : zero(T)
    end
    return remove_ghost_bits(BitString{B,N,T}(SVector{N,T}(s)))
end
function BitString{B}(i::BigInt) where {B}
    N, T = num_chunks(Val(B))
    s = zero(SVector{N,T})
    j = N
    while i ≠ 0
        chunk = i & typemax(T) % T
        i >>>= 64 # Can use 64 here, as only 1-chunk addresses can be smaller
        s = setindex(s, chunk, j)
        j -= 1
    end
    return remove_ghost_bits(BitString{B,N,T}(s))
end

function Base.zero(S::Type{<:BitString{B}}) where {B}
    N, T = num_chunks(Val(B))
    BitString{B,N,T}(zero(SVector{N,T}))
end
Base.zero(s::BitString) = zero(typeof(s))

function Base.show(io::IO, s::BitString{B,N}) where {B,N}
    str = join(map(i -> repr(i)[3:end], s.chunks), '_')

    print(io, "BitString{$B}(big\"0x", str, "\")")
end
Base.bitstring(s::BitString{B}) where {B} = join(bitstring.(s.chunks))[(end - B + 1):end]

###
### Operations on BitStrings
###
for op in (:⊻, :&, :|)
    @eval (Base.$op)(l::S, r::S) where S<:BitString = S($op.(l.chunks, r.chunks))
end
Base.:~(s::S) where S<:BitString = remove_ghost_bits(S(.~(s.chunks)))

Base.count_ones(s::BitString) = sum(count_ones, s.chunks)
Base.count_zeros(s::BitString) = num_bits(s) - count_ones(s)

function _trailing(f, s::BitString)
    result = 0
    i = 0
    # Idea: if all whole chunk is the same digit, you have to look at the next one.
    # This gets compiled away if N=1
    for i in num_chunks(s):-1:1
        r = f(s.chunks[i])
        result += r
        r == chunk_bits(s, i) || break
    end
    # If top chunk occupies the whole integer, result will always be smaller or equal to B.
    if f ≢ trailing_ones && top_chunk_bits(s) ≠ 64
        return min(num_bits(s), result)
    else
        return result
    end
end

function _leading(f, s::BitString)
    N = sizeof(chunk_type(s)) * 8
    # First chunk is a special case - we have to ignore the empty space before the string.
    result = min(f(s.chunks[1] << (N - top_chunk_bits(s))), top_chunk_bits(s))

    # This gets compiled away if N=1
    if num_chunks(s) > 1 && result == top_chunk_bits(s)
        for i in 2:num_chunks(s)
            r = f(s.chunks[i])
            result += r
            r == 64 || break
        end
    end
    return result
end

Base.trailing_ones(s::BitString) = _trailing(trailing_ones, s)
Base.trailing_zeros(s::BitString) = _trailing(trailing_zeros, s)
Base.leading_ones(s::BitString) = _leading(leading_ones, s)
Base.leading_zeros(s::BitString) = _leading(leading_zeros, s)

@generated function _right_shift(s::S, k) where {S<:BitString}
    N = num_chunks(S)
    quote
        $(Expr(:meta, :inline))
        # equivalent to d, r = divrem(k, 64)
        d = k >>> 0x6
        r = k & 63
        ri = 64 - r
        mask = ~zero(UInt64) >>> ri # 2^r-1 # 0b0...01...1 with `r` 1s
        c = chunks(s)

        @nif $(N + 1) l -> (d < l) l -> (
            S(SVector((@ntuple l - 1 k -> zero(UInt64))... ,c[1] >>> r,
                      (@ntuple $N-l q -> (c[q + 1] >>> r | ((c[q] & mask) << ri)))...
                      ))
        ) l -> (
            return zero(S)
        )
    end
end

function _left_shift(s::S, k) where {S<:BitString}
    result = zeros(MVector{num_chunks(S),UInt64})
    # d, r = divrem(k, 64)
    d = k >>> 0x6
    r = k & 63

    shift = s.chunks .<< (r % UInt64)
    carry = s.chunks .>>> ((64 - r) % UInt64)

    for i in d + 1:length(result)
        @inbounds result[i - d] = shift[i] | get(carry, i + 1, zero(UInt64))
    end
    # This bit removes ghost bits.
    result[1] &= ghost_bit_mask(S)
    return S(SVector(result))
end

Base.:>>(s::BitString, k) = k ≥ 0 ? _right_shift(s, k) : _left_shift(s, -k)
Base.:<<(s::BitString, k) = k > 0 ? _left_shift(s, k) : _right_shift(s, -k)
Base.:>>>(s::BitString, k) = s >> k

# remove ghost bits must be applied to both because k might be negative.
Base.:>>(s::S, k) where S<:BitString{<:Any,1} = remove_ghost_bits(S(s.chunks .>> k))
Base.:>>(s::S, k::Unsigned) where S<:BitString{<:Any,1} = S(s.chunks .>> k)
Base.:<<(s::S, k) where S<:BitString{<:Any,1} = remove_ghost_bits(S(s.chunks .<< k))

function Base.isless(s1::B, s2::B) where {B<:BitString}
    for i in 1:num_chunks(B)
        if chunks(s1)[i] ≠ chunks(s2)[i]
            return chunks(s1)[i] < chunks(s2)[i]
        end
    end
    return false
end
Base.isodd(s::BitString) = isodd(chunks(s)[end])
Base.iseven(s::BitString) = iseven(chunks(s)[end])

# For compatibility. Changing any of the hashes will slightly change results and make the
# tests fail.
Base.hash(b::BitString{<:Any,1}, h::UInt) = hash(b.chunks[1], h)
Base.hash(b::BitString, h::UInt) = hash(b.chunks.data, h)

"""
    partial_left_shift(bs::BitString, i, j)

Shift a part of the bitstring left by one place with boundaries `i < j`.
In a `BoseFS` bitstring, it moves a particle at offset `i` to the position at
offset `j`.

See also: [`excitation`](@ref), [`partial_right_shift`](@ref).
"""
function partial_left_shift(chunk::T, i, j) where {T<:Unsigned}
    # Mask of one spanning from i to j
    mask = (T(1) << T(j - i + 1) - T(1)) << T(i)
    # Shift the part of the string that needs to be shifted, ensure a one is added at the end
    # swap shift to move in other direction
    #println(bitstring(mask))
    shifted_part = ((chunk & mask) << 0x1) & mask
    # Leave the rest intact
    intact_part = chunk & ~mask

    return shifted_part | intact_part | T(1) << T(i)
end

"""
    partial_right_shift(bs::BitString, i, j)

Shift a part of the bitstring right by one place with boundaries `i < j`.
In a `BoseFS` bitstring, it moves a particle at offset `j` to the position at
offset `i`.

See also: [`partial_left_shift`](@ref), [`excitation`](@ref).
"""
function partial_right_shift(chunk::T, i, j) where {T<:Unsigned}
    # Mask of one spanning from i to j
    mask = (T(1) << T(j - i + 1) - T(1)) << T(i)
    # Shift the part of the string that needs to be shifted, ensure a one is added at the end
    # swap shift to move in other direction
    shifted_part = ((chunk & mask) >> 0x1) & mask
    # Leave the rest intact
    intact_part = chunk & ~mask
    #println(lpad("↑" * " "^j, length(bitstring(chunk))))

    return shifted_part | intact_part | T(1) << T(j)
end

function partial_left_shift(bs::S, i, j) where {S<:BitString{<:Any,1}}
    return S(partial_left_shift(bs.chunks[1], i, j))
end

function partial_right_shift(bs::S, i, j) where {S<:BitString{<:Any,1}}
    return S(partial_right_shift(bs.chunks[1], i, j))
end

function partial_left_shift(bs::S, i, j) where {N,S<:BitString{<:Any,N}}
    result = MVector(bs.chunks)
    lo_idx = N - (i >>> 0x6)
    hi_idx = N - (j >>> 0x6)
    lo_off = i & 63
    hi_off = j & 63
    @inbounds if hi_idx == lo_idx
        result[hi_idx] = partial_left_shift(result[hi_idx], lo_off, hi_off)
    else
        # Top part first.
        chunk = result[hi_idx]
        chunk = partial_left_shift(chunk, 0, hi_off)
        # Carry bit.
        chunk &= -UInt(1) << 0x1
        chunk |= result[hi_idx + 1] >> 63
        result[hi_idx] = chunk

        idx = hi_idx + 1
        while idx < lo_idx
            chunk = result[idx]
            chunk <<= 0x1
            chunk |= result[idx + 1] >> 63
            result[idx] = chunk
            idx += 1
        end

        # Bottom part.
        chunk = result[lo_idx]
        chunk = partial_left_shift(chunk, lo_off, 64)
        result[lo_idx] = chunk
    end
    return S(SVector(result))
end

function partial_right_shift(bs::S, i, j) where {N,S<:BitString{<:Any,N}}
    result = MVector(bs.chunks)
    lo_idx = N - (i >>> 0x6)
    hi_idx = N - (j >>> 0x6)
    lo_off = i & 63
    hi_off = j & 63
    @inbounds if hi_idx == lo_idx
        result[hi_idx] = partial_right_shift(result[hi_idx], lo_off, hi_off)
    else
        # Bottom first
        chunk = result[lo_idx]
        chunk = partial_right_shift(chunk, lo_off, 64)
        # Carry bit.
        chunk &= -UInt(1) >> 0x1
        chunk |= result[lo_idx - 1] << 63
        result[lo_idx] = chunk

        idx = lo_idx - 1
        while idx > hi_idx
            chunk = result[idx]
            chunk >>= 0x1
            chunk |= result[idx - 1] << 63
            result[idx] = chunk
            idx -= 1
        end

        # Top part.
        chunk = result[hi_idx]
        chunk = partial_right_shift(chunk, 0, hi_off)
        result[hi_idx] = chunk
    end
    return S(SVector(result))
end

function bit_mask(bs::BitString, i, j)

end


text: @testset "prior" begin end


text:
import UnPack

opencsv(path) = open(read, path, CSV_ENC)

# FIXED: Переделанная версия '@unpack' с возможностью добавления суффикса к названию переменных
macro unpack(args)
    args.head != :(=) && error("Expression needs to be of form `a, b = c`")
    items, suitecase = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suffix, items... = isa(first(items), AbstractString) ? items : ("", items...)
    suitecase_instance = gensym()
    kd = map(items) do key
        var = Symbol(string(key) * suffix)
        :( $var = $UnPack.unpack($suitecase_instance, Val{$(Expr(:quote, key))}()) )
    end
    kdblock = Expr(:block, kd...)
    expr = quote
        local $suitecase_instance = $suitecase # handles if suitecase is not a variable but an expression
        $kdblock
        $suitecase_instance # return RHS of `=` as standard in Julia
    end
    esc(expr)
end

x_abschange(x, x_previous) = maxdiff(x, x_previous)
maxdiff(x, y) = mapreduce((a, b) -> abs(a - b), max, x, y)

function assess_convergence(x, x_prev, f_x, x_tol, f_tol)
    x_converged, f_converged = false, false
    if x_abschange(x, x_prev) ≤ x_tol * maximum(abs, x)
        x_converged = true
    end
    if maximum(abs, f_x) ≤ f_tol
        f_converged = true
    end
    return x_converged, f_converged
end

getparams(df::AbstractDataFrame, ::Val{:tanks}) = @view df[df.Parameter .∉ Ref((:Gw, :Jinj, :Jp)), :]
getparams(df::AbstractDataFrame, ::Val{:wells}) = @view df[df.Parameter .∈ Ref((:Gw, :Jinj, :Jp)), :]
getparams(df::AbstractDataFrame, ::Val{:tanks}, ::Val{:var}) = @view df[(df.Parameter .∉ Ref((:Gw, :Jinj, :Jp))) .& .!df.Const, :]
getparams(df::AbstractDataFrame, ::Val{:tanks}, ::Val{:const}) = @view df[(df.Parameter .∉ Ref((:Gw, :Jinj, :Jp))) .& df.Const, :]
getparams(df::AbstractDataFrame, ::Val{:wells}, ::Val{:var}) = @view df[(df.Parameter .∈ Ref((:Gw, :Jinj, :Jp))) .& .!df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{:wells}, ::Val{:const}) = @view df[(df.Parameter .∈ Ref((:Gw, :Jinj, :Jp))) .& df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{S}) where {S} = @view df[df.Parameter .=== S, :]
getparams(df::AbstractDataFrame, ::Val{S}, ::Val{:var}) where {S} = @view df[(df.Parameter .=== S) .& .!df.Const, :]
getparams(df::AbstractDataFrame, ::Val{S}, ::Val{:const}) where {S} = @view df[(df.Parameter .=== S) .& df.Const, :]
getparams(df::AbstractDataFrame, ::Val{:Gw}, ::Val{:var}) = @view df[(df.Parameter .=== :Gw) .& .!df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{:Jinj}, ::Val{:var}) = @view df[(df.Parameter .=== :Jinj) .& .!df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{:Jp}, ::Val{:var}) = @view df[(df.Parameter .=== :Jp) .& .!df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{:Gw}, ::Val{:const}) = @view df[(df.Parameter .=== :Gw) .& df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{:Jinj}, ::Val{:const}) = @view df[(df.Parameter .=== :Jinj) .& df.Ignore, :]
getparams(df::AbstractDataFrame, ::Val{:Jp}, ::Val{:const}) = @view df[(df.Parameter .=== :Jp) .& df.Ignore, :]


text:
using Caesar, Gadfly, YAML, DelimitedFiles
using Cairo, Fontconfig

#Test with sines
xin = 0:0.001:2*pi
raw = sin.(100*collect(xin)) + sin.(10*collect(xin)) + sin.(collect(xin));
f1 = 80; f2=120;
cztlen=1024;
Fs = 1000;
m = nextpow(2,length(raw))
signal = zeros(m);
signal[1:length(raw)] = raw;

#Test with raw data
#cfgd = YAML.load(open(joinpath(dirname(pathof(Caesar)),"config","SAS2D.yaml")))
#dataFile = joinpath(ENV["HOME"],"data","sas","sample_data","waveform1600.csv");
#raw = readdlm(dataFile,',',Float64,'\n')
#f1 = cfgd["freq_lower_cutoff"]; f2=cfgd["freq_upper_cutoff"];
#Fs = cfgd["freq_sampling"];
#cztlen = cfgd["nfft_czt"];
#m = nextpow(2,size(raw,2)
#signal = zeros(m);
#signal[1:size(raw,2)] = raw[1,:];

w = exp(-2im*pi*(f2-f1)/m/Fs)
myczt = Caesar.prepCZTFilter(m,1,w,cztlen)
out = zeros(Complex{Float64}, cztlen)
@time myczt(collect(signal),out)

pl = Gadfly.plot(
 Gadfly.layer(y=real(signal), Geom.line),
 Gadfly.Coord.cartesian(xmin=0, xmax=m),
 Guide.ylabel(nothing),Guide.xlabel("Samples")
 ) ; pl |> PDF("/tmp/test1.pdf")

 pl = Gadfly.plot(
  Gadfly.layer(y=real(out),Geom.line, Theme(default_color="red")),
  Gadfly.Coord.cartesian(xmin=0, xmax=1024),
  Guide.ylabel(nothing),Guide.xlabel("Samples")
  ) ; pl |> PDF("/tmp/test2.pdf")


text:
@static if VERSION >= v"0.7.0-"
    using REPL
    using REPL: LineEdit
else
    using Base: REPL, LineEdit
end

"""
    afterreplinit(f)

Like `atreplinit` but triggers `f` even after REPL is initialized when
it is called.
"""
function afterreplinit(f)
    # See: https://github.com/JuliaLang/Pkg.jl/blob/v1.0.2/src/Pkg.jl#L338
    function wrapper(repl)
        if isinteractive() && repl isa REPL.LineEditREPL
            f(repl)
        end
    end
    if isdefined(Base, :active_repl)
        wrapper(Base.active_repl)
    else
        atreplinit() do repl
            @async begin
                wait_repl_interface(repl)
                wrapper(repl)
            end
        end
    end
end

function wait_repl_interface(repl)
    for _ in 1:20
        try
            repl.interface.modes[1].keymap_dict
            return
        catch
        end
        sleep(0.05)
    end
end

# Register keybind '.' in Julia REPL:

function on_dot_press(s, _...)
    if isempty(s) || position(LineEdit.buffer(s)) == 0
        @static if VERSION >= v"0.7.0-"
            start_ipython()
        else
            # Force current_module() inside IPython to be Main:
            Base.eval(Main, :($start_ipython()))
        end
        println()
        LineEdit.refresh_line(s)
    else
        LineEdit.edit_insert(s, '.')
    end
end

function init_repl(repl)
    ipy_prompt_keymap = Dict{Any,Any}('.' => on_dot_press)

    main_mode = repl.interface.modes[1]
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict,
                                                  ipy_prompt_keymap)
end
# See: https://github.com/JuliaInterop/RCall.jl/blob/master/src/RPrompt.jl


text:
"""
    li3(z::ComplexF64)::ComplexF64

Returns the complex trilogarithm of a complex number `z` of type `ComplexF64`.

Author: Alexander Voigt

License: MIT

# Example
```julia
li3(1.0 + 1.0im)
```
"""
function li3(z::ComplexF64)::ComplexF64
    function clog(z)
        az::Float64 = angle(z)
        return 0.5*log(abs2(z)) + (imag(z) == 0.0 && az < 0.0 ? -az : az)*1.0im
    end

    z2::Float64 = 1.6449340668482264
    z3::Float64 = 1.2020569031595943
    bf = (
        1.0, -3.0/8.0, 17.0/216.0, -5.0/576.0,
        1.2962962962962963e-04,  8.1018518518518519e-05,
       -3.4193571608537595e-06, -1.3286564625850340e-06,
        8.6608717561098513e-08,  2.5260875955320400e-08,
       -2.1446944683640648e-09, -5.1401106220129789e-10,
        5.2495821146008294e-11,  1.0887754406636318e-11,
       -1.2779396094493695e-12, -2.3698241773087452e-13,
        3.1043578879654623e-14,  5.2617586299125061e-15
    )

    if imag(z) == 0.0
        if real(z) == 0.0
            return 0.0 + 0.0im
        end
        if real(z) == 1.0
            return z3 + 0.0im
        end
        if real(z) == -1.0
            return -0.75*z3 + 0.0im
        end
        if real(z) == 0.5
            return 0.53721319360804020 + 0.0im
        end
    end

    nz::Float64  = abs2(z)
    pz::Float64  = angle(z)
    lnz::Float64 = 0.5*log(nz)

    if lnz*lnz + pz*pz < 1.0 # |log(z)| < 1
        v::ComplexF64  = lnz + pz*im
        v2::ComplexF64 = v*v
        v4::ComplexF64 = v2*v2
        v8::ComplexF64 = v4*v4
        c0::ComplexF64 = z3 + v*(z2 - v2/12.0)
        c1::ComplexF64 = 0.25 * (3.0 - 2.0*clog(-v))

        cs = (
            -3.4722222222222222e-03, 1.1574074074074074e-05,
            -9.8418997228521038e-08, 1.1482216343327454e-09,
            -1.5815724990809166e-11, 2.4195009792525152e-13,
            -3.9828977769894877e-15
        )

        return c0 +
            c1*v2 +
            v4*(cs[1] + v2*cs[2]) +
            v8*(cs[3] + v2*cs[4] + v4*(cs[5] + v2*cs[6])) +
            v8*v8*cs[7]
    end

    (u::ComplexF64, rest::ComplexF64) = if nz <= 1.0
        (-clog(1.0 - z), 0.0 + 0.0im)
    else # nz > 1.0
        arg::Float64 = pz > 0.0 ? pz - pi : pz + pi
        lmz::ComplexF64 = lnz + arg*im # clog(z)
        (-clog(1.0 - 1.0/z), -lmz*(lmz*lmz/6.0 + z2))
    end

    u2::ComplexF64 = u*u
    u4::ComplexF64 = u2*u2
    u8::ComplexF64 = u4*u4

    rest +
    u*bf[1] +
    u2*(bf[2] + u*bf[3]) +
    u4*(bf[4] + u*bf[5] + u2*(bf[6] + u*bf[7])) +
    u8*(bf[8] + u*bf[9] + u2*(bf[10] + u*bf[11]) +
        u4*(bf[12] + u*bf[13] + u2*(bf[14] + u*bf[15]))) +
    u8*u8*(bf[16] + u*bf[17] + u2*bf[18])
end


text:
# 3 main data structures: Node, Element, Structure

# defining all nodes: Node([position], type, boundary conditions)
node1 = Node([72.0, 0.0, 0.0], :truss, :yfixed)
node2 = Node([0.0, 36.0, 0.0], :truss, :fixed)
node3 = Node([0.0, 36.0, 72.0], :truss, :fixed)
node4 = Node([0.0, 0.0, -48.0], :truss, :fixed)

# The following is equivalent: Node([position], [DOFs])
node1 = Node([72.0, 0.0, 0.0], [true, false, true])
node2 = Node([0.0, 36.0, 0.0], [false, false, false])
node3 = Node([0.0, 36.0, 72.0], [false, false, false])
node4 = Node([0.0, 0.0, -48.0], [false, false, false])

# group nodes together
nodes = [node1, node2, node3, node4]

# Material properties
E = 1.2e6

# defining elements
e1 = Element(nodes, [1,2], E, 0.302)
e2 = Element(nodes, [1,3], E, 0.729)
e3 = Element(nodes, [1,4], E, 0.187)

# group elements
elements = [e1, e2, e3]

# loads
l1 = Load(nodes, node1.position, [0.0,0.0,-1000.0])

# group loads (in this case just 1)
loads = [l1]

# assemble structure
structure = Structure(nodes, elements, loads)

# analyze structure
analyze(structure)

# relevant information is now stored in the Structure datatype, e.g.
displacements = structure.U

# the Geometry type stores things like displaced positions for visualization
geo = Geometry(structure)

# structurePlot is used to plot the Geometry type
fig = structurePlot(geo)

max_stars_repo_name: keithjlee/Asap

max_issues_repo_name: keithjlee/Asap

max_forks_repo_name: keithjlee/Asap


text:
# This file contains code to add extra records using annotation.tbl and
# a FASTA file.
#=
"""
    NCBIInfluenzaData.Extra

This module contains functionality to create new `SegmentData` objects from existing
sequences not from NCBI, such that influenza sequences from other sources can be
parsed to the same data structure as those from NCBI.

It requires two files: A FASTA file where the header is of the following format:
"IDENTIFIER|HOST|SEGMENT|SEROTYPE|YEAR|ISOLATE".
* Identifier can be any string, not including `|` or trailing/leading whitespace
* Host can be "Human", "Swine", "Avian", or "Other", or in lowercase.
* Segment must be parsable into InfluenzaCore.Segment
* Serotype can be the empty string or parsable into InfluenzaCore.SeroType
* Year must be 4 characters in 0-9
* Isolate must be the standardized e.g. "A/Denmark/1/2021"

The second file is an annotation.tbl file produced by The NCBI Influenza Virus
Sequence Annotation Tool. The headers in the annotation.tbl file must match the
headers in the FASTA file.
"""

function try_parse_fasta_header(
    ::Type{IncompleteSegmentData},
    fields::Vector{SubString{String}},
    line::Union{String, SubString{String}}
)::Option{IncompleteSegmentData}
    # Check there are 6 fields
    (@? try_split!(fields, line, UInt8('|'))) == 6 || return none
    id, s_host, s_segment, s_serotype, s_year, isolate = fields
    # Check ID is not empty
    isempty(id) && return none
    # Check host
    host = parse(Host, s_host)
    host in (Hosts.human, Hosts.avian, Hosts.swine) || return none
    # Check segment
    segment::Segment = let
        s = tryparse(Segment, s_segment)
        s === nothing && return none
        s
    end
    # Check serotype
    serotype::Option{SeroType} = let
        if isempty(s_serotype)
            none(SeroType)
        else
            s = tryparse(SeroType, s_serotype)
            s === nothing && return none
            some(s)
        end
    end
    # Check year
    year = let
        s = tryparse(Int16, s_year)
        s === nothing && return none
        s
    end
    # Check isolate format
    validate_isolate_format(isolate) || return none
    return some(IncompleteSegmentData(
        id, host, segment, serotype, year, isolate, ProteinORF[], none(LongDNASeq)
    ))
end

function try_parse_fasta_header(::Type{IncompleteSegmentData}, line::Union{String, SubString{String}})
    try_parse_fasta_header(IncompleteSegmentData, Vector{SubString{String}}(undef, 6), line)
end

struct Unbounded end
const FeatureRange = Tuple{Union{Unbounded, UInt}, Union{Unbounded, UInt}}
struct Feature
    key::String
    ranges::Vector{FeatureRange}
    qualifiers::Dict{String, String}
end

# This format is absolutely horrendous. The parser here is not a generic parser for
# genbank formats.
function try_parse_genbank_feature_table(lines::Vector{SubString{String}})::Option{Vector{Feature}}
    # First split it into individual features
    # The features consists of N "headers" where the first field is not a tab
    # followed by M lines where the field IS a tab.
    fieldbuffer = Vector{SubString{String}}(undef, 5)
    result = Vector{Feature}()
    feature_groups = Vector{Vector{SubString{String}}}()
    filter!(!isempty, lines)
    start = 1
    last_started_tab = false
    for (i, line) in enumerate(lines)
        this_starts_tab = codeunit(line, 1) === UInt8('\t')
        if last_started_tab && !this_starts_tab
            push!(feature_groups, lines[start:i-1])
            start = i
        end
        last_started_tab = this_starts_tab
    end
    push!(feature_groups, lines[start:end])

    # Then parse each of the features
    for feature in feature_groups
        push!(result, @? try_parse_genbank_feature(fieldbuffer, feature))
    end

    return some(result)
end

# Example:
# <3522    3572    CDS
# 3706    >4197
#                         product        Yip2p
#                         prot_desc      similar to human polyposis locus protein 1 (YPD)
function try_parse_genbank_feature(
    v::Vector{SubString{String}},
    lines::Vector{SubString{String}}
)::Option{Feature}
    length(lines) < 2 && return none

    ranges = FeatureRange[]
    is_first_line = true
    line_num = 1
    key = nothing
    n_fields = (@? try_split!(v, lines[line_num], UInt8('\t')))

    # While the first field is not empty, we are still parsing the feature range
    while !isempty(@inbounds v[1])
        line_num += 1
        line_num > length(lines) && return none
        if is_first_line
            n_fields < 3 && return none
            key = @inbounds v[3]
            isempty(key) && return none
            is_first_line = false
        else
            n_fields > 2 && !isempty(v[3]) && return none
        end
        fst = @? try_parse_feature_number(@inbounds(v[1]), true)
        lst = @? try_parse_feature_number(@inbounds(v[2]), false)
        push!(ranges, (fst, lst))
        n_fields = (@? try_split!(v, lines[line_num], UInt8('\t')))
    end

    # Parse the qualifier lines
    qualifiers = Dict{String, String}()
    while line_num ≤ length(lines)
        n_fields = (@? try_split!(v, lines[line_num], UInt8('\t')))
        n_fields == 5 || return none
        all(isempty(v[i]) for i in 1:3) || return none
        qualkey, qualval = (@inbounds v[4]), (@inbounds v[5])
        (isempty(qualkey) || isempty(qualval)) && return none
        qualifiers[qualkey] = qualval
        line_num += 1
    end

    return some(Feature(key, ranges, qualifiers))
end

function try_parse_feature_number(
    s::Union{SubString{String}, String},
    isfirst::Bool
)::Option{Union{Unbounded, UInt}}
    isempty(s) && return none
    firstbyte = codeunit(s, 1)
    if (isfirst  & (firstbyte === UInt8('<'))) ||
       (!isfirst & (firstbyte === UInt8('>')))
        return some(Unbounded())
    end
    n = tryparse(UInt, s)
    n === nothing && return none
    return some(n::UInt)
end

function fill_orfs_extra_records!(records::Dict{String, IncompleteSegmentData}, annotpath::String)
    chunks = maybe_gzip(annotpath) do file
        partition_chunks(file)
    end
    filled_ids = Set{String}()
    for chunk in chunks
        parse_annot_segment_data!(records, filled_ids, chunk)
    end

    # If any were missing, error here
    if length(filled_ids) < length(records)
        missing_ids = setdiff(keys(records), filled_ids)
        error("Missing header from annot.tbl: \"$(first(missing_ids))\"")
    end
    return records
end

# Please have a look at the files in annotation.tbl to understand this
# it partitions each individual record which is separated by equal signs.
function partition_chunks(io::IO)::Vector{Vector{String}}
    lines = collect(eachline(io))
    chunks = Vector{Vector{String}}()
    chunk = String[]
    for line in lines
        if startswith(line, "=========")
            push!(chunks, copy(chunk))
            empty!(chunk)
        else
            push!(chunk, line)
        end
    end
    push!(chunks, chunk)
end

function parse_annot_segment_data!(
    records::Dict{String, SegmentData},
    filled_ids::Set{String},
    chunk::Vector{String}
)
    @assert startswith(first(chunk), ">Feature ")
    header = first(chunk)[10:end]
    id = first(split(header, '|'))
    @assert haskey(records, id)
    push!(filled_ids, id)
    data = records[id]
    orfs = UnitRange{UInt16}[]
    reading_cds = false
    for line in @view chunk[2:end]
        isempty(strip(line)) && continue
        fields = split(line)
        if isnumeric(line[1]) && length(fields) > 2 && fields[3] == "CDS"
            reading_cds = true
        end
        if isnumeric(line[1]) && reading_cds
            push!(orfs, parse(UInt16, fields[1]) : parse(UInt16, fields[2]))
        end
        if reading_cds && fields[1] == "gene"
            variant = expect(parse(Protein, fields[2]), "Error on chunk $header")
            reading_cds = false
            push!(data.proteins, ProteinORF(variant, copy(orfs)))
            empty!(orfs)
        end
    end
    records
end

function add_extra_records!(
    segment_data::Dict{String, SegmentData},
    annotpath::String,
    fastapath::String,
    deduplicate::Bool
)
    # id => data from the FASTA file
    records = get_extra_records(fastapath, deduplicate)

    # Add info from annotated ORFS
    fill_orfs_extra_records!(records, annotpath)

    # Merge with the original segment_data dict
    @assert isdisjoint(keys(segment_data), keys(records))
    merge!(segment_data, records)
end

function get_extra_records(fastapath::String, deduplicate::Bool)
    records = Dict{String, SegmentData}()
    maybe_gzip(fastapath) do io
        reader = FASTA.Reader(io)
        record = FASTA.Record()
        while !eof(reader)
            read!(reader, record)
            data = parse_data_from_header(FASTA.header(record)::String, deduplicate)
            data.seq[] = some(FASTA.sequence(LongDNASeq, record))
            @assert !haskey(records, data.id) "Duplicate key $(data.id)"
            records[data.id] = data
        end
    end
    records
end

function parse_data_from_header(header::String, deduplicate::Bool)::SegmentData
    # Looks like this "EPI1843298|avian|NS|H1N1|2021|A/Denmark/1/2021|CLADE"
    # some of the fields may be empty
    fields = split(header, '|')
    @assert length(fields) == 7 header
    id, s_species, s_segment, s_subtype, s_year, name, clade = fields
    @assert !isempty(id)
    @assert !isempty(name)
    species = parse(Species, s_species)
    if species == other
        println(s_species)
    end
    @assert species != other
    return SegmentData(
        id,
        species,
        unwrap(parse(Segment, s_segment)),
        unwrap(parse(SubType, s_subtype)),
        deduplicate,
        isempty(clade) ? none(String) : some(String(clade)),
        parse(Int16, s_year),
        name,
        ProteinORF[],
        RefValue(none(LongDNASeq))
    )
end
=#


text:
module NedelecFESpacesTests
##
using Test
using Gridap
using Gridap.MapApply
using Gridap.CellValuesGallery
using LinearAlgebra

D = 3
order = 2

dom = fill(1,D)
part = Tuple(fill(5,D))
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0,0.0,1.0), partition=part)
# model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=part)

p = Polytope(fill(HEX_AXIS,D)...)
reffe = NedelecRefFE(p,Float64,order)

grid = Grid(model,D)
trian = Triangulation(grid)
graph = GridGraph(model)

labels = FaceLabels(model)
# tags = [5,6,7,8]
tags = Int[]
fun(x) = VectorValue(x[2],x[1],x[2]*x[1])
# fun(x) = VectorValue(1.0,1.0)
# tags = Int[]

fesp = ConformingFESpace(reffe,trian,graph,labels,tags)
dfesp = TrialFESpace(fesp,fun)

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

uh = interpolate(dfesp,fun)
uh.free_dofs
uh.diri_dofs

e = fun - uh
uh.free_dofs
uh.diri_dofs

el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
@test el2 < 1.0e-10

##

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Endpoint represents a single logical \"backend\" implementing a service.

    IoK8sApiDiscoveryV1beta1Endpoint(;
        addresses=nothing,
        conditions=nothing,
        hostname=nothing,
        targetRef=nothing,
        topology=nothing,
    )

    - addresses::Vector{String} : addresses of this endpoint. The contents of this field are interpreted according to the corresponding EndpointSlice addressType field. Consumers must handle different types of addresses in the context of their own capabilities. This must contain at least one address but no more than 100.
    - conditions::IoK8sApiDiscoveryV1beta1EndpointConditions : conditions contains information about the current status of the endpoint.
    - hostname::String : hostname of this endpoint. This field may be used by consumers of endpoints to distinguish endpoints from each other (e.g. in DNS names). Multiple endpoints which use the same hostname should be considered fungible (e.g. multiple A values in DNS). Must pass DNS Label (RFC 1123) validation.
    - targetRef::IoK8sApiCoreV1ObjectReference : targetRef is a reference to a Kubernetes object that represents this endpoint.
    - topology::Dict{String, String} : topology contains arbitrary topology information associated with the endpoint. These key/value pairs must conform with the label format. https://kubernetes.io/docs/concepts/overview/working-with-objects/labels Topology may include a maximum of 16 key/value pairs. This includes, but is not limited to the following well known keys: * kubernetes.io/hostname: the value indicates the hostname of the node   where the endpoint is located. This should match the corresponding   node label. * topology.kubernetes.io/zone: the value indicates the zone where the   endpoint is located. This should match the corresponding node label. * topology.kubernetes.io/region: the value indicates the region where the   endpoint is located. This should match the corresponding node label.
"""
mutable struct IoK8sApiDiscoveryV1beta1Endpoint <: SwaggerModel
    addresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: addresses
    conditions::Any # spec type: Union{ Nothing, IoK8sApiDiscoveryV1beta1EndpointConditions } # spec name: conditions
    hostname::Any # spec type: Union{ Nothing, String } # spec name: hostname
    targetRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: targetRef
    topology::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: topology

    function IoK8sApiDiscoveryV1beta1Endpoint(;addresses=nothing, conditions=nothing, hostname=nothing, targetRef=nothing, topology=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("addresses"), addresses)
        setfield!(o, Symbol("addresses"), addresses)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("hostname"), hostname)
        setfield!(o, Symbol("hostname"), hostname)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("targetRef"), targetRef)
        setfield!(o, Symbol("targetRef"), targetRef)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("topology"), topology)
        setfield!(o, Symbol("topology"), topology)
        o
    end
end # type IoK8sApiDiscoveryV1beta1Endpoint

const _property_map_IoK8sApiDiscoveryV1beta1Endpoint = Dict{Symbol,Symbol}(Symbol("addresses")=>Symbol("addresses"), Symbol("conditions")=>Symbol("conditions"), Symbol("hostname")=>Symbol("hostname"), Symbol("targetRef")=>Symbol("targetRef"), Symbol("topology")=>Symbol("topology"))
const _property_types_IoK8sApiDiscoveryV1beta1Endpoint = Dict{Symbol,String}(Symbol("addresses")=>"Vector{String}", Symbol("conditions")=>"IoK8sApiDiscoveryV1beta1EndpointConditions", Symbol("hostname")=>"String", Symbol("targetRef")=>"IoK8sApiCoreV1ObjectReference", Symbol("topology")=>"Dict{String, String}")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1Endpoint))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1Endpoint[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1Endpoint[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1Endpoint)
    (getproperty(o, Symbol("addresses")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }, name::Symbol, val)
end


text:
# This file is a part of BATBase.jl, licensed under the MIT License (MIT).

using BATBase

using Distributions, ValueShapes

struct MyMvDistributionDensity{D<:Distribution{Multivariate,Continuous}} <: AbstractDistributionDensity
    d::D
end

Base.convert(::Type{Distribution}, density::MyMvDistributionDensity) = density.d

BATBase.eval_logval_unchecked(density::MyMvDistributionDensity, v::AbstractVector{<:Real}) =
    logpdf(convert(Distribution, density), v)

ValueShapes.varshape(density::MyMvDistributionDensity) = varshape(density.d)

Distributions.sampler(density::MyMvDistributionDensity) = sampler(convert(Distribution, density))

# Note: Only for testing. Normally implemented generically in BAT.jl, not in BATBase.jl:
Base.convert(::Type{AbstractDistributionDensity}, d::Distribution{Multivariate,Continuous}) =
    MyMvDistributionDensity(d)

# Note: Only for testing. Normally implemented generically in BAT.jl, not in BATBase.jl:
BATBase.logdensityof(density::MyMvDistributionDensity, v::AbstractVector{<:Real}) =
    BATBase.eval_logval_unchecked(density, v)

# Note: Only for testing. Normally implemented generically in BAT.jl, not in BATBase.jl:
BATBase.logdensityof(density::MyMvDistributionDensity) = v -> logdensityof(density, v)


text:
#
# InfiniteMPS
#

# TODO: store the cell 1 as an MPS
# Implement `getcell(::InfiniteMPS, n::Integer) -> MPS`
mutable struct InfiniteMPS <: AbstractInfiniteMPS
  data::CelledVector{ITensor}
  llim::Int #RealInfinity
  rlim::Int #RealInfinity
  reverse::Bool
end

#
# InfiniteCanonicalMPS
#

# L and R are the orthogonalized MPS tensors.
# C are the center bond matrices (singular values of C
# are the singular values of the MPS at the specified
# cut)
struct InfiniteCanonicalMPS <: AbstractInfiniteMPS
  AL::InfiniteMPS
  C::InfiniteMPS
  AR::InfiniteMPS
end

# TODO: check if `isempty(ψ.AL)`, if so use `ψ.AR`
nsites(ψ::InfiniteCanonicalMPS) = nsites(ψ.AL)
isreversed(ψ::InfiniteCanonicalMPS) = isreversed(ψ.AL)
#ITensors.data(ψ::InfiniteCanonicalMPS) = data(ψ.AL)
ITensors.data(ψ::InfiniteCanonicalMPS) = ψ.AL.data

Base.copy(ψ::InfiniteCanonicalMPS) = InfiniteCanonicalMPS(copy(ψ.AL), copy(ψ.C), copy(ψ.AR))

function fmap(f, ψ::InfiniteCanonicalMPS)
  return InfiniteCanonicalMPS(f(ψ.AL), f(ψ.C), f(ψ.AR))
end

function ITensors.prime(ψ::InfiniteCanonicalMPS, args...; kwargs...)
  return fmap(x -> prime(x, args...; kwargs...), ψ)
end

function ITensors.prime(f::typeof(linkinds), ψ::InfiniteCanonicalMPS, args...; kwargs...)
  # Doesn't prime ψ.C
  #return fmap(x -> prime(f, x, args...; kwargs...), ψ)
  AL = prime(linkinds, ψ.AL, args...; kwargs...)
  C = prime(ψ.C, args...; kwargs...)
  AR = prime(linkinds, ψ.AR, args...; kwargs...)
  return InfiniteCanonicalMPS(AL, C, AR)
end

function ITensors.dag(ψ::InfiniteCanonicalMPS, args...; kwargs...)
  return fmap(x -> dag(x, args...; kwargs...), ψ)
end

ITensors.siteinds(f::typeof(only), ψ::InfiniteCanonicalMPS) = siteinds(f, ψ.AL)

# For now, only represents nearest neighbor interactions
# on a linear chain
struct InfiniteITensorSum
  data::CelledVector{ITensor}
end
InfiniteITensorSum(N::Int) = InfiniteITensorSum(Vector{ITensor}(undef, N))
InfiniteITensorSum(data::Vector{ITensor}) = InfiniteITensorSum(CelledVector(data))
function Base.getindex(l::InfiniteITensorSum, n1n2::Tuple{Int,Int})
  n1, n2 = n1n2
  @assert n2 == n1 + 1
  return l.data[n1]
end
nsites(h::InfiniteITensorSum) = length(l.data)

## HDF5 support for the InfiniteCanonicalMPS type

function HDF5.write(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ψ::InfiniteCanonicalMPS
)
  g = create_group(parent, name)
  attributes(g)["type"] = "InfiniteCanonicalMPS"
  attributes(g)["version"] = 1
  N = nsites(ψ)
  write(g, "length", N)
  for n in 1:N
    write(g, "AL[$(n)]", ψ.AL[n])
    write(g, "AR[$(n)]", ψ.AR[n])
    write(g, "C[$(n)]", ψ.C[n])
  end
end

function HDF5.read(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{InfiniteCanonicalMPS}
)
  g = open_group(parent, name)
  if read(attributes(g)["type"]) != "InfiniteCanonicalMPS"
    error("HDF5 group or file does not contain InfiniteCanonicalMPS data")
  end
  N = read(g, "length")
  AL = InfiniteMPS([read(g, "AL[$(i)]", ITensor) for i in 1:N])
  AR = InfiniteMPS([read(g, "AR[$(i)]", ITensor) for i in 1:N])
  C = InfiniteMPS([read(g, "C[$(i)]", ITensor) for i in 1:N])
  return InfiniteCanonicalMPS(AL, C, AR)
end


text:
using Base.Test
using Clang.wrap_c
include("strip_line_numbers.jl")

wc = wrap_c.init(;  headers = [joinpath(dirname(@__FILE__), "cxx/cbasic.h")])
hdr_tunits = wrap_c.parse_c_headers(wc)

# test func1
tunit = hdr_tunits[joinpath(dirname(@__FILE__), "cxx/cbasic.h")]
func1 = search(tunit, "func1")[1]
buf = Any[]

wrap_c.wrap(wc, buf, func1, "test")

exc = :(
    function func1(a::Cint,b::Cdouble,c,d)
        ccall((:func1,test),Cint,(Cint,Cdouble,Ptr{Cdouble},Ptr{Void}),a,b,c,d)
    end)

strip_line_numbers!(exc)

@test buf[1] == exc


text:
const WFLOAT = Sys.WORD_SIZE == 64 ? Float64 : Float32

#=========
Abstraction
==========#

abstract type RobustVariance end
abstract type HAC{G} <: RobustVariance end
abstract type HC <: RobustVariance end
abstract type CRHC{V} <: RobustVariance end

#=========
HAC Types
=========#

abstract type BandwidthType{G} end
abstract type OptimalBandwidth end

struct NeweyWest<:OptimalBandwidth end
struct Andrews<:OptimalBandwidth end

struct Fixed<:BandwidthType{G where G} end
struct Optimal{G<:OptimalBandwidth}<:BandwidthType{G where G<:OptimalBandwidth} end

struct Prewhitened end
struct Unwhitened end

struct TruncatedKernel{G<:BandwidthType, F}<:HAC{G}
  bwtype::G
  bw::Vector{F}
  weights::Vector{F}
  prewhiten::Bool
end

struct BartlettKernel{G<:BandwidthType, F}<:HAC{G}
    bwtype::G
    bw::Vector{F}
    weights::Vector{F}
    prewhiten::Bool
end

struct ParzenKernel{G<:BandwidthType, F}<:HAC{G}
    bwtype::G
    bw::Vector{F}
    weights::Vector{F}
    prewhiten::Bool
end

struct TukeyHanningKernel{G<:BandwidthType, F}<:HAC{G}
    bwtype::G
    bw::Vector{F}
    weights::Vector{F}
    prewhiten::Bool
end

struct QuadraticSpectralKernel{G<:BandwidthType, F}<:HAC{G}
    bwtype::G
    bw::Vector{F}
    weights::Vector{F}
    prewhiten::Bool
end

const TRK=TruncatedKernel
const BTK=BartlettKernel
const PRK=ParzenKernel
const THK=TukeyHanningKernel
const QSK=QuadraticSpectralKernel

#=========
HC Types
=========#
struct HC0  <: HC end
struct HC1  <: HC end
struct HC2  <: HC end
struct HC3  <: HC end
struct HC4  <: HC end
struct HC4m <: HC end
struct HC5  <: HC end

#const CLVector{T<:Integer} = DenseArray{T,1}

mutable struct CRHC0{V<:AbstractVector}  <: CRHC{V}
    cl::V
end

mutable struct CRHC1{V<:AbstractVector}  <: CRHC{V}
    cl::V
end

mutable struct CRHC2{V<:AbstractVector}  <: CRHC{V}
    cl::V
end

mutable struct CRHC3{V<:AbstractVector}  <: CRHC{V}
    cl::V
end

struct CovarianceMatrix{T2<:Factorization, T3<:CovarianceMatrices.RobustVariance, F1, T1<:AbstractMatrix{F1}} <: AbstractMatrix{F1}
    F::T2       ## Factorization
    K::T3       ## RobustVariance, e.g. HC0()
    V::T1       ## The Variance Covariance
end

#=======
Caches
=======#

abstract type AbstractCache end

struct HACCache{TYPE, F<:AbstractMatrix, V<:AbstractVector} <: AbstractCache
    prew::TYPE
    X_demean::F ## Should call this q n_origin x p
    YY::F
    XX::F
    Y_lagged::F
    X_lagged::F
    μ::F        ## p x 1
    Q::F        ## p x p
    V::F        ## p x p
    D::F       ## p x p
    U::V
    ρ::V
    σ⁴::V
    u::F
end


text:
export generate_mesh

function generate_mesh(;dim::Int64 = 2,
                        order::Int64 = 1,
                        faces::String = "Triangle",
                        opt_iters::Int64 = 0,
                        force_quads::Bool = false)
    @info "Generating a mesh of dimension $dim with $faces faces of order $order"
    gmsh.option.set_number("Mesh.SecondOrderIncomplete", 1)
    if dim != 2
        error("Only 2D meshes are currently supported.")
    end
    if faces == "Triangle"
        # Delaunay (5) handles large element size gradients better
        gmsh.option.set_number("Mesh.Algorithm", 5)
        if order == 1
            gmsh.model.mesh.generate(2)
            for _ in 1:opt_iters
                gmsh.model.mesh.optimize("Laplace2D")
            end
        elseif order == 2
            gmsh.option.set_number("Mesh.HighOrderOptimize", 2) # elastic + optimization
            gmsh.model.mesh.generate(2)
            gmsh.model.mesh.set_order(2)
            for _ in 1:opt_iters
                gmsh.model.mesh.optimize("HighOrderElastic")
                gmsh.model.mesh.optimize("Relocate2D")
                gmsh.model.mesh.optimize("HighOrderElastic")
            end
        else
            error("Mesh order must be 1 or 2.")
        end
    elseif faces == "Quadrilateral"
        gmsh.option.set_number("Mesh.RecombineAll", 1)
        gmsh.option.set_number("Mesh.Algorithm", 8) # Frontal-Delaunay for quads.
        if force_quads
            gmsh.option.set_number("Mesh.RecombinationAlgorithm", 2) # simple full-quad
            gmsh.option.set_number("Mesh.SubdivisionAlgorithm", 1) # All quads
        end
        if order == 1
            gmsh.model.mesh.generate(2)
            for _ in 1:opt_iters
                gmsh.model.mesh.optimize("Laplace2D")
                gmsh.model.mesh.optimize("Relocate2D")
                gmsh.model.mesh.optimize("Laplace2D")
            end
        elseif order == 2
            gmsh.option.set_number("Mesh.HighOrderOptimize", 2)
            gmsh.model.mesh.generate(2)
            gmsh.model.mesh.set_order(2)
            for _ in 1:opt_iters
                gmsh.model.mesh.optimize("HighOrderElastic")
                gmsh.model.mesh.optimize("Relocate2D")
                gmsh.model.mesh.optimize("HighOrderElastic")
            end
        else
            error("Mesh order must be 1 or 2.")
        end
    else
        error("Could not identify face type.")
    end

    return nothing
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using MatrixCompletion
import LinearAlgebra

function unit_test_lpspace(p)
    tc = LpSpace(p)
    test_vec = rand(100)
    @test tc.p==p && tc.norm(test_vec) == LinearAlgebra.norm(test_vec,tc.p)
end

@testset "$(format("Concepts: LpSpace[Construction]"))" begin
    [unit_test_lpspace(i) for i in 1:10]
end

@testset "$(format("Concepts: Type Convertsion[Symbol->Exponential Family]"))" begin
    @test typeof(convert(ExponentialFamily,:Poisson)) == typeof(Poisson())
end

@testset "$(format("Concepts: Comparator[construction]"))" begin
    @test typeof(Comparator{Int64}()) == Comparator{Int64}
    @test typeof(Comparator(MGF(:Gaussian)))  == Comparator{MGF{Gaussian}}
    @test typeof(Comparator(MGF)) == Comparator{MGF}
    @test typeof(Comparator(:MGF)) == Comparator{MGF}
    let
        tc = Comparator{MGF}(MGF(),eval_at = [1,2,3])
        @test tc.field[:eval_at] == [1,2,3]
    end
end


text:
"""
Solve MVDA problem with cyclic coordinate descent.
"""
struct CyclicVDA <: AbstractMMAlg end

# Initialize data structures.
function __mm_init__(::CyclicVDA, problem, ::Nothing)
    return nothing
end

# Check for data structure allocations; otherwise initialize.
function __mm_init__(::CyclicVDA, problem, extras)
    if :apply_projection in keys(extras) && :buffer in keys(extras) # TODO
        return extras
    else
        __mm_init__(CyclicVDA(), problem, nothing)
    end
end

# Apply one update in regularized version.
function __mm_iterate__(::CyclicVDA, problem, ϵ, δ, λ₁, λ₂, extras)
    @unpack X, intercept, coeff, res = problem
    MAX_NEWTON_STEPS = 100
    MAX_BACKTRACK_STEPS = 5

    B = coeff.all
    R = res.main.all

    for j in axes(B, 1) # loop over features
        for k in axes(B, 2) # loop over dimensions in vertex space
            objective, d1objective, d2objective = fetch_objective(problem, j, k, ϵ, δ, λ₁, λ₂)
            new_objective = objective

            if abs(d1objective) ≤ 0 continue end

            for _ in 1:MAX_NEWTON_STEPS # minimize via Newton's method
                β = B[j,k]
                ∇ = -d1objective / d2objective

                # Evaluate the new objective. Use step-halving strategy for line search.
                for _ in 1:MAX_BACKTRACK_STEPS
                    A = β + ∇ - B[j,k]
                    B[j,k] = β + ∇

                    # Update residuals, then reevaluate objective to check for descent property.
                    for i in axes(X, 1)
                        R[i,k] = R[i,k] - A * X[i,j]
                    end
                    new_objective, d1objective, d2objective = fetch_objective(problem, j, k, ϵ, δ, λ₁, λ₂)
                    descent = new_objective ≤ objective

                    if descent break end
                    ∇ = 1//2 * ∇
                end # step-halving
                converged = abs(objective - new_objective) ≤ 1e-4

                if converged break end
            end # Newton

            # Reset parameters close to 0 to 0.
            β = B[j,k]
            is_near_zero = abs(β) ≤ 1e-8
            B[j,k] = ifelse(is_near_zero, zero(β), β)
        end # dimensions
    end # features

    return nothing
end

function fetch_objective(problem, j, k, ϵ, δ, λ₁, λ₂) # feature j, dimension k
    @unpack X, coeff, res = problem
    cases, features, _ = probdims(problem)

    # Constants.
    C1 = ϵ - δ
    C2 = ϵ + δ
    C3 = 1 / (16 * δ^3)
    C4 = 1 / (4 * δ^3)
    C5 = 3 / (4 * δ^3)
    μ₁ = cases * λ₁
    μ₂ = cases * λ₂

    # Initialize outputs.
    objective = 0.0
    d1objective = 0.0
    d2objective = 0.0

    R = res.main.all
    B = coeff.all

    # Loop over all cases i.
    for i in axes(R, 1)
        rᵢ = view(R, i, :)
        norm_rᵢ = norm(rᵢ)
        norm_rᵢ² = norm_rᵢ^2

        # Skip the loss associated with case i if its predicted value is within
        # a distance ϵ-δ of its associated vertex.
        if norm_rᵢ < C1 continue end

        x = X[i,j]
        rx = rᵢ[k] * x

        # Compute the loss for case i if it lies beyond a distance ϵ+δ of its associated vertex.
        if norm_rᵢ > C2
            objective = objective + norm_rᵢ - ϵ
            s = 1 / norm_rᵢ
            d1objective = d1objective - s*rx
            d2objective = d2objective + s*(x*x - rx*rx/norm_rᵢ²)
        else
            s = norm_rᵢ - C1
            t = C2 - norm_rᵢ
            d²p = C5 * s * t
            s² = s*s
            t = t + δ
            d¹p = C4 * s² * t
            t = t + δ
            s³ = s²*s
            p = C3 * s³ * t
            objective = objective + p
            s = d¹p / norm_rᵢ
            d1objective = d1objective - s*rx
            t = d²p / norm_rᵢ²
            rx² = rx * rx
            d2objective = d2objective + t*rx² + s*(x*x - rx²/norm_rᵢ²)
        end
    end

    # Compute the penalty for the current parameter.
    update_penalty = problem.intercept ? (j != features + 1) : true
    if update_penalty
        βⱼ = view(B, j, :)
        normβⱼ = norm(βⱼ, 2)
        normβⱼ² = normβⱼ^2
        objective = objective + μ₂*normβⱼ + μ₁* norm(βⱼ, 1)

        # If the parameter is parked at 0, then reset the first derivative of the objective
        # so that it either remains at 0 or moves in the correct direction.
        if normβⱼ ≤ 0
            if abs(d1objective) ≤ μ₂ + μ₁
                d1objective = zero(d1objective)
            elseif d1objective + μ₂ + μ₁ ≤ 0
                d1objective = d1objective + μ₂ + μ₁
            else
                d1objective = d1objective - μ₂ - μ₁
            end
        else # Deal with the Euclidean penalty away from the origin.
            s = B[j,k]
            t = s / normβⱼ
            d1objective = d1objective + μ₂*t
            d2objective = d2objective + μ₂*(1-t*t)/normβⱼ

            if abs(d1objective) ≤ μ₁
                d1objective = zero(d1objective)
            elseif d1objective + μ₁ ≤ 0
                d1objective = d1objective + μ₁
            else
                d1objective = d1objective - μ₁
            end
        end
    end

    return objective, d1objective, d2objective
end


text:
@sumProductRule(:node_type     => Wishart,
                :outbound_type => Message{Wishart},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPWishartOutNPP)

@naiveVariationalRule(:node_type     => Wishart,
                 	  :outbound_type => Message{Wishart},
                 	  :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                 	  :name          => VBWishartOut)


text:
mutable struct GoToTargetDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Target1, Target2}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    target1_pos::CartesianIndex{2}
    target2_pos::CartesianIndex{2}
    target::Union{Target1, Target2}
    terminal_reward::T
    terminal_penalty::T
    done::Bool
end

@generate_getters(GoToTargetDirected)
@generate_setters(GoToTargetDirected)

mutable struct GoToTargetUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Target1, Target2}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    target1_pos::CartesianIndex{2}
    target2_pos::CartesianIndex{2}
    target::Union{Target1, Target2}
    terminal_reward::T
    terminal_penalty::T
    done::Bool
end

@generate_getters(GoToTargetUndirected)
@generate_setters(GoToTargetUndirected)

#####
# Directed
#####

function GoToTargetDirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, TARGET1, TARGET2)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    target1_pos = CartesianIndex(height - 1, width - 1)
    world[TARGET1, target1_pos] = true
    target2_pos = CartesianIndex(height - 1, width - 2)
    world[TARGET2, target2_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    agent_dir = RIGHT
    reward = zero(T)
    target = rand(rng, (TARGET1, TARGET2))
    terminal_reward = one(T)
    terminal_penalty = -one(T)
    done = false

    env = GoToTargetDirected(world, agent_pos, agent_dir, reward, rng, target1_pos, target2_pos, target, terminal_reward, terminal_penalty, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::GoToTargetDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const GO_TO_TARGET_DIRECTED_LAYERS = SA.SVector(2, 3, 4)
RLBase.state(env::GoToTargetDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), GO_TO_TARGET_DIRECTED_LAYERS)

RLBase.state_space(env::GoToTargetDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::GoToTargetDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::GoToTargetDirected, ::RLBase.DefaultPlayer) = DIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::GoToTargetDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::GoToTargetDirected) = get_done(env)

function RLBase.reset!(env::GoToTargetDirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[TARGET1, get_target1_pos(env)] = false
    world[TARGET2, get_target2_pos(env)] = false

    new_target1_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_target1_pos!(env, new_target1_pos)
    world[TARGET1, new_target1_pos] = true

    new_target2_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_target2_pos!(env, new_target2_pos)
    world[TARGET2, new_target2_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_target!(env, rand(rng, (TARGET1, TARGET2)))
    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::GoToTargetDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)
    target = get_target(env)
    agent_pos = get_agent_pos(env)

    if world[TARGET1, agent_pos] || world[TARGET2, agent_pos]
        set_done!(env, true)
        if world[target, agent_pos]
            set_reward!(env, get_terminal_reward(env))
        else
            set_reward!(env, get_terminal_penalty(env))
        end
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::GoToTargetDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)
    target = get_target(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    agent_pos = get_agent_pos(env)
    if world[TARGET1, agent_pos] || world[TARGET2, agent_pos]
        set_done!(env, true)
        if world[target, agent_pos]
            set_reward!(env, get_terminal_reward(env))
        else
            set_reward!(env, get_terminal_penalty(env))
        end
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Undirected
#####

function GoToTargetUndirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, TARGET1, TARGET2)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    target1_pos = CartesianIndex(height - 1, width - 1)
    world[TARGET1, target1_pos] = true
    target2_pos = CartesianIndex(height - 1, width - 2)
    world[TARGET2, target2_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    reward = zero(T)
    target = rand(rng, (TARGET1, TARGET2))
    terminal_reward = one(T)
    terminal_penalty = -one(T)
    done = false

    env = GoToTargetUndirected(world, agent_pos, reward, rng, target1_pos, target2_pos, target, terminal_reward, terminal_penalty, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::GoToTargetUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const GO_TO_TARGET_UNDIRECTED_LAYERS = SA.SVector(2, 3, 4)
RLBase.state(env::GoToTargetUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), GO_TO_TARGET_UNDIRECTED_LAYERS)

RLBase.state_space(env::GoToTargetUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::GoToTargetUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::GoToTargetUndirected, player::RLBase.DefaultPlayer) = UNDIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::GoToTargetUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::GoToTargetUndirected) = get_done(env)

function RLBase.reset!(env::GoToTargetUndirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[TARGET1, get_target1_pos(env)] = false
    world[TARGET2, get_target2_pos(env)] = false

    new_target1_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_target1_pos!(env, new_target1_pos)
    world[TARGET1, new_target1_pos] = true

    new_target2_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_target2_pos!(env, new_target2_pos)
    world[TARGET2, new_target2_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_target!(env, rand(rng, (TARGET1, TARGET2)))
    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::GoToTargetUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)
    target = get_target(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    agent_pos = get_agent_pos(env)
    if world[TARGET1, agent_pos] || world[TARGET2, agent_pos]
        set_done!(env, true)
        if world[target, agent_pos]
            set_reward!(env, get_terminal_reward(env))
        else
            set_reward!(env, get_terminal_penalty(env))
        end
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end


text:
# ---
# title: 1202. Smallest String With Swaps
# id: problem1202
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Union Find
# link: <https://leetcode.com/problems/smallest-string-with-swaps/description/>
# hidden: true
# ---
#
# You are given a string `s`, and an array of pairs of indices in the string
# `pairs` where `pairs[i] = [a, b]` indicates 2 indices(0-indexed) of the
# string.
#
# You can swap the characters at any pair of indices in the given `pairs`  **any
# number of times**.
#
# Return the lexicographically smallest string that `s` can be changed to after
# using the swaps.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "dcab", pairs = [[0,3],[1,2]]
#     Output: "bacd"
#     **Explaination:**
#     Swap s[0] and s[3], s = "bcad"
#     Swap s[1] and s[2], s = "bacd"
#
#
# **Example 2:**
#
#
#
#     Input: s = "dcab", pairs = [[0,3],[1,2],[0,2]]
#     Output: "abcd"
#     **Explaination:**
#     Swap s[0] and s[3], s = "bcad"
#     Swap s[0] and s[2], s = "acbd"
#     Swap s[1] and s[2], s = "abcd"
#
# **Example 3:**
#
#
#
#     Input: s = "cba", pairs = [[0,1],[1,2]]
#     Output: "abc"
#     **Explaination:**
#     Swap s[0] and s[1], s = "bca"
#     Swap s[1] and s[2], s = "bac"
#     Swap s[0] and s[1], s = "abc"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `0 <= pairs.length <= 10^5`
#   * `0 <= pairs[i][0], pairs[i][1] < s.length`
#   * `s` only contains lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
###############################################################################
#
# Temperature effects on surface tension and viscosity
#
###############################################################################
"""
    temperature_effects!(
                hs::LeafHydraulics{FT}
    ) where {FT<:AbstractFloat}
    temperature_effects!(
                hs::AbstractHydraulicOrgan{FT}
    ) where {FT<:AbstractFloat}
    temperature_effects!(
                hs::AbstractHydraulicOrgan{FT},
                T::FT
    ) where {FT<:AbstractFloat}
    temperature_effects!(
                tree::TreeSimple{FT}
    ) where {FT<:AbstractFloat}

Update temperature effetcs on hydralic system, given
- `hs` [`AbstractHydraulicOrgan`](@ref) type struct
- `T` Given temperature
- `tree` [`TreeSimple`](@ref) type struct
"""
function temperature_effects!(hs::LeafHydraulics)
    if hs.T_sap != hs.T_old
        hs.f_st  = relative_surface_tension(hs.T_sap);
        hs.f_vis = relative_viscosity(hs.T_sap);
        hs.p_crt = xylem_p_crit(hs.vc);
        hs.T_old = hs.T_sap;
    end

    return nothing
end

function temperature_effects!(hs::Union{RootHydraulics,StemHydraulics})
    if hs.T_sap != hs.T_old
        hs.f_st  = relative_surface_tension(hs.T_sap);
        hs.f_vis = relative_viscosity(hs.T_sap);
        hs.T_old = hs.T_sap;
    end

    return nothing
end

function temperature_effects!(
            hs::Union{LeafHydraulics{FT},
                      RootHydraulics{FT},
                      StemHydraulics{FT}},
            T::FT
) where {FT<:AbstractFloat}
    if T != hs.T_sap
        hs.T_sap = T;
        temperature_effects!(hs);
    end

    return nothing
end

function temperature_effects!(tree::GrassLikeOrganism)
    for root in tree.roots
        temperature_effects!(root);
    end
    for leaf in tree.leaves
        temperature_effects!(leaf);
    end

    return nothing
end

function temperature_effects!(tree::PalmLikeOrganism)
    for root in tree.roots
        temperature_effects!(root);
    end
    temperature_effects!(tree.trunk);
    for leaf in tree.leaves
        temperature_effects!(leaf);
    end

    return nothing
end

function temperature_effects!(tree::TreeLikeOrganism)
    for root in tree.roots
        temperature_effects!(root);
    end
    temperature_effects!(tree.trunk);
    for stem in tree.branch
        temperature_effects!(stem);
    end
    for leaf in tree.leaves
        temperature_effects!(leaf);
    end

    return nothing
end

function temperature_effects!(tree::TreeSimple)
    temperature_effects!(tree.root);
    temperature_effects!(tree.stem);
    temperature_effects!(tree.leaf);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
function eye(m::Int64)
    Matrix(1.0I, m, m)
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics