text:
using ..CUDASupport
using ..TopOpt: @init_cuda
@init_cuda()

using ..GPUUtils
import ..TopOpt: whichdevice

whichdevice(s::AbstractDisplacementSolver) = whichdevice(s.u)

for T in (:DirectDisplacementSolver, :PCGDisplacementSolver)
    @eval @inline CuArrays.cu(p::$T) = error("$T does not support the GPU.")
end

function update_f!(f::CuVector{T}, values, prescribed_dofs, applyzero, dof_cells, cell_dofs, black, white, Kes, xmin, penalty, vars, varind, M) where {T}
    args = (f, values, prescribed_dofs, applyzero, dof_cells.offsets, dof_cells.values, cell_dofs, black, white, Kes, xmin, penalty, vars, varind, M)
    callkernel(dev, bc_kernel, args)
    CUDAdrv.synchronize(ctx)
    return
end

function bc_kernel(f::AbstractVector{T}, values, prescribed_dofs, applyzero, dof_cells_offsets, dof_cells_values, cell_dofs, black, white, Kes, xmin, penalty, vars, varind, M) where {T}

    ind = @thread_global_index()
    offset = @total_threads()
    while ind <= length(values)
        d = prescribed_dofs[ind]
        v = values[ind]
        m = size(Kes[ind], 1)

        r = dof_cells_offsets[d] : dof_cells_offsets[d+1]-1
        if !applyzero && v != 0
            for idx in r
                (i, j) = dof_cells_values[idx]
                if PENALTY_BEFORE_INTERPOLATION
                    px = ifelse(black[i], one(T),
                    ifelse(white[i], xmin,
                    density(penalty(vars[varind[i]]), xmin)))
                else
                    px = ifelse(black[i], one(T),
                    ifelse(white[i], xmin,
                    penalty(density(vars[varind[i]], xmin))))
                end
                if eltype(Kes) <: Symmetric
                    Ke = Kes[i].data
                else
                    Ke = Kes[i]
                end
                for row in 1:m
                    f[cell_dofs[row,i]] -= px * v * Ke[row,j]
                end
            end
        end
        f[d] = M*v
        ind += offset
    end
    return
end

@define_cu(StaticMatrixFreeDisplacementSolver, :f, :problem, :vars, :cg_statevars, :elementinfo, :penalty, :prev_penalty, :u, :fixed_dofs, :free_dofs, :xes, :lhs, :rhs)

whichdevice(m::MatrixFreeOperator) = whichdevice(m.vars)

function mul!(y::TV, A::MatrixFreeOperator, x::TV) where {TV <: CuArrays.CuVector}
    T = eltype(y)
    nels = length(A.elementinfo.Kes)
    ndofs = length(A.elementinfo.fixedload)
    meandiag = A.meandiag

    @unpack Kes, metadata, black, white, varind = A.elementinfo
    @unpack cell_dofs, dof_cells = metadata
    @unpack penalty, xmin, vars, fixed_dofs, free_dofs, xes = A

    args1 = (xes, x, black, white, vars, varind, cell_dofs, Kes, xmin, penalty, nels)
    callkernel(dev, mul_kernel1, args1)
    CUDAdrv.synchronize(ctx)

    args2 = (y, x, dof_cells.offsets, dof_cells.values, xes, fixed_dofs, free_dofs, meandiag)
    callkernel(dev, mul_kernel2, args2)
    CUDAdrv.synchronize(ctx)

    return y
end

# CUDA kernels
function mul_kernel1(xes::AbstractVector{TV}, x, black, white, vars, varind, cell_dofs, Kes, xmin, penalty, nels) where {N, T, TV<:SVector{N, T}}
    i = @thread_global_index()
    offset = @total_threads()
    while i <= nels
        if PENALTY_BEFORE_INTERPOLATION
            px = ifelse(black[i], one(T),
            ifelse(white[i], xmin,
                    density(penalty(vars[varind[i]]), xmin)
                )
            )
        else
            px = ifelse(black[i], one(T),
            ifelse(white[i], xmin,
                    penalty(density(vars[varind[i]], xmin))
                )
            )
        end
        xe = xes[i]
        for j in 1:N
            xe = @set xe[j] = x[cell_dofs[j, i]]
        end
        if eltype(Kes) <: Symmetric
            xe = SVector{1, T}((px,)) .* (bcmatrix(Kes[i]).data * xe)
        else
            xe = SVector{1, T}((px,)) .* (bcmatrix(Kes[i]) * xe)
        end

        xes[i] = xe
        i += offset
    end

    return
end

function mul_kernel2(y, x, dof_cells_offsets, dof_cells_values, xes, fixed_dofs, free_dofs, meandiag)
    T = eltype(y)
    offset = @total_threads()
    n_fixeddofs = length(fixed_dofs)
    ndofs = length(y)

    i = @thread_global_index()
    while i <= n_fixeddofs
        dof = fixed_dofs[i]
        y[dof] = meandiag * x[dof]
        i += offset
    end

    i = @thread_global_index()
    while n_fixeddofs < i <= ndofs
        dof = free_dofs[i - n_fixeddofs]
        yi = zero(T)
        r = dof_cells_offsets[dof] : dof_cells_offsets[dof+1]-1
        for ind in r
            k, m = dof_cells_values[ind]
            yi += xes[k][m]
        end
        y[dof] = yi
        i += offset
    end
    return
end

whichdevice(r::LinearElasticityResult) = whichdevice(r.u)

function fill_vars!(vars::GPUArray, topology; round)
    if round
        vars .= round.(typeof(vars)(topology))
    else
        copyto!(vars, topology)
    end
    return vars
end


text:
mutable struct image_sync_t <: LCMType
    utime::Int64
end

@lcmtypesetup(image_sync_t)

ext: jl


text:
using PrettyTables

struct PolicyTestResults
    results::Vector{PolicyTestResult}
end

function PolicyTestResults(x::MuZeroConfig)
    getpolicytests()
    d = policy_testset(x.game_type)
    tests = []
    for (idx, kv) in enumerate(d)
        game, actions = kv
        push!(tests, PolicyTest(idx, game, actions))
    end
    PolicyTestResults(tests)
end

struct PolicyTestResultsSVG
    tests::PolicyTestResults
end

using Random

svg(x::PolicyTestResults) = PolicyTestResultsSVG(x)
Base.rand(x::PolicyTestResults) = rand(x.results)
Base.getindex(x::PolicyTestResults, idx::Integer) = x.results[idx]
Base.getindex(x::PolicyTestResults, idx::Vector{<:Integer}) = x.results[idx]
Base.length(x::PolicyTestResults) = length(x.results)
Base.lastindex(x::PolicyTestResults) = length(x)
Base.iterate(x::PolicyTestResults, state=1) = state > length(x) ? nothing : (x[state], state+1)

function df(x::PolicyTestResults; range=1:length(x), limit=3)
    test_num = "Test #"
    test_pf = "P/F"
    test_result = "Pred (π*)"
    optimal_action = "Actual (π*)"
    d = Dict(
        test_num=>[],
        test_pf=>[],
        test_result=>[],
        optimal_action=>[],
    )
    for result in x.results[range]
        idx = result.test.idx
        optimal = result.test.actions
        pred = result.prediction
        _, greedy_action = findmax(pred)
        pass = in(greedy_action, optimal)
        from = 1
        to = min(limit, length(pred))
        actions = reverse(sortperm(pred))[from:to]
        prs = pred[actions]
        action_prs = map(x->begin
            action, pr = x
            pr_s = @sprintf "%.2f" pr
            "$action ($(pr_s[2:end]))"
        end, zip(actions, prs))
        result_policy = join(action_prs, ", ")
        push!(d[test_num], idx)
        push!(d[test_pf], pass ? "Pass" : "Fail")
        push!(d[test_result], result_policy)
        push!(d[optimal_action], join(optimal, ","))
    end
    DataFrame(
        test_num=>d[test_num],
        test_pf=>d[test_pf],
        test_result=>d[test_result],
        optimal_action=>d[optimal_action],
    )
end

function Base.show(io::IO, m::MIME"text/plain", x::PolicyTestResults)
    len = length(x.results)
    println()
    pretty_results(
        io,
        df(x);
        title="Policy Test Results (all $len)"
    )
end

function pretty_results(io::IO, df::DataFrame; title)
    pf_hl = Highlighter((data, i, j)->(j==2) && data[i,j] == "Pass", crayon"green");
    pretty_table(
        io,
        df;
        show_row_number=false,
        nosubheader=true,
        title=title,
        highlighters=(pf_hl),
    )
end

function Base.getindex(x::PolicyTestResults, r::UnitRange{<:Integer})
    println()
    pretty_results(
        stdout,
        df(x; range=r),
        title="Policy Test Results ($(r))"
    )
end


text:
function get_geometry_free_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getGeometryFreeCombination", GeometryFreeCombination, (SatelliteSystem,), arg0)
end

function get_graphic_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getGRAPHICCombination", GRAPHICCombination, (SatelliteSystem,), arg0)
end

function get_ionosphere_free_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getIonosphereFreeCombination", IonosphereFreeCombination, (SatelliteSystem,), arg0)
end

function get_melbourne_wubbena_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getMelbourneWubbenaCombination", MelbourneWubbenaCombination, (SatelliteSystem,), arg0)
end

function get_narrow_lane_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getNarrowLaneCombination", NarrowLaneCombination, (SatelliteSystem,), arg0)
end

function get_phase_minus_code_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getPhaseMinusCodeCombination", PhaseMinusCodeCombination, (SatelliteSystem,), arg0)
end

function get_wide_lane_combination(::Type{MeasurementCombinationFactory}, arg0::SatelliteSystem)
    return jcall(MeasurementCombinationFactory, "getWideLaneCombination", WideLaneCombination, (SatelliteSystem,), arg0)
end


text:
module Bobby

    export Bitboard
    using Printf
    using Crayons

    include("structs.jl")
    include("constants.jl")
    include("cvt.jl")
    include("bitboard.jl")
    include("king.jl")
    include("knight.jl")
    include("rook.jl")
    include("bishop.jl")
    include("magic.jl")
    include("pawn.jl")
    include("move.jl")
    include("perft.jl")
end


text:
"""
    padding(n::Integer)

Determines the total number of bytes needed to store `n` bytes with padding.
Note that the Arrow standard requires buffers to be aligned to 8-byte boundaries.
"""
padding(n::Integer) = ((n + ALIGNMENT - 1) ÷ ALIGNMENT) * ALIGNMENT

paddinglength(n::Integer) = padding(n) - n

function writezeros(io::IO, n::Integer)
    s = 0
    for i ∈ 1:n
        s += Base.write(io, 0x00)
    end
    s
end

# efficient writing of arrays
function writearray(io::IO, ::Type{T}, col) where {T}
    if col isa Vector{T}
        n = Base.write(io, col)
    elseif isbitstype(T) && (col isa Vector{Union{T, Missing}} || col isa SentinelVector{T, T, Missing, Vector{T}})
        # need to write the non-selector bytes of isbits Union Arrays
        n = Base.unsafe_write(io, pointer(col), sizeof(T) * length(col))
    else
        n = 0
        for x in col
            n += Base.write(io, coalesce(x, default(T)))
        end
    end
    return n
end

"""
    getbit

This deliberately elides bounds checking.
"""
getbit(v::UInt8, n::Integer) = Bool((v & 0x02^(n - 1)) >> (n - 1))

"""
    setbit

This also deliberately elides bounds checking.
"""
function setbit(v::UInt8, b::Bool, n::Integer)
    if b
        v | 0x02^(n - 1)
    else
        v & (0xff ⊻ 0x02^(n - 1))
    end
end

"""
    bitpackedbytes(n[, pad=true])

Determines the number of bytes used by `n` bits, optionally with padding.
"""
function bitpackedbytes(n::Integer, pad::Bool=true)
    a, b = divrem(n, 8)
    ℓ = a + (b > 0)
    pad && (ℓ += paddinglength(ℓ))
    return ℓ
end

# count # of missing elements in an iterable
nullcount(col) = count(ismissing, col)

# like startswith/endswith for strings, but on byte buffers
function _startswith(a::AbstractVector{UInt8}, pos::Integer, b::AbstractVector{UInt8})
    for i = 1:length(b)
        @inbounds check = a[pos + i - 1] == b[i]
        check || return false
    end
    return true
end

function _endswith(a::AbstractVector{UInt8}, endpos::Integer, b::AbstractVector{UInt8})
    aoff = endpos - length(b) + 1
    for i = 1:length(b)
        @inbounds check = a[aoff] == b[i]
        check || return false
        aoff += 1
    end
    return true
end

# read a single element from a byte vector
# copied from read(::IOBuffer, T) in Base
function readbuffer(t::AbstractVector{UInt8}, pos::Integer, ::Type{T}) where {T}
    GC.@preserve t begin
        ptr::Ptr{T} = pointer(t, pos)
        x = unsafe_load(ptr)
    end
end

flatten(x) = Iterators.flatten(x)

# we need to treat missing specially w/ length of flattened array
_length(x) = length(x)
_length(g::Base.Generator) = _length(g.iter)
_length(x::Iterators.Flatten) = sum(i -> i === missing ? 1 : _length(i), x)

# argh me mateys, don't nobody else go pirating this method
# this here be me own booty!
if !applicable(iterate, missing)
Base.iterate(::Missing, st=1) = st === nothing ? nothing : (missing, nothing)
end

# need a custom representation of Union types since arrow unions
# are ordered, and possibly indirected via separate typeIds array
# here, T is Meta.UnionMode.Dense or Meta.UnionMode.Sparse,
# typeIds is a NTuple{N, Int32}, and U is a Tuple{...} of the
# unioned types
struct UnionT{T, typeIds, U}
end

unionmode(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = T
typeids(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = typeIds
Base.eltype(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = U

# convenience wrappers for signaling that an array shoudld be written
# as with dense/sparse union arrow buffers
struct DenseUnionVector{T, U} <: AbstractVector{UnionT{Meta.UnionMode.Dense, nothing, U}}
    itr::T
end

DenseUnionVector(x::T) where {T} = DenseUnionVector{T, Tuple{eachunion(eltype(x))...}}(x)
Base.IndexStyle(::Type{<:DenseUnionVector}) = Base.IndexLinear()
Base.size(x::DenseUnionVector) = (length(x.itr),)
Base.eltype(x::DenseUnionVector{T, U}) where {T, U} = UnionT{Meta.UnionMode.Dense, nothing, U}
Base.iterate(x::DenseUnionVector, st...) = iterate(x.itr, st...)
Base.getindex(x::DenseUnionVector, i::Int) = getindex(x.itr, i)

struct SparseUnionVector{T, U} <: AbstractVector{UnionT{Meta.UnionMode.Sparse, nothing, U}}
    itr::T
end

SparseUnionVector(x::T) where {T} = SparseUnionVector{T, Tuple{eachunion(eltype(x))...}}(x)
Base.IndexStyle(::Type{<:SparseUnionVector}) = Base.IndexLinear()
Base.size(x::SparseUnionVector) = (length(x.itr),)
Base.eltype(x::SparseUnionVector{T, U}) where {T, U} = UnionT{Meta.UnionMode.Sparse, nothing, U}
Base.iterate(x::SparseUnionVector, st...) = iterate(x.itr, st...)
Base.getindex(x::SparseUnionVector, i::Int) = getindex(x.itr, i)

# iterate a Julia Union{...} type, producing an array of unioned types
function eachunion(U::Union, elems=nothing)
    if elems === nothing
        return eachunion(U.b, Type[U.a])
    else
        push!(elems, U.a)
        return eachunion(U.b, elems)
    end
end

function eachunion(T, elems)
    push!(elems, T)
    return elems
end

# dense union child array producer
# for dense union children, we split the parent array
# into N children arrays, each with one of the unioned types
# only the 1st child needs to include missing values
struct Filtered{T, I}
    itr::I
    len::Int64
end

function filtered(::Type{T}, itr::I) where {T, I}
    len = Int64(0)
    for x in itr
        len += x isa T
    end
    return Filtered{T, I}(itr, len)
end

Base.length(f::Filtered) = f.len
Base.eltype(f::Filtered{T}) where {T} = T

function Base.iterate(f::Filtered{T}, st=()) where {T}
    st = iterate(f.itr, st...)
    st === nothing && return nothing
    x, state = st
    while !(x isa T)
        st = iterate(f.itr, state)
        st === nothing && return nothing
        x, state = st
    end
    return x, (state,)
end

# sparse union child array producer
# for sparse unions, we split the parent array into
# N children arrays, each having the same length as the parent
# but with one child array per unioned type; each child
# should include the elements from parent of its type
# and other elements can be missing/default
struct Replaced{T, I}
    itr::I
end

replaced(::Type{T}, itr::I) where {T, I} = Replaced{T, I}(itr)

Base.length(r::Replaced) = _length(r.itr)
Base.eltype(r::Replaced{T}) where {T} = Union{T, Missing}

function Base.iterate(r::Replaced{T}, st=()) where {T}
    st = iterate(r.itr, st...)
    st === nothing && return nothing
    x, state = st
    return ifelse(x isa T, x, missing), (state,)
end

# convenience wrapper to signal that an input column should be
# dict encoded when written to the arrow format
# note that only top-level columns are supported for dict encoding
# currently; (i.e. no nested dict encoding)
struct DictEncode{T, A <: AbstractVector} <: AbstractVector{T}
    data::A
end

DictEncode(x::A) where {A} = DictEncode{eltype(A), A}(x)
Base.IndexStyle(::Type{<:DictEncode}) = Base.IndexLinear()
Base.size(x::DictEncode) = (length(x.data),)
Base.eltype(x::DictEncode{T, A}) where {T, A} = T
Base.iterate(x::DictEncode, st...) = iterate(x.data, st...)
Base.getindex(x::DictEncode, i::Int) = getindex(x.data, i)

encodingtype(n) = n < div(typemax(Int8), 2) ? Int8 : n < div(typemax(Int16), 2) ? Int16 : n < div(typemax(Int32), 2) ? Int32 : Int64

struct Converter{T, A} <: AbstractVector{T}
    data::A
end

converter(::Type{T}, x::A) where {T, A} = Converter{eltype(A) >: Missing ? Union{T, Missing} : T, A}(x)

Base.IndexStyle(::Type{<:Converter}) = Base.IndexLinear()
Base.size(x::Converter) = (length(x.data),)
Base.eltype(x::Converter{T, A}) where {T, A} = T
Base.getindex(x::Converter{T}, i::Int) where {T} = convert(T, getindex(x.data, i))

maybemissing(::Type{T}) where {T} = T === Missing ? Missing : Base.nonmissingtype(T)

function getfooter(filebytes)
    len = readbuffer(filebytes, length(filebytes) - 9, Int32)
    FlatBuffers.getrootas(Meta.Footer, filebytes[end-(9 + len):end-10], 0)
end

function getrb(filebytes)
    f = getfooter(filebytes)
    rb = f.recordBatches[1]
    return filebytes[rb.offset+1:(rb.offset+1+rb.metaDataLength)]
    # FlatBuffers.getrootas(Meta.Message, filebytes, rb.offset)
end

function readmessage(filebytes, off=9)
    @assert readbuffer(filebytes, off, UInt32) === 0xFFFFFFFF
    len = readbuffer(filebytes, off + 4, Int32)
    @show len
    FlatBuffers.getrootas(Meta.Message, filebytes, off + 8)
end


text:
# ------------------------------------------------------------------------------------------
# # Julia is fast
#
# Very often, benchmarks are used to compare languages.  These benchmarks can lead to long
# discussions, first as to exactly what is being benchmarked and secondly what explains the
# differences.  These simple questions can sometimes get more complicated than you at first
# might imagine.
#
# The purpose of this notebook is for you to see a simple benchmark for yourself.  One can
# read the notebook and see what happened on the author's Macbook Pro with a 4-core Intel
# Core I7, or run the notebook yourself.
#
# (This material began life as a wonderful lecture by Steven Johnson at MIT:
# https://github.com/stevengj/18S096-iap17/blob/master/lecture1/Boxes-and-registers.ipynb.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # `sum`: An easy enough function to understand
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Consider the  **sum** function `sum(a)`, which computes
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i,
# $$
# where $n$ is the length of `a`.
# ------------------------------------------------------------------------------------------

a = rand(10^7) # 1D vector of random numbers, uniform on [0,1)

sum(a)

# ------------------------------------------------------------------------------------------
# The expected result is 0.5 * 10^7, since the mean of each entry is 0.5
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Benchmarking a few ways in a few languages
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia has a `BenchmarkTools.jl` package for easy and accurate benchmarking:
# ------------------------------------------------------------------------------------------

#Pkg.add("BenchmarkTools")

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# #  1. The C language
#
# C is often considered the gold standard: difficult on the human, nice for the machine.
# Getting within a factor of 2 of C is often satisfying. Nonetheless, even within C, there
# are many kinds of optimizations possible that a naive C writer may or may not get the
# advantage of.
#
# The current author does not speak C, so he does not read the cell below, but is happy to
# know that you can put C code in a Julia session, compile it, and run it. Note that the
# `"""` wrap a multi-line string.
# ------------------------------------------------------------------------------------------

C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # make a temporary file

# compile to a shared library by piping C_code to gcc
# (works only if you have gcc installed):

open(`gcc  -fPIC -O3 -ffast-math -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# define a Julia function that calls the C function:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a) # type \approx and then <TAB> to get the ≈ symbolb

≈  # alias for the `isapprox` function

?isapprox

# ------------------------------------------------------------------------------------------
# We can now benchmark the C code directly from Julia:
# ------------------------------------------------------------------------------------------

c_bench = @benchmark c_sum($a)

println("C: Fastest time was $(minimum(c_bench.times) / 1e6) msec")

#d = Dict()  # a "dictionary", i.e. an associative array
d["C"] = minimum(c_bench.times) / 1e6  # in milliseconds
d

using Plots
gr()

t = c_bench.times / 1e6 # times in milliseconds
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ------------------------------------------------------------------------------------------
# # 2. Python's built in `sum`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The `PyCall` package provides a Julia interface to Python:
# ------------------------------------------------------------------------------------------

#Pkg.add("PyCall")

using PyCall

# Call a low-level PyCall function to get a Python list, because
# by default PyCall will convert to a NumPy array instead (we benchmark NumPy below):

apy_list = PyCall.array2py(a, 1, 1)

# get the Python built-in "sum" function:
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($apy_list)

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 3. Python: `numpy`
#
# ## Takes advantage of hardware "SIMD", but only works when it works.
#
# `numpy` is an optimized C library, callable from Python.
# It may be installed within Julia as follows:
# ------------------------------------------------------------------------------------------

using Conda
#Conda.add("numpy")

numpy_sum = pyimport("numpy")["sum"]
apy_numpy = PyObject(a) # converts to a numpy array by default

py_numpy_bench = @benchmark $numpy_sum($apy_numpy)

numpy_sum(apy_list) # python thing

numpy_sum(apy_list) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 4. Python, hand-written
# ------------------------------------------------------------------------------------------

py"""
def py_sum(a):
    s = 0.0
    for x in a:
        s = s + x
    return s
"""

sum_py = py"py_sum"

py_hand = @benchmark $sum_py($apy_list)

sum_py(apy_list)

sum_py(apy_list) ≈ sum(a)

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 5. Julia (built-in)
#
# ## Written directly in Julia, not in C!
# ------------------------------------------------------------------------------------------

@which sum(a)

j_bench = @benchmark sum($a)

d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 6. Julia (hand-written)
# ------------------------------------------------------------------------------------------

function mysum(A)
    s = 0.0  # s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end

j_bench_hand = @benchmark mysum($a)

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # Summary
# ------------------------------------------------------------------------------------------

for (key, value) in sort(collect(d))
    println(rpad(key, 20, "."), lpad(round(value, 1), 8, "."))
end

for (key, value) in sort(collect(d), by=x->x[2])
    println(rpad(key, 20, "."), lpad(round(value, 2), 10, "."))
end

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/22. Julia is fast.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
import Pangolin
using ModernGL

## create window

width = 640;#//1280;
height = 480;#//980;
panel = 0;#//205;
name = "Main2"

Pangolin.create_window(name, width, height)

# glEnable(GL_DEPTH_TEST);

## create s_cam

s_cam = Pangolin.OpenGlRenderState()
handler = Pangolin.Handler3D(s_cam)
## create d_cam
# void* pangolin_create_display(void* renderstate_ptr)
d_cam = Pangolin.View(handler)

while !Pangolin.should_quit()

      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
      # sleep(0.1)

      Pangolin.activate(d_cam, s_cam)
      println(Pangolin.get_projection_view_mat(s_cam))

      Pangolin.glDrawColouredCube()

      Pangolin.finish_frame()

end

Pangolin.destroy_window(name)


text:
using ..TopOpt: ElementFEAInfo, ElementMatrix
using ..TopOpt.TopOptProblems: make_cload, convert

"""
    ElementFEAInfo(sp, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}) where {mat_type}

Constructs an instance of `ElementFEAInfo` from a stiffness **truss** problem `sp` using a Gaussian quadrature order of `quad_order`. The element matrix and vector types will be:
1. `SMatrix` and `SVector` if `mat_type` is `:SMatrix` or `:Static`, the default,
2. `MMatrix` and `MVector` if `mat_type` is `:MMatrix`, or
3. `Matrix` and `Vector` otherwise.

The static matrices and vectors are more performant and GPU-compatible therefore they are used by default.
"""
function ElementFEAInfo(sp::TrussProblem, quad_order = 1, ::Type{Val{mat_type}} = Val{:Static},) where {mat_type}
    # weights: self-weight element load vectors, all zeros now
    Kes, weights, cellvalues, facevalues = make_Kes_and_fes(sp, quad_order, Val{mat_type})
    element_Kes = convert(Vector{<:ElementMatrix}, Kes;
        bc_dofs = sp.ch.prescribed_dofs, dof_cells = sp.metadata.dof_cells)

    # * concentrated load
    # ? why convert a sparse vector back to a Vector?
    fixedload = Vector(make_cload(sp))
    # * distributed load (if any)
    # assemble_f!(fixedload, sp, dloads)

    cellvolumes = get_cell_volumes(sp, cellvalues)
    cells = sp.ch.dh.grid.cells
    ElementFEAInfo(element_Kes, weights, fixedload, cellvolumes,
        cellvalues, facevalues,
        sp.metadata, sp.black, sp.white, sp.varind, cells)
end

####################################

function get_cell_volumes(sp::TrussProblem{xdim, T}, cellvalues) where {xdim, T}
    dh = sp.ch.dh
    As = getA(sp)
    cellvolumes = zeros(T, getncells(dh.grid))
    for (i, cell) in enumerate(CellIterator(dh))
        truss_reinit!(cellvalues, cell, As[i])
        cellvolumes[i] = sum(Ferrite.getdetJdV(cellvalues, q_point) for q_point in 1:Ferrite.getnquadpoints(cellvalues))
    end
    return cellvolumes
end


text:
##  ----------------------------------------------------------------------------
# pol vox volumen
let
    model0 = iJR.load_model("max_model")

    bins = 50
    Hd_rxns_map = iJR.load_rxns_map()
    offsetf = 1.1

    exglcidx = ChU.rxnindex(model0, iJR.EX_GLC_IDER)
    biomidx = ChU.rxnindex(model0, iJR.BIOMASS_IDER)
    exglcL, exglcU = ChU.bounds(model0, exglcidx)

    maxD = maximum(Hd.val(:D)) * offsetf
    max_cgD_X = -maximum(Hd.ciD_X(:GLC)) * offsetf

    Ds = range(0.01, maxD; length = bins)
    cgD_Xs = range(max_cgD_X, exglcU; length = bins)

    box_vols = zeros(bins, bins)

    Hd_ids = ["GLC", "AC", "NH4"]
    model_ids = [Hd_rxns_map[Hd_id] for Hd_id in Hd_ids]
    model_idxs = [ChU.rxnindex(model0, model_id) for model_id in model_ids]

    # feeding task
    nths = nthreads()
    Ch = Channel(nths) do Ch_
        @showprogress for (Di, D) in enumerate(Ds)
            for (cgD_Xi, cgD_X) in enumerate(cgD_Xs)
                put!(Ch_, (Di, D, cgD_Xi, cgD_X))
            end
        end
    end

    # compute volume map
    @threads for _ in 1:nths
        thid = threadid()
        thid == 1 && nths > 1 && continue
        model = deepcopy(model0)
        for (Di, D, cgD_Xi, cgD_X) in Ch

            # Reduce Pol
            ChU.lb!(model, exglcidx, cgD_X)
            ChU.bounds!(model, biomidx, D, D)
            try
                L, U = ChLP.fva(model, model_idxs)
                vol = prod(abs.(L .- U))
                box_vols[Di, cgD_Xi] = max(0.0, log10(vol + 1e-50))
            catch
                box_vols[Di, cgD_Xi] = NaN
            end
        end
    end

    # vol map
    p = heatmap(Ds, cgD_Xs, box_vols';
        title = "Polytope Box volume", label = "",
        xlabel = "D", ylabel = "cgD/ X"
    )

    # exp vals
    exp_Ds = [Hd.val(:D, exp) for exp in EXPS]
    exp_cgD_Xs = [-Hd.ciD_X(:GLC, exp) for exp in EXPS]
    scatter!(p, exp_Ds, exp_cgD_Xs;
        label = "exp", color = :white,
        m = 8
    )

    mysavefig(p, "pol_box_volume"; bins)
end


text:
export quantile_residuals, pearson_residuals

function quantile_residuals(obs::Vector{T}, gas::Model{D, T};
                            initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}

    n = length(obs)
    len_ini_par = length(initial_params)
    quant_res = Vector{T}(undef, n - len_ini_par)
    params_fitted = score_driven_recursion(gas, obs; initial_params = initial_params)

    for t in axes(params_fitted[len_ini_par + 1:end - 1, :], 1)
        dist = update_dist(D, params_fitted[len_ini_par + 1:end - 1, :], t)
        # Continuous distributions PIT
        if supertype(D).parameters[2] == Continuous
            prob = cdf(dist, obs[t + len_ini_par])
            quant_res[t] = quantile(Normal(0, 1), prob)
        else # Discrete distributions PIT
            lb = cdf(dist, obs[t + len_ini_par] - 1)
            ub = cdf(dist, obs[t + len_ini_par])
            prob = rand(Uniform(lb, ub))
            quant_res[t] = quantile(Normal(0, 1), prob)
        end

        # Treat possible infinity values
        if quant_res[t] == Inf
            quant_res[t] = 1e5
        elseif quant_res[t] == -Inf
            quant_res[t] = -1e5
        end
    end

    return quant_res
end

function pearson_residuals(obs::Vector{T}, gas::Model{D, T};
                             initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}

    n = length(obs)
    len_ini_par = length(initial_params)
    pearson = Vector{T}(undef, n - len_ini_par)
    params_fitted = score_driven_recursion(gas, obs; initial_params = initial_params)

    for t in axes(params_fitted[len_ini_par + 1:end - 1, :], 1)
        dist = update_dist(D, params_fitted[len_ini_par + 1:end - 1, :], t)
        pearson[t] = (obs[t + len_ini_par] - mean(dist))/sqrt(var(dist))
    end

    return pearson
end


text:
function init(datafile, ::Type{ExaPF.ReducedSpaceEvaluator})
    constraints = Function[
        ExaPF.voltage_magnitude_constraints,
        ExaPF.reactive_power_constraints,
        ExaPF.active_power_constraints,
        # ExaPF.flow_constraints,
    ]
    return ExaPF.ReducedSpaceEvaluator(datafile; constraints=constraints)
end
function init(datafile, ::Type{ExaPF.ProxALEvaluator})
    nlp = ExaPF.ReducedSpaceEvaluator(datafile)
    time = ExaPF.Normal
    return ExaPF.ProxALEvaluator(nlp, time)
end
function init(datafile, ::Type{ExaPF.SlackEvaluator})
    return ExaPF.SlackEvaluator(datafile)
end

@testset "AugLagEvaluator with $Evaluator backend" for Evaluator in [
    ExaPF.ReducedSpaceEvaluator,
    ExaPF.ProxALEvaluator,
    ExaPF.SlackEvaluator,
]
    @testset "Inactive constraints" begin
        datafile = joinpath(INSTANCES_DIR, "case9.m")
        # Build reference evaluator
        nlp = init(datafile, Evaluator)
        u0 = ExaPF.initial(nlp)
        # Build penalty evaluator
        pen = ExaPF.AugLagEvaluator(nlp, u0)

        u = u0
        # Update nlp to stay on manifold
        ExaPF.update!(pen, u)
        # Compute objective
        c = ExaPF.objective(pen, u)
        c_ref = ExaPF.inner_objective(pen, u)
        @test isa(c, Real)
        @test c >= c_ref

        # Compute gradient of objective
        g = similar(u)
        g_ref = similar(u)
        fill!(g, 0)
        ExaPF.gradient!(pen, g, u)
        fill!(g_ref, 0)
        function reduced_cost(u_)
            ExaPF.update!(pen, u_)
            return ExaPF.objective(pen, u_)
        end
        g_ref = FiniteDiff.finite_difference_gradient(reduced_cost, u)
        @test isapprox(g_ref, g, rtol=1e-5)
        # Update penalty weigth
        ExaPF.update_penalty!(pen)
        # Utils
        inf_pr1 = ExaPF.primal_infeasibility(nlp, u)
        @test inf_pr1 >= 0.0
        # Test reset
        ExaPF.reset!(pen)
    end
    @testset "Active constraints" begin
        datafile = joinpath(INSTANCES_DIR, "case57.m")
        # Build reference evaluator
        nlp = init(datafile, Evaluator)
        u0 = ExaPF.initial(nlp)
        w♭, w♯ = ExaPF.bounds(nlp, ExaPF.Variables())
        # Build penalty evaluator
        for scaling in [true, false]
            ExaPF.reset!(nlp)
            pen = ExaPF.AugLagEvaluator(nlp, u0; scale=scaling)
            u = w♭
            # Update nlp to stay on manifold
            ExaPF.update!(pen, u)
            # Compute objective
            c = ExaPF.objective(pen, u)
            c_ref = ExaPF.inner_objective(pen, u)
            @test isa(c, Real)
            # For case57.m some constraints are active, so penalty are >= 0
            @test c > c_ref
            inf_pr2 = ExaPF.primal_infeasibility(pen, u)
            @test inf_pr2 > 0.0

            pen.ρ = 10.0
            # Update penalty weigth with a large factor to have
            # a meaningful derivative check
            ExaPF.update_penalty!(pen, η=1e3)
            ExaPF.update_multipliers!(pen)
            ExaPF.update!(pen, u)
            obj = ExaPF.objective(pen, u)
            # Compute gradient of objective
            g = similar(u)
            fill!(g, 0)
            ExaPF.gradient!(pen, g, u)
            # Compare with finite differences
            function reduced_cost(u_)
                ExaPF.update!(pen, u_)
                return ExaPF.objective(pen, u_)
            end
            grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, u)
            @test isapprox(grad_fd, g, rtol=1e-6)

            # Test Hessian only on ReducedSpaceEvaluator and SlackEvaluator
            if isa(nlp, ExaPF.ReducedSpaceEvaluator) || isa(nlp, ExaPF.SlackEvaluator)
                n = length(u)
                ExaPF.update!(pen, u)
                hv = similar(u) ; fill!(hv, 0)
                w = similar(u) ; fill!(w, 0)
                w[1] = 1
                ExaPF.hessprod!(pen, hv, u, w)
                H = similar(u, n, n) ; fill!(H, 0)
                ExaPF.hessian!(pen, H, u)
                # Is Hessian vector product relevant?
                @test H * w ≈ hv
                # Is Hessian correct?
                hess_fd = FiniteDiff.finite_difference_hessian(reduced_cost, u)
                @test isapprox(H, hess_fd, rtol=1e-6)
            end
            # Test estimation of multipliers only on SlackEvaluator
            if isa(nlp, ExaPF.SlackEvaluator)
                λ = ExaPF.estimate_multipliers(pen, u)
            end
        end
    end
end


text:
#printIO
function show(io::IO, b::BasisBlade)
	if GIdLayout[3]=="GA" || GIdLayout[3]=="Projective"
		if b.index>0
			if b.index == 1
				print(io, "id")
			else
				print(io, "e")
				for i=1:GIdLayout[1]+GIdLayout[2]
					if LogicalBasisBlade[b.index][i]==true
						print(io, "$i" )
					end
				end
			end
		else
			if b.index == -1
				print(io, "-id")
			else
				if b.index==0
					print(io,"0")
				else
					print(io, "-e")
					for i=1:GIdLayout[1]+GIdLayout[2]
						if LogicalBasisBlade[-b.index][i]==true
							print(io, "$i" )
						end
					end
				end
			end
		end
	end
	if GIdLayout[3]=="Conformal"
		if b.index>0
			if b.index == 1
				print(io, "id")
			else
				print(io, "e")
				for i=1:GIdLayout[1]+GIdLayout[2]
					if LogicalBasisBlade[b.index][i]==true
						if i<GIdLayout[1]+GIdLayout[2]-1
							print(io, "$i" )
						else
							if i==GIdLayout[1]+GIdLayout[2]-1
								print(io,"+")
							else
								print(io,"-")
							end
						end
					end
				end
			end
		else
			if b.index == -1
				print(io, "-id")
			else
				if b.index==0
					print(io,"0")
				else
					print(io, "-e")
					for i=1:GIdLayout[1]+GIdLayout[2]
						if LogicalBasisBlade[-b.index][i]==true
							if i<GIdLayout[1]+GIdLayout[2]-1
								print(io, "$i" )
							else
								if i==GIdLayout[1]+GIdLayout[2]-1
									print(io,"+")
								else
									print(io,"-")
								end
							end
						end
					end
				end
			end
		end
	end
end
function show(io::IO,M::MultiVector)
	#len=length(M.comp)
	#if maximum(abs.(M.comp))<=LigaPrecision
	#	print(io,"0.0")
	#else
	#	if M.comp[1]!=0.0
	#		if M.comp[1]>0.0
	#			print(io,"( $(M.comp[1])$(BasisBlade(1)) )")
	#		else
	#	end
	#	for i=2:len
	#		if M.comp[i]!=0.0
	#			if M.comp[i]>0.0
	#				print(io," + ( $(M.comp[i])$(BasisBlade(i)) )")
	#			else
	#				print(io," - ( $(abs(M.comp[i]))$(BasisBlade(i)) )")
	#			end
	#		end
	#	end
	#end

		ind,val=findnz(M.comp)
		if maximum(abs.(M.comp))<=LigaPrecision
			print(io,"0.0")
		else
			k=true
			for i in ind
				if abs(M.comp[i])>LigaPrecision
					if k==true
						if M.comp[i]>0.0
							print(io,"( $(round(M.comp[i],digits=4)) $(BasisBlade(i)) )")
						else
							print(io," - ( $(round(abs(M.comp[i]),digits=4)) $(BasisBlade(i)) )")
						end
						k=false
					else
						if M.comp[i]>0.0
							print(io," + ( $(round(M.comp[i],digits=4)) $(BasisBlade(i)) )")
						else
							print(io," - ( $(round(abs(M.comp[i]),digits=4)) $(BasisBlade(i)) )")
						end
					end
				end
			end
		end
end

function show(io::IO,B::Blade)
	print(io,"($(B.scalar)) * $(B.vectors[1])")
	for i=2:B.grade
		print(io," ∧ $(B.vectors[i])")
	end
end


text:
module Autoerror

include("dualnumbers.jl")
include("basic_extensions.jl")

function uncert(f::Dual, σ::Vector{<:Real})
  x = f.grad .* σ
  (val = f.val, err = √(x'x))
end

end

max_stars_repo_name: kaosmicadei/AutoError

max_issues_repo_name: kaosmicadei/AutoError

max_forks_repo_name: kaosmicadei/AutoError


text:
function install_projs(root, deep, dry_run = false)

    ## ---------------------------------------------------------------------------------------------
    # Looking for projects
    cd(root)
    println("\n", "-" ^ 45)
    @info("Looking for projects", root)
    projs = find_projects(root)

    ## ---------------------------------------------------------------------------------------------
    # Install source projects
    pkgs_pool = Dict()
    let
        println("\n", "-" ^ 45)
        @info("Installing source projects", root)

        for proj_file in projs
            try
                # Install and update
                println("\n", "-" ^ 45)
                proj_file = relpath(proj_file, root)
                @info("Found on disk", proj_file)
                proj_file = joinpath(root, proj_file)

                # Collect devs
                pkgs = extract_pkgdeps(proj_file)
                merge!(pkgs_pool, pkgs)

                # Installing
                if !dry_run

                    # Temp env
                    proj_dir = tempenv(dirname(proj_file))
                    try
                        proj_file = joinpath(proj_dir, basename(proj_file))
                        Pkg.activate(proj_file)

                        for manf_name in Base.manifest_names
                            manf_file = joinpath(proj_dir, manf_name)
                            isfile(manf_file) && Pkg.instantiate()
                        end
                        Pkg.resolve()
                        Pkg.update()
                        Pkg.build()

                    finally
                        rm(proj_dir; force = true, recursive = true)
                    end
                end

            catch err
                @warn("ERROR", proj_file, err)
            end
        end
    end

    ## ---------------------------------------------------------------------------------------------
    # Install version range
    let
        pkgs_count = length(pkgs_pool)
        println("\n", "-" ^ 45)
        @info("Installing pkg versions", pkgs_count, deep)

        deep <= 0 && return
        for (uuidpkg, pkgname) in pkgs_pool
            try
                versions = VersionNumber[]
                for pkg_regpath in findin_regs(pkgname, uuidpkg)
                    uuidfile = extract_reg_pkguuid(pkg_regpath)
                    uuidpkg != uuidfile && continue
                    push!(versions, extract_reg_pkgversions.(pkg_regpath)...)
                end
                sort!(unique!(versions), rev = true)

                println("\n" ^ 3, "-" ^ 45)
                @info("Found in register", pkgname, uuidpkg, length(versions))

                # Installing version range
                c = 0
                for version in versions

                    c >= deep && break
                    c += 1

                    println("\n", "-" ^ 45)
                    @info("Installing", pkgname, uuidpkg, version)

                    if !dry_run
                        proj_dir = tempenv()
                        mkpath(proj_dir)
                        try
                            Pkg.activate(proj_dir)
                            pkg = Pkg.Types.PackageSpec(pkgname, Base.UUID(uuidpkg), version)
                            Pkg.add(pkg)
                            Pkg.build()

                        catch err
                            @warn("ERROR", pkgname, uuidpkg, version, err)
                        finally
                            rm(proj_dir; force = true, recursive = true)
                        end
                    end
                end
            catch err
                @warn("ERROR", pkgname, uuidpkg, err)
            end
        end
    end
end


text:
## Diagonostic code
##
import LinearAlgebra: diagm

## create a rotation matrix
function rotm(a::T,b, i, N) where {T}
    r = matrix{T}(I, N,  N)
    r[i:i+1, i:i+1] = [a -conj(b); b conj(a)]
    r
end

## make a rotator into a full matrix
function as_full(a::AbstractRotator{T}, N::Int) where {T}
    c,s = vals(a)
    i = idx(a)
    i < N || error("$i >= $N; too big")
    A = Matrix{Complex{T}}(I, N, N)
    A[i:i+1, i:i+1] = [c -conj(s); s conj(c)]
    A
end

# make zeros zeros
function zero_out!(A::Array{T}, tol=1e-12) where {T}
    A[norm.(A) .<= tol] = zero(T)
end

function zero_out!(A::Array{Complex{T}}, tol=1e-12) where {T}
    for i in eachindex(A)
        c = A[i]
        cr, ci = real(c), imag(c)
        if abs(cr) < tol
            cr = zero(T)
        end
        if abs(ci) < tol
            ci = zero(T)
        end
        A[i] = complex(cr, ci)
    end
end

## diagnostic

## create Full matrix from state object. For diagnostic purposes.
# we may or may not have a diagonal matrix to keep track or
D_matrix(state::FactorizationType{T, Val{:SingleShift}, P, Tw}) where {T,P,Tw} = diagm(0 => state.D)
D_matrix(state::FactorizationType) = diagm(0 => ones(state.N+1))#I

## We can compute yt from the decomposition of R into: R = Ct *(B + yt e1) by multiplying on left by e_n+1^T...
function _compute_yt(Cts, Bs, N, T)
    Ct = as_full(Cts[1], N+1); for i in 2:N Ct =  as_full(Cts[i],N+1)*Ct end
    B = as_full(Bs[1],N+1); for i in 2:N B = B * as_full(Bs[i],N+1) end
    e1 = zeros(T, N+1); e1[1]=one(T)
    en = zeros(T, N+1); en[N] = one(T)
    en1 = zeros(T, N+1); en1[N+1] = one(T)

    rho = transpose(en1) * Ct * e1  # scalar
    yt = -1/rho * transpose(en1)  * Ct * B
    # clean
    for i in eachindex(yt)
        if norm(yt[i]) < 1e-12
            yt[i] = 0
        end
    end

    Ct, B, yt
end

function _compute_R(Ct, B, yt, N, T)
    ## we have R = Z + x = Ct * (B  + e1 * yt)
    Z = Ct * B
    zero_out!(Z)

    e1 = zeros(T, N+1); e1[1]=one(T)
    x = Ct * e1 * yt

    R = (Z + x)
    zero_out!(R)
    R
end

function _compute_R(Ct, B, yt, D, N, T)
    ## we have R = Z + x = Ct * (B * D + e1 * yt)
    Z = Ct * B
    zero_out!(Z)

    e1 = zeros(T, N+1); e1[1]=one(T)
    x = Ct * e1 * yt
    R = Z + x
    R = diagm(0 => diag(D)) * R
    zero_out!(R)
    R
end

function compute_R(Cts, Bs, N, T)
    Ct, B, yt = _compute_yt(Cts, Bs, N, T)
    R = _compute_R(Ct, B, yt, N, T)
    R
end

function compute_R(Cts, Bs, D, N, T)
    Ct, B, yt = _compute_yt(Cts, Bs, N, T)
    R = _compute_R(Ct, B, yt, D, N, T)
    R
end

# Go from efficiently stored state to full matrix
function full_matrix(state::FactorizationType{T, St, Val{:NoPencil}, Val{:NotTwisted}}, what=:A) where {T, St}
    N = state.N
    Q = as_full(state.Q[1],N+1); for i in 2:N Q = Q * as_full(state.Q[i],N+1) end
    D = D_matrix(state)

    Ct, B, yt = _compute_yt(state.Ct, state.B, N, T)
    R = _compute_R(Ct, B, yt, D, N, T)

    what == :R && return R

    A = Q * R
    zero_out!(A)
    A
end

function full_matrix(state::FactorizationType{T, St, Val{:HasPencil}, Val{:NotTwisted}}, what=:A) where {T, St}
    n = state.N

    Q = (prod(as_full.(state.Q, n+1)))[1:n, 1:n]
    P = diagm(0 => ones(n+1))[:, 1:n]
    if St == Val{:SingleShift}
        V = P' * compute_R(state.Ct, state.B, state.D, state.N, T) * P
        W = P' * compute_R(state.Ct1, state.B1, state.D1, state.N, T) * P
    else
        V = P' * compute_R(state.Ct, state.B, state.N, T) * P
        W = P' * compute_R(state.Ct1, state.B1, state.N, T) * P
    end

    Q * V * inv(W)
end

# simple graphic to show march of algorithm
function show_status(state::FactorizationType)
    qs = [norm(u.s) for u in state.Q[state.ctrs.start_index:state.ctrs.stop_index]]
    minq = length(qs) > 0 ?  minimum(qs) : 0.0

    x = fill(".", state.N+2)
    x[state.ctrs.zero_index+1] = "α"
    x[state.ctrs.start_index+1] = "x"
    x[state.ctrs.stop_index+2] = "Δ"
    println(join(x, ""), " ($minq)")
end


text:
module Reference

import CFTime: isleap, DateTimeJulian, DateTimeProlepticGregorian

# Adapted
# from https://github.com/Unidata/cftime/blob/dc75368cd02bbcd1352dbecfef10404a58683f94/src/cftime/_cftime.pyx
# Licence MIT
# by Jeff Whitaker (https://github.com/jswhit)

@inline function month_lengths(::Type{T}, year::Integer, has_year_zero) where T
    if isleap(T, year, has_year_zero)
        return (31,29,31,30,31,30,31,31,30,31,30,31)
    else
        return (31,28,31,30,31,30,31,31,30,31,30,31)
    end
end

@inline function datetuple_ymd(::Type{T}, delta_days, julian_gregorian_mixed, has_year_zero) where T
    year = 1858
    month = 11
    day = 17

    month_length = month_lengths(T, year, has_year_zero)

    n_invalid_dates =
        if julian_gregorian_mixed
            10
        else
            0
        end

    @inbounds while delta_days < 0
        if (year == 1582) && (month == 10) && (day > 14) && (day + delta_days < 15)
            delta_days -= n_invalid_dates    # skip over invalid dates
        end

        if day + delta_days < 1
            delta_days += day
            # decrement month
            month -= 1
            if month < 1
                month = 12
                year -= 1
                if (year == 0) && !has_year_zero
                    year = -1
                end
                month_length = month_lengths(T, year, has_year_zero)
            end

            day = month_length[month]
        else
            day += delta_days
            delta_days = 0
        end
    end

    @inbounds while delta_days > 0
        if (year == 1582) && (month == 10) && (day < 5) && (day + delta_days > 4)
            delta_days += n_invalid_dates    # skip over invalid dates
        end

        if day + delta_days > month_length[month]
            delta_days -= month_length[month] - (day - 1)
            # increment month
            month += 1
            if month > 12
                month = 1
                year += 1
                if (year == 0) && !has_year_zero
                    year = 1
                end
                month_length = month_lengths(T, year, has_year_zero)
            end
            day = 1
        else
            day += delta_days
            delta_days = 0
        end
    end

    return year,month,day
end

function datetuple_ymd(::Type{DateTimeProlepticGregorian},Z)
    has_year_zero = false
    julian_gregorian_mixed = false
    return datetuple_ymd(DateTimeProlepticGregorian, Z, julian_gregorian_mixed, has_year_zero)
end

end


text:
function get_incompatible_devices(devices_template::Dict)
    incompatible_device_types = Vector{DataType}()
    for model in values(devices_template)
        formulation = get_formulation(model)
        if formulation == FixedOutput
            if !isempty(get_services(model))
                @info "$(formulation) for $(get_component_type(model)) is not compatible with the provision of reserve services"
            end
            push!(incompatible_device_types, get_component_type(model))
        end
    end
    return incompatible_device_types
end

function construct_services!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    services_template::ServicesModelContainer,
    devices_template::DevicesModelContainer,
)
    isempty(services_template) && return
    incompatible_device_types = get_incompatible_devices(devices_template)

    function _construct_valid_services!(service_model::ServiceModel)
        @debug "Building $(service_model.component_type) with $(service_model.formulation) formulation"
        services = service_model.component_type[]
        if validate_services!(
            service_model.component_type,
            services,
            incompatible_device_types,
            sys,
        )
            construct_service!(
                optimization_container,
                services,
                sys,
                service_model,
                devices_template,
                incompatible_device_types,
            )
        end
    end

    groupservice = nothing

    for (key, service_model) in services_template
        if service_model.formulation === GroupReserve  # group service needs to be constructed last
            groupservice = key
            continue
        end
        _construct_valid_services!(service_model)
    end
    groupservice === nothing || _construct_valid_services!(services_template[groupservice])
    return
end

function construct_service!(
    optimization_container::OptimizationContainer,
    services::Vector{SR},
    sys::PSY.System,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Vector{<:DataType},
) where {SR <: PSY.Reserve}
    services_mapping = PSY.get_contributing_device_mapping(sys)
    time_steps = model_time_steps(optimization_container)
    names = [PSY.get_name(s) for s in services]

    if model_has_parameters(optimization_container)
        add_param_container!(
            optimization_container,
            UpdateRef{SR}("service_requirement", "requirement"),
            names,
            time_steps,
        )
    end

    add_cons_container!(
        optimization_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_devices =
            services_mapping[(type = SR, name = PSY.get_name(service))].contributing_devices
        if !isempty(incompatible_device_types)
            contributing_devices = [
                d for d in contributing_devices if
                typeof(d) ∉ incompatible_device_types && PSY.get_available(d)
            ]
        end
        # Services without contributing devices should have been filtered out in the validation
        @assert !isempty(contributing_devices)
        # Variables
        add_variables!(
            optimization_container,
            ActiveServiceVariable,
            service,
            contributing_devices,
            RangeReserve(),
        )
        # Constraints
        service_requirement_constraint!(optimization_container, service, model)
        modify_device_model!(devices_template, model, contributing_devices)

        # Cost Function
        cost_function!(optimization_container, service, model)
    end

    if get_feedforward(model) !== nothing
        feedforward!(optimization_container, PSY.Device[], model, get_feedforward(model))
    end

    return
end

function construct_service!(
    optimization_container::OptimizationContainer,
    services::Vector{SR},
    sys::PSY.System,
    model::ServiceModel{SR, StepwiseCostReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Vector{<:DataType},
) where {SR <: PSY.Reserve}
    services_mapping = PSY.get_contributing_device_mapping(sys)
    time_steps = model_time_steps(optimization_container)
    names = [PSY.get_name(s) for s in services]
    # Does not use the standard implementation of add_variable!()
    add_variable!(
        optimization_container,
        ServiceRequirementVariable(),
        services,
        StepwiseCostReserve(),
    )
    add_cons_container!(
        optimization_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_devices =
            services_mapping[(
                type = typeof(service),
                name = PSY.get_name(service),
            )].contributing_devices
        if !isempty(incompatible_device_types)
            contributing_devices = [
                d for d in contributing_devices if
                typeof(d) ∉ incompatible_device_types && PSY.get_available(d)
            ]
        end
        # Variables
        add_variables!(
            optimization_container,
            ActiveServiceVariable,
            service,
            contributing_devices,
            StepwiseCostReserve(),
        )
        # Constraints
        service_requirement_constraint!(optimization_container, service, model)
        modify_device_model!(devices_template, model, contributing_devices)
        # Cost Function
        cost_function!(optimization_container, service, model)
    end

    if get_feedforward(model) !== nothing
        feedforward!(optimization_container, PSY.Device[], model, get_feedforward(model))
    end

    return
end

function construct_service!(
    optimization_container::OptimizationContainer,
    services::Vector{PSY.AGC},
    sys::PSY.System,
    ::ServiceModel{PSY.AGC, T},
    devices_template::Dict{Symbol, DeviceModel},
    ::Vector{<:DataType},
) where {T <: AbstractAGCFormulation}
    # Order is important in the addition of these variables
    # for device_model in devices_template
    # TODO: make a check for the devices' models
    #end
    agc_areas = [PSY.get_area(agc) for agc in services]
    areas = PSY.get_components(PSY.Area, sys)
    for area in areas
        if area ∉ agc_areas
            #    throw(IS.ConflictingInputsError("All area most have an AGC service assigned in order to model the System's Frequency regulation"))
        end
    end
    add_variables!(optimization_container, SteadyStateFrequencyDeviation)
    add_variables!(optimization_container, AreaMismatchVariable, areas, T())
    add_variables!(optimization_container, SmoothACE, areas, T())
    add_variables!(optimization_container, LiftVariable, areas, T())
    add_variables!(optimization_container, ActivePowerVariable, areas, T())
    add_variables!(optimization_container, DeltaActivePowerUpVariable, areas, T())
    add_variables!(optimization_container, DeltaActivePowerDownVariable, areas, T())
    # add_variables!(optimization_container, AdditionalDeltaActivePowerUpVariable, areas)
    # add_variables!(optimization_container, AdditionalDeltaActivePowerDownVariable, areas)
    balancing_auxiliary_variables!(optimization_container, sys)

    absolute_value_lift(optimization_container, areas)
    frequency_response_constraint!(optimization_container, sys)
    area_control_initial_condition!(optimization_container, services, T())
    smooth_ace_pid!(optimization_container, services)
    aux_constraints!(optimization_container, sys)
end

"""
    Constructs a service for StaticReserveGroup.
"""
function construct_service!(
    optimization_container::OptimizationContainer,
    services::Vector{SR},
    ::PSY.System,
    model::ServiceModel{SR, GroupReserve},
    ::Dict{Symbol, DeviceModel},
    ::Vector{<:DataType},
) where {SR <: PSY.StaticReserveGroup}
    time_steps = model_time_steps(optimization_container)
    names = [PSY.get_name(s) for s in services]

    if model_has_parameters(optimization_container)
        add_param_container!(
            optimization_container,
            UpdateRef{SR}("service_requirement", "requirement"),
            names,
            time_steps,
        )
    end

    add_cons_container!(
        optimization_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_services = PSY.get_contributing_services(service)

        # check if variables exist
        check_activeservice_variables(optimization_container, contributing_services)
        # Constraints
        service_requirement_constraint!(
            optimization_container,
            service,
            model,
            contributing_services,
        )
    end
    return
end

function construct_service!(
    optimization_container::OptimizationContainer,
    services::Vector{SR},
    sys::PSY.System,
    model::ServiceModel{SR, RampReserve},
    devices_template::Dict{Symbol, DeviceModel},
    incompatible_device_types::Vector{<:DataType},
) where {SR <: PSY.Reserve}
    services_mapping = PSY.get_contributing_device_mapping(sys)
    time_steps = model_time_steps(optimization_container)
    names = [PSY.get_name(s) for s in services]

    if model_has_parameters(optimization_container)
        add_param_container!(
            optimization_container,
            UpdateRef{SR}("service_requirement", "requirement"),
            names,
            time_steps,
        )
    end

    add_cons_container!(
        optimization_container,
        make_constraint_name(REQUIREMENT, SR),
        names,
        time_steps,
    )

    for service in services
        contributing_devices =
            services_mapping[(
                type = typeof(service),
                name = PSY.get_name(service),
            )].contributing_devices
        if !isempty(incompatible_device_types)
            contributing_devices = [
                d for d in contributing_devices if
                typeof(d) ∉ incompatible_device_types && PSY.get_available(d)
            ]
        end
        # Variables
        add_variables!(
            optimization_container,
            ActiveServiceVariable,
            service,
            contributing_devices,
            RampReserve(),
        )
        # Constraints
        service_requirement_constraint!(optimization_container, service, model)
        ramp_constraints!(optimization_container, service, contributing_devices, model)
        modify_device_model!(devices_template, model, contributing_devices)

        # Cost Function
        cost_function!(optimization_container, service, model)
    end
    return
end


text:
using JuliennedArrays
using Test
using Documenter: doctest

if VERSION >= v"1.6"
    doctest(JuliennedArrays)
end

@testset "JuliennedArrays.jl" begin
    @testset "Align" begin
        Xs = [rand(2, 3) for _ = 1:4]
        X = @inferred Align(Xs, True(), False(), True())
        @test size(X) == (2, 4, 3)
        @test permutedims(cat(Xs...; dims = 3), (1, 3, 2)) == X
        @test X[1] == Xs[1][1] # test linear indexing
        @test X[1, 2, 3] == Xs[2][1, 3] # test cartesian indexing

        Xs = reshape(Xs, 1, 4)
        X = @inferred Align(Xs, True(), False(), True(), False())
        @test size(X) == (2, 1, 3, 4)
        @test permutedims(X, (1, 3, 2, 4))[:] == cat(Xs...; dims = 3)[:]

        # type is inferrable for constant integer alongs
        align_1_3(x) = Align(x, 1, 3)
        @inferred align_1_3(Xs)

        @test_throws DimensionMismatch Align([rand(2, 3) for _ in 1:4], 1)
        @test_throws MethodError Align(ones(2, 3, 4), 1, 2, 3)

        @test size(Align(Slices(Slices(randn(3, 4, 5, 2), 3), 3), 3)) == (3, 4, 2)
    end

    @testset "Slice" begin
        X = rand(2, 3, 4, 5)
        Xs = @inferred Slices(X, True(), False(), False(), False())

        Xs = Slices(X, 1)
        @test size(Xs) == (3, 4, 5)
        @test Xs[1, 1, 1] == X[:, 1, 1, 1]

        Xs = Slices(X, 2)
        @test size(Xs) == (2, 4, 5)
        @test Xs[1, 1, 1] == X[1, :, 1, 1]

        Xs = Slices(X, 1, 3)
        @test size(Xs) == (3, 5)
        @test Xs[1, 2] == X[:, 1, :, 2]
        @test Align(Xs, 1, 3) == X # Slices is the inverse of Align
        @test Xs[1] == X[:, 1, :, 1] # test linear indexing
        @test Xs[1, 3] == X[:, 1, :, 3] # test cartesian indexing

        # type is inferrable for constant integer alongs
        slices_1_3(x) = Slices(x, 1, 3)
        @inferred slices_1_3(X)

        X = rand(2, 3, 4, 5)
        @test_throws MethodError Slices(X, True())
        @test_throws MethodError Slices(X, True(), False(), False(), False(), False())
        @test_throws ArgumentError("All alongs values (5,) should be less than or equal to 4") Slices(X, 5)
    end
end


text:
# module fiberDataProcLib
# export  CalDatainFile
using Statistics
using Printf
using Dates
include("../../RawFileofflineProc/dataStruct.jl")

function getRelPos(CumSumDataobj::spaceCumSumData ,vPos::Int )
    return vPos- CumSumDataobj.startID +1

end
function getAbsPos(CumSumDataobj::spaceCumSumData ,vPos::Int )
    return vPos- CumSumDataobj.startID +1

end
function getSpaceLoopRng(spaceCumSumDataObj::spaceCumSumData)
    # println("RNG:  $(spaceCumSumDataObj.startID) : $(spaceCumSumDataObj.endID) " )
    return spaceCumSumDataObj.startID : spaceCumSumDataObj.endID
end

function InitObj(spaceCumSumDataObj::spaceCumSumData , AlrmProcStartID1::Int,AlrmProcEndtID1::Int , startID1::Int ,endID1::Int,maxTimeCnt1::Int= 100 )
    spaceCumSumDataObj.startID= startID1
    spaceCumSumDataObj.endID =endID1
    SegdataNum = endID1-startID1+1

    spaceCumSumDataObj.AlrmProcStartID=AlrmProcStartID1
    spaceCumSumDataObj.AlrmProcEndtID=AlrmProcEndtID1

    spaceCumSumDataObj.maxTimeCnt=maxTimeCnt1

    spaceCumSumDataObj.defaultSpaceCnt =10
    spaceCumSumDataObj.defaultTimeCnt = 1000

    spaceCumSumDataObj.m1 = Array{UInt32}(undef,spaceCumSumDataObj.defaultTimeCnt,SegdataNum)
    spaceCumSumDataObj.m1 .= 0
    spaceCumSumDataObj.rectSum = Array{UInt64}(undef,SegdataNum)
    spaceCumSumDataObj.curRow = 0
    spaceCumSumDataObj.dt= ""
    spaceCumSumDataObj.isInit = false

end

########################
#计算当前的行号
###########################
function getCurRowCnt(  temobj::spaceCumSumData)
    idx = temobj.curRow +1
    dataCnt =  length( @view temobj.m1[:,1] )
    temobj.curRow = idx %dataCnt
    # @printf( "====%d === %d==============\n" , idx ,dataCnt  )
    return idx
 end

###########  计算做累加计算需要的 矩阵 的时间累加时  的  下标，
###因为 矩阵的数据结构为 循环填写有效数据，所以一共有两对   时间下标，没对包含起点、终点
####################
function getTimeRng(tcMax,curtc,vTcNum ) # time is Row
        if curtc < 1 ;curtc = tcMax ;end
        if curtc > tcMax ;curtc = tcMax ;end
        if vTcNum > tcMax; vTcNum =tcMax;end

        tctop1=0;tcbottom1 =0
        tctop2=0;tcbottom2 = 0

        tctop1= curtc
        tcbottom1  = tctop1-  vTcNum+1
        if tcbottom1  <  1;
            tcbottom1 =1;
            cnt =  ( tctop1- tcbottom1)+1
            if  cnt < vTcNum;
                tctop2= tcMax
                tcbottom2 = tctop2- (vTcNum-cnt-1)
            end
            # @printf( "%d,%d,%d----,%d,%d,%d---%d\n",  tctop2,tcbottom2,(tc-cnt),  tctop1,tcbottom1, cnt  ,tc )
        end

      (  (  tctop1,tcbottom1), (tctop2,tcbottom2,) )
end

#############################
# 根据光缆空间位置的 中心位置 vspacePos ,以及计算累加的总空间点数vsc,计算出  累加计算的两个下标位置。
# 当 vsc为偶数时，调整为单数。
##############

function getSpaceRng( cumsumObj::spaceCumSumData, vspacePos,vsc)
        vsideMax =cumsumObj.endID
        vsideMin = cumsumObj.startID

        if vspacePos <  vsideMin || vspacePos  >  vsideMax ;
            error("vspacePos 越界  ")
        end

        if vsc % 2 ==0
            vsc+=1
        end
        left = vsideMin ;
        right = vsideMax;
        halfsc = trunc(Int,vsc/2)

        left = vspacePos-halfsc-1 ;
        right =vspacePos+halfsc;
        if left < vsideMin ;
            left=0;
        end
        if  right >vsideMax ;
            right =vsideMax;
        end
       return  (left ,right)
end

function getRelativeSpaceRng( cumsumObj::spaceCumSumData, vposIndex,vsc)
    vsideMax =cumsumObj.endID
    vsideMin = cumsumObj.startID

    if vposIndex <  0 || vposIndex  >  vsideMax ;
        error("vspacePos 越界  ")
    end

    if vsc % 2 ==0
        vsc+=1
    end
    left = 0 ;
    right = vsideMax;
    halfsc = trunc(Int,vsc/2)

    left = vposIndex-halfsc-1 ;
    right =vposIndex+halfsc;
    if left < 0 ;
        left=0;
    end
    if  right >vsideMax  ;
        right =vsideMax;
    end
   return (left ,right)
end

######################
# vsc  计算均值时的 空间点颗粒 的 个数
# vTcNum  计算均值时的 时间点的 个数
# vPos  当前计算点的空间 位置，空间索引点 ，对应本结构下的 rectSum的下标位置,不是光缆的绝对参考 颗粒号
###############
function calRectSum( CumSumDataobj::spaceCumSumData, vPos ,vsc ,vTcNum   )
    temobj = CumSumDataobj
    rectsumObj =temobj.rectSum

    vsideMax  = length( rectsumObj )  #？

    timeRngMax1 = length(　temobj.m1[:,1])　 #？

    curtcIdx = temobj.curRow
    # (left,right )=getSpaceRng( vPos, vsideMax,vsc)
    (absleft1,absright1 )=getSpaceRng( temobj,vPos,vsc)
    a = getTimeRng(timeRngMax1,curtcIdx,vTcNum )
    # println(  "time rng : $a " )
    right = getRelPos(CumSumDataobj ,absright1 )
    if absleft1 != 0;
        left = getRelPos(CumSumDataobj ,absleft1 )
    else
        left=0
    end

    numbercnt =0
    rows=0;
    cols=0
    if   left == 0 ;
        cols =  right
    else
        cols = right-left
    end

    # left = absleft1###########
    # right= absright1 ##########
    # println(size( temobj.m1 ))
    sumv =0
    for ( top1,bottom1) in a;  ## time range loop
        # println( "absleft1,absright1,left,right,bottom1,top1::$absleft1,$absright1,$left,$right, $bottom1, $top1 ")
            # println( "VPos,absleft1,absright1,left,right,bottom1,top1::$vPos  :$absleft1,$absright1,$left,$right, $bottom1, $top1 ")
            if top1 == 0;
                continue;
            end

            rows += top1 - bottom1+1

            v1 = @view temobj.m1[ bottom1:top1,right  ]
            sum1= sum( v1)
            # println( sum1 )

            sum0=0 ;
            if left != 0;
                v0 =@view temobj.m1[bottom1:top1,left  ]
                sum0 = sum(  v0)
            end

            sumv += sum1-sum0

    end   # for ( top1,bottom1) in a;

        numbercnt =rows*cols
        if vTcNum*vsc != numbercnt
            # error(  "Logic Error vTcNum != numbercnt " )
            # println( "Logic Error vTcNum*vsc != numbercnt :$vTcNum * $vsc  != $numbercnt" )
        end
        # println( "(left,right):($left,$right)")
        curtcIdx = getRelPos(CumSumDataobj ,vPos )
        # numbercnt =1
        avg = round( UInt,sumv /numbercnt)
        rectsumObj[ curtcIdx ]  = avg
        # println( "curtcIdx:$curtcIdx,sumv:$sumv ,avg:$avg ")
        # println( "numbercnt,sumv ,sumv /numbercnt : $numbercnt,$sumv ,$avg:  "  )
        # rectsumObj[ vPos  ] =  round( UInt,sumv /numbercnt)
end

function AlarmCal( CumSumDataobj::spaceCumSumData,threshold= 4000)
    println("Nothing to do  in AlarmCal")
    return
end

########### 计算空间距离上的累加数据，并保存到对应的数组中。
function cumSumProc(spaceCumSumDataobj::spaceCumSumData ,featuredata::Array{UInt16,1} )
    temobj = spaceCumSumDataobj
    idx = getCurRowCnt(  temobj)
    va = @view temobj.m1[idx , : ]  #
    pieceData = @view featuredata[ temobj.startID : temobj.endID ]
    # println( "length(pieceData) ,length(va ) ，length(featuredata）: $(length(pieceData)) ,$(length(va )),$(length(featuredata) )" )
    va .=   pieceData
    va .= cumsum(  va )

    # for d in pieceData; print(d,"," );end
    # println("Cumsum,")
    # for d in va; print(d,"," );end ;println()
end

function setDTValue(cumObj::spaceCumSumData,dt  )
    cumObj.dt =  dt
end

function  postProcFeatureDataRT(spaceCumSumDataobj::spaceCumSumData ,data::Array{UInt16,1} )
    cumSumProc( spaceCumSumDataobj,  data  )
    # println("     ????                               ")
    vsc = spaceCumSumDataobj.defaultSpaceCnt
    vTcNum =spaceCumSumDataobj.defaultTimeCnt
    ######### calRectSum( CumSumDataobj::spaceCumSumData, vPos ,vsc ,vTcNum   )
    for vpos = getSpaceLoopRng( spaceCumSumDataobj );
        # println("vpos ,vsc ,vTcNum :$vpos ,$vsc ,$vTcNum   "  )
        calRectSum( spaceCumSumDataobj, vpos ,vsc ,vTcNum   )
    end
    # println()
    AlarmCal(spaceCumSumDataobj)
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# include(joinpath(dirname(@__FILE__), "datafiles.jl"))
#
# if !isdir(datadir)
#   mkdir(datadir)
# end
#
# info("Downloading test datasets:")
#
# for f in datafiles
#   fn = joinpath(datadir, f)
#   if !isfile(fn)
#     file_url = publicexamples*f*"?raw=true"
#     info("Downloading $(file_url)")
#     download(file_url, fn)
#   end
# end
#
# info("Download complete.")


text:
# This file includes the MonteCarlo simulation of FSK modulation scheme and compares
# the numerical results with the theoretical results.

using DigitalCommunications
using Plots

# Settings
pulse = RectangularPulse()          # Modulation pulse
fs = 10                             # Samling frequency
ts = 1 / fs                         # Sampling period
tb = pulse.duration                 # Pulse duration
k = 3                               # Bits per symbol
M = 2^k                             # Costellation size
nsymbols = Int(1e6)                 # Number of symbols
ebno = collect(-2 : 8)              # Energy per bit to noise power spectral density ratio
esno = ebno .+ 10 * log10(k)        # Energy per symbol to noise power spectral densit ratio.

# Communcation system components
gen = SymbolGenerator(nsymbols, M)
modulator = Modulator(FSK(M), RectangularPulse(), ts)
channel = AWGNChannel(1., ts, tb)
detector = Detector(modulator(1:M))

# Monte Carlo simulation
message = gen.symbols
symerr = zeros(length(esno))
for i in 1 : length(symerr)
    channel.esno = esno[i]  # Update channel noise level
    mbar = message |> modulator |> channel |> detector  # Run communication system
    symerr[i] = sum(mbar .!= message) / length(message)  # Symbol error rate
end

# Plots
plt = plot(title="$M-PSK", xlabel="ebno [dB]", ylabel="Pe")
plot!(ebno, berfsk.(esno, M), marker=:circle, yscale=:log10, label="theoretical")
plot!(ebno, symerr, marker=:circle, yscale=:log10, label="montecarlo")


text:
"""

$(SIGNATURES)

Helper infrastructure to compile and sample models using `cmdstan`.
"""
module StanDiagnose

using StanBase
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF

import StanBase: stan_sample, get_cmdstan_home
import StanBase: cmdline, stan_summary, read_summary
import StanBase: RandomSeed, Init, Output, StanModelError

include("stanmodel/diagnose_types.jl")
include("stanmodel/DiagnoseModel.jl")
include("stanrun/cmdline.jl")
include("stansamples/read_diagnose.jl")

stan_diagnose = stan_sample

export
  DiagnoseModel,
  StanModelError,
  stan_diagnose,
  read_diagnose,
  read_summary,
  stan_summary

end # module


text:
using Test, ToggleableAsserts

function foo(u, v)
    @toggled_assert length(u) == length(v)
    1
end

@test_throws AssertionError foo([1, 2], [1])
@test_throws AssertionError("3 is an odd number!") @toggled_assert iseven(3) "3 is an odd number!"

toggle(false) # toggle off assertions

@test foo([1, 2], [1]) == 1

Threads.@threads for i in 1:20
    rand(Bool) ? toggle(true) : toggle(false)
end


text:
"""
relaxedSequencingDD(
    model::SequencingModel,
    maxWidth::U,
    splitOrder::Vector{U};

    dd::Union{Vector{Vector{RelaxedSequencingNode}},Vector{Vector{RelaxedStateSequencingNode}},Nothing}=nothing,
    startingRoot::Union{RelaxedSequencingNode,RelaxedStateSequencingNode,Nothing}=nothing,
    bestKnownSolution::Union{Vector{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, V}=nothing,
    timeElapsed::W=0.0,
    startTime::X=time(),
    widthOne::Bool=false,
    debugOn::Bool=false
)where{T<:Real,U<:Int}

Generate a relaxed bound for a sequencing problem using a relaxed decision diagram
Returns dd (Vector{Vector{RelaxedSequencingNode}}), bound (Real), relaxedSol (Vector{Int}), bestKnownSolution (Vector{Int}), bestKnownValue (Real), exact (Bool)
'exact' is true if the diagram exactly represents the problem (or subproblem) it was given

# Arguments
- `model::SequencingModel`: The sequencing problem
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `splitOrder::Vector{U}`: A heuristic ordering of nodes by importance
- `dd::Union{Vector{Vector{RelaxedSequencingNode}},Vector{Vector{RelaxedStateSequencingNode}},Nothing}`: An optional diagram to start from
- `startingRoot::Union{RelaxedSequencingNode,RelaxedStateSequencingNode,Nothing}`: An optional root node to start from
- `bestKnownSolution::Union{Vector{Int},Nothing}`: An optional parameter that gives the DD a starting solution
- `bestKnownValue::Union{T,Nothing}`: An optional parameter, The value of the best known feasible solution
- `rrbSetting::Union{RRBSetting,Nothing}`: An optional parameter that chooses an rough relaxed bounding method
- `rrbPacket`: Optional data to be passed to the rrbFunction
- `timeLimit::Union{Nothing, V}`: An optional parameter that adds a time limit
- `timeElapsed::W`: An optional parameter denoting how much time passed already
- `startTime::X`: An optional parameter denoting when the countdown started if there is a timer
- `widthOne::Bool`: An optional parameter to indicate that a widdth one diagram is the starting diagram
- `rootDomain::BitVector`: the domain of the rootNode if using widthOne
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
"""
function relaxedSequencingDD(
    model::SequencingModel,
    maxWidth::U,
    splitOrder::Vector{U};
    dd::Union{Vector{Vector{RelaxedSequencingNode}},Vector{Vector{RelaxedStateSequencingNode}},Nothing}=nothing,
    startingRoot::Union{RelaxedSequencingNode,RelaxedStateSequencingNode,Nothing}=nothing,
    bestKnownSolution::Union{Vector{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, V}=nothing,
    timeElapsed::W=0.0,
    startTime::X=time(),
    widthOne::Bool=false,
    rootDomain::BitVector=BitVector(),
    debugOn::Bool=false
)where{T<:Real,U<:Int, V<:Real, W<:Real, X<:Real}
    timedOut() = isTimedOut(timeLimit,startTime, timeElapsed)
    #create initial DD from initial constraints and domains
    if isnothing(dd)
        if isnothing(startingRoot)
            if widthOne
                dd = createWidthOneDDFromNode(model::SequencingModel)
            else
                dd = createSquareStartDD(model)
            end
        else
            if widthOne
                dd = createWidthOneDDFromNode(model,rootNode=startingRoot,masterDomain=rootDomain)
            else
                dd = createSquareStartDD(model, rootNode=startingRoot)
            end
        end
    end
    #for each layer
    for i in 1:length(dd)-1
        #only run this loop if the dd isnt too wide (or its the exact layer)
        if length(dd[i])<maxWidth || i==2
            #for each job in the ordering
            for assignment in splitOrder
                queue = selectNodes(dd[i],assignment)
                for node in queue
                    #split node into two new nodes
                    newNode = splitNode!(node, assignment, model)

                    #trim useless arcs
                    filterArcs!(node.arcsOut,model;bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)
                    filterArcs!(newNode.arcsOut,model;bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)

                    #if newNode has remaining out arcs, add it to the dd
                    if !isempty(newNode.arcsOut)
                        push!(dd[i], newNode)
                    else #delete in arcs so parent nodes may be removed in subsequent passes
                        deleteArcs!(newNode.arcsIn)
                    end
                    #if old node has no out arcs, remove it from the dd
                    if isempty(node.arcsOut)
                        deleteArcs!(node.arcsIn)
                        filter!(x->x!=node, dd[i])
                    end
                    #continue until max width of layer is reached except on the first layer
                    if length(dd[i])>=maxWidth && i!= 2
                        break
                    end
                    #timeout check
                    if timedOut()
                        return Vector{Vector{RelaxedSequencingNode}}(), 0, bestKnownSolution,bestKnownSolution,bestKnownValue, false
                    end
                end#queue loop
                #continue until max width of layer is reached except on the first layer
                if length(dd[i])>=maxWidth && i!= 2
                    break
                end
                #timeout check
                if timedOut()
                    return Vector{Vector{RelaxedSequencingNode}}(), 0, bestKnownSolution,bestKnownSolution,bestKnownValue, false
                end
            end#end assignment iteration
        end#if statement

        #timeout check
        if timedOut()
            return Vector{Vector{RelaxedSequencingNode}}(), 0, bestKnownSolution,bestKnownSolution,bestKnownValue, false
        end
        if i<length(dd)
            for node in dd[i+1]
                if isempty(node.arcsIn)
                    deleteArcs!(node.arcsOut)
                end
            end
            filter!(x->length(x.arcsIn)!=0, dd[i+1])
        end
        if dd isa Vector{Vector{RelaxedStateSequencingNode}}
            for node in dd[i]
                updateOutgoingArcValues!(model, node)
            end
        end
        #prep the next layer
        for node in dd[i+1]
            updateNodeDownVariables!(node,model)
            if node isa RelaxedStateSequencingNode
                updateOutgoingArcValues!(model, node)
            end
            filterArcs!(node.arcsOut,model;bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)

            if isempty(node.arcsOut) && i+1 != length(dd)
                deleteArcs!(node.arcsIn)
            end
        end
        filter!(x->length(x.arcsIn)!=0, dd[i+1])
        #timeout check
        if timedOut()
            return Vector{Vector{RelaxedSequencingNode}}(), 0, bestKnownSolution,bestKnownSolution,bestKnownValue, false
        end
    end#end layer iteration

    #check if empty
    if isempty(last(dd))
        return Vector{Vector{RelaxedSequencingNode}}(), bestKnownValue, bestKnownSolution,bestKnownSolution,bestKnownValue, true
    end

    #timeout check
    if timedOut()
        return Vector{Vector{RelaxedSequencingNode}}(), 0, bestKnownSolution,bestKnownSolution,bestKnownValue, false
    end

    #iteratively delete nodes going up
    bottomUpUpdate(model, dd, bestKnownValue)

    if widthOne
        return dd, last(dd)[1].lengthToRoot, [], bestKnownSolution, bestKnownValue, last(dd)[1].exact
    else
        #update best value
        bestKnownSolution, bestKnownValue = checkFeasibleTerminalSolutions!(model, last(dd)[1],bestKnownValue=bestKnownValue,bestKnownSolution=bestKnownSolution)
        return dd, last(dd)[1].lengthToRoot, last(dd)[1].bestSeq, bestKnownSolution, bestKnownValue, last(dd)[1].exact
    end
end

"""
    getFrontierNode(dd::Vector{Vector{RelaxedSequencingNode}}, relaxedSolution::Vector{T})where{T<:Int}

Get the deepest layer indexed exact node from the best path through the diagram

# Arguments
- `dd::Vector{Vector{RelaxedSequencingNode}}`: The dd to search
- `relaxedSolution::Vector{T}`: The problem solution
"""
function getFrontierNode(dd::Vector{Vector{RelaxedSequencingNode}}, relaxedSolution::Vector{T})where{T<:Int}
    currentNode = dd[1][1]
    for i in length(currentNode):length(relaxedSolution)
        for arc in currentNode.arcsOut
            if relaxedSolution[i+1] == arc.assignment
                if arc.destination.exact
                    currentNode = arc.destination
                    break
                else
                    return currentNode, (i - length(dd[1][1]) + 1)
                end
            end
        end
    end
end

"""
    getLastExactNode(dd::Vector{Vector{RelaxedSequencingNode}}, relaxedSolution::Vector{T})where{T<:Int}

Get the first node from the best path through the diagram that is a parent to a non-exact node

# Arguments
- `dd::Vector{Vector{RelaxedSequencingNode}}`: The dd to search
- `relaxedSolution::Vector{T}`: The problem solution
"""
function getLastExactNode(dd::Vector{Vector{RelaxedSequencingNode}}, relaxedSolution::Vector{T})where{T<:Int}
    currentNode = dd[1][1]
    for i in length(currentNode):length(relaxedSolution)
        nextNode = nothing
        for arc in currentNode.arcsOut
            if relaxedSolution[i+1] == arc.assignment
                if arc.destination.exact
                    nextNode = arc.destination
                else
                    return currentNode, (i - length(dd[1][1]) + 1)
                end
            elseif !arc.destination.exact
                return currentNode, (i - length(dd[1][1]) + 1)
            end
        end
        currentNode = nextNode
    end
end

"""
    getLastExactLayer(dd::Vector{Vector{RelaxedStateSequencingNode}})

Get the last exact layer of a dd

# Arguments
- `dd::Vector{Vector{RelaxedStateSequencingNode}}`: The dd to search
"""
function getLastExactLayer(dd::Vector{Vector{RelaxedStateSequencingNode}})
    currentLayer = dd[1]
    for layer in dd
        for node in layer
            if !node.exact
                return currentLayer
            end
        end
        currentLayer = layer
    end
end


text:
module NonconvexNLopt

export NLoptAlg, NLoptOptions

using Reexport, Parameters, SparseArrays, Zygote
@reexport using NonconvexCore
using NonconvexCore: @params, VecModel, AbstractResult
using NonconvexCore: AbstractOptimizer, CountingFunction
import NonconvexCore: optimize!, Workspace
import NLopt

include("nlopt.jl")

end


text:
#-----------------------------------------------------------------------# ModelingType
abstract type ModelingType{T} <: OnlineStat{T} end
nobs(o::ModelingType) = nobs(o.stat)
_fit!(o::ModelingType, y) = _fit!(o.stat, y)

#### Numerical
struct Numerical <: ModelingType{Number}
    stat::Variance
end
Numerical() = Numerical(Variance())
function Base.show(io::IO, o::Numerical)
    print(io, "📈  Numerical: μ=", mean(o), ", σ=", std(o))
end
Statistics.mean(o::Numerical) = mean(o.stat)
Statistics.std(o::Numerical) = std(o.stat)
width(o::Numerical) = 1
transform(o::Numerical, xi, j) = (xi - mean(o)) / std(o)

#### Categorical
struct Categorical{T} <: ModelingType{T}
    stat::CountMap{T, OrderedDict{T, Int}}
end
Categorical(T) = Categorical(CountMap(T))
function Base.show(io::IO, o::Categorical{T}) where {T}
    print(io, "📊  Categorical{$T}: nclasses = ", length(o.stat.value))
end
width(o::Categorical) = max(0, length(o.stat.value) - 1)
function transform(o::Categorical, xi, j)
    out = 0
    for (k, ky) in enumerate(keys(o))
        j == k && xi == ky && (out = 1)
    end
    out
end

#-----------------------------------------------------------------------# DataPreprocessor
# S = standardize?
# G = Group of Numerical/Categorical
struct DataPreprocessor{S, G} <: OnlineStat{VectorOb}
    group::G
    x::Vector{Float64}
end
function DataPreprocessor(g::Group; standardize::Bool=true)
    DataPreprocessor{standardize, typeof(g)}(g, zeros(sum(width, g.stats)))
end
function DataPreprocessor(row; kw...)
    stats = []
    for item in row
        push!(stats, make_modeling_type(item))
    end
    DataPreprocessor(Group(stats...); kw...)
end
nobs(o::DataPreprocessor) = nobs(o.group)
Base.show(io::IO, o::DataPreprocessor) = print(io, "DataPreprocessor:\n", o.group)

make_modeling_type(x::Number) = Numerical()
make_modeling_type(x::T) where {T<:Union{AbstractString, Bool, Char}} = Categorical(T)

_fit!(o::DataPreprocessor, y) = _fit!(o.group, y)

function preprocess(itr, hints::Pair...; kw...)
    row = first(itr)
    p = Pair.(_keys(row), make_modeling_type.(values(row)))
    d = OrderedDict{Any, Any}(p...)
    for (k,v) in hints
        d[k] = v
    end
    g = fit!(Group(collect(values(d))...), itr)
    DataPreprocessor(g; kw...)
end
_keys(o) = keys(o)
_keys(o::Tuple) = 1:length(o)

# function width(o::DataPreprocessor, lhs::Expr)
#     out = 0
#     for ex in lhs.args[2:end]
#         if ex isa Symbol

#         end
#     end
# end

# #-----------------------------------------------------------------------# Terms
# abstract type ModelTerm end

# struct SingleTerm{T <: ModelingType} <: ModelTerm
#     a::T
# end

# struct Interaction{A<:ModelingType, B <: ModelingType} <: ModelTerm
#     a::A
#     b::B
# end
# Interaction(a::Numerical, b::Categorical) = Interaction(b, a)

# #-----------------------------------------------------------------------# Formula
# struct ModelFormula{L,R,D<:DataPreprocessor}
#     lhs
#     rhs::Vector{ModelTerm}
# end


text:
using Documenter
using BlockOpt

makedocs(
    sitename = "BlockOpt",
    format = Documenter.HTML(
        # See https://github.com/JuliaDocs/Documenter.jl/issues/868
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.MathJax2(),
        collapselevel = 1,
    ),
    modules = [BlockOpt],
    pages = [
        "Introduction" => "index.md",
        "Overview" => "overview/design.md",
        "Manual" => [
            "manual/model.md",
            "manual/driver.md",
            "manual/options.md",
            "manual/simulation.md",
        ],
        "Tutorials" => ["tutorials/installation.md", "tutorials/simple.md"],
    ],
)

deploydocs(
    devbranch = "main",
    push_preview = true,
    repo = "github.com/danphenderson/BlockOpt.jl.git",
)


text:
import Pkg
println("# Activate"); flush(stdout)
Pkg.activate(".")
println("# using Test"); flush(stdout)
using Test
println("# using CuArrays"); flush(stdout)
using CuArrays

println("# Try cufill"); flush(stdout)
N = 128
x_d = cufill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = cufill(2.0f0, N)  # a vector stored on the GPU filled with 2.0

println("# Try y_d .+= x_d"); flush(stdout)
y_d .+= x_d

println("# Test"); flush(stdout)
@test all(Array(y_d) .== 3.0f0)
println("# Completed"); flush(stdout)


text:
module TopicModelsVB

using SpecialFunctions
using Distributions
using OpenCL

export VectorList, MatrixList
export Document, Corpus
export TopicModel, GPUTopicModel, BaseTopicModel
export AbstractLDA, AbstractfLDA, AbstractCTM, AbstractfCTM, AbstractDTM, AbstractCTPF
export LDA, fLDA, CTM, fCTM, CTPF, DTM
export gpuLDA, gpuCTM, gpuCTPF
export isnegative, ispositive, logsumexp, addlogistic, partition
export checkdoc, checkcorp, readcorp, writecorp, abridgecorp!, trimcorp!, compactcorp!, padcorp!, cullcorp!, fixcorp!, showdocs, getlex, getusers
export train!, @gpu
export fixmodel!, gendoc, gencorp, showtopics, showlibs, showdrecs, showurecs

include("macros.jl")
include("utils.jl")
include("Corpus.jl")

include("TopicModel.jl")
include("GPUTopicModel.jl")

include("LDA.jl")
include("gpuLDA.jl")
include("AbstractLDA.jl")

include("fLDA.jl")
include("AbstractfLDA.jl")

include("CTM.jl")
include("gpuCTM.jl")
include("AbstractCTM.jl")

include("fCTM.jl")
include("AbstractfCTM.jl")

include("BaseTopicModel.jl")

include("DTM.jl")
include("AbstractDTM.jl")

include("CTPF.jl")
include("gpuCTPF.jl")
include("AbstractCTPF.jl")

include("modelutils.jl")

end


text:
#i64o-bitpatterns.jl

#bitpattern detection of unum integer numbers.  Certain patterns of numbers are
#well suited for detection for certain fraction tasks.  In particular, all-zeros
#is necessary for deciding if number is an exact power of two or zero.  And
#is_top helps detect the subnormal one.

doc"""
`is_all_zero` outputs whether or not all of the bits in an UInt64 or an `ArrayNum`
are zero.  Used for checking zero at the Unum level.
"""
is_all_zero(n::UInt64) = (n == 0)
function is_all_zero{FSS}(n::ArrayNum{FSS})
  for idx = 1:__cell_length(FSS)
    @inbounds (n.a[idx] != 0) && return false   #accumulate bits.
  end
  return true
end

doc"""
`is_not_zero` outputs whether or not any of the bits in an UInt64 or an `ArrayNum`
are one.  Used for checking not_zero at the Unum level.
"""
is_not_zero(n::UInt64) = (n != 0)
function is_not_zero{FSS}(n::ArrayNum{FSS})
  for idx = 1:__cell_length(FSS)
    @inbounds (n.a[idx] != 0) && return true
  end
  return false
end

doc"""
`is_all_ones(n)` outputs whether or not all of the bits in an UInt64 or an `ArrayNum`
are one.
`is_all_ones(n, fsize)` does the same, but constrained to a certain fsize.
"""
function is_all_ones{FSS}(n::ArrayNum{FSS})
  for idx = 1:__cell_length(FSS)
    @inbounds (n.a[idx] != f64) && return false
  end
  return true
end
function is_all_ones(n::UInt64, fsize::UInt16)
  n == mask_top(fsize)
end
function is_all_ones{FSS}(n::ArrayNum{FSS}, fsize::UInt16)
  cell_index = (fsize ÷ 0x0040) + o16
  bit_index = fsize % 0x0040
  for idx = 1:cell_index - 1
    @inbounds (n.a[idx] != f64) && return false
  end
  @inbounds return is_all_ones(n.a[cell_index], bit_index)
end

doc"""
`is_not_ones` outputs whether or not any of the bits in an UInt64 or an `ArrayNum`
is zero
"""
function is_not_ones{FSS}(n::ArrayNum{FSS})
  for idx = 1:__cell_length(FSS)
    @inbounds (n.a[idx] != f64) && return true
  end
  return false
end

doc"""
`is_top` outputs whether or not the most significant bit of an UInt64 or an
`ArrayNum` is a one, while the rest are zero.  Used for checking the subnormal
one at the Unum level.
"""
is_top(a::UInt64) = (a == t64)
function is_top{FSS}(n::ArrayNum{FSS})
  (n.a[1] != t64) && return false
  for idx = 2:__cell_length(FSS)
    @inbounds (n.a[idx] != 0) && return false
  end
  return true
end

doc"""
`top_bit` outputs whether or not the top bit is one.
"""
top_bit(a::UInt64) = (a & t64) != z64
top_bit{FSS}(a::ArrayNum{FSS}) = @inbounds top_bit(a.a[1])
@universal frac_top_bit(a::Unum) = top_bit(a.fraction)

doc"""
`is_not_top` outputs if the most significant bit of an UInt64 or an
`ArrayNum` isn't one, or if any of the other bits are one.  Used for checking
the subnormal one at the Unum level.
"""
is_not_top(a::UInt64) = (a != t64)
function is_not_top{FSS}(n::ArrayNum{FSS})
  (n.a[1] != t64) && return true
  for idx = 2:__cell_length(FSS)
    @inbounds (n.a[idx] != 0) && return true
  end
  return false
end

doc"""
`is_mmr_frac(::ArrayNum)` checks if the fraction looks like mmr.
`is_mmr_frac(::UInt64, ::Type{Val{FSS}})` also checks if the fraction looks like mmr.
"""
function is_mmr_frac{FSS}(n::UInt64, ::Type{Val{FSS}})
  if FSS == 0
    n == z64
  else
    n == mask_top(max_fsize(FSS) - 0x0001)
  end
end
function is_mmr_frac{FSS}(n::ArrayNum{FSS})
  l = __cell_length(FSS)
  for idx = 1:(l - 1)
    @inbounds (n.a[idx] != f64) && return false
  end
  @inbounds return (n.a[l] == 0xFFFF_FFFF_FFFF_FFFE)
end

bool_bottom_bit{FSS}(fraction::UInt64, ::Type{Val{FSS}}) = (bottom_bit(FSS) & fraction) != 0
bool_bottom_bit{FSS}(n::ArrayNum{FSS}) = (n.a[__cell_length(FSS)] & o64) != 0
bool_bottom_bit{ESS,FSS}(x::UnumSmall{ESS,FSS}) = bool_bottom_bit(x.fraction, Val{FSS})
bool_bottom_bit{ESS,FSS}(x::UnumLarge{ESS,FSS}) = bool_bottom_bit(x.fraction)

doc"""
  `Unums.bool_indexed_bit(fraction, index)`
"""
bool_indexed_bit(fraction::UInt64, index::UInt16) = ((t64 >> index) & fraction) != 0
function bool_indexed_bit{FSS}(fraction::ArrayNum{FSS}, index::UInt16)
  cell_index = (index ÷ 0x0040) + o16
  index = index % 0x0040
  bool_indexed_bit(fraction.a[cell_index], index)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Tokens
using Test

function testflytoken()

    t1 = T_TEXT"hello world"
    @test_throws BoundsError codeunit(t1,0)
    @test_throws BoundsError codeunit(t1,12)
    @test t1 == "hello world"
    @test category(t1) == T_TEXT
    t1a = DToken(T_TEXT,"hello")
    t1b = DToken(T_TEXT,"world")
    @test t1 == t1a * " " * t1b
    @test t1 > t1a
    @test t1a < t1b
    @test t1 < t1b

    t2 = D1"123456"
    @test_throws BoundsError codeunit(t2,7)
    @test_throws BoundsError codeunit(t2,0)
    @test t2 == "123456"
    @test category(t2) == T_INT
    s1 = t1 * t2
    @test s1 == "hello world123456"

    d1 = D1"direct."
    d2 = D2"direct."
    d1identical = D1"direct."
    @test usize(d1)==7
    @test d1==d2
    @test d1identical === d1
end