text:
# Level-2 LinearAlgebra.BLAS interface.
#
# TODO: support gbmv!, sbmv! and hbmv!.
#       hpmv! and spmv! is not possible within current frame.

# Level2 BLAS size checker.
bli_check_lv2(trans::BliTrans,
              m ::Integer, n ::Integer,
              n_::Integer, m_::Integer) = begin

    ((trans.enum & BLIS_TRANS_BIT) != 0) && ((m, n) = (n, m))
    (m == m_) || throw(DimensionMismatch("Target buffer size mismatch."))
    (n == n_) || throw(DimensionMismatch("Contracted size mismatch."))

    nothing
end

macro blis_interface_linalg_lv2_gemv(Tc1, T1, T2, Tc2, T3, targetfunc, bliname)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(tA::AbstractChar,
                           α::$Tc1,
                           A::StridedMatrix{<:$T1},
                           x::StridedVector{<:$T2},
                           β::$Tc2,
                           y::StridedVector{<:$T3}) = begin

            bli_tA = char_to_trans[tA]

            bli_check_lv2(bli_tA,
                          size(A)...,
                          length(x),
                          length(y))

            oα = BliObj(α)
            oA = BliObj(A)
            ox = BliObj(x)
            oβ = BliObj(β)
            oy = BliObj(y)

            ObjectBackend.bli_obj_set_onlytrans!(bli_tA, oA.obj)
            $blifunc(oα, oA, ox, oβ, oy)
            y

        end
    end
end

@doc """
        gemv!(tA, α, A, x, β, y)
    BLIS-based GEMV with strides support & mixed-precision.
    """ gemv!
@blis_interface_linalg_lv2_gemv(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                gemv!, gemv!)
@blis_interface_linalg_lv2_gemv Float32    Float32    Float32    Float32    Float32    gemv! gemv!
@blis_interface_linalg_lv2_gemv Float64    Float64    Float64    Float64    Float64    gemv! gemv!
@blis_interface_linalg_lv2_gemv ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 gemv! gemv!
@blis_interface_linalg_lv2_gemv ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 gemv! gemv!

macro blis_interface_linalg_lv2_hemv(Tc1, T1, T2, Tc2, T3, targetfunc, bliname, bli_struc)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(ul::AbstractChar,
                           α::$Tc1,
                           A::StridedMatrix{<:$T1},
                           x::StridedVector{<:$T2},
                           β::$Tc2,
                           y::StridedVector{<:$T3}) = begin

            bli_ul = char_to_uplo[ul]

            bli_check_lv2(BLIS_NO_TRANSPOSE,
                          size(A)...,
                          length(x),
                          length(y))

            oα = BliObj(α)
            oA = BliObj(A)
            ox = BliObj(x)
            oβ = BliObj(β)
            oy = BliObj(y)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oA)
            ObjectBackend.bli_obj_set_struc!($bli_struc, oA.obj)
            $blifunc(oα, oA, ox, oβ, oy)
            y

        end
    end
end

@doc """
        hemv!(ul, α, A, x, β, y)
    BLIS-based HEMV with strides support & mixed-precision.
    """ hemv!
@blis_interface_linalg_lv2_hemv(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                hemv!, hemv!,
                                BLIS_HERMITIAN)
@blis_interface_linalg_lv2_hemv Float32    Float32    Float32    Float32    Float32    hemv! hemv! BLIS_HERMITIAN
@blis_interface_linalg_lv2_hemv Float64    Float64    Float64    Float64    Float64    hemv! hemv! BLIS_HERMITIAN
@blis_interface_linalg_lv2_hemv ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 hemv! hemv! BLIS_HERMITIAN
@blis_interface_linalg_lv2_hemv ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 hemv! hemv! BLIS_HERMITIAN

@doc """
        symv!(ul, α, A, x, β, y)
    BLIS-based SYMV with strides support & mixed-precision.
    """ symv!
@blis_interface_linalg_lv2_hemv(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                symv!, symv!,
                                BLIS_SYMMETRIC)
@blis_interface_linalg_lv2_hemv Float32    Float32    Float32    Float32    Float32    symv! symv! BLIS_SYMMETRIC
@blis_interface_linalg_lv2_hemv Float64    Float64    Float64    Float64    Float64    symv! symv! BLIS_SYMMETRIC
@blis_interface_linalg_lv2_hemv ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 symv! symv! BLIS_SYMMETRIC
@blis_interface_linalg_lv2_hemv ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 symv! symv! BLIS_SYMMETRIC

macro blis_interface_linalg_lv2_trmv(T1, T2, targetfunc, bliname)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(ul::AbstractChar,
                           tA::AbstractChar,
                           dA::AbstractChar,
                           A::StridedMatrix{<:$T1},
                           b::StridedVector{<:$T2}) = begin

            bli_ul = char_to_uplo[ul]
            bli_dA = char_to_diag[dA]
            bli_tA = char_to_trans[tA]

            bli_check_lv2(bli_tA,
                          size(A)...,
                          length(b),
                          length(b))

            oα = BliObj($T1(1.0))
            oA = BliObj(A)
            ob = BliObj(b)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oA)
            ObjectBackend.bli_obj_set_diag!(bli_dA, oA)
            ObjectBackend.bli_obj_set_onlytrans!(bli_tA, oA)
            $blifunc(oα, oA, ob)
            b

        end
    end
end

@doc """
        trmv!(ul, tA, dA, A, b)
    BLIS-based TRMV with strides support & mixed-precision.
    """
@blis_interface_linalg_lv2_trmv(BliCompatibleType,
                                BliCompatibleType,
                                trmv!, trmv!)
@blis_interface_linalg_lv2_trmv Float32    Float32    trmv! trmv!
@blis_interface_linalg_lv2_trmv Float64    Float64    trmv! trmv!
@blis_interface_linalg_lv2_trmv ComplexF32 ComplexF32 trmv! trmv!
@blis_interface_linalg_lv2_trmv ComplexF64 ComplexF64 trmv! trmv!

@doc """
        trsv!(ul, tA, dA, A, b)
    BLIS-based TRSV with strides support & mixed-precision.
    """
@blis_interface_linalg_lv2_trmv(BliCompatibleType,
                                BliCompatibleType,
                                trsv!, trsv!)
@blis_interface_linalg_lv2_trmv Float32    Float32    trsv! trsv!
@blis_interface_linalg_lv2_trmv Float64    Float64    trsv! trsv!
@blis_interface_linalg_lv2_trmv ComplexF32 ComplexF32 trsv! trsv!
@blis_interface_linalg_lv2_trmv ComplexF64 ComplexF64 trsv! trsv!

macro blis_interface_linalg_lv2_ger(Tc1, T1, T2, T3, targetfunc, bliname)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(α::$Tc1,
                           x::StridedVector{<:$T1},
                           y::StridedVector{<:$T2},
                           A::StridedMatrix{<:$T3}) = begin

            bli_check_lv2(BLIS_NO_TRANSPOSE,
                          size(A)...,
                          length(y),
                          length(x))

            oα = BliObj(α)
            ox = BliObj(x)
            oy = BliObj(y)
            oA = BliObj(A)

            $blifunc(oα, ox, oy, oA)
            A

        end
    end
end

@doc """
        ger!(α, x, y, A)
    BLIS-based GER with strides support & mixed-precision.
    """ ger!
@blis_interface_linalg_lv2_ger(BliCompatibleType,
                               BliCompatibleType,
                               BliCompatibleType,
                               BliCompatibleType,
                               ger!, ger!)
@blis_interface_linalg_lv2_ger Float32    Float32    Float32    Float32    ger! ger!
@blis_interface_linalg_lv2_ger Float64    Float64    Float64    Float64    ger! ger!
@blis_interface_linalg_lv2_ger ComplexF32 ComplexF32 ComplexF32 ComplexF32 ger! ger!
@blis_interface_linalg_lv2_ger ComplexF64 ComplexF64 ComplexF64 ComplexF64 ger! ger!

macro blis_interface_linalg_lv2_her(Tc1, T1, T2, targetfunc, bliname, bli_struc)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(ul::AbstractChar,
                           α::$Tc1,
                           x::StridedVector{<:$T1},
                           A::StridedMatrix{<:$T2}) = begin

            bli_ul = char_to_uplo[ul]

            bli_check_lv2(BLIS_NO_TRANSPOSE,
                          size(A)...,
                          length(x),
                          length(x))

            oα = BliObj(α)
            ox = BliObj(x)
            oA = BliObj(A)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oA)
            ObjectBackend.bli_obj_set_struc!($bli_struc, oA.obj)
            $blifunc(oα, ox, oA)
            A

        end
    end
end

@doc """
        her!(uplo, α, x, A)
    BLIS-based HER with strides support & mixed-precision.
    """ her!
@blis_interface_linalg_lv2_her(BliCompatibleType,
                               BliCompatibleType,
                               BliCompatibleType,
                               her!, her!,
                               BLIS_HERMITIAN)
@blis_interface_linalg_lv2_her Float32    Float32    Float32    her! her! BLIS_HERMITIAN
@blis_interface_linalg_lv2_her Float64    Float64    Float64    her! her! BLIS_HERMITIAN
@blis_interface_linalg_lv2_her Float32    ComplexF32 ComplexF32 her! her! BLIS_HERMITIAN
@blis_interface_linalg_lv2_her Float64    ComplexF64 ComplexF64 her! her! BLIS_HERMITIAN

@doc """
        syr!(uplo, α, x, A)
    BLIS-based SYR with strides support & mixed-precision.
    """ syr!
@blis_interface_linalg_lv2_her(BliCompatibleType,
                               BliCompatibleType,
                               BliCompatibleType,
                               syr!, syr!,
                               BLIS_SYMMETRIC)
@blis_interface_linalg_lv2_her Float32    Float32    Float32    syr! syr! BLIS_SYMMETRIC
@blis_interface_linalg_lv2_her Float64    Float64    Float64    syr! syr! BLIS_SYMMETRIC
@blis_interface_linalg_lv2_her ComplexF32 ComplexF32 ComplexF32 syr! syr! BLIS_SYMMETRIC
@blis_interface_linalg_lv2_her ComplexF64 ComplexF64 ComplexF64 syr! syr! BLIS_SYMMETRIC


text:
using Caesar, IncrementalInference, RoME, HDF5, JLD
using CloudGraphs, Neo4j

# # Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true)
# session = addrdict["session"]

# interactive operation
# session = "SESSROX"
# Nparticles = 100
# include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

# get all the data

d = jldopen("test/data/VicPrkBasic.jld", "r") do file
  read(file, "dBasic")
end

f = jldopen("test/data/VicPrkBasic.jld", "r") do file
  read(file, "fBasic")
end
# MM = jldopen("test/data/VicPrkBasic.jld", "r") do file
#   read(file, "MMBasic")
# end
include("test/data/Vic120MM.jl")
MMr = MMwithRedirects(MM);

isamdict = jldopen("test/data/fgWithISAMrefERR01.jld", "r") do file
  read(file, "isamdict")
end
fgu = jldopen("test/data/fgWithISAMrefERR01.jld", "r") do file
  read(file, "fgu")
end

# do inference

T=100
fg = emptyFactorGraph();
fg.cg = cloudGraph

idx = appendFactorGraph!(fg, d, f, toT=T, lcmode=:mmodal, MM=MMr);

# tree = prepBatchTree!(fg, ordering=:qr,drawpdf=true);
# inferOverTree!(fg,tree, N=100);
#
# tree = wipeBuildNewTree!(fg)
# removeGenericMarginals!(conn)
# inferOverTree!(fg, tree)

ext: jl


text:
using LinearAlgebra

struct TimeGF <: AbstractArray{ComplexF64, 2}
  data::Array{ComplexF64, 2}
  grid::TimeGrid
end

function TimeGF(grid::TimeGrid)
  N = length(grid)
  TimeGF(zeros(ComplexF64, N, N), grid)
end

function TimeGF(f::Function, grid::TimeGrid; lower = false, time_invariant = false)
  if time_invariant
    δ = minimum(abs.(grid.step)) / 10

    make_key = (t1, t2) -> begin
      theta = θ(t1.val, t2.val)
      Δt = t1.val.val - t2.val.val
      return (Complex{Int}(round(Δt / δ)), theta)
    end

    T = typeof(f(grid[1], grid[1])) # extra evaluation to get type for cache

    cache = Dict{Tuple{Complex{Int}, Bool}, T}()

    g = (t1, t2) -> begin
      key = make_key(t1, t2)
      key ∈ keys(cache) ? cache[key] : cache[key] = f(t1, t2)
    end

    return TimeGF(g, grid, lower=lower, time_invariant=false)
  else
    N = length(grid)
    gf = TimeGF(grid)

    for t1 in grid
      for t2 in grid
        lower && t1.idx < t2.idx && continue
        gf[t1, t2] = f(t1, t2)
      end
    end
    return gf
  end
end

### AbstractArray Interface ###
IndexStyle(::Type{<:TimeGF}) = IndexLinear()
size(gf::TimeGF) = size(gf.data)
Base.@propagate_inbounds getindex(gf::TimeGF, i::Int) = gf.data[i]
Base.@propagate_inbounds setindex!(gf::TimeGF, v, i::Int) = gf.data[i] = v
function similar(gf::TimeGF, ::Type{S}) where S
  data = similar(gf.data, S)
  TimeGF(data, gf.grid) # TODO copy grid?
end

### BroadCasting ###
Base.BroadcastStyle(::Type{<:TimeGF}) = Broadcast.ArrayStyle{TimeGF}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{TimeGF}}, ::Type{ElType}) where ElType
    # Scan the inputs for first TimeGF
    gf = find_time_gf(bc)
    # copy grid to output
    TimeGF(similar(Array{ElType}, axes(bc)), gf.grid) # TODO copy grid?
end

"`A = find_time_gf(As)` returns the first TimeGF among the arguments."
find_time_gf(bc::Base.Broadcast.Broadcasted) = find_time_gf(bc.args)
find_time_gf(args::Tuple) = find_time_gf(find_time_gf(args[1]), Base.tail(args))
find_time_gf(x) = x
find_time_gf(a::TimeGF, rest) = a
find_time_gf(::Any, rest) = find_time_gf(rest)

# indexing with TimeGridPoint
Base.@propagate_inbounds function getindex(gf::TimeGF, t1::TimeGridPoint, t2::TimeGridPoint, gtr=true)
  val = gf[t1.idx, t2.idx]
  (!gtr && t1.idx == t2.idx) && (val += jump(gf))
  return val
end

Base.@propagate_inbounds setindex!(gf::TimeGF, v, t1::TimeGridPoint, t2::TimeGridPoint) = gf[t1.idx, t2.idx] = v

function jump(gf::TimeGF)
  t0_plus = branch_bounds(gf.grid, forward_branch)[1]
  t0_minus = branch_bounds(gf.grid, backward_branch)[2]
  return gf[t0_plus, t0_minus] - gf[t0_plus, t0_plus]
end

"""
Transposed view of gf taking into account diagonal discontinuity
"""
struct TimeGFTranspose <: AbstractArray{ComplexF64, 2}
  gf::TimeGF
end

adjoint(gf::TimeGF) = TimeGFTranspose(gf)

### AbstractArray Interface ###
IndexStyle(::Type{<:TimeGFTranspose}) = IndexCartesian()
size(gfa::TimeGFTranspose) = size(gfa.gf)
Base.@propagate_inbounds getindex(gfa::TimeGFTranspose, i::Int, j::Int) = i == j ? gfa.gf[i,i] + jump(gfa.gf) : gfa.gf[j,i]
Base.@propagate_inbounds setindex!(gfa::TimeGFTranspose, v, i::Int, j::Int) = gfa.gf[j,i] = v

# indexing with TimeGridPoint
Base.@propagate_inbounds getindex(gfa::TimeGFTranspose, t1::TimeGridPoint, t2::TimeGridPoint) = gfa[t1.idx, t2.idx]
Base.@propagate_inbounds setindex!(gfa::TimeGFTranspose, v::ComplexF64, t1::TimeGridPoint, t2::TimeGridPoint) = gfa[t1.idx, t2.idx] = v

function getindex(gf::TimeGF, b1::BranchEnum, b2::BranchEnum)
  grid = gf.grid
  @assert b1 ∈ grid.contour && b2 ∈ grid.contour
  x = b1 == backward_branch ? reverse(grid[b1]) : grid[b1]
  y = b2 == backward_branch ? reverse(grid[b2]) : grid[b2]
  return [gf[t1, t2] for t1 in x, t2 in y]
end

function getindex(gf::TimeGF, component::Symbol)
  grid = gf.grid
  if component == :greater
    return gf[backward_branch, forward_branch]
  elseif component == :lesser
    return gf[forward_branch, backward_branch]
  elseif component == :matsubara
    gf[imaginary_branch, imaginary_branch][:,1]
  elseif component == :retarded
    gfʳ = gf[:greater] - gf[:lesser]
    θ = LowerTriangular(ones(size(gfʳ)...))
    return θ .* gfʳ
  else
    throw(ArgumentError("component $component not recognized"))
  end
end

"""
    density(gf::TimeGF)

n(t) = <c†(t)c(t)> = -i G<(t,t)
"""
function density(gf::TimeGF)
  return -1.0im * diag(gf[:lesser])
end

"""
    equilibrium_spectrum(gf::TimeGF, ω)

A(ω) = Im -1/π ∫dt Gʳ(t, 0) exp(iωt)
"""
function equilibrium_spectrum(gf::TimeGF, ω)
  grid = gf.grid
  t = realtimes(grid)
  trapz = ft -> sum((ft[i] + ft[i+1]) * (t[i+1] - t[i]) / 2 for i in 1:(length(t) - 1))
  gfʳ = gf[:retarded][:,1]
  Aω = [imag((-1.0 / π) * trapz(gfʳ .* exp.(1.0im .* ωi .* t))) for ωi in ω]
end

function aux_spectrum(gf::TimeGF, ω::Number)
  grid = gf.grid

  delta_f = (t1, t2) -> -1.0im *(θ(t1.val, t2.val) - 1) * exp(-1.0im * (t1.val.val - t2.val.val) * ω)
  delta_e = (t1, t2) -> -1.0im *(θ(t1.val, t2.val) - 0) * exp(-1.0im * (t1.val.val - t2.val.val) * ω)

  At = map(zip(grid[forward_branch], reverse(grid[backward_branch]))) do (t⁺, t⁻)
    Iω_f = -2.0 * integrate(t -> @inbounds(gf[t⁺, t] * delta_f(t, t⁻)), grid, t⁺, t⁻)
    Iω_e = 2.0 * integrate(t -> @inbounds(delta_e(t⁺, t) * gf[t, t⁻]), grid, t⁺, t⁻)
    (1 / (2π)) * (Iω_e - Iω_f)
  end
end

"""
    aux_spectrum(gf::TimeGF, ω)

Compute auxiliary current spectral function of gf for frequencies ω
See Eq. 24 from [1]

[1] Cohen, Guy, David R. Reichman, Andrew J. Millis, and Emanuel Gull. “Green’s
Functions from Real-Time Bold-Line Monte Carlo.” Physical Review B 89, no. 11
(March 31, 2014): 115139. https://doi.org/10.1103/PhysRevB.89.115139.
"""
function aux_spectrum(gf::TimeGF, ω)
  reduce(vcat, (aux_spectrum(gf, ωi)' for ωi in ω))
end


text:
# LSAC JUNE 2007 Form 8LSN75
# https://www.lsac.org/sites/default/files/legacy/docs/default-source/jd-docs/sampleptjune.pdf

# Questions 6-10

"""
Exactly three films—Greed, Harvest, and Limelight—are
shown during a film club’s festival held on Thursday, Friday,
and Saturday. Each film is shown at least once during the
festival but never more than once on a given day. On each day
at least one film is shown. Films are shown one at a time. The
following conditions apply:
"""
julia> using AbstractLogic
Start the repl in command prompt by typing `=`.

# Three films Greed, Harvest, and Limelight as well as 0 no film
# t1 is the first film, t2 the second, and t3 are the last film on Thursday
# f1, f2, f3 and s1, s2, s3 are similarly coded for friday and saturday
t.1, t.2, t.3, f.1, f.2, f.3, s.1, s.2, s.3  ∈  _, Greed, Harvest, Limelight [clear]

# On each day at least one film is shown. Since 3 is last slot then .3 can't be empty.
{{j}}.3 != '_'

# Let's impose the sttructural contrainst that the highest number 3 is the
# last film shown. Thus if t2|t1 != 0 => t3 != 0
# t1 != 0 ==> t2 !=0", logicset)
{{j}}.1 != '_' ==> {{j}}.2 != '_'

# Each film is shown at least once during the festival
{{i}} == 'Greed' {{1,}}
{{i}} == 'Harvest' {{1,}}
{{i}} == 'Limelight' {{1,}}

# Never more than once on a given day. If t1 does not equal '_' then t1 != t2.
# Since we know that t3 can't be '_' then we can assert that t2 != t3
{{j}}.1 != '_' ==> {{j}}.1 != {{j}}.2 &&& {{j}}.1 != {{j}}.3 &&& {{j}}.2 != {{j}}.3

# On Thursday Harvest is shown, and no film is shown after it on that day
t.3 = 'Harvest'

# On Friday either Greed or Limelight, but not both, and no film is
# shown after it on that day.
f.3 = 'Greed'|'Limelight'; f.2,f.1 != 'Greed','Limelight'

# On Saturday either Greed or Harvest, but not both, is shown, and no film is
# shown after it on that day.
s.3 = 'Greed'|'Harvest'; s.1,s.2 != 'Greed','Harvest'

export lsat3

#6. Which one of the following could be a complete and accurate description of
# the order in which the films are shown at the festival?
#(A) Thursday: Limelight, then Harvest; Friday: Limelight; Saturday: Harvest
check t.2 = 'Limelight' ; t.3 = 'Harvest' ; f.3 = 'Limelight' ; s.3 = 'Harvest' ; t.1,f.2,f.1,s.2,s.1='_'
# false 0 out of 64 possible combinations true

#(B) Thursday: Harvest; Friday: Greed, then Limelight; Saturday: Limelight, then Greed
check t.3 = 'Harvest' ; f.2 = 'Greed' ; f.3 = 'Limelight' ; s.2 = 'Limelight' ; s.3 = 'Greed'
# false 0 out of 64 possible combinations true

#(C) Thursday: Harvest; Friday: Limelight; Saturday: Limelight, then Greed
check t.3 = 'Harvest' ; f.3 = 'Limelight' ; s.2 = 'Limelight' ; s.3 = 'Greed'; t.1, t.2, f.1, f.2, s.1='_'
# true 1 out of 64 possible combinations true

#(D) Thursday: Greed, then Harvest, then Limelight; Friday: Limelight; Saturday: Greed
check t.1 = 'Greed' ; t.2 = 'Harvest' ; t.3 = 'Limelight' ; f.3 = 'Limelight' ; s.3 = 'Greed'
# false, 0 out of 64 possible combinations 'true'.

#(E) Thursday: Greed, then Harvest; Friday: Limelight, then Harvest; Saturday: Harvest
check t.2 = 'Greed' ; t.3 = 'Harvest' ; f.2 = 'Limelight' ; f.3 = 'Harvest' ; s.3 = 'Harvest'
# false, 0 out of 64 possible combinations 'true'.

#7. Which one of the following CANNOT be true?
#(A) Harvest is the last film shown on each day of the festival.
check t.3,f.3,s.3 = 'Harvest'
# false, 0 out of 64 possible combinations 'true'

# (B) Limelight is shown on each day of the festival.
any {{i}} == 'Limelight' {{3,}}
# true,  10 out of 72 possible combinations 'true'.

#(C) Greed is shown second on each day of the festival.
any f.1,f.2,f.3 |= 'Greed'
# true,  32 out of 72 possible combinations 'true'

# (D) A different film is shown first on each day of the festival.
## This will be a little trickier. Define .0 as first film

t.0, f.0, s.0 ∈ Greed, Harvest, Limelight
{{j}}.0 == {{j}}.1 if {{j}}.1 != '_'
{{j}}.1 == '_' &&& {{j}}.2 != '_' ==> {{j}}.0 == {{j}}.2
{{j}}.2 == '_' ==> {{j}}.0 == {{j}}.3
check t.0 != f.0, s.0 ; f.0 != s.0

import lsat3

# possible,  19 out of 64 possible combinations 'true'.

# (E) A different film is shown last on each day of the festival.
check t.3, f.3 != s.3; t.3 != f.3
# possible,  24 out of 74 possible combinations 'true'.

# If Limelight is never shown again during the festival once Greed is shown,
{{i}} = 'Greed' ==> {{>i}} != 'Limelight'

# which one of the following is the maximum number of film showings that could
# occur during the festival?
#(A) three
any: {{i}} != '_' {{3}}
# true

#(B) four
any: {{i}} != '_' {{4}}
# true

#(C) five
any: {{i}} != '_' {{5}}
# true

# (D) six
any: {{i}} != '_' {{6}}
# true
# * Answer

# (E) seven
any: {{i}} != '_' {{7}}
# false

# If Greed is shown exactly three times, Harvest is shown exactly twice, and
# Limelight is shown exactly once, then which one of the following must be true?
import lsat3

{{i}} = 'Greed' {{3}}; {{i}} = 'Harvest' {{2}}; {{i}} = 'Limelight' {{1}}

#(A) All three films are shown on Thursday.
prove: t.1 != '_'
# false, 2 out of 3 possible combinations 'true'.

#(B) Exactly two films are shown on Saturday.
prove: s.1 = '_'
# true, 3 out of 3 possible combinations 'true'.

#(C) Limelight and Harvest are both shown on Thursday.
prove: t.1,t.2,t.3 |= 'Limelight'; t.1,t.2,t.3 |= 'Harvest'
# false, 2 out of 3 possible combinations 'true'.

#(D) Greed is the only film shown on Saturday.
prove: s.1,s.2 != '_'; s.3 = 'Greed'
# false, 0 out of 3 possible combinations 'true'.

#(E) Harvest and Greed are both shown on Friday
prove: s.1,s.2,s.3 |= 'Harvest'; s.1,s.2,s.3 |= 'Greed'
# false, 0 out of 3 possible combinations 'true'.

# If Limelight is shown exactly three times, Harvest is shown exactly twice, and
# Greed is shown exactly once...
import lsat3
{{i}} = 'Limelight' {{3}}; {{i}} = 'Harvest' {{2}}; {{i}} = 'Greed' {{1}}

# Reintroduce first film
t.0, f.0, s.0 ∈ Greed, Harvest, Limelight
{{j}}.0 == {{j}}.1 if {{j}}.1 != '_'
{{j}}.2 != '_' &&& {{j}}.1 == '_' ==> {{j}}.0 == {{j}}.2
{{j}}.2 == '_' ==> {{j}}.0 == {{j}}.3
# Order of operations helps here as the &&& block is evaluated before the if operator

# then which one of the following is a complete and accurate list of the films
# that could be the first film shown on Thursday?
range t.1
# Greed, Limelight

# (A) Harvest
# (B) Limelight
# (C) Greed, Harvest
# (D) Greed, Limelight
# (E) Greed, Harvest, Limelight


text:
export  excelTransPlot

function excelTransPlot( mpSim::MPsim, filename::AbstractString,
    sheetname::AbstractString="Transition flux plots" )

    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    if !ispath( filename )
        @warn string( "File '", filename,
            "' not found. Can't generate plots/reports." )
        return
    end  # if !ispath( filename )

    plotsToMake = Dict{Tuple, Vector{MP.TransitionType}}()

    XLSX.openxlsx( filename ) do xf
        if !XLSX.hassheet( xf, sheetname )
            @warn string( "Sheet '", sheetname,
                "' not found in file. Can't generate plots/reports." )
            return
        end  # if !XLSX.hassheet( xf, sheetname )

        xs = xf[sheetname]
        nToProcess = xs["B4"]

        for ii in (1:nToProcess) .+ 6
            # Read transition type, if given.
            transName = xs[ii, 1]
            sourceNode, targetNode = xs[ii, 2], xs[ii, 3]
            isNamed = !( transName isa Missing )
            isTargetted = !( ( sourceNode isa Missing ) &&
                ( targetNode isa Missing ) )

            if isNamed &&
                !( haskey( mpSim.transitionsByName, string( transName ) ) ||
                haskey( mpSim.recruitmentByName, string( transName ) ) )
                continue
            end  # if isNamed && ...

            if !isNamed && !isTargetted
                continue
            end  # if !isNamed && ...

            transName = transName isa Missing ? "" : string( transName )
            sourceNode = sourceNode isa Missing ? "" : string( sourceNode )
            sourceNode = sourceNode ∈ MP.specialNodes ? "" : sourceNode
            targetNode = targetNode isa Missing ? "" : string( targetNode )
            targetNode = targetNode ∈ MP.specialNodes ? "" : targetNode

            # Check if source/target pair is valid.
            isSourceValid = haskey( mpSim.baseNodeList, sourceNode ) ||
                ( lowercase( sourceNode ) ∈ vcat( "", "active" ) )
            isTargetValid = haskey( mpSim.baseNodeList, targetNode ) ||
                ( lowercase( targetNode ) ∈ vcat( "", "active" ) )

            if !( isSourceValid && isTargetValid ) ||
                ( isTargetted && ( sourceNode == targetNode ) )
                continue
            end  # if !( isSourceValid && isTargetValid ) || ...

            # The time resolution.
            timeRes = xs[ii, 4]

            if !( timeRes isa Real ) || ( timeRes <= 0 )
                continue
            end  # if !( timeRes isa Real ) || ...

            # Show plot in browser.
            showPlot = lowercase( string( xs[ii, 5] ) ) == "yes"

            # Save plot.
            savePlot = lowercase( string( xs[ii, 6] ) ) == "yes"

            # Generate report.
            makeReport = lowercase( string( xs[ii, 7] ) ) == "yes"

            plotkey = (timeRes, showPlot, savePlot, makeReport)

            trans = isNamed ?
                ( isTargetted ? (transName, sourceNode, targetNode) :
                transName ) : (sourceNode, targetNode)

            if haskey( plotsToMake, plotkey )
                push!( plotsToMake[plotkey], trans )
            else
                plotsToMake[plotkey] = [trans]
            end  # if haskey( plotsToMake, plotkey )
        end  # for ii in (1:nToProcess) .+ 6
    end  # XLSX.openxlsx( filename ) do xf

    plotdir = filename[1:(end-5)]

    if !ispath( plotdir )
        mkpath( plotdir )
    end  # if !ispath( plotdir )

    dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
    isFirst = true
    reportname = string( plotdir, "/transition flux ", dateStr )

    for plotkey in keys( plotsToMake )
        if !( any( [plotkey[2], plotkey[3], plotkey[4]] ) )
            continue
        end  # if !( any( [plotkey[2], plotkey[3], plotkey[4]] ) )

        # Generate the report.
        tStart = now()
        fluxReport = transitionFluxReport( mpSim, plotkey[1],
            plotsToMake[plotkey]... )
        reportGenerationTime = ( now() - tStart ).value / 1000.0

        if isempty( fluxReport )
            continue
        end  # if isempty( fluxReport )

        if plotkey[2] || plotkey[3]
            dateStr = string( " (", Dates.format( now(), "yyyymmdd HHMMSS" ),
                ";", plotkey[1], ")" )
            plotname = string( plotdir, "/transition flux", dateStr )
            generateTransFluxPlot( fluxReport, 12.0, plotkey[2], plotkey[3],
                plotname )
        end  # if plotkey[2] || ...

        if plotkey[4]
            writeTransFluxData( fluxReport, mpSim, 12.0, reportGenerationTime,
                reportname, isFirst )
            isFirst = false
        end  # if plotkey[4]
    end  # for plotkey in keys( plotsToMake )

end  # excelTransPlot( mpSim, filename, sheetname )


text:
module BioinformaticsBISC195

using Base: stop_reading
export normalizeDNA
export composition
export gc_content
export complement
export reverse_complement
export parse_fasta
export uniqueKmers
export kmerdist
export kmerloc
export kmertime
export kmertimes
export lengthcount
export minMax
export pairdist
export distsort

# ### NormalizeDNA Function
"""
    normalizeDNA(::AbstractString)

Function that takes a sequence and returns the uppercase and normalized version of it.

Returns an error if invalid bases are encountered.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> normalizeDNA("acatg")
    "ACATG"

    julia> normalizeDNA("LAC")
    ERROR: key 'L' not found

    julia> normalizeDNA(7)
    ERROR: MethodError: no method matching uppercase(::Int64)

    julia> normalizeDNA('G')
    "G"

    julia> s = normalizeDNA('C');

    julia> typeof(s)
    String
"""
function normalizeDNA(sequence)
    sequence= uppercase(sequence)
    rep = Dict('R' => 'N', 'Y' => 'N', 'S' => 'N', 'K' => 'N', 'M' => 'N', 'B' => 'N', 'D' => 'N', 'H' => 'N', 'W' => 'N','V' => 'N', '.' => 'N', '-' => 'N', 'A' => 'A', 'G' => 'G', 'C' => 'C', 'T' => 'T', 'N' => 'N')
    return join([rep[c] for c in sequence])
end

# ### Composition Function
"""
    composition(sequence)

Takes a DNA sequence and counts the number of each type of base in the DNA sequence.

Returns a Dictionary of how many of each exist in the order A, C, G, T.

Converts any ambiguous bases not in AGCT to "N".

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> composition("ACGGG")
    'A' => 1
    'G' => 3
    'C' => 1

    julia> composition("ACC")
    'A' => 1
    'C' => 2

    julia> A,C,T = composition("acctt")
    'A' => 1
    'T' => 2
    'C' => 2

    julia> composition("BAGGGRRRR")
    'A' => 1
    'G' => 3
    'N' => 5
"""
function composition(sequence)
    sequence = normalizeDNA(sequence) #make uppercase string, check invalid bases
    comp = Dict()
    for i in 1:length(sequence)
        pos = sequence[i]
        if haskey(comp, pos) == true
            comp[pos] = comp[pos] + 1
        else
            comp[pos] = 1
        end
    end
    return comp
end

# ### GC_Content Function
"""
    gc_content(sequence)

Calculates the GC ratio of a DNA sequence.

The GC ratio is the total number of G and C bases divided by the total number of G's, C's, A's, and T's within the sequence.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> gc_content("CATA")
    0.25

    julia> gc_content("cccg") * 20
    20.0

    julia> gc_content("AaTt")
    0.0

    julia> gc_content("ATty")
    0.0

    julia> gc_content("GGCCNNNN")
    1.0
"""
function gc_content(sequence)
    comp = composition(sequence)
    ng = get(comp, 'G', 0) #searches for each of these bases using the composition() function above
    nc = get(comp, 'C', 0)
    na= get(comp, 'A', 0)
    nt= get(comp, 'T', 0)
    return (ng + nc) / (ng + nc + nt + na) #divided by sum of AGCT w/ composition function
end

# ### Complement Function
"""
    complement(base)

Get the DNA complement of the provided base:

    A <-> T
    G <-> C

Accepts uppercase or lowercase `String` or `Char`,
but always returns an uppercase `String` (orig says Char).

If a valid base is not provided, the function returns an error.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> complement("A")
    "T"

    julia> complement("GCT")
    "CGA"

    julia> complement("AGCTACC")
    "TCGATGG"

    julia> complement("ZAG")
    ERROR: KeyError: key 'Z' not found

"""
function complement(sequence)
    sequence = normalizeDNA(sequence)
    ret = "" #initialize empty string to place the complements
    for i in 1:length(sequence)
        if sequence[i] == 'A'
            ret = ret * 'T' #if position equals A return "T" in that new position as the complement
        elseif sequence[i] == 'T'
            ret = ret * 'A'
        elseif sequence[i] == 'G'
            ret = ret * 'C'
        elseif sequence[i] == 'C'
            ret = ret * 'G'
        elseif sequence[i] == 'N'
            ret = ret * 'N'
        end
    end
    return ret
end

# ### Reverse Complement Function
"""
    reverse_complement(sequence)

Takes a DNA sequence and returns the reverse complement
of that sequence.

Takes lowercase or uppercase sequences, but always returns uppercase.

    Examples
    ≡≡≡≡≡≡≡≡≡≡
        julia> reverse_complement("ACCTTT")
        "AAAGGT"

        julia> reverse_complement("CCGTAGTA")
        "TACTACGG"

        julia> reverse_complement("ttacg")
        "CGTAA"

        julia> rc = reverse_complement("ACAG");

        julia> println(rc)
        CTGT
    """
function reverse_complement(sequence)
    st= "" #initialize empty string to place the reverse complements
    sequence =reverse(sequence)
    for i in 1:length(sequence)
        st = st*string(complement(sequence[i]))
    end
    return st
end

# ### Parse Function

"""
    function parse_fasta(path)

Reads a fasta-formated file and returns 2 vectors,
one containing parsed headers,
the other containing the entire sequence as a `String`.

Note: function does not validate DNA sequences for correctness.

Example
≡≡≡≡≡≡≡≡≡
    julia> ex1 = parse_fasta("data/ex1.fasta");

    julia> ex1 isa Tuple
    true

    julia> ex1[1]
    2-element Array{Tuple{String,String},1}:
      ("ex1.1", "easy")
      ("ex1.2", "multiline")

    julia> ex1[2]
    2-element Array{String,1}:
    "AATTATAGC"
    "CGCCCCCCAGTCGGATT"

    julia> ex2 = parse_fasta("data/ex2.fasta");

    julia> ex2[1]
    4-element Array{Tuple{String,String},1}:
      ("ex2.1", "oneper")
      ("ex2.2", "wrong")
      ("ex2.3", "wronger")
      ("ex2.4", "wrongest")

    julia> ex2[2]
    4-element Array{String,1}:
      "ATCCGT"
      "ATCGTGGaact"
      "ATCGTGGaact"
      "this isn't a dna string,but parse it anyway"
"""
function parse_fasta(path)
    headers= String[] #initialize empty array to hold the headers
    sequences= [] #initialize empty array to hold sequences
    tmp= ""
    for line in eachline(path)
        if startswith(line, '>')
            if !isempty(tmp)
                push!(sequences, tmp) #if it is empty, put the sequences within tmp
            end
            tmp= ""
            push!(headers, SubString(line, 2))
        else
            line = normalizeDNA(line)
            tmp= tmp*line
        end
    end
        push!(sequences, tmp)
    return (headers, sequences)
end

# ### Mean Lengths and Counts
"""
    function lengthcount(path)

Takes data and returns the mean and standard deviation of sequence lengths and of gc content within one tuple.

Returns a KeyError if sequence is not within "AGCT".

    Example
    ≡≡≡≡≡≡≡≡≡≡
        julia> using Statistics

        julia> lengthcount("data/ex1.fasta")
        (13.0, 5.656854249492381, 0.46405228758169936, 0.34199935822094457)

        julia> lengthcount("data/ex2.fasta")
        ERROR: KeyError: key 'I' not found
    """
function lengthcount(path)
    lengths= []
    counts= []
    data= parse_fasta(path) #parsing the data in order to separate headers and sequences
    for i in data[2]
        push!(lengths, length(i))
        gs = count(==('G'), i) #counts the number of g's
        cs = count(==('C'), i) #counts the number of c's
        GCcontent= (gs + cs)/length(i) #finds the GC content by dividing by length of sequence
        push!(counts, GCcontent)
    end
    return (mean(lengths), std(lengths), mean(counts), std(counts)) #using Statistics, takes mean and SD of each calculation
end

# ### Minimum and Maximum Function
"""
    minMax(path)

Takes data and returns a Tuple of the minimum and maximum of sequence lengths.

    Examples
    ≡≡≡≡≡≡≡≡≡≡
        julia> minMax("data/ex1.fasta")
        (9, 17)

        julia> minMax("data/datatry.fasta")
        (3, 64)

        julia> minMax("data/datasort.fasta")
        (21, 51)

        julia> typeof(minMax("data/datasort.fasta"))
        Tuple{Int64, Int64}

        julia> minMax("data/refined_data.fasta")
        (29833, 29881)
    """
function minMax(path)
    ret = parse_fasta(path)[2]
    retmin = length(minimum(ret))
    retmax = length(maximum(ret))
    return (retmin, retmax)
end

# ### Sequence Lengths Function
"""
        seqlength(path)

Takes data and returns a vector of the sequence lengths within the data.

    Examples
    ≡≡≡≡≡≡≡≡≡≡
        julia> seqlength("data/ex1.fasta")
        2-element Vector{Any}:
            9
            17

        julia> seqlength("data/datatry.fasta")
        3-element Vector{Any}:
            3
            64
            4

    """
function seqlength(path)
    lengths = []
    data = parse_fasta(path)
    for i in data[2]
         push!(lengths, length(i))
    end
    return lengths
end

# ### Unique Kmer Function
"""
    uniqueKmers(sequence, k)

Takes a sequence and a kmer length (k)
and returns a set of strings of the unique kmers that appear within the DNA sequence.

Converts invalid bases to "N" and converts lowercase sequences to uppercase sequences.

Example
≡≡≡≡≡≡≡≡≡
    julia> uniqueKmers("ATGCGATG", 4)
    Set{Any} with 5 elements:
    "TGCG"
    "ATGC"
    "GATG"
    "CGAT"
    "GCGA"

    julia> uniqueKmers("ATGCN", 2)
    Set{String} with 4 elements:
    "CN"
    "AT"
    "GC"
    "TG"

    julia> uniqueKmers("ACT", 2)
    Set{Any} with 2 elements:
    "AC"
    "CT"

    julia>uniqueKmers("RAcGcN", 4)
    Set{String} with 3 elements:
    "NACG"
    "CGCN"
    "ACGC"

    julia> typeof(uniqueKmers("RAcGcN", 4))
    Set{String}
"""
function uniqueKmers(sequence, k)
    sequence= normalizeDNA(sequence)
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = Set{String}()
    stopindex = length(sequence) - k + 1
    for i in 1:stopindex
        kmer= sequence[i:i+k-1]
        kmer = uppercase(kmer)
        push!(kmers, kmer)
    end
    return kmers
end

# ### Sorting Sequences to 5%
"""
    function sortingseq(path)

Takes a path (data set) and returns vectors:
one with the sequence lengths greater than 29500 and one with the corresponding headers.

Example
≡≡≡≡≡≡≡≡≡
    julia> sorting("data/data_sorting.fasta")
    (Any[29903, 29838, 29771, 29782], ["NC_045512.2 |China|Homo sapiens", "MT614596.1 |Egypt|Homo sapiens", "MT669321.1 |India|Homo sapiens", "MT670022.1 |Chile|Homo sapiens"])

    julia> g= sorting("data/data_sorting.fasta");

    julia> typeof(g)
    Vector{Any} (alias for Array{Any, 1})

"""
function sortingseq(path)
    sequences = seqlength(path)
    headers = parse_fasta(path)[1]
    check = findall(x->x<29800, sequences)
    deleteat!(sequences, check)
    deleteat!(headers, check)
    return (sequences, headers)
end

# ### KmerDistance Function
"""
    function kmerdist(set1, set2)

Takes two kmer sets and returns the distance between the kmer sets as a Float64.
Must return a positive number between 0 and 1.

Returns 0.0 for identical sets.
Returns 1.0 for sets with no similarity.

Example
≡≡≡≡≡≡≡≡≡
    julia> kmerdist(uniqueKmers("ATCGATG",2), uniqueKmers("GCATACC",2))
    0.9

    julia> kmerdist(uniqueKmers("GCGCAT",2), uniqueKmers("ATAT",2))
    0.8

    julia> kmerdist(uniqueKmers("AATA", 4), uniqueKmers("CGGCCCG", 4))
    1.0

    julia> typeof(kmerdist(uniqueKmers("ATCGATG",2), uniqueKmers("GCATACC",2)))
    Float64
"""
function kmerdist(set1, set2)
    return 1 - (length(intersect(set1, set2))/length(union(set1,set2))) #length of intersection of the sets divided by length of the union
end

# ### Kmertime Function
"""
    function kmertime(headers, sequences, k)

Takes a dataset and returns three vectors:
of the unique kmers along with corresponding headers for early, middle, and late time periods.

In this data: (and my project: "data/refined_data.fasta"), early: 2019, middle: 2020, late: 2021.

Example
≡≡≡≡≡≡≡≡≡
    julia> headers, sequences = parse_fasta("data/datatry.fasta");

    julia> kmertime(headers, sequences, 3)
    (Any[Set(["ACT"])], Any[Set(["TGT", "CTT", "GAC", "AGT", "GAA", "TGC", "TTC", "ACC", "ACA", "TAG"  …  "CCA", "TTT", "AGA", "TGG", "ATA", "AAG", "CAC", "AAA", "TTA", "CCG"])], Any[Set(["GGT", "AGG"])])

    julia> j= kmertime(headers, sequences, 3);

    julia> typeof(j)
    Tuple{Vector{Any}, Vector{Any}, Vector{Any}}
"""
function kmertime(headers, sequences, k=3) #whatto set k to?
    early= []
    middle= []
    late= [] #Initializing empty arrays for 3 time periods: early(2019), middle(2020), and late(2021)
    for i in 1:length(sequences)
        if occursin("2019", headers[i]) #Findall occurrences of 2019 per header
            push!(early, uniqueKmers(sequences[i], k)) #If the header contains the date "2019", the kmer will be pushed into the "early" array, calling the uniqueKmer function to process how many unique kmers exist in the sequence
        end
        if occursin("2020", headers[i])
            push!(middle, uniqueKmers(sequences[i], k)) #If the header contains the date "2020", the kmer will be pushed into the "middle" array
        end
        if occursin("2021", headers[i])
            push!(late, uniqueKmers(sequences[i], k)) #If the header contains the date "2021", the kmer will be pushed into the "late" array
        end
    end
    return (early, middle, late)
end

# ### Kmertimes Function
"""
    function kmertimes(path)

Takes a dataset and returns an Tuple of integers for
the number of unique kmers within early, middle, and late time periods.

Example
≡≡≡≡≡≡≡≡≡

    julia> headers, sequences = parse_fasta("data/datatry.fasta");

    julia> kmertimes(headers, sequences)

    julia> (1, 36, 2)

    julia> typeof(kmertimes(headers, sequences))
    Tuple{Int64, Int64, Int64}
"""
function kmertimes(headers, sequences, k=3)
    early, middle, late = kmertime(headers, sequences, k) #calling kmertime to separate the sequences by time period: 2019, 2020, 2021
    early_kmers = length(union(early...)) #calculating length of union of the sets within each time period
    middle_kmers = length(union(middle...))
    late_kmers = length(union(late...))
    return(early_kmers, middle_kmers, late_kmers) #returns the number of unique kmers in each time period
end

# ### Pairwise Distance Function
"""
    pairdist(path)

Takes a dataset and initiates a matrix for sequences and
returns the distance between pairs of sequences.

Example
≡≡≡≡≡≡≡≡≡
    julia> ("data/refined_data.fasta")[35][1]
    0.189873417721519

    julia> pairdist("data/refined_data.fasta")[70][1]
    0.24705882352941178

    julia> typeof(pairdist("data/refined_data.fasta"))
    Matrix{Float64} (alias for Array{Float64, 2})
"""

function pairdist(path)
    h = kmertime(parse_fasta(path)[1], parse_fasta(path)[2], 3)[1] #takes kmertime of headers and sequences and stores them with these letters
    j =  kmertime(parse_fasta(path)[1], parse_fasta(path)[2], 3)[2]
    k =  kmertime(parse_fasta(path)[1], parse_fasta(path)[2], 3)[3]
    mesh = vcat(h,j,k) #concatenates h, j, k
    ret = zeros(36, 36) #initiates empty matrix
    for i in 1:36
        for j in 1:36
            i <= j && continue #skips upper triangle, focuses bottom half of matrix
            d = kmerdist(mesh[i], mesh[j]) #finds distance between the kmers
            ret[i, j] = d
        end
    end
    return ret #returns the matrix with the added distance between kmers of sequences
end

# ### Distance Sort
"""
    distsort(mat)

Takes a matrix (of my "data/refined_data.fasta") and returns a Tuple of Vectors
containing the distances between unique kmers of the sequences, organized by time
period.

Example
≡≡≡≡≡≡≡≡≡
    julia> distsort(pairdist("data/refined_data.fasta"))[3][142]
    0.16883116883116878

    julia> typeof(distsort(pairdist("data/refined_data.fasta")))
    Tuple{Vector{Any}, Vector{Any}, Vector{Any}}
"""

function distsort(mat) #takes matrix as an argument
    early_dist = [] #3 separate one dimensional arrays initialized for each time period
    middle_dist = []
    late_dist = []
    for i in 1:12
        for j in 1:12 #goes through the early time periods (1:12 are the first 12 in my refined_data)
            push!(early_dist, mat[i,j])
        end
    end
    for i in 13:24 #goes thru the middle time periods
        for j in 13:24
            push!(middle_dist, mat[i,j])
        end
    end
    for i in 25:36 #goes thru the late time periods
        for j in 25:36
            push!(late_dist, mat[i,j])
        end
    end
    return (early_dist, middle_dist, late_dist)
end
end


text:
"""
    QR factorization and solve methods
"""
abstract type QRSolver{n,T} <: LinearSolver end

function triangularize(rs::AbstractVector, n::Int)
    R = zeros(n, n)
    off = 1
    for j = 1:n
        for k = 1:j
            R[k, j] = rs[off][1]
            off += 1
        end
    end
    return UpperTriangular(R)
end

function triu_perm(k::Int, j::Int)
    # TODO would be better to precompute those indices and store them in the QRSolver structure.
    ind = Int((j - 1) * j / 2) + k
    return ind::Int
end

function qr_solve!(gs_solver::QRSolver, x, b)
    qr_solve!(gs_solver.qs, gs_solver.rs, x, b)
end

function qr_solve!(qs, rs, x, b)
    n = length(b)
    for j = 1:n
        x[j] = transpose(qs[j]) * b
    end
    for j = n:-1:1
        for k = j+1:n
            x[j] -= rs[triu_perm(j, k)][1] * x[k]
        end
        x[j] /= rs[triu_perm(j, j)][1]
    end
    return nothing
end

function qr_matrix_solve!(gs_solver::QRSolver, X::Matrix{T}, B::Matrix{T}) where {T}
    qr_matrix_solve!(gs_solver.qs, gs_solver.rs, X, B)
end

function qr_matrix_solve!(qs, rs, X::Matrix{T}, B::Matrix{T}) where {T}
    n,m = size(B)
    for j in eachindex((1:n))
        for k in eachindex((1:m))
            @inbounds @views X[j,k] = transpose(qs[j]) * B[:,k]
        end
    end
    for j = n:-1:1
        for k = j+1:n
            @inbounds @views @. X[j,:] -= rs[triu_perm(j, k)][1] * X[k,:]
        end
        @inbounds @views @. X[j,:] /= rs[triu_perm(j, j)][1]
    end
    return nothing
end

"""
    Static Dense Modified Gram-Schmidt
"""
mutable struct SDMGSSolver{n,T} <: QRSolver{n,T}
    as::Vector{SVector{n,T}}
    qs::Vector{SVector{n,T}}
    rs::Vector{T}
    xv::Vector{T}
    xs::SVector{n,T}
end

function SDMGSSolver(n::Int; T::DataType=Float64)
    as = Vector{SVector{n,T}}([zeros(SVector{n,T}) for i=1:n])
    qs = Vector{SVector{n,T}}([zeros(SVector{n,T}) for i=1:n])
    rs = zeros(T,Int((n + 1) * n / 2))
    xv = zeros(T,n)
    xs = zeros(SVector{n,T})
    return SDMGSSolver{n,T}(as, qs, rs, xv, xs)
end

function SDMGSSolver(A::AbstractMatrix{T}) where T
    n, m = size(A)
    @assert n == m
    gs_solver = SDMGSSolver(n; T=T)
    factorize!(gs_solver, A)
    return gs_solver
end

"""
    Gram-Schmidt algorithm applied to A.
"""
function factorize!(gs_solver::SDMGSSolver{n,T}, A::AbstractMatrix{T}) where {n,T}
    for j in eachindex(1:n)
        @inbounds @views gs_solver.as[j] = A[:,j]
    end
    factorize!(gs_solver)
    return nothing
end

"""
    Gram-Schmidt algorithm applied to a.
"""
function factorize!(gs_solver::SDMGSSolver{n,T}, a::Vector{SVector{n,T}}) where {n,T}
    gs_solver.as .= a
    factorize!(gs_solver)
    return nothing
end

"""
    Gram-Schmidt algorithm applied to gs_solver.a.
"""
function factorize!(gs_solver::SDMGSSolver{n,T}) where {n,T}
    # Unpack
    as = gs_solver.as
    qs = gs_solver.qs
    rs = gs_solver.rs

    off = 1
    for j in eachindex(1:n)
        # qi
        qs[j] = as[j]
        for k in eachindex(1:j-1)
            # rk
            @inbounds rs[off] = dot(qs[j], qs[k])
            # qu
            @inbounds qs[j] -= qs[k] .* rs[off]
            off += 1
        end
        # re
        rs[off] = norm(qs[j],2)
        # q
        qs[j] /= rs[off]
        off += 1
    end
    return nothing
end

"""
    QR Back substitution after Gram-Schmidt.
"""
function qr_solve!(gs_solver::SDMGSSolver{n,T}, b) where {n,T}
    qs = gs_solver.qs
    rs = gs_solver.rs
    xv = gs_solver.xv

    for j in eachindex(1:n)
        xv[j] = transpose(qs[j]) * b
    end
    for j = n:-1:1
        for k = j+1:n
            xv[j] -= rs[triu_perm(j, k)] * xv[k]
        end
        xv[j] /= rs[triu_perm(j, j)]
    end
    gs_solver.xs = xv
    return nothing
end

"""
    QR solver
"""

function linear_solve!(solver::QRSolver, x::Vector{T}, A::SparseMatrixCSC{T,Int}, b::Vector{T}) where T
    factorize!(solver, A)
    qr_solve!(solver, x, b)
end

function linear_solve!(solver::QRSolver, X::Matrix{T}, A::AbstractMatrix{T}, B::Matrix{T}) where T
    factorize!(solver, A)
    qr_matrix_solve!(solver, X, B)
end

function sdmgs_solver(A::Array{T, 2}) where T
    SDMGSSolver(A)
end


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

name = "spectra"
version = v"1.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/yixuan/spectra/archive/refs/tags/v$(version).tar.gz", "45228b7d77b916b5384245eb13aa24bc994f3b0375013a8ba6b85adfd2dafd67")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/spectra-*
mkdir build && cd build/

cmake .. \
-DCMAKE_INSTALL_PREFIX=$prefix \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
#header only, so no products
products = Product[]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="Eigen_jll", uuid="bc6bbf8a-a594-5541-9c57-10b0d0312c70"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Clang
using Darknet_jll

# LIBDARKNET_HEADERS are those headers to be wrapped.
const LIBDARKNET_INCLUDE = joinpath(dirname(Darknet_jll.libdarknet_path), "..", "include") |> normpath
const LIBDARKNET_HEADERS = [joinpath(LIBDARKNET_INCLUDE, header) for header in readdir(LIBDARKNET_INCLUDE) if endswith(header, ".h")]

wc = init(; headers = LIBDARKNET_HEADERS,
            output_file = joinpath(@__DIR__, "libdarknet_api.jl"),
            common_file = joinpath(@__DIR__, "libdarknet_common.jl"),
            clang_includes = vcat(LIBDARKNET_INCLUDE, CLANG_INCLUDE),
            clang_args = ["-I", joinpath(LIBDARKNET_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"libdarknet",
            clang_diagnostics = true,
            )

run(wc)


text:
export MtlCompileOptions

const MTLCompileOptions = Ptr{MtCompileOptions}

mutable struct MtlCompileOptions
    handle::MTLCompileOptions
end

Base.unsafe_convert(::Type{MTLCompileOptions}, opts::MtlCompileOptions) = opts.handle

Base.:(==)(a::MtlCompileOptions, b::MtlCompileOptions) = a.handle == b.handle
Base.hash(opts::MtlCompileOptions, h::UInt) = hash(opts.handle, h)

function MtlCompileOptions()
    handle = mtNewCompileOpts()
    obj = MtlCompileOptions(handle)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(opts::MtlCompileOptions)
    mtRelease(opts.handle)
end

## properties

Base.propertynames(::MtlCompileOptions) = (:fastMathEnabled, :languageVersion)

const language_versions = Dict(
    MtLanguageVersion1_0 => v"1.0",
    MtLanguageVersion1_1 => v"1.1",
    MtLanguageVersion1_2 => v"1.2",
    MtLanguageVersion2_0 => v"2.0",
    MtLanguageVersion2_1 => v"2.1",
    MtLanguageVersion2_2 => v"2.2",
    MtLanguageVersion2_3 => v"2.3",
    MtLanguageVersion2_4 => v"2.4",
)

function Base.getproperty(opts::MtlCompileOptions, f::Symbol)
    if f === :fastMathEnabled
        mtCompileOptsFastMath(opts)
    elseif f === :languageVersion
        ver = mtCompileOptsLanguageVersion(opts)
        haskey(language_versions, ver) || error("Unknown language version $ver; please file an issue.")
        language_versions[ver]
    else
        getfield(opts, f)
    end
end

function Base.setproperty!(opts::MtlCompileOptions, f::Symbol, val)
    if f === :fastMathEnabled
        mtCompileOptsFastMathSet(opts, val)
    elseif f === :languageVersion
        isa(val, VersionNumber) ||
            throw(ArgumentError("languageVersion property should be a version number"))
        for (enum,ver) in language_versions
            if ver === val
                mtCompileOptsLanguageVersionSet(opts, enum)
                return
            end
        end
        error("Unknown language version $val")
    else
        setfield!(opts, f, val)
    end
end

## display

function Base.show(io::IO, opts::MtlCompileOptions)
    print(io, "CompileOptions(…)")
end

function Base.show(io::IO, ::MIME"text/plain", opts::MtlCompileOptions)
    println(io, "CompileOptions:")
    println(io, " fast math:        ", opts.fastMathEnabled)
    print(io,   " language version: ", opts.languageVersion)
end


text:
@testset "couple_flux_matrix.jl" begin

end


text:
#--------------------------------------------------------------------------------------------------
# Linear mapping tests
#--------------------------------------------------------------------------------------------------
using Test
using Hensel

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

@test Hensel.to01(0.5, 0, 1)    == 0.5
@test Hensel.to01(1.5, 1, 2)    == 0.5
@test Hensel.to01(1.3, 1, 2)    ≈  0.3
@test Hensel.to01(1.7, 1, 2)    ≈  0.7
@test Hensel.to01(-1.3, -2, -1) ≈  0.7
@test Hensel.to01(-1.7, -2, -1) ≈  0.3

@test Hensel.to01(3//2, 1, 2) == 1//2

@test Hensel.to01(0.5, (0, 1))    == 0.5
@test Hensel.to01(1.5, (1, 2))    == 0.5
@test Hensel.to01(1.3, (1, 2))    ≈  0.3
@test Hensel.to01(1.7, (1, 2))    ≈  0.7
@test Hensel.to01(-1.3, (-2, -1)) ≈  0.7
@test Hensel.to01(-1.7, (-2, -1)) ≈  0.3

@test Hensel.from01(1//2, 1, 2) == 3//2

@test Hensel.from01(0.5, 0, 1)   == 0.5
@test Hensel.from01(0.5, 1, 2)   == 1.5
@test Hensel.from01(0.3, 1, 2)   ≈  1.3
@test Hensel.from01(0.7, 1, 2)   ≈  1.7
@test Hensel.from01(0.7, -2, -1) ≈  -1.3
@test Hensel.from01(0.3, -2, -1) ≈  -1.7

@test Hensel.from01(0.5, (0, 1))   == 0.5
@test Hensel.from01(0.5, (1, 2))   == 1.5
@test Hensel.from01(0.3, (1, 2))   ≈  1.3
@test Hensel.from01(0.7, (1, 2))   ≈  1.7
@test Hensel.from01(0.7, (-2, -1)) ≈  -1.3
@test Hensel.from01(0.3, (-2, -1)) ≈  -1.7

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

max_stars_repo_name: ayarodionov/Hensel

max_issues_repo_name: ayarodionov/Hensel

max_forks_repo_name: ayarodionov/Hensel


text:
include("TBIO.jl")
module ProgramTunnel_fs

    using Formatting
    using ..TBIO

    export ProgramTunnelInfo, sendData, recvData!

    mutable struct Timing
        chk_freq             :: Float64
        timeout              :: Float64
        buffer               :: Float64
        recv_first_sleep     :: Float64
        recv_first_sleep_max :: Float64

        timeout_limit_cnt    :: Integer
        buffer_cnt           :: Integer
        recv_first_cnt       :: Integer

        function Timing(;
            chk_freq             :: Float64,
            timeout              :: Float64,
            buffer               :: Float64,
            recv_first_sleep_max :: Float64,
            recv_first_sleep     :: Float64,
        )

            T = new(
                chk_freq,
                timeout,
                buffer,
                recv_first_sleep,
                recv_first_sleep_max,
                0,
                0,
                0,
            )

            updateCount!(T)

            return T
        end

    end

    function updateCount!(T::Timing)
        T.timeout_limit_cnt = ceil(T.timeout          / T.chk_freq)
        T.buffer_cnt        = ceil(T.buffer           / T.chk_freq)
        T.recv_first_cnt    = ceil(T.recv_first_sleep / T.chk_freq)
    end

    mutable struct ProgramTunnelInfo

        nchars       :: Integer

        recv_fn      :: AbstractString
        send_fn      :: AbstractString

        recv_channels :: AbstractArray{Timing}

        rotate       :: Integer
        recv_trackno :: Integer
        send_trackno :: Integer

        path         :: AbstractString
        error_sleep  :: Float64

        history_len  :: Integer

        recv_fns     :: AbstractArray
        send_fns     :: AbstractArray

        function ProgramTunnelInfo(;
            path                 :: Union{AbstractString, Nothing} = "x_tmp",
            nchars               :: Integer            = 512,
            recv_fn              :: AbstractString     = "Y2X",
            send_fn              :: AbstractString     = "X2Y",
            recv_channels        :: Integer            = 1,
            chk_freq             :: AbstractFloat      = 0.05,
            timeout              :: AbstractFloat      = 60.0 * 30,
            buffer               :: AbstractFloat      = 0.1,
            recv_first_sleep_max :: AbstractFloat      = 5.00,
            recv_first_sleep     :: AbstractFloat      = 0.0,
            reverse_role         :: Bool               = false,
            rotate               :: Integer            = 100,
            error_sleep          :: Float64            = 0.05,
            history_len          :: Integer            = 5,
        )

            if chk_freq <= 0.0
                ErrorException("chk_freq must be positive.") |> throw
            end

            recv_channels = [
                Timing(
                    chk_freq             = chk_freq,
                    timeout              = timeout,
                    buffer               = buffer,
                    recv_first_sleep_max = recv_first_sleep_max,
                    recv_first_sleep     = recv_first_sleep,
                )
                for i = 1:recv_channels
            ]

            PTI = new(
                nchars,
                recv_fn,
                send_fn,
                recv_channels,
                rotate, 1, 1,
                path,
                1.0, history_len,
                [], [],
            )

            if reverse_role
                reverseRole!(PTI)
            end

            mkpath(path)
            updateFiles!(PTI)

            return PTI
        end
    end

    function updateFiles!(PTI::ProgramTunnelInfo)
        recv_fns = []
        send_fns = []
        for i = 1:PTI.rotate
            push!(recv_fns, joinpath(PTI.path, format("{:s}_{:03d}.tb", PTI.recv_fn, i)))
            push!(send_fns, joinpath(PTI.path, format("{:s}_{:03d}.tb", PTI.send_fn, i)))
        end

        PTI.recv_fns = recv_fns
        PTI.send_fns = send_fns
    end

    function cleanHistory(PTI::ProgramTunnelInfo)
        cleanDir(PTI, ((PTI.send_trackno-1 - PTI.history_length -1) % PTI.rotate) + 1)
    end

    function reverseRole!(PTI::ProgramTunnelInfo)
        PTI.recv_fn, PTI.send_fn = PTI.send_fn, PTI.recv_fn
    end

    function incTrackno(PTI::ProgramTunnelInfo, which::Symbol)
        if which == :recv
            rm( PTI.recv_fns[mod( PTI.recv_trackno - 1 - PTI.history_len, PTI.rotate) + 1], force=true)
            PTI.recv_trackno += 1

        elseif which == :send
            rm( PTI.send_fns[mod( PTI.send_trackno - 1 - PTI.history_len, PTI.rotate) + 1], force=true)
            PTI.send_trackno += 1

        end
    end

    function sendData(
        PTI  :: ProgramTunnelInfo,
        msg  :: AbstractString,
        arrs :: AbstractArray,
    )

        send_fn = PTI.send_fns[mod(PTI.send_trackno - 1, PTI.rotate) + 1]

        while true
            try

                writeTB(
                    send_fn,
                    format("{:d}#{:s}", PTI.send_trackno, msg),
                    PTI.nchars,
                    arrs,
                )
            catch ex
                throw(ex)
                println(string(ex))
                println("keep sending msg.")
                sleep(PTI.error_sleep)
            end
            break
        end

        incTrackno(PTI, :send)

    end

    function recvData!(
        PTI     :: ProgramTunnelInfo,
        arrs    :: AbstractArray;
        which   :: Integer = 1,
    )

        channel = PTI.recv_channels[which]

        recv_fn = PTI.recv_fns[mod(PTI.recv_trackno - 1, PTI.rotate) + 1]

        println(format("[recvData!] [{:d}] Expecting file: {:s}. First sleep: {:.2f} ", which, recv_fn, channel.recv_first_sleep))
        get_through = false
        sleep(channel.recv_first_sleep)
        if isfile(recv_fn)
            channel.recv_first_sleep -= channel.chk_freq
            channel.recv_first_sleep = max(0.0, channel.recv_first_sleep)

            get_through = true
            println(format("[recvData!] [{:d}] Message is already there. Adjust recv_first_sleep to : {:.2f}", which, channel.recv_first_sleep))
        else
            for cnt in 1:(channel.timeout_limit_cnt - channel.recv_first_cnt)

                sleep(channel.chk_freq)

                if isfile(recv_fn)
                    get_through = true

                    if cnt <= channel.buffer_cnt

                        println(format("[recvData!] [{:d}] Good guess of the recv_first_sleep : {:.2f}", which, channel.recv_first_sleep))

                    elseif channel.recv_first_sleep < channel.recv_first_sleep_max

                        # Out of buffer, need to adjust: increase PTI.recv_first_sleep
                        channel.recv_first_sleep += channel.chk_freq
                        channel.recv_first_sleep = min(channel.recv_first_sleep_max, channel.recv_first_sleep)
                        println(format("[recvData!] [{:d}] Out of buffer. recv_first_sleep to : {:.2f} ", which, channel.recv_first_sleep))

                    else
                        println(format("[recvData!] [{:d}] Out of buffer. But reach to recv_first_sleep_max : {:.2f}", which, channel.recv_first_sleep))
                    end

                    break

                end

            end
        end

        updateCount!(channel)

        if ! get_through
            ErrorException(format("[recvData!] [{:d}] No further incoming message within timeout.", which)) |> throw
        end

        local msg

        while true
            try
                msg = readTB!(
                    recv_fn,
                    PTI.nchars,
                    arrs,
                )

                if msg == nothing

                    println("File does not exist / not the expected size / checksum failed. Keep receiving...")
                    sleep(PTI.error_sleep)
                    continue

                end

                recv_no, msg = split(msg, "#")
                if parse(Int, recv_no) != PTI.recv_trackno
                    println(format("Recive file trackno does not match. Expect {:d} but got {:s}. Keep receiving...", PTI.recv_trackno, recv_no))
                    sleep(PTI.error_sleep)
                    continue
                end

            catch ex
                println(string(ex))
                println("Keep receiving...")
                sleep(PTI.error_sleep)
                continue
            end

            break
        end

        incTrackno(PTI, :recv)

        return msg
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# This file autogenerated through https://github.com/timholy/SnoopCompile.jl.  Precompile every day!
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Nettle.call, (Type{Nettle.Decryptor}, ASCIIString, Array{UInt8, 1},))
    precompile(Nettle.call, (Type{Nettle.HMACState}, ASCIIString, Array{UInt8, 1},))
    precompile(Nettle.call, (Type{Nettle.CipherType}, Nettle.NettleCipher,))
    precompile(Nettle.call, (Type{Nettle.HashType}, Nettle.NettleHash, Ptr{Void},))
    precompile(Nettle.call, (Type{Nettle.Encryptor}, ASCIIString, Array{UInt8, 1},))
    precompile(Nettle.call, (Type{Nettle.Hasher}, ASCIIString,))
    precompile(Nettle.call, (Type{Nettle.HMACState}, ASCIIString, ASCIIString,))
    precompile(Nettle.encrypt!, (Nettle.Encryptor, Array{UInt8, 1}, Array{UInt8, 1},))
    precompile(Nettle.decrypt!, (Nettle.Decryptor, Array{UInt8, 1}, Array{UInt8, 1},))
    precompile(Nettle.digest!, (Nettle.HMACState,))
    precompile(Nettle.get_cipher_types, ())
    precompile(Nettle.get_hash_types, ())
    precompile(Nettle.digest!, (Nettle.Hasher,))
    precompile(Nettle.call, (Type{Nettle.HashType}, ASCIIString, UInt32, UInt32, UInt32, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void},))
    precompile(Nettle.get_libnettle_version, ())
    precompile(Nettle.call, (Type{Nettle.CipherType}, ASCIIString, UInt32, UInt32, UInt32, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void},))
    precompile(Nettle.update!, (Nettle.HMACState, Array{UInt8, 1},))
    precompile(Nettle.update!, (Nettle.Hasher, ASCIIString,))
    precompile(Nettle.decrypt, (ASCIIString, Array{UInt8, 1}, Array{UInt8, 1},))
    precompile(Nettle.update!, (Nettle.HMACState, ASCIIString,))
    precompile(Nettle.encrypt, (ASCIIString, Array{UInt8, 1}, Array{UInt8, 1},))
    precompile(Nettle.hexdigest!, (Nettle.Hasher,))
    precompile(Nettle.__init__, ())
end


text:
@testset "30.substring-with-concatenation-of-all-words.jl" begin
    @test find_word_concatenation("barfoothefoobarman", ["foo","bar"]) == [1, 10]
    @test find_word_concatenation("wordgoodgoodgoodbestword", ["word","good","best","word"]) == Int[]
    @test find_word_concatenation("barfoofoobarthefoobarman", ["bar","foo","the"]) == [7, 10, 13]
end


text:
using Base: batch_size_err_str
@testset "d4rl_pybullet" begin
    ds = dataset(
        "hopper-bullet-mixed-v0";
        repo="d4rl-pybullet",
        style = style,
        rng = rng,
        is_shuffle = true,
        batch_size = batch_size
    )

    n_s = 15
    n_a = 3

    N_samples = 59345

    data_dict = ds.dataset

    @test size(data_dict[:state]) == (n_s, N_samples)
    @test size(data_dict[:action]) == (n_a, N_samples)
    @test size(data_dict[:reward]) == (1, N_samples)
    @test size(data_dict[:terminal]) == (1, N_samples)

    for sample in Iterators.take(ds, 3)
        @test typeof(sample) <: NamedTuple{SARTS}
        @test size(sample[:state]) ==  (n_s, batch_size)
        @test size(sample[:action]) ==  (n_a, batch_size)
        @test size(sample[:reward]) ==  (1, batch_size) || size(sample[:reward]) == (batch_size,)
        @test size(sample[:terminal]) ==  (1, batch_size) || size(sample[:terminal]) == (batch_size,)
    end

end


text:
"""
    conv_aux(conv_otf, rec, otf)

Calculate the convolution between `rec` and `otf`.
The used convolution function is `conv_otf`.
`conv_otf` can be exchanged to be a rfft, fft or plan_fft based routine.

This function is just defined to speed up automatic differentiation
and it's custom defined adjoint.
"""
function conv_aux(conv_otf, rec, otf)
    return conv_otf(rec, otf)
end

 # define custom adjoint for conv_aux
function ChainRulesCore.rrule(::typeof(conv_aux), conv, rec, otf)
    Y = conv_aux(conv, rec, otf)
    function conv_aux_pullback(barx)
        return zero(eltype(rec)), zero(eltype(rec)), conv(barx, conj(otf)), zero(eltype(rec))
    end
    return Y, conv_aux_pullback
end


text:
#
# Function that obtains the best fit between a SIRD simulation and trajectory
# data obtained from simulations
#

using Optim

# The two possible functions to be fited: all data or one of the curves, defined by ipop

function fit(md_input :: MDInput, md_traj :: MDTraj;
             ipop :: Int64 = -1,
             bounded :: Bool = true, lower :: Vector{Float64} = zeros(4), upper :: Vector{Float64} = [10.,10.,10.,10.])

  sird_input = SIRDInput(tmax=md_input.tmax,nsave=md_traj.nsteps,
                         Si=md_input.Si,
                         print=false)
  sird_traj = SIRDTraj(sird_input)

  k0 = [ rand()*sird_input.kc,
         rand()*sird_input.kd,
         rand()*sird_input.ki,
         rand()sird_input.kiu ]

  fall(k0 :: Vector{Float64}) = fiteval!( k0, md_traj, sird_input, sird_traj )
  fpop(k0 :: Vector{Float64}) = fiteval!( k0, md_traj, sird_input, sird_traj, ipop )

  if bounded
    if ipop == -1
      x = Optim.optimize( fall, lower, upper, k0, Fminbox(NelderMead()) )
    else
      x = Optim.optimize( fpop, lower, upper, k0, Fminbox(NelderMead()) )
    end
  else
    if ipop == -1
      x = Optim.optimize( fall, k0, NelderMead() )
    else
      x = Optim.optimize( fpop, k0, NelderMead() )
    end
  end

  # Perform simulation with optimal parameters

  sird_input = SIRDInput( kc = x.minimizer[1],
                          kd = x.minimizer[2],
                          ki = x.minimizer[3],
                          kiu = x.minimizer[4],
                          tmax = md_input.tmax,
                          nsave = md_traj.nsteps )
  if ipop == -1
    fitresult = fiteval!(x.minimizer,md_traj,sird_input,sird_traj)
  else
    fitresult = fiteval!(x.minimizer,md_traj,sird_input,sird_traj,ipop)
  end

  println("-------------------------------------------------------")
  println(" Fitting result: ")
  println(" Error: ", fitresult )
  println(" Best ks found: ")
  println(" kc = ", sird_input.kc)
  println(" kd = ", sird_input.kd)
  println(" ki = ", sird_input.ki)
  println(" kiu = ", sird_input.kiu)
  println("-------------------------------------------------------")

  return sird_input, sird_traj

end


text:
module LBNumber

include("./../../constraints/geometric/constants.jl")
include("./1_if.jl")
include("./2_calc_number.jl")
using JuMP
using .Constants
using .IfSimpleLoadBearing
using .CalcNumberLoadBearing

export cons_lb_number_load_bearing

function cons_lb_number_load_bearing(m)

    m = if_simple_lb(m)
    m = calc_number_lb(m)

    return m
end

end


text:
using FastGaussQuadrature

abstract type SauterSchwabStrategy end

struct CommonFace{A}       <: SauterSchwabStrategy qps::A end
struct CommonEdge{A}       <: SauterSchwabStrategy qps::A end
struct CommonVertex{A}     <: SauterSchwabStrategy qps::A end
struct PositiveDistance{A} <: SauterSchwabStrategy qps::A end

function _legendre(n,a,b)
    x, w = FastGaussQuadrature.gausslegendre(n)
    w .*= (b-a)/2
    x = (x.+1)/2*(b-a).+a
    collect(zip(x,w))
end

function sauterschwab_parameterized(integrand, method::CommonFace)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_cf(integrand, η1, η2, η3, ξ)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ, w4) in qps)
end

function sauterschwab_parameterized(integrand, method::CommonEdge)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_ce(integrand, η1, η2, η3, ξ)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ, w4) in qps)
end

function sauterschwab_parameterized(integrand, method::CommonVertex)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_cv(integrand,η1, η2, η3, ξ)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ, w4) in qps)
end

function sauterschwab_parameterized(integrand, method::PositiveDistance)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_pd(integrand, η1, η2, η3, ξ)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ, w4) in qps)
end


text:
"""
    SGHMC(n_iters::Int, learning_rate::Float64, momentum_decay::Float64)

Stochastic Gradient Hamiltonian Monte Carlo sampler.

Usage:

```julia
SGHMC(1000, 0.01, 0.1)
```

Arguments:

- `n_iters::Int` : Number of samples to pull.
- `learning_rate::Float64` : The learning rate.
- `momentum_decay::Float64` : Momentum decay variable.

Example:

```julia
@model example begin
  ...
end

sample(example, SGHMC(1000, 0.01, 0.1))
```
"""
mutable struct SGHMC{AD, T} <: StaticHamiltonian{AD}
    n_iters::Int       # number of samples
    learning_rate::Float64   # learning rate
    momentum_decay::Float64   # momentum decay
    space::Set{T}    # sampling space, emtpy means all
    gid::Int
end
SGHMC(args...) = SGHMC{ADBackend()}(args...)
function SGHMC{AD}(learning_rate::Float64, momentum_decay::Float64, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return SGHMC{AD, eltype(_space)}(1, learning_rate, momentum_decay, _space, 0)
end
function SGHMC{AD}(n_iters, learning_rate, momentum_decay) where AD
    return SGHMC{AD, Any}(n_iters, learning_rate, momentum_decay, Set(), 0)
end
function SGHMC{AD}(n_iters, learning_rate, momentum_decay, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return SGHMC{AD, eltype(_space)}(n_iters, learning_rate, momentum_decay, _space, 0)
end
function SGHMC{AD1}(alg::SGHMC{AD2, T}, new_gid::Int) where {AD1, AD2, T}
    return SGHMC{AD1, T}(alg.n_iters, alg.learning_rate, alg.momentum_decay, alg.space, new_gid)
end
function SGHMC{AD, T}(alg::SGHMC, new_gid::Int) where {AD, T}
    return SGHMC{AD, T}(alg.n_iters, alg.learning_rate, alg.momentum_decay, alg.space, new_gid)
end

function step(model, spl::Sampler{<:SGHMC}, vi::VarInfo, is_first::Val{true})
    spl.alg.gid != 0 && link!(vi, spl)

    # Initialize velocity
    v = zeros(Float64, size(vi[spl]))
    spl.info[:v] = v

    spl.alg.gid != 0 && invlink!(vi, spl)
    return vi, true
end

function step(model, spl::Sampler{<:SGHMC}, vi::VarInfo, is_first::Val{false})
    # Set parameters
    η, α = spl.alg.learning_rate, spl.alg.momentum_decay

    @debug "X-> R..."
    if spl.alg.gid != 0
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    @debug "recording old variables..."
    θ, v = vi[spl], spl.info[:v]
    _, grad = gradient(θ, vi, model, spl)
    verifygrad(grad)

    # Implements the update equations from (15) of Chen et al. (2014).
    @debug "update latent variables and velocity..."
    θ .+= v
    v .= (1 - α) .* v .- η .* grad .+ rand.(Normal.(zeros(length(θ)), sqrt(2 * η * α)))

    @debug "saving new latent variables..."
    vi[spl] = θ

    @debug "R -> X..."
    spl.alg.gid != 0 && invlink!(vi, spl)

    @debug "always accept..."
    return vi, true
end


text:
#######################################################################
# System
#######################################################################
struct System{LA<:Lattice,HA<:Hamiltonian}
    lattice::LA
    hamiltonian::HA
end

struct ParametricSystem{LA<:Lattice,HA<:ParametricHamiltonian}
    lattice::LA
    hamiltonian::HA
end

displaycompatible(lat, ham) = iscompatible(lat, ham) ? "compatible" : "incompatible"

function Base.show(io::IO, sys::System)
    ioindent = IOContext(io, :indent => "  ")
    print(io,
        "System : bundled Lattice and Hamiltonian ($(displaycompatible(
            sys.lattice, sys.hamiltonian)))",
        "\n")
    print(ioindent, sys.lattice, "\n")
    print(ioindent, sys.hamiltonian)
end

function Base.show(io::IO, sys::ParametricSystem)
    ioindent = IOContext(io, :indent => "  ")
    print(io,
        "Parametric System : bundled Lattice and Hamiltonian ($(displaycompatible(
            sys.lattice, sys.hamiltonian.hamiltonian)))",
        "\n")
    print(ioindent, sys.lattice, "\n")
    print(ioindent, sys.hamiltonian.hamiltonian)
end

(ps::ParametricSystem)(;kw...) = System(ps.lattice, ps.hamiltonian(;kw...))

# External API #

system(lat::Lattice, t::AbstractTightbindingModel...; kw...) =
    System(lat, hamiltonian(lat, t...; kw...))
system(lat::Lattice, f::Function, t::AbstractTightbindingModel...; kw...) =
    ParametricSystem(lat, hamiltonian(lat, f, t...; kw...))

superlattice(sys::System, v...; kw...) =
    System(superlattice(sys.lattice, v...; kw...), sys.hamiltonian)
hamiltonian(sys::System) =
    System(sys.lattice, hamiltonian(sys.lattice, sys.hamiltonian))
hamiltonian(sys::System, t::AbstractTightbindingModel...; kw...) =
    System(sys.lattice, hamiltonian(sys.lattice, t...; kw...))
lattice(sys::System) =
    System(lattice(sys.lattice), sys.hamiltonian)

# """
#     System(sublats::Sublat... [, br::Bravais[, model::Model]]; dim::Val, postype, hamtype)

# Build a `System{E,L,T,Tv}` of `L` dimensions in `E`-dimensional embedding
# space and composed of `T`-typed sites and `Tv`-typed Hamiltonian. See
# `Sublat`, `Bravais` and `Model` for syntax.  To override the embedding dimension `E`, use
# keyword `dim = Val(E)`. Similarly, override types `T` and `Tv` with `ptype = T` and
# `htype = Tv`.

#     System(presetname::$(NameType)[, model]; kw...)

# Build a system from a given preset, and optionally add a `model` to it. Preset
# names are defined in the dictionary `Elsa.systempresets`, together with their
# allowed keyword arguments `kw`. These typically include `norbitals` for sublattices.

#     System(sys::System, model::Model)

# Build a new system with the same sublattices and bravias vectors as `sys`
# but replacing its Hamiltonian with a new `model`.

# # Examples
# ```jldoctest
# julia> System(Sublat((1,0), (0,1); name = :C, norbitals = 2),
#               Sublat((0.5,0.5); name = :D),
#               Bravais((1,0), (0,2)),
#               Model(Hopping((r,dr) -> @SMatrix[0.1; dr[1]], sublats = (:C, :D), range = 1)))
# System{2,2,Float64,Complex{Float64}} : 2D system in 2D space
#   Bravais vectors     : ((1.0, 0.0), (0.0, 2.0))
#   Sublattice names    : (:C, :D)
#   Sublattice orbitals : (2, 1)
#   Total sites         : 3 [Float64]
#   Total hoppings      : 8 [Complex{Float64}]
#   Coordination        : 2.6666666666666665

# julia> System(LatticePresets.honeycomb(), Model(Hopping(@SMatrix[1 2; 0 1], sublats = (1,2))),
#               dim = Val(3), htype = Float32, norbitals = 2)
# System{3,2,Float64,Float32} : 2D system in 3D space
#   Bravais vectors     : ((0.5, 0.866025, 0.0), (-0.5, 0.866025, 0.0))
#   Sublattice names    : (:A, :B)
#   Sublattice orbitals : (2, 2)
#   Total sites         : 2 [Float64]
#   Total hoppings      : 6 [Float32]
#   Coordination        : 3.0

# julia> Tuple(keys(Elsa.systempresets))
# (LatticePresets.bcc(), LatticePresets.cubic(), LatticePresets.honeycomb(), :linear, :graphene_bilayer, :square, :triangular)
# ```

# # See also:

#     `Model`, `Onsite`, `Hopping`, `grow`, `transform`, `combine`
# """
# struct System{E,L,T,Tv,S<:SystemInfo,EL}
#     lattice::Lattice{E,L,T,EL}
#     hamiltonian::Operator{Tv,L}
#     velocity::Operator{Tv,L}
#     sysinfo::S
# end

# System(s1::Sublat, s2...; kw...) = System(_collectlattice((s1,), s2...)...; kw...)

# _collectlattice(ss::Tuple, s1::Sublat, s2...) = _collectlattice((ss..., s1), s2...)
# _collectlattice(ss::Tuple, b::Bravais, m::Model = Model()) = ((b, ss...), m)
# _collectlattice(ss::Tuple, m::Model = Model()) = (ss, m)

# System(latparts::Tuple, m::Model; kw...) = System(Lattice(latparts...; kw...), m; kw...)

# function System(lat::Lattice{E,L,T}, model::Model{Tv} = Model();
#                 htype::Type{Tv2} = Tv, kw...) where {E,L,T,Tv,Tv2}
#     actualmodel = convert(Model{Tv2}, model)
#     sysinfo = SystemInfo(lat, actualmodel)
#     hamiltonian = Operator(lat, sysinfo)
#     velocity = boundaryoperator(hamiltonian)
#     return System(lat, hamiltonian, velocity, sysinfo)
# end

# # System(name::NameType; kw...) = systempresets[name](; kw...)
# # System(name::NameType, model::Model; kw...) = combine(System(name; kw...), model)

# System(sys::System{E,L,T,Tv}, model::Model; kw...) where {E,L,T,Tv} =
#     System(sys.lattice, convert(Model{Tv}, model); kw...)

# Base.show(io::IO, sys::System{E,L,T,Tv}) where {E,L,T,Tv} = print(io,
# "System{$E,$L,$T,$Tv} : $(L)D system in $(E)D space
#   Bravais vectors     : $(displayvectors(sys.lattice.bravais))
#   Sublattice names    : $((sublatnames(sys.lattice)... ,))
#   Sublattice orbitals : $((norbitals(sys)... ,))
#   Total sites         : $(nsites(sys)) [$T]
#   Total hoppings      : $(nlinks(sys)) [$Tv]
#   Coordination        : $(coordination(sys))")

# # Treat System as scalar
# Broadcast.broadcastable(sys::System) = Ref(sys)

# #######################################################################
# # System internal API
# #######################################################################

# nsublats(sys::System) = length(sys.lattice.sublats)
# nsublats(lat::Lattice) = length(lat.sublats)

# norbitals(sys::System) = sys.sysinfo.norbitals

# dim(s::System) = isempty(s.sublats) ? 0 : sum(dim, s.sublats)

# nsites(sys::System) = sum(sys.sysinfo.nsites)

# nlinks(sys::System) = nlinks(sys.hamiltonian)

# isunlinked(sys::System) = nlinks(sys) == 0

# coordination(sys::System) = nlinks(sys.hamiltonian)/nsites(sys)

# site(sys::System, s::Int, i::Int) = sys.lattice.sublats[s].sites[i]

# function boundingbox(sys::System{E,L,T}) where {E,L,T}
#     bmin = zero(MVector{E,T})
#     bmax = zero(MVector{E,T})
#     foreach(sl -> foreach(s -> _boundingbox!(bmin, bmax, s), sl.sites), sys.lattice.sublats)
#     return (SVector(bmin), SVector(bmax))
# end

# function _boundingbox!(bmin, bmax, site)
#     bmin .= min.(bmin, site)
#     bmax .= max.(bmax, site)
#     return nothing
# end

# #######################################################################
# # External API
# #######################################################################
# """
#     transform!(s::Union{Lattice,System}, f::Function; sublats)

# Change `s` in-place by moving positions `r` of sites in sublattices specified by `sublats`
# (all by default) to `f(r)`. Bravais vectors are also updated, but the system Hamiltonian is
# unchanged.

#     s |> transform!(f::Function; kw...)

# Functional syntax, equivalent to `transform!(s, f; kw...)`

# # Examples
# ```jldoctest
# julia> transform!(LatticePresets.honeycomb(dim = Val(3)), r -> 2r + SVector(0,0,1))
# Lattice{3,2,Float64} : 2D lattice in 3D space
#   Bravais vectors     : ((1.0, 1.732051, 0.0), (-1.0, 1.732051, 0.0))
#   Sublattice names    : (:A, :B)
#   Sublattice orbitals : ((:noname), (:noname))
# ```
# """
# transform!

# """
#     transform(s::Union{Lattice,System}, f::Function; sublats)

# Perform a `transform!` on a `deepcopy` of `s`, so that the result is completely decoupled
# from the original. See `transform!` for more information.
# ```
# """
# transform

# """
#     combine(systems::System... [, model::Model])

# Combine sublattices of a set of systems (giving them new names if needed) and optionally
# apply a new `model` that describes their coupling.

#     system |> combine(model)

# Functional syntax, equivalent to `combine(system, model)`, which adds `model` to an existing
# system.

# # Examples
# ```jldoctest
# julia> combine(System(LatticePresets.honeycomb()), System(LatticePresets.honeycomb()), Model(Hopping(1, sublats = (:A, 4))))
# System{2,2,Float64,Complex{Float64}} : 2D system in 2D space
#   Bravais vectors     : ((0.5, 0.866025), (-0.5, 0.866025))
#   Sublattice names    : (:A, :B, 3, 4)
#   Sublattice orbitals : (1, 0, 0, 1)
#   Total sites         : 4 [Float64]
#   Total hoppings      : 6 [Complex{Float64}]
#   Coordination        : 1.5
# ```
# """
# combine

# """
#     grow(system::System{E,L}; supercell = SMatrix{L,0,Int}(), region = r -> true)

# Transform an `L`-dimensional `system` into another `L2`-dimensional system with a different
# supercell, so that the new Bravais matrix is `br2 = br * supercell`, and only sites with
# `region(r) == true` in the unit cell are included.

# `supercell` can be given as an integer matrix `s::SMatrix{L,L2,Int}`, a single integer
# `s::Int` (`supercell = s * I`), a single `NTuple{L,Int}` (`supercell` diagonal), or a tuple
# of  `NTuple{L,Int}`s (`supercell` columns). Note that if the new system dimension `L2` is
# smaller than the original, a bounded `region` function should be provided to determine the
# extension of the remaining dimensions.

#     system |> grow(region = f, supercell = s)

# Functional syntax, equivalent to `grow(system; region = f, supercell = s)

# # Examples
# ```jldoctest
# julia> grow(System(:triangular, Model(Hopping(1))), supercell = (3, -3), region = r-> 0 < r[2] < 12)
# System{2,2,Float64,Complex{Float64}} : 2D system in 2D space
#   Bravais vectors     : ((1.5, 2.598076), (1.5, -2.598076))
#   Sublattice names    : (1,)
#   Sublattice orbitals : (1,)
#   Total sites         : 3 [Float64]
#   Total hoppings      : 6 [Complex{Float64}]
#   Coordination        : 2.0
# ```

# # See also

#     'Region`
# """
# grow

# """
#     bound(system::System{E,L,T}; except = ()))

# Remove the periodicity of the `system`'s lattice along the specified `axes`, dropping
# Bloch Hamiltonian harmonics accordingly.

#     system |> bound(; kw...)

# Functional syntax, equivalent to `bound(system; kw...)``

# # Examples
# ```jldoctest
# julia> bound(System(LatticePresets.cubic()), except = (1, 3))
# System{3,2,Float64,Complex{Float64}} : 2D system in 3D space
#   Bravais vectors     : ((1.0, 0.0, 0.0), (0.0, 0.0, 1.0))
#   Sublattice names    : (1,)
#   Sublattice orbitals : (0,)
#   Total sites         : 1 [Float64]
#   Total hoppings      : 0 [Complex{Float64}]
#   Coordination        : 0.0

# julia> System(:triangular) |> bound()
# System{2,0,Float64,Complex{Float64}} : 0D system in 2D space
#   Bravais vectors     : ()
#   Sublattice names    : (1,)
#   Sublattice orbitals : (0,)
#   Total sites         : 1 [Float64]
#   Total hoppings      : 0 [Complex{Float64}]
#   Coordination        : 0.0
# ```
# """
# bound

# """
#     hamiltonian(system; k, ϕn, avoidcopy = false)

# Return the Bloch Hamiltonian of an `L`-dimensional `system` in `E`-dimensional space at
# a given `E`-dimensional Bloch momentum `k`, or alternatively `L`-dimensional normalised
# Bloch phases `ϕn = k*B/2π`, where `B` is the system's Bravais matrix.
# By default the Hamiltonian at zero momentum (Gamma point) is returned. For `0`-dimensional
# systems, the Bloch Hamiltonian is simply the Hamiltonian of the system. `avoidcopy = true`
# is a performance optimization that allows returning the internal system hamiltonian matrix
# instead of a copy (can lead to confusion/bugs - a second call to `hamiltonian` will
# overwrite the result of the first when `avoidcopy = true`).

# # Examples
# ```jldoctest
# julia> hamiltonian(System(LatticePresets.honeycomb(), Model(Hopping(1))), ϕn = (0,0.5))
# 2×2 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
#   [1, 1]  =  -2.0+0.0im
#   [2, 1]  =  1.0-1.22465e-16im
#   [1, 2]  =  1.0+1.22465e-16im
#   [2, 2]  =  -2.0+0.0im
# ```

# # See also

#     'velocity`
# """
# hamiltonian

# """
#     velocity(system::System{E,L}, axis = missing; k, ϕn, avoidcopy = false)

# Return the velocity operator `∂H(k)` along the specified crystallographic `axis`, which can
# be an integer from 1 to `L` (unit axis), an `SVector{L}` or an `NTuple{L}`. If no `axis` is
# given, a tuple of all velocity operators along each of the `L` axis is returned. See
# `hamiltonian` for details on the `k` and `ϕn` keywords. `avoidcopy = true`
# is a performance optimization that allows returning the internal system velocity matrix
# instead of a copy (can lead to confusion/bugs - a second call to `velocity` will
# overwrite the result of the first when `avoidcopy = true`).

# # Examples
# ```jldoctest
# julia> velocity(System(LatticePresets.honeycomb(), Model(Hopping(1))), 1; ϕn = (0,0.5))
# 2×2 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
#   [1, 1]  =  -2.0+0.0im
#   [2, 1]  =  1.0-1.22465e-16im
#   [1, 2]  =  1.0+1.22465e-16im
#   [2, 2]  =  -2.0+0.0im
# ```

# # See also

#     'hamiltonian`
# """
# velocity

# """
#     bravaismatrix(s)

# Return the Bravais matrix of `s`. Its columns are the Bravais vectors.
# """
# bravaismatrix(sys::System) = sys.lattice.bravais.matrix
# bravaismatrix(lat::Lattice) = lat.bravais.matrix
# bravaismatrix(br::Bravais) = br.matrix

# """
#     sitepositions(system[, sublat::Union{Integer, $NameType}[, siteindex::Int]])
#     sitepositions(system, (sublat, siteindex))

# Return a vector with all the site positions of `system`, or only those in `sublat`
# or with a given `siteindex`, if specified.
# """
# sitepositions(sys::System, s, i) = sys.lattice.sublats[sublatindex(sys, s)].sites[i]
# sitepositions(sys::System, (s, i)::Tuple) = sitepositions(sys, s, i)
# sitepositions(sys::System, s::Union{NameType,Integer}) = sys.lattice.sublats[sublatindex(sys, s)].sites
# sitepositions(sys::System) = [sl.sites for sl in sys.lattice.sublats]

# """
#     neighbors(system, sublat, siteindex; targetsublat = missing)
#     neighbors(system, (sublat, siteindex); targetsublat = missing)

# Return a vector with `(targetsublat, targetsiteindex)` of all the neighbors of input site
# (with `sublatindex` in sublattice `sublat`). A non-missing `targetsublat` restricts the search
# to a specific sublat.

# # Examples
# ```jldoctest
# julia> sys = System(LatticePresets.honeycomb(), Model(Hopping(1))) |> grow(region = Region(:square, 2))
# System{2,0,Float64,Complex{Float64}} : 0D system in 2D space
#   Bravais vectors     : ()
#   Sublattice names    : (:A, :B)
#   Sublattice orbitals : (1, 1)
#   Total sites         : 10 [Float64]
#   Total hoppings      : 50 [Complex{Float64}]
#   Coordination        : 5.0

# julia> neighbors(sys, 1, 4)
# 4-element Array{Tuple{Int64,Int64},1}:
#  (1, 2)
#  (1, 5)
#  (2, 3)
#  (2, 5)

# julia> neighbors(sys, 1, 4, targetsublat = :A)
# 2-element Array{Tuple{Int64,Int64},1}:
#  (1, 2)
#  (1, 5)

# ```
# """
# neighbors(sys, (sublat, siteindex)::Tuple; kw...) = neighbors(sys, sublat, siteindex; kw...)

# function neighbors(sys, sublat, siteindex; targetsublat = missing)
#     h = sys.hamiltonian.matrix
#     rows = rowvals(h)
#     source = torow(siteindex, sublat, sys.sysinfo)
#     targets = Tuple{Int,Int}[]
#     tsublatindex = targetsublat === missing ? 0 : sublatindex(sys, targetsublat)
#     for ptr in nzrange(h, source)
#         target, offset, tsublat = tosite(rows[ptr], sys.sysinfo)
#         offset == 0 && (targetsublat === missing || tsublatindex === tsublat) &&
#             (tsublat != sublat || target != source) && push!(targets, (tsublat, target))
#     end
#     return sort!(targets)
# end


text:
# Tests based on the VRT format tutorial found at http://www.gdal.org/gdal_vrttut.html

utmsmall = joinpath(@__DIR__, "data/utmsmall.tif")
utmsmall_vrt = joinpath(@__DIR__, "tmp/utmsmall.vrt")

# To create a VRT dataset that is a clone of an existing dataset use the CreateCopy() method.
driver = GDAL.gdalgetdriverbyname("VRT")
src_dataset = GDAL.gdalopenshared(utmsmall, GDAL.GA_ReadOnly)
progressfunc = convert(Ptr{GDAL.GDALProgressFunc}, C_NULL)
vrt_dataset =
    GDAL.gdalcreatecopy(driver, utmsmall_vrt, src_dataset, 0, C_NULL, progressfunc, C_NULL)

GDAL.gdalsetmetadataitem(vrt_dataset, "SourceAgency", "JuliaGeo", "")

GDAL.gdalclose(src_dataset)
GDAL.gdalclose(vrt_dataset)

# checking the contents of the VRT XML
vrt_xml = read(utmsmall_vrt, String)
@test occursin("<VRTDataset rasterXSize=\"100\" rasterYSize=\"100\">", vrt_xml)
@test occursin("<MDI key=\"SourceAgency\">JuliaGeo</MDI>", vrt_xml)
@test occursin("<SourceFilename relativeToVRT=\"0\">$utmsmall</SourceFilename>", vrt_xml)
rm(utmsmall_vrt)


text:
#debug fullLocalGraphCopy!
using Caesar, RoME, IncrementalInference
using CloudGraphs, Neo4j
# using Distributions
using JSON

using KernelDensityEstimate

using Base.Test

## Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true)
# session = addrdict["session"]
# Nparticles = parse(Int,addrdict["num particles"])

# interactive operation
session = "SESSROX"
Nparticles = 100
include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

fg = Caesar.initfg(sessionname=session, cloudgraph=cloudGraph)

# setDBAllReady!(fg)
setBackendWorkingSet!(session)

fullLocalGraphCopy!(fg)

IDs = getAllExVertexNeoIDs(conn, sessionname=fg.sessionname)

tosee = Dict()
tosee2 = Dict()
for i in 1:length(IDs)
  vid, nid = IDs[i]
  tosee[vid] = nid
  tosee2[nid] = vid
end

tosee2[149088]

writeGraphPdf(fg)
# @async run(`evince fg.pdf`)

tree = wipeBuildNewTree!(fg,drawpdf=true)
# @async run(`evince bt.pdf`)

spyCliqMat(tree, :l200050)

inferOverTreeR!(fg, tree, N=N, dbg=true)

# Some inspection methods to see what is going on
spyCliqMat(tree, :l200050)

plotMCMC(tree, :l200050)

plotPriorsAtCliq(tree, :l200050, :l200050, dims=[1;2], levels=1)

plotUpMsgsAtCliq(tree, :l200050, :l200050, marg=[1;2])

# basic plotting of Pose2D and landmarks
drawPosesLandms(fg ,spscale=0.5)

lbll=:l200050
cliq = whichCliq(tree, string(lbll))
cliqdbg = cliq.attributes["debug"]
vidx = 1
for lb in cliqdbg.mcmc[1].lbls
  if lb == lbll
    break;
  else
    vidx += 1
  end
end
ppp = kde!(cliqdbg.mcmc[3].prods[vidx].product)
plotKDE(ppp)

vc = startdefaultvisualization(draworigin=true)

drawposepoints!(vc, fg, :x4)

visualizeallposes!(vc, fg, drawlandms=false, drawtype=:fit)

drawPoses(fg,from=4,to=4)

# prior on :x1 issue

fc = getVert(fg, 100007)

vid = fg.IDs[:l200050]
vid = fg.IDs[:x4]

fc

Profile.clear()

@profile pts = evalFactor2(fg, fc, vid)
plotKDE(kde!(pts[1:2,:]))

getSample( getData(fc).fnc.usrfnc!, 100 )

plotKDE(getData(fc).fnc.usrfnc!.range)

#


text:
let
    d = SparseCat([:a, :b, :d], [0.4, 0.5, 0.1])
    c = collect(weighted_iterator(d))
    @test c == [:a=>0.4, :b=>0.5, :d=>0.1]
    @test pdf(d, :c) == 0.0
    @test pdf(d, :a) == 0.4
    @test mode(d) == :b
    @test Random.gentype(d) == Symbol
    @test Random.gentype(typeof(d)) == Symbol
    @inferred rand(Random.GLOBAL_RNG, d)

    dt = SparseCat((:a, :b, :d), (0.4, 0.5, 0.1))
    c = collect(weighted_iterator(dt))
    @test c == [:a=>0.4, :b=>0.5, :d=>0.1]
    @test pdf(dt, :c) == 0.0
    @test pdf(dt, :a) == 0.4
    @test mode(dt) == :b
    @test Random.gentype(dt) == Symbol
    @test Random.gentype(typeof(dt)) == Symbol
    @inferred rand(Random.GLOBAL_RNG, dt)

    # rand(::SparseCat)
    samples = Symbol[]
    N = 100_000
    @time for i in 1:N
        push!(samples, rand(d))
    end
    @test isapprox(count(samples.==:a)/N, pdf(d,:a), atol=0.005)
    @test isapprox(count(samples.==:b)/N, pdf(d,:b), atol=0.005)
    @test isapprox(count(samples.==:c)/N, pdf(d,:c), atol=0.005)
    @test isapprox(count(samples.==:d)/N, pdf(d,:d), atol=0.005)

    # rand(rng, ::SparseCat)
    rng = MersenneTwister(14)
    samples = Symbol[]
    N = 100_000
    @time for i in 1:N
        push!(samples, rand(rng, d))
    end
    @test isapprox(count(samples.==:a)/N, pdf(d,:a), atol=0.005)
    @test isapprox(count(samples.==:b)/N, pdf(d,:b), atol=0.005)
    @test isapprox(count(samples.==:c)/N, pdf(d,:c), atol=0.005)
    @test isapprox(count(samples.==:d)/N, pdf(d,:d), atol=0.005)

    @test_throws ErrorException rand(Random.GLOBAL_RNG, SparseCat([1], [0.0]))

    @test sprint((io,d)->show(io,MIME("text/plain"),d), d) == sprint((io,d)->showdistribution(io,d,title="SparseCat distribution"), d)
end