text:
using Pkg
Pkg.activate(pwd())
include("COAST.jl")
using .COAST

progress_test = [0.0]
try
    COAST.App.launchServer(progress_test,parse(Int, ARGS[1]))
catch e
    print("Not on server . . . attempting to start locally \n")
    COAST.App.launchServer(progress_test,8000)
end


text:
using Documenter, MendelIHT

makedocs(
    doctest = false, #will set to true later
    format = Documenter.HTML(),
    sitename = "MendelIHT",
    authors = "Benjamin Chu, Kevin Keys",
    clean = true,
    pages = [
        "Home"            => "index.md",
        "Getting Started" => "man/getting_started.md",
        "Examples"        => "man/examples.md",
        "Contributing"    => "man/contributing.md",
    ]
)

deploydocs(
    repo   = "github.com/biona001/MendelIHT.jl.git",
    target = "build",
  	deps   = nothing,
  	make   = nothing,
)

ext: jl


text:
module TestDLCRQ

using Base.Experimental: @sync
using ConcurrentCollections
using ConcurrentCollections.Implementations:
    MPCRQSlot, DATA, ANTIDATA, denqueue!, MPCRQ_ENQUEUED
using ProgressLogging: @logprogress, @withprogress
using Test
using ..Utils: ⊏

function test_MPCRQSlot()
    for index in [111, 222],
        safe in [false, true],
        polarity in [DATA, ANTIDATA],
        storage in UInt32[0xaaa, 0xbbb]

        @test NamedTuple(MPCRQSlot(; index, safe, polarity, storage)) ==
              (; index, safe, polarity, storage)
    end
end

function test_push_pop_once_int()
    q = DualLinkedConcurrentRingQueue{Int}()
    push!(q, 111)
    @test popfirst!(q) == 111
end

function test_push_pop_once_any()
    q = DualLinkedConcurrentRingQueue()
    push!(q, 111)
    @test popfirst!(q) == 111
end

function test_error()
    q = DualLinkedConcurrentRingQueue{Int}()
    t = @task popfirst!(q)
    yield(t)
    msg = "Interrupting DLCRQ @$(time_ns())"
    schedule(t, ErrorException(msg); error = true)
    err = try
        Some(wait(t))
    catch e
        e
    end
    @test err isa TaskFailedException
    @test occursin(msg, sprint(showerror, err))
    push!(q, 111)  # this dequeues but ignores the interrupted waiter
    @test popfirst!(q) == 111
end

function test_push_pop_100()
    n = 100
    q = DualLinkedConcurrentRingQueue{Int}(; log2ringsize = 3)
    foldl(push!, 1:n; init = q)
    ys = Int[]
    for _ in 1:n
        y = popfirst!(q)
        push!(ys, y)
    end
    @test ys == 1:n
end

function test_push_pop_100_wait_first()
    n = 100
    q = DualLinkedConcurrentRingQueue{Int}(; log2ringsize = 3)
    task = Threads.@spawn begin
        ys = Int[]
        for _ in 1:n
            y = popfirst!(q)
            push!(ys, y)
        end
        return ys
    end
    sleep(0.01)
    foldl(push!, 1:n; init = q)
    ys = fetch(task)
    @test ys == 1:n
end

function test_push_pop_100_inline()
    n = 100
    q = DualLinkedConcurrentRingQueue{Int16}(; log2ringsize = 3)
    @test q.data.data === nothing
    foldl(push!, 1:n; init = q)
    ys = Int[]
    for _ in 1:n
        y = popfirst!(q)
        push!(ys, y)
    end
    @test ys == 1:n
end

function unfair_sleep(seconds::Real)
    t0 = time_ns()
    ns = seconds * 1e9
    while time_ns() - t0 < ns
        GC.safepoint()
    end
end

function concurrent_push_pop!(q, nitems::Integer, nsend::Integer, nrecv::Integer)
    received = Vector{Int}[]
    tasks = Task[]
    global TASKS = tasks
    activesenders = Threads.Atomic{Int}(nsend)
    @sync begin
        for offset in 1:nsend
            t = Threads.@spawn begin
                for i in offset:nsend:nitems
                    push!(q, i)
                end
                if Threads.atomic_sub!(activesenders, 1) == 1
                    for _ in 1:nrecv
                        push!(q, -1)
                    end
                end
            end
            push!(tasks, t)
        end
        for _ in 1:nrecv
            ys = Int[]
            push!(received, ys)
            t = Threads.@spawn begin
                while true
                    y = popfirst!(q)
                    y == -1 && break
                    push!(ys, y)
                end
            end
            push!(tasks, t)
        end
    end
    return received
end

function check_consecutive(xs)
    notfound = Int[]
    dups = Int[]
    pre = xs[begin] - 1
    for x in xs
        e = pre + 1
        append!(notfound, e:x-1)
        append!(dups, x:e-1)
        pre = x
    end
    return (; notfound, dups)
end

function test_concurrent_push_pop(ntrials = 100)
    @withprogress name = "concurrent push-pop" begin
        @testset for trial in 1:ntrials
            @logprogress (trial - 1) / ntrials
            q = DualLinkedConcurrentRingQueue{Int}(; log2ringsize = 5)
            # q = Channel{Int}(Inf)
            check_concurrent_push_pop!(q)
        end
    end
end

function check_concurrent_push_pop!(q; nitems = 2^20)
    nsend = cld(Threads.nthreads(), 2)
    nrecv = max(1, Threads.nthreads() - nsend)
    received = concurrent_push_pop!(q, nitems, nsend, nrecv)
    allreceived = reduce(vcat, received)
    sort!(allreceived)
    (; notfound, dups) = check_consecutive(allreceived)
    @test length(allreceived) == nitems
    @test notfound == []
    @test dups == []
    @test allreceived == 1:nitems
end

function test_print()
    q = DualLinkedConcurrentRingQueue{Int}()
    str = sprint(show, "text/plain", q)
    @test "Dual LCRQ: " ⊏ str
    @test "0 data item" ⊏ str

    push!(q, 333)
    str = sprint(show, "text/plain", q)
    @test "Dual LCRQ: " ⊏ str
    @test "1 data item" ⊏ str

    popfirst!(q)
    t = @task popfirst!(q)
    yield(t)
    str = sprint(show, "text/plain", q)
    push!(q, 444)
    wait(t)
    @test "Dual LCRQ: " ⊏ str
    @test "1 waiter" ⊏ str
end

end  # module


text:
# QuEST.jl/src/_Tomog.jl
# (License stuff)

@doc raw"

Module `_Tomog`

Hosts the code for quantum process tomography.

## Exports
* Iterator type `StateTomog_Iter`
* Function `distance()` between two iterators

# Extended Help

## State tomography

State tomography is based on iterators.  That avoids having to store the whole matrix.
Julia iterators involve 3 concepts: Items, Iters, and State. Here:
* Items: Float64
* Iters: struct `StateTomog_Iter`
* State: UInt64 number representing the Pauli observable (4 bit per qubit)

## Process tomography
...
"
module _Tomog

export StateTomog_Iter
export distance

import ..QuEST
using ..QuEST: Qureg, numQubits, isDensityMatrix
using ..QuEST: QuESTEnv, createCloneQureg, destroyQureg
using ..QuEST: calcExpecPauliProd

import Base: iterate, eltype, length              # Iterator interface

import LinearAlgebra

using  Quadmath: Float128
const LngDbl = Float128

############################################################################
#                                                                          #
#   State tomography                                                       #
#                                                                          #
############################################################################
@doc raw"
Struct `StateTomog_Iter`

Iterator type for state tomograpy.
Constructor function:
    StateTomog_Iter( ::Qureg, ::QuESTEnv)
"
struct StateTomog_Iter
    qr        ::Qureg
    questenv  ::QuESTEnv
    qubits    ::Vector{Cint}   # = [0,1,2, ..., n-1] -- this will not change!

    # constructor
    function StateTomog_Iter(qr::Qureg, questenv::QuESTEnv)
        @assert numQubits(qr) ≤ 31

        local n = numQubits(qr)
        qubits = Cint[ j   for j =  0 : n-1 ]
        new(qr,questenv,qubits)
    end
end

struct StateTomog_State
    σ             ::Vector{Cint}
    count         ::UInt64
    workspace     ::Qureg
end
function StateTomog_State(iter::StateTomog_Iter) ::StateTomog_State
    workspace = createCloneQureg(iter.qr,iter.questenv)
    σ = [ Cint(3) for j = 1 : numQubits(iter.qr) ]
    return StateTomog_State(σ,0,workspace)
end

"""
Helper function `_inc(σ)`

Moves to the next multi-qubit Pauli (in lexicographic ordering)
"""
function _inc!(σ ::Vector{Cint}) ::Vector{Cint}
    @assert 0 < length(σ) ≤ 31
    local n = length(σ)

    for idx = 1:n
        @assert σ[idx] ∈ 0:3
        σ[idx] = ( σ[idx] + 1 ) % 4
        if σ[idx] != 0    break   end
    end
    return σ
end

function iterate(iter  ::StateTomog_Iter,
                 state ::StateTomog_State =StateTomog_State(iter)) ::Union{ Nothing, Tuple{Float64, StateTomog_State} }
    local n ::UInt64  = numQubits(iter.qr)
    if state.count ≥ (  1<<(2*n)  )
        destroyQureg(state.workspace,iter.questenv)
        return nothing
    else
        _inc!(state.σ)
        local expval = calcExpecPauliProd(iter.qr, iter.qubits,
                                          state.σ,
                                          state.workspace)
        return (  expval,
                  StateTomog_State(state.σ,
                                   state.count + 1,
                                   state.workspace)  )
    end
end

eltype(::Type{StateTomog_Iter})      =    Float64
length(iter ::StateTomog_Iter) ::Int =    1 << (2*length(iter.qubits))

@doc raw"
Function `distance(::StateTomog_Iter, ::StateTomog_Iter ; 𝑝=1)`

Computes the 𝑝-norm distance between the Pauli-basis vector representations of the two states given as state-tomography iterators.
"
function distance(ψ ::StateTomog_Iter,
                  ϕ ::StateTomog_Iter
                  ;
                  p ::Real = 1)             ::Float64
    @assert numQubits(ψ.qr) == numQubits(ϕ.qr)
    @assert 0 < p

    if p == Inf
        local maxval ::Float64 = 0
        for (e1,e2) ∈ zip(ψ,ϕ)
            maxval = max(maxval,
                         abs(e1 - e2) )
        end
        return maxval
    else
        local q   ::Float64 = 1/p
        local sum ::LngDbl  = 0
        for (e1,e2) ∈ zip(ψ,ϕ)
            sum += abs( LngDbl(e1)-LngDbl(e2) )^q
        end
        return Float64( sum^p )
    end
end

@doc raw"
Function `distance(ψ ::StateTomog_Iter, ϕ::Function ; 𝑝=1)`

Computes the 𝑝-norm distance between the Pauli-basis vector representations
of the state ψ and the vector defined by the function ϕ.

The function ϕ allow as argument a vector containing a Pauli code for each qubit.
It must return the expectation value of that observable for whatever state it represents.
"
function distance(ψ ::StateTomog_Iter,
                  ϕ ::Function
                  ;
                  p ::Real = 1)
    @assert 0 < p

    if p == Inf
        local maxval ::Float64 = 0
        local it = iterate(ψ)
        while it !== nothing
            (v,s) = it
            maxval = max(maxval,
                         abs(v - ϕ(s.σ) ))
            it = iterate(ψ,s)
        end
        return maxval
    else
        local q   ::Float64 = 1/p
        local sum ::LngDbl  = 0
        local it = iterate(ψ)
        while it !== nothing
            (v,s) = it
            sum += abs( LngDbl(v)-LngDbl(ϕ(s.σ)) )^q
            it = iterate(ψ,s)
        end
        return Float64( sum^p )
    end
end

############################################################################
#                                                                          #
#   Process tomography                                                     #
#                                                                          #
############################################################################

@doc raw"
Function `prepare!(qureg ::Qureg, 𝑥 ::UInt16) ::Qureg`

Prepares the given quantum register in the state |𝑥⟩.

### Requires
* 0 ≤ 𝑥 < 2ⁿ, where 𝑛 is the number of qubits in `qureg`
"
function prepare!(qureg ::Qureg, x ::UInt16) ::Qureg
    @assert 0 ≤ x < numQubits(qureg)

    @error "Not yet implemented"
end

@doc raw"
Function `prepareX!(qureg ::Qureg, 𝑥 ::UInt16, 𝑦 ::UInt16, minus::Bool) ::Qureg`

Prepares the given quantum register in the state ( |𝑥⟩ ± |𝑦⟩ )/√2, where the sign is determined by `minus`.

### Requires
* 0 ≤ 𝑥,y < 2ⁿ, where 𝑛 is the number of qubits in `qureg`
* 𝑥 ≠ 𝑦
"
function prepare!(qureg ::Qureg, x ::UInt16, y::UInt16, minus::Bool) ::Qureg
    @assert 0 ≤ x < numQubits(qureg)
    @assert 0 ≤ y < numQubits(qureg)
    @assert x ≠ y

    @error "Not yet implemented"
end

struct ProcTomog_Iter
    things
end

# function iterate()
# function eltype()
# function isequal()

end #^ module _Tomog
#EOF


text:
function GetShots(in)

	filename_h = join([filename ".seish"])
	stream = open(filename_h)
	sx = Array(Any)
	sy = Array(Any)
	key = ["sx","sy"]
	curr = zeros(length(key),1)
	prev = copy(curr)
	ntraces = int(filesize(stream)/(4*length(names(Header))))
	nsx = 1
	h1 = GrabHeader(stream,1)
	for ikey=1:length(key)
		curr[ikey] = getfield(h1,Symbol(key[ikey]))
	end
	prev = copy(curr)
	sx = push!(sx,h1.sx)
	sy = push!(sy,h1.sy)
	for j=2:ntraces
		h1 = GrabHeader(stream,j)
		for ikey=1:length(key)
			curr[ikey] = getfield(h1,Symbol(key[ikey]))
		end
		if curr != prev
			nx = j - itrace
			sx = push!(sx,h1.sx)
			sy = push!(sy,h1.sy)
		end
		prev = copy(curr)
	end
	close(stream)
	return sx,sy

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
module InlineTest

@static if isdefined(Base, :Experimental) &&
             isdefined(Base.Experimental, Symbol("@optlevel"))
    Base.Experimental.@optlevel 0
end

export @testset

# from Test:
export Test,
    @test, @test_throws, @test_broken, @test_skip,
    @test_warn, @test_nowarn,
    @test_logs, @test_deprecated,
    @inferred,
    detect_ambiguities, detect_unbound_args

using Test: Test,
    @test, @test_throws, @test_broken, @test_skip,
    @test_warn, @test_nowarn,
    @test_logs, @test_deprecated,
    @inferred,
    detect_ambiguities, detect_unbound_args

# it seems difficult to use gensym here, as other gensym calls (from other
# modules) can lead to collisions or mismatches, whether we statically
# initialize INLINE_TEST or at __init__ time; so let's use a long enough
# static random hexadecimal string
const INLINE_TEST = Symbol("##InlineTest-01b48f5c342f65df7fcd07f28f0d2cacbb09f0a0")

const TESTED_MODULES = Union{Module,Nothing}[]

get_tests(m::Module) = getfield(m, INLINE_TEST).tests
register(m::Module) = push!(TESTED_MODULES, m)

function retest end

"""
    @testset args...

Similar to `Test.@testset args...`, but the contained tests are not run
immediately, and are instead stored for later execution, triggered by
[`retest()`](@ref) or `runtests()`.
Invocations of `@testset` can be nested, but qualified invocations of
`ReTest.@testset` can't.
Internally, `@testset` expressions are converted to an equivalent of
`Test.@testset` at execution time.
"""
macro testset(args...)
    let mod = __module__
        while !isdefined(mod, :runtests)
            @eval mod begin
                """
                    $($mod).runtests(pattern...; kwargs...)

                Equivalent to `ReTest.retest($($mod), pattern...; kwargs...)`.
                This function is defined automatically in any module containing
                a `@testset`, possibly nested within submodules.
                """
                function runtests(specs...; kwargs...)
                    $retest($mod, specs...; kwargs...)
                end
            end
            mod = parentmodule(mod)
        end
    end
    if !isdefined(__module__, INLINE_TEST)
        # Credit to Takafumi Arakaki for the idea of creating a submodule
        # within modules using InlineTest in order to be able to define
        # an `__init__` function which does the registration
        # (as we can't directly create/modify `__init__` for the given module).
        @eval __module__ module $INLINE_TEST
            const tests = (tests=[], news=[], map=Dict{Union{String,Expr},Int}())
            __init__() = $register($__module__)
        end
    end
    # this must take effect at compile/run time rather than parse time, e.g.
    # if the @testset if in a `if false` branch
    # TODO: test that
    quote
        push!(get_tests($__module__).news, (ts=$args, source=$(QuoteNode(__source__))))
        nothing
    end
end

end # InlineTest


text:
##############################################################
# Distributed Multinomial Regression (DMR)
##############################################################

"Abstract Distributed Counts Regression (DCR) returned object"
abstract type DCR <: RegressionModel end

"Abstract DMR returned object"
abstract type DMR <: DCR end

"""
Relatively heavy object used to return DMR results when we care about the regulatrization paths.
"""
struct DMRPaths <: DMR
  nlpaths::Vector{Union{Missing,GammaLassoPath}} # independent Poisson GammaLassoPath for each phrase
  intercept::Bool               # whether to include an intercept in each Poisson regression
                                # (only kept with remote cluster, not with local cluster)
  n::Int                      # number of observations. May be lower than provided after removing all zero obs.
  d::Int                      # number of categories (terms/words/phrases)
  p::Int                      # number of covariates
end

const defsegselect = MinAICc()

"""
Relatively light object used to return DMR results when we only care about estimated coefficients.
"""
struct DMRCoefs{T<:AbstractMatrix, S<:SegSelect} <: DMR
  coefs::T                    # model coefficients
  intercept::Bool             # whether to include an intercept in each Poisson regression
  n::Int                      # number of observations. May be lower than provided after removing all zero obs.
  d::Int                      # number of categories (terms/words/phrases)
  p::Int                      # number of covariates
  select::S                   # path segment selector
end

function DMRCoefs(m::DMRPaths, select::SegSelect=defsegselect)
  coefs = coef(m, select)
  DMRCoefs(coefs, m.intercept, m.n, m.d, m.p, select)
end

"""
    fit(DMR,covars,counts; <keyword arguments>)
    dmr(covars,counts; <keyword arguments>)

Fit a Distributed Multinomial Regression (DMR) of counts on covars.

DMR fits independent poisson gamma lasso regressions to each column of counts to
approximate a multinomial, picks a segement of each path, and
returns a coefficient matrix (wrapped in DMRCoefs) representing point estimates
for the entire multinomial (includes the intercept if one was included).

# Example:
```julia
  m = fit(DMR,covars,counts)
```
# Arguments
- `covars` n-by-p matrix of covariates
- `counts` n-by-d matrix of counts (usually sparse)

# Keywords
- `intercept::Bool=false` include an intercept in each poisson
- `parallel::Bool=true` parallelize the poisson fits
- `local_cluster::Bool=true` use local_cluster mode that shares memory across
    parallel workers that is appropriate on a single multicore machine, or
    remote cluster mode that is more appropriate when distributing across machines
    for which sharing memory is costly.
- `verbose::Bool=true`
- `showwarnings::Bool=false`
- `select::SegSelect=MinAICc()` which path segment to pick
- `kwargs...` additional keyword arguments passed along to fit(GammaLassoPath,...)
"""
function StatsBase.fit(::Type{D}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
  kwargs...) where {T<:AbstractFloat, D<:DMR}

  dmr(covars, counts; kwargs...)
end

"""
    fit(DMRPaths,covars,counts; <keyword arguments>)
    dmrpaths(covars,counts; <keyword arguments>)

Fit a Distributed Multinomial Regression (DMR) of counts on covars, and returns
the entire regulatrization paths, which may be useful for plotting or picking
coefficients other than the AICc optimal ones. Same arguments as [`fit(::DMR)`](@ref).
"""
function StatsBase.fit(::Type{DMRPaths}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
  local_cluster=false, # ignored. will always assume remote_cluster
  kwargs...) where {T<:AbstractFloat}

  dmrpaths(covars, counts; kwargs...)
end

"""
    fit(DMR,@model(c ~ x1*x2),df,counts; <keyword arguments>)

Fits a DMR but takes a model formula and dataframe instead of the covars matrix.
See also [`fit(::DMR)`](@ref).

`c` must be specified on the lhs to indicate the model for counts.
"""
function StatsBase.fit(::Type{T}, m::Model, df, counts::AbstractMatrix;
  contrasts::Dict{Symbol,<:Any} = Dict{Symbol,Any}(), kwargs...) where {T<:DMR}

  # parse and merge rhs terms
  trms = getrhsterms(m, :c)

  # create model matrix
  covars, counts, as = modelcols(trms, df, counts; model=T, contrasts=contrasts)

  # fit and wrap in TableCountsRegressionModel
  TableCountsRegressionModel(fit(T, covars, counts; kwargs...), df, counts, m, as)
end

"""
    coef(m::DMRCoefs)

Returns the coefficient matrices fitted with DMR using
the segment selected during fit (MinAICc by default).

# Example:
```julia
  m = fit(DMR,covars,counts)
  coef(m)
```
"""
StatsBase.coef(m::DMRCoefs) = m.coefs

coefspace(p, d, nλ, select::SegSelect) = zeros(p,d)
coefspace(p, d, nλ, select::AllSeg) = zeros(nλ,p,d)

coeffill!(coefs, path::Missing, p, j, select::SegSelect) = nothing

function coeffill!(coefs, path::RegularizationPath, p, j, select::SegSelect)
  cj = coef(path, select)
  coeffill!(coefs, cj, p, j, select)
end

function coeffill!(coefs, cj::AbstractArray, p, j, select::SegSelect)
  for i=1:p
    coefs[i,j] = cj[i]
  end
end
function coeffill!(coefs, cj::AbstractArray, p, j, select::AllSeg)
  for i=1:p
    for s=1:size(cj,2)
      coefs[s,i,j] = cj[i,s]
    end
  end
end

"""
    coef(m::DMRPaths, select::SegSelect=MinAICc())

Returns all or selected coefficients matrix fitted with DMR.

# Example:
```julia
  m = fit(DMRPaths,covars,counts)
  coef(m, MinCVKfold{MinCVmse}(5))
```
"""
# StatsBase.coef(m::DMRPaths; select=defsegselect) = coef(m, select)
function StatsBase.coef(m::DMRPaths, select::SegSelect=defsegselect)
  # get dims
  d = length(m.nlpaths)
  d < 1 && return nothing

  # drop missing paths
  nonmsngpaths = skipmissing(m.nlpaths)

  # get number of coefs from paths object
  p = ncoefs(m)

  # establish maximum path lengths
  nλ = 0
  if !isempty(nonmsngpaths)
    nλ = maximum(size(nlpath)[2] for nlpath in nonmsngpaths)
  end

  # allocate space
  coefs = coefspace(p, d, nλ, select)

  # iterate over paths
  for j=1:d
    path = m.nlpaths[j]
    coeffill!(coefs, path, p, j, select)
  end

  coefs
end

"Whether the model includes an intercept in each independent counts (e.g. hurdle) regression"
hasintercept(m::DCR) = m.intercept

"Number of observations used. May be lower than provided after removing all zero obs."
StatsBase.nobs(m::DCR) = m.n

"Number of categories (terms/words/phrases) used for DMR estimation"
Distributions.ncategories(m::DCR) = m.d

"Number of covariates used for DMR estimation"
ncovars(m::DMR) = m.p

"Number of coefficient potentially including intercept used for each independent Poisson regression"
ncoefs(m::DMR) = ncovars(m) + (hasintercept(m) ? 1 : 0)

# some helpers for converting to SharedArray
Base.convert(::Type{SharedArray}, A::SubArray) = (S = SharedArray{eltype(A)}(size(A)); copyto!(S, A))
function Base.convert(::Type{SharedArray}, A::SparseMatrixCSC{T,N}) where {T,N}
  S = SharedArray{T}(size(A))
  fill!(S,zero(T))
  rows = rowvals(A)
  vals = nonzeros(A)
  n, m = size(A)
  for j = 1:m
     for i in nzrange(A, j)
        row = rows[i]
        val = vals[i]
        S[row,j] = val
     end
  end
  S
end

"convert counts matrix elements to Float64 if necessary"
fpcounts(counts::M) where {V, N, M<:SparseMatrixCSC{V,N}} = convert(SparseMatrixCSC{Float64,N}, counts)
fpcounts(counts::M) where {V, M<:AbstractMatrix{V}} = Matrix{Float64}(counts)
fpcounts(counts::M) where {V<:GLM.FP, N, M<:SparseMatrixCSC{V,N}} = counts
fpcounts(counts::M) where {V<:GLM.FP, M<:AbstractMatrix{V}} = counts

totalcounts(counts, prespecifiedm::Nothing) = vec(sum(counts, dims=2))
totalcounts(counts, prespecifiedm::AbstractVector) = convert(Vector{GLM.FP}, prespecifiedm)

"""
Computes DMR shifters (μ=log(m)) and removes all zero observations.
Optionally uses a prespecifed total counts `m`, to allow computation in batches.
"""
function shifters(::Type{DMR}, covars::AbstractMatrix, counts::AbstractMatrix{C}, showwarnings::Bool,
  prespecifiedm::Union{Nothing, AbstractVector}) where C

  # standardize counts matrix to conform to GLM.FP
  counts = fpcounts(counts)

  m = totalcounts(counts, prespecifiedm)

  if any(iszero,m)
      # omit observations with no counts
      ixposm = findall(x->x!=zero(C), m)
      showwarnings && @warn("omitting $(length(m)-length(ixposm)) observations with no counts")
      m = m[ixposm]
      counts = counts[ixposm,:]
      covars = covars[ixposm,:]
  end

  μ = log.(m)

  n = length(m)

  covars, counts, μ, n
end

"Destandardize coefficents estimated using a potentially standardized covars matrix"
function destandardize!(coefs::AbstractMatrix{T}, covarsnorm::AbstractVector{T},
    standardize, intercept) where T

  if standardize
    if intercept
      covarsnorm = [one(T); covarsnorm]
    end

    # destandardize all coefs only once
    lmul!(Diagonal(covarsnorm), coefs)
  end
  nothing
end

"Destandardize path.coefs estimated using a potentially standardized covars matrix"
function destandardize!(path::RegularizationPath{S,T}, covarsnorm::AbstractVector{T},
    standardize) where {S,T}

  if standardize && isdefined(path,:coefs)
    # destandardize all coefs only once
    lmul!(Diagonal(covarsnorm), path.coefs)
    # not really used but set just in case it is in the future
    path.Xnorm = covarsnorm
  end

  path
end

"""
This version is built for local clusters and shares memory used by both inputs and outputs if run in parallel mode.
"""
function dmr_local_cluster(covars::AbstractMatrix{T},counts::AbstractMatrix{V},
          parallel,verbose,showwarnings,intercept; select=defsegselect, m=nothing,
          standardize=true, kwargs...) where {T<:AbstractFloat,V}
  # get dimensions
  n, d = size(counts)
  n1,p = size(covars)
  @assert n==n1 "counts and covars should have the same number of observations"
  verbose && @info("fitting $n observations on $d categories, $p covariates ")

  # add one coef for intercept
  ncoef = p + (intercept ? 1 : 0)

  covars, counts, μ, n = shifters(DMR, covars, counts, showwarnings, m)

  # standardize covars only once if needed
  covars, covarsnorm = Lasso.standardizeX(covars, standardize)

  # fit separate GammaLassoPath's to each dimension of counts j=1:d and pick its min AICc segment
  if parallel
    verbose && @info("distributed poisson run on local cluster with $(nworkers()) nodes")
    scounts = convert(SharedArray,counts)
    scoefs = SharedMatrix{T}(ncoef,d)
    scovars = convert(SharedArray,covars)
    # μ = convert(SharedArray,μ) incompatible with GLM

    @sync @distributed for j=1:d
      tryfitgl!(scoefs, j, scovars, scounts; offset=μ, verbose=false,
        showwarnings=showwarnings, intercept=intercept, select=select,
        standardize=false, kwargs...)
    end

    coefs = convert(Matrix{T}, scoefs)
  else
    verbose && @info("serial poisson run on a single node")
    coefs = Matrix{T}(undef,ncoef,d)
    for j=1:d
      tryfitgl!(coefs, j, covars, counts; offset=μ, verbose=false,
        showwarnings=showwarnings, intercept=intercept, select=select,
        standardize=false, kwargs...)
    end
  end

  # destandardize coefs only once if needed
  destandardize!(coefs, covarsnorm, standardize, intercept)

  DMRCoefs(coefs, intercept, n, d, p, select)
end

"This version does not share memory across workers, so may be more efficient for small problems, or on remote clusters."
function dmr_remote_cluster(covars::AbstractMatrix{T},counts::AbstractMatrix{V},
          parallel,verbose,showwarnings,intercept; select=defsegselect, kwargs...) where {T<:AbstractFloat,V}
  paths = dmrpaths(covars, counts; parallel=parallel, verbose=verbose, showwarnings=showwarnings, intercept=intercept, kwargs...)
  DMRCoefs(paths, select)
end

"Shorthand for fit(DMRPaths,covars,counts). See also [`fit(::DMRPaths)`](@ref)"
function dmrpaths(covars::AbstractMatrix{T},counts::AbstractMatrix;
      intercept=true,
      parallel=true,
      verbose=true, showwarnings=false,
      m=nothing,
      standardize=true,
      kwargs...) where {T<:AbstractFloat}
  # get dimensions
  n, d = size(counts)
  n1,p = size(covars)
  @assert n==n1 "counts and covars should have the same number of observations"
  verbose && @info("fitting $n observations on $d categories, $p covariates ")

  covars, counts, μ, n = shifters(DMR, covars, counts, showwarnings, m)

  # standardize covars only once if needed
  covars, covarsnorm = Lasso.standardizeX(covars, standardize)

  function tryfitgl(countsj::AbstractVector)
    try
      # we make it dense remotely to reduce communication costs
      path = fit(GammaLassoPath,covars,Vector(countsj),Poisson(),LogLink(); offset=μ, standardize=false, verbose=false, kwargs...)
      destandardize!(path, covarsnorm, standardize)
    catch e
      showwarnings && @warn("fitgl failed for countsj with frequencies $(countmapstr(countsj)) and will return missing path ($e)")
      missing
    end
  end

  # counts generator
  countscols = (counts[:,j] for j=1:d)

  if parallel
    verbose && @info("distributed poisson run on remote cluster with $(nworkers()) nodes")
    mapfn = pmap
  else
    verbose && @info("serial poisson run on a single node")
    mapfn = map
  end

  nlpaths = allowmissing(mapfn(tryfitgl,countscols))

  DMRPaths(nlpaths, intercept, n, d, p)
end

"Fits a regularized poisson regression counts[:,j] ~ covars saving the coefficients in coefs[:,j]"
function poisson_regression!(coefs::AbstractMatrix{T}, j::Int, covars::AbstractMatrix{T},counts::AbstractMatrix{V};
  select=defsegselect, kwargs...) where {T<:AbstractFloat,V}
  cj = Vector(counts[:,j])
  path = fit(GammaLassoPath,covars,cj,Poisson(),LogLink(); kwargs...)
  coefs[:,j] = coef(path, select)
  nothing
end

"Wrapper for poisson_regression! that catches exceptions in which case it sets coefs to zero"
function tryfitgl!(coefs::AbstractMatrix{T}, j::Int, covars::AbstractMatrix{T},counts::AbstractMatrix{V};
            showwarnings = false,
            kwargs...) where {T<:AbstractFloat,V}
  try
    poisson_regression!(coefs, j, covars, counts; kwargs...)
  catch e
    showwarnings && @warn("fitgl! failed on count dimension $j with frequencies $(countmapstr(counts[:,j])) and will return zero coefs ($e)")
    # zero out coef column because coefs array is uninitialized at first
    for i=1:size(coefs,1)
      coefs[i,j] = zero(T)
    end
  end
end

"Shorthand for fit(DMR,covars,counts). See also [`fit(::DMR)`](@ref)"
function dmr(covars::AbstractMatrix{T},counts::AbstractMatrix{V};
          intercept=true,
          parallel=true, local_cluster=true,
          verbose=true, showwarnings=false,
          kwargs...) where {T<:AbstractFloat,V}
  if local_cluster || !parallel
    dmr_local_cluster(covars,counts,parallel,verbose,showwarnings,intercept; kwargs...)
  else
    dmr_remote_cluster(covars,counts,parallel,verbose,showwarnings,intercept; kwargs...)
  end
end

# We take care of the intercept ourselves, without relying on StatsModels, because
# it is unregulated, so we drop it from formula
StatsModels.drop_intercept(::Type{T}) where {T<:DCR} = true

# delegate f(m::TableCountsRegressionModel,...) to f(m.model,...)
StatsModels.@delegate TableCountsRegressionModel.model [hasintercept,
  Distributions.ncategories, ncovars, ncoefs, ncovarszero, ncovarspos,
  ncoefszero, ncoefspos,
  StatsBase.coef, StatsBase.confint,
  StatsBase.deviance, StatsBase.nulldeviance,
  StatsBase.loglikelihood, StatsBase.nullloglikelihood,
  StatsBase.dof, StatsBase.dof_residual, StatsBase.nobs,
  StatsBase.stderror, StatsBase.vcov,
  StatsBase.residuals, StatsBase.response,
  StatsBase.predict, StatsBase.predict!]

"""
    predict(m,newcovars; <keyword arguments>)

Predict counts using a fitted DMRPaths object and given newcovars.

# Example:
```julia
  m = fit(DMRPaths,covars,counts)
  newcovars = covars[1:10,:]
  countshat = predict(m, newcovars; select=MinAICc())
```

# Arguments
- `m::DMRPaths` fitted DMRPaths model (DMRCoefs currently not supported)
- `newcovars` n-by-p matrix of covariates of same dimensions used to fit m.

# Keywords
- `select=MinAICc()` See [`coef(::RegularizationPath)`](@ref).
- `kwargs...` additional keyword arguments passed along to predict() for each
  category j=1..size(counts,2)
"""
function StatsBase.predict(m::DMRPaths, newcovars::AbstractMatrix{T};
  select=defsegselect, kwargs...) where {T<:AbstractFloat}

  _predict(m,newcovars;select=select,kwargs...)
end

# internal mothod used by both dmr and hdmr
function _predict(m, newcovars::AbstractMatrix{T};
  select=defsegselect, kwargs...) where {T<:AbstractFloat}

  @assert !isa(select, AllSeg) "select cannot be AllSeg and must choose a particular segment (e.g. select=MinAICc())"

  # dimensions
  newn = size(newcovars,1)

  # offset does not matter here because we rescale in the end
  newoffset = zeros(T,size(newcovars,1))

  η = zeros(T,newn,m.d)
  for j=1:m.d
    path = m.nlpaths[j]
    if !ismissing(path)
      η[:,j] = predict(path, newcovars;offset=newoffset, select=select, kwargs...)
    end
  end

  lmul!(Diagonal(one(T)./vec(sum(η, dims=2))),η)

  η
end

function StatsBase.predict(m::M, newcovars::AbstractMatrix{T};
  select=defsegselect, kwargs...) where {T<:AbstractFloat, M<:DMR}

  error("predict(m::DMR,...) can currently only be evaluated for DMRPaths structs returned from fit(DMRPaths,...)")
end


text:
module TestAmbiguity
using Test

@testset "method ambiguity" begin
    # Ambiguity test is run inside a clean process.
    # https://github.com/JuliaLang/julia/issues/28804
    script = joinpath(@__DIR__, "detect_ambiguities.jl")
    code = """
    $(Base.load_path_setup_code())
    include($(repr(script)))
    """
    cmd = Base.julia_cmd()
    if Base.JLOptions().color == 1
        cmd = `$cmd --color=yes`
    end
    cmd = `$cmd --startup-file=no -e $code`
    @test success(pipeline(cmd; stdout=stdout, stderr=stderr))
end

end  # module


text:
"""
    struct DiscreteSpace <: AbstractDiscreteSpace
        n::Int
    end

The elements in a `DiscreteSpace` is `1:n`
"""
struct DiscreteSpace <: AbstractDiscreteSpace
    n::Int
end

size(d::DiscreteSpace) = (d.n,)
eltype(d::DiscreteSpace) = Int
sample(d::DiscreteSpace) = rand(1 : d.n)
in(x::Int, d::DiscreteSpace) = 1 ≤ x ≤ d.n


text:
module TapFilters

# 1nd order
p = [0.037659  0.249153  0.426375  0.249153  0.037659]
d1 =[0.109604  0.276691  0.000000 -0.276691 -0.109604]
dx = p'*d1
dy = d1'*p

# 2nd order
p  = [0.030320  0.249724  0.439911  0.249724  0.030320]
d1 = [0.104550  0.292315  0.000000 -0.292315 -0.104550]
d2 = [0.232905  0.002668 -0.471147  0.002668  0.232905]
dxx = p'*d2
dyy = d2'*p
dxy = p'*d1*d1'*p

export dx, dy, dxx, dyy, dxy

end


text:
using Mamba
include("BayesTest.jl")
import Base.mean, Base.median

function ppv(samples::Mamba.ModelChains)
    predict(samples) # TODO: implement posterior predictive p values
end


text:
using Test
using NumberTheory

@testset "Mersenne Numbers" begin
    # Tests for the isPowerOfTwoFunction
    @test isPowerOfTwo(4) == true
    @test isPowerOfTwo(256) == true
    @test isPowerOfTwo(1024) == true
    @test isPowerOfTwo(1073741824) == true
    @test isPowerOfTwo(1073741825) == false
    @test isPowerOfTwo(1073741823) == false
    @test isPowerOfTwo(1) == true
    @test isPowerOfTwo(0) == false

    # Tests for the isMersenneNumber
    @test isMersenneNumber(2) == false
    @test isMersenneNumber(1073741824) == false
    @test isMersenneNumber(1073741823) == true
    @test isMersenneNumber(131071) == true

    # Tests for the getMersenneNumber
    @test getMersenneNumber(17) == 131071
    @test getMersenneNumber(1) == 1
    @test getMersenneNumber(5) == 31
end

max_stars_repo_name: EllianCarlos/number-theory-but-in-julia

max_issues_repo_name: EllianCarlos/number-theory-but-in-julia

max_forks_repo_name: EllianCarlos/number-theory-but-in-julia


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())
 cd("NeidSolarScripts")
 using Pkg
 Pkg.activate(".")
end

using CSV, DataFrames, Query
 using JLD2, FileIO
 using Statistics, Dates
 using RvSpectMLBase, RvSpectML
 using Plots

data_to_read = ["times", "rvs_ccf_espresso", "σ_rvs_espresso",
                "times_binned", "rvs_binned", "rv_rms",
                "rv_0", "rv_mean", "rv_slope",
                "order_snr", "order_rvs_g",
                "v_grid", "ccfs_espresso", "ccf_vars_espresso"]
 df = DataFrame(:filename=>readdir("output")) |> @filter(endswith(_.filename,".jld2")) |> DataFrame

 for field in data_to_read
   df[!,Symbol(field)] = map(i->load(joinpath("output",i.filename),field),eachrow(df))
 end
 df[!,"jd"] = map(i->datetime2julian(DateTime(parse(Int,SubString(df.filename[i],7:10)), parse(Int,SubString(df.filename[i],11:12)),
               parse(Int,SubString(df.filename[i],13:14)))), 1:size(df,1) )
 #df = df |> @filter( ! (2459213.4 < _.jd < 2459213.6) ) |> DataFrame
 #df = df |> @filter( ! (2459224.4 < _.jd < 2459224.6) ) |> DataFrame
 t0 = datetime2julian(DateTime(2021,1,22))

plt = scatter(df[!,"jd"].-t0, first.(df[!,"rv_0"]), xlabel="Day", ylabel="RV fit at noon (m/s)", label=:none)
 savefig("rv_solar_noon_vs_t_all.png")

df = df |> @filter( ! (2459213.4 < _.jd < 2459213.6) ) |> DataFrame
  df = df |> @filter( ! (2459224.4 < _.jd < 2459224.6) ) |> DataFrame
  plt = scatter(df[!,"jd"].-t0, first.(df[!,"rv_0"]), xlabel="Day", ylabel="RV fit at noon (m/s)", label=:none)
  first_obs_to_fit = 3
  x = [df[!,"jd"][first_obs_to_fit:end].-t0  ones(length(df[!,"jd"][first_obs_to_fit:end])) ]
  y = first.(df[!,"rv_0"])[first_obs_to_fit:end]
  linfit = (x'*x)\x'*y
  y_fit = (x*linfit)
  slope = linfit[1]
  rms = std(y.-y_fit)
  plot!(plt,x[:,1],y_fit, label=:none)
  title!(plt,"slope = " * string(round(slope,digits=2)) * " m/s/day  rms = " * string(round(rms,digits=2)) * " m/s")
  savefig("rv_solar_noon_vs_t.png")
  #ylims!(plt, 132,143)

plt = scatter(df[!,"jd"].-t0, first.(df[!,"rv_slope"]), xlabel="Day", ylabel="RV slope (m/s/hr)", label=:none)
 savefig("rv_slope_vs_t.png")

plt1 = plot()
 plt2 = plot()
 for i in 1:size(df,1)
    plot!(plt1,df[i,"times"].*24, df[i,"order_snr"][60,:], label=string(df[i,"jd"]-t0), color=i, legend=:none)
    plot!(plt2,df[i,"times"].*24, df[i,"order_snr"][20,:]./df[i,"order_snr"][80,:] .- maximum(df[i,"order_snr"][20,:]./df[i,"order_snr"][80,:]), color=i, label=string(df[i,"jd"]-t0), legend=:none)
 end
 ylims!(plt1,2.8e4,3.2e4)
 plt = plot(plt1,plt2, layout = (2,1))
 xlabel!(plt,"Time (hr - solar noon)")
 ylabel!(plt1,"SNR in Order 60")
 display(plt1)
 savefig("snr_vs_t_day.png")
 ylabel!(plt2,"Δ (SNR_20/SNR_80)")
 display(plt)
 savefig("snr_color_vs_t_day.png")

daily_ccfs = zeros(size(df[1,"ccfs_espresso"],1),size(df,1))
for i in 1:size(daily_ccfs,2)
  idx_to_use = -1.5 .<= df[i,"times"]*24.0 .<= 1.5
  daily_ccfs[:,i] .= vec(mean(df[i,"ccfs_espresso"][:,idx_to_use],dims=2))
end

obsid_start = 3
obsid_stop = size(daily_ccfs,2)-1
 daily_ccfs_norm = daily_ccfs[:,obsid_start:obsid_stop]./maximum(daily_ccfs[:,obsid_start:obsid_stop],dims=1)
 ccf_template = vec(mean(daily_ccfs_norm,dims=2))
 plot(daily_ccfs_norm.-ccf_template)
 heatmap(df[1,"v_grid"],1:length(df[!,"jd"][obsid_start:obsid_stop]),daily_ccfs_norm'.-ccf_template')
 title!("Residuals of Mid-day mean CCF")
 xlabel!("v (m/s)")
 ylabel!("Obs ID")
 xlims!(-30e3,30e3)
 savefig("daily_ccf_resid.png")


text:
# BioSequence
# ===========
#
# A general purpose biological sequence representation.
#
# It is cheap to create a subsequence from a sequence because sequences can
# share underlying data: creating a subsequence is just allocating a new
# sequence object defining a part of the original sequence without copy.
# Destructive operations create a copy of the underlying data if and only if it
# is shared between (sub)sequences. This is often called as copy-on-write
# strategy in computer science and should be transparent to the user.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# About Internals
# ---------------
#
# The `data` field of a `BioSequence{A}` object contains binary representation
# of a biological character sequence. Each character is encoded with an encoder
# corresponding to the alphabet `A` and compactly packed into `data`. To extract
# a character from a sequence, you should decode this binary sequence with a
# decoder that is a pair of the encoder. The length of encoded binary bits is
# fixed, and hence a character at arbitrary position can be extracted in a
# constant time. To know the exact location of a character at a position, you
# can use the `bitsof(seq, i)` function, which returns a pair of element's index
# containing binary bits and bits' offset. As a whole, character extraction
# `seq[i]` can be written as:
#
#     j = bitindex(seq, i)
#     decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
#
#  index :        index(j) - 1       index(j)       index(j) + 1
#   data : ....|xxxxx...........|xxXxxxxxxxxxxxxx|............xxxx|....
# offset :                          |<-offset(j)-|
#  width :      |<---- 64 ---->| |<---- 64 ---->| |<---- 64 ---->|
#
#  * '.' : unused (4 bits/char)
#  * 'x' : used
#  * 'X' : used and pointed by index `i`

"""
Biological sequence data structure indexed by an alphabet type `A`.
"""
type BioSequence{A<:Alphabet} <: Sequence
    data::Vector{UInt64}  # encoded character sequence data
    part::UnitRange{Int}  # interval within `data` defining the (sub)sequence
    shared::Bool          # true if and only if `data` is shared between sequences

    function BioSequence(data::Vector{UInt64}, part::UnitRange{Int}, shared::Bool)
        return new(data, part, shared)
    end
end

typealias DNASequence       BioSequence{DNAAlphabet{4}}
typealias RNASequence       BioSequence{RNAAlphabet{4}}
typealias AminoAcidSequence BioSequence{AminoAcidAlphabet}
typealias CharSequence      BioSequence{CharAlphabet}

# String Decorators
# -----------------

remove_newlines(s) = replace(s, r"\r|\n", "")

macro dna_str(seq, flags...)
    return DNASequence(remove_newlines(seq))
end

macro rna_str(seq, flags...)
    return RNASequence(remove_newlines(seq))
end

macro aa_str(seq, flags...)
    return AminoAcidSequence(remove_newlines(seq))
end

macro char_str(seq, flags...)
    return CharSequence(remove_newlines(seq))
end

# Constructors
# ------------

function seq_data_len{A}(::Type{A}, len::Integer)
    return cld(len, div(64, bitsof(A)))
end

function (::Type{BioSequence{A}}){A<:Alphabet}(len::Integer)
    return BioSequence{A}(Vector{UInt64}(seq_data_len(A, len)), 1:len, false)
end

BioSequence(::Type{DNANucleotide}) = DNASequence()
BioSequence(::Type{RNANucleotide}) = RNASequence()
BioSequence(::Type{AminoAcid}) = AminoAcidSequence()
BioSequence(::Type{Char}) = CharSequence()

function BioSequence()
    return BioSequence{VoidAlphabet}(Vector{UInt64}(0), 0:-1, false)
end

function (::Type{BioSequence{A}}){A<:Alphabet}(
        src::Union{AbstractString,AbstractVector},
        startpos::Integer=1,
        stoppos::Integer=length(src))
    len = stoppos - startpos + 1
    seq = BioSequence{A}(len)
    return encode_copy!(seq, 1, src, startpos, len)
end

# create a subsequence
function BioSequence{A,T<:Integer}(other::BioSequence{A}, part::UnitRange{T})
    checkbounds(other, part)
    start = other.part.start + part.start - 1
    stop = start + length(part) - 1
    subseq = BioSequence{A}(other.data, start:stop, true)
    other.shared = true
    return subseq
end

function (::Type{BioSequence{A}}){A}(other::BioSequence{A}, part::UnitRange)
    return BioSequence(other, part)
end

# concatenate chunks
function (::Type{BioSequence{A}}){A}(chunks::BioSequence{A}...)
    len = 0
    for chunk in chunks
        len += length(chunk)
    end
    seq = BioSequence{A}(len)
    offset = 1
    for chunk in chunks
        copy!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

function Base.repeat{A}(chunk::BioSequence{A}, n::Integer)
    seq = BioSequence{A}(length(chunk) * n)
    offset = 1
    for _ in 1:n
        copy!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

# operators for concat and repeat
Base.:*{A}(chunk::BioSequence{A}, chunks::BioSequence{A}...) =
    BioSequence{A}(chunk, chunks...)
Base.:^(chunk::BioSequence, n::Integer) = repeat(chunk, n)

# conversion between different alphabet size
for A in [DNAAlphabet, RNAAlphabet]
    # 4 bits => 2 bits
    @eval function Base.convert(::Type{BioSequence{$(A{2})}}, seq::BioSequence{$(A{4})})
        # TODO: make it faster with bit-parallel algorithm
        newseq = BioSequence{$(A{2})}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end

    # 2 bits => 4 bits
    @eval function Base.convert(::Type{BioSequence{$(A{4})}}, seq::BioSequence{$(A{2})})
        newseq = BioSequence{$(A{4})}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end
end

# conversion between DNA and RNA
for (A1, A2) in [(DNAAlphabet, RNAAlphabet), (RNAAlphabet, DNAAlphabet)], n in (2, 4)
    # NOTE: assumes that binary representation is identical between DNA and RNA
    @eval function Base.convert(::Type{BioSequence{$(A1{n})}},
                                seq::BioSequence{$(A2{n})})
        newseq = BioSequence{$(A1{n})}(seq.data, seq.part, true)
        seq.shared = true
        return newseq
    end
end

# from a vector
function Base.convert{A<:DNAAlphabet}(::Type{BioSequence{A}},
                                      seq::AbstractVector{DNANucleotide})
    return BioSequence{A}(seq, 1, endof(seq))
end
function Base.convert{A<:RNAAlphabet}(::Type{BioSequence{A}},
                                      seq::AbstractVector{RNANucleotide})
    return BioSequence{A}(seq, 1, endof(seq))
end
function Base.convert(::Type{AminoAcidSequence}, seq::AbstractVector{AminoAcid})
    return AminoAcidSequence(seq, 1, endof(seq))
end

# to a vector
Base.convert(::Type{Vector}, seq::BioSequence) = collect(seq)
function Base.convert{A<:DNAAlphabet}(::Type{Vector{DNANucleotide}},
                                      seq::BioSequence{A})
    return collect(seq)
end
function Base.convert{A<:RNAAlphabet}(::Type{Vector{RNANucleotide}},
                                      seq::BioSequence{A})
    return collect(seq)
end
Base.convert(::Type{Vector{AminoAcid}}, seq::AminoAcidSequence) = collect(seq)

# from/to a string
function Base.convert{S<:AbstractString}(::Type{S}, seq::BioSequence)
    return convert(S, [Char(x) for x in seq])
end
Base.convert{S<:AbstractString,A}(::Type{BioSequence{A}}, seq::S) = BioSequence{A}(seq)

# Summaries
# ---------

Base.summary{A<:DNAAlphabet}(seq::BioSequence{A}) = string(length(seq), "nt ", "DNA Sequence")
Base.summary{A<:RNAAlphabet}(seq::BioSequence{A}) = string(length(seq), "nt ", "RNA Sequence")
Base.summary(seq::AminoAcidSequence) = string(length(seq), "aa ", "Amino Acid Sequence")
Base.summary(seq::CharSequence) = string(length(seq), "char ", "Char Sequence")

# Basic Operators
# ---------------

alphabet{A}(::Type{BioSequence{A}}) = alphabet(A)

Base.length(seq::BioSequence) = length(seq.part)
Base.eltype{A}(::Type{BioSequence{A}}) = eltype(A)

function Base.count(f::Function, seq::BioSequence)
    n = 0
    for x in seq
        if f(x)
            n += 1
        end
    end
    return n
end

function Base.map(f::Function, seq::BioSequence)
    return map!(f, copy(seq))
end

function Base.filter(f::Function, seq::BioSequence)
    return filter!(f, copy(seq))
end

function Base.checkbounds(seq::BioSequence, range::UnitRange)
    if 1 ≤ range.start && range.stop ≤ endof(seq)
        return true
    end
    throw(BoundsError(seq, range))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector{Bool})
    if length(seq) == length(locs)
        return true
    end
    throw(BoundsError(seq, locs))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector)
    for i in locs
        checkbounds(seq, i)
    end
    return true
end

function checkdimension(from::Integer, to::Integer)
    if from == to
        return true
    end
    throw(DimensionMismatch(string(
        "attempt to assign ",
        from, " elements to ",
        to,   " elements")))
end

function checkdimension(seq::BioSequence, locs::AbstractVector)
    return checkdimension(length(seq), length(locs))
end

function checkdimension(seq::BioSequence, locs::AbstractVector{Bool})
    return checkdimension(length(seq), sum(locs))
end

# creates a bit mask for given number of bits `n`
mask(n::Integer) = (UInt64(1) << n) - 1
mask{A<:Alphabet}(::Type{A}) = mask(bitsof(A))

# assumes `i` is positive and `bitsof(A)` is a power of 2
@inline function bitindex{A}(seq::BioSequence{A}, i::Integer)
    return BitIndex((i + first(seq.part) - 2) << trailing_zeros(bitsof(A)))
end

@inline function inbounds_getindex{A}(seq::BioSequence{A}, i::Integer)
    j = bitindex(seq, i)
    @inbounds return decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
end

Base.getindex(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)
Base.view(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A},
                                      other::BioSequence{A},
                                      locs::AbstractVector{T})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    for (i, x) in zip(locs, other)
        unsafe_setindex!(seq, x, i)
    end
    return seq
end

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A},
                                      other::BioSequence{A},
                                      locs::UnitRange{T})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    return copy!(seq, locs.start, other, 1)
end

function Base.setindex!{A}(seq::BioSequence{A},
                           other::BioSequence{A},
                           locs::AbstractVector{Bool})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    i = j = 0
    while (i = findnext(locs, i + 1)) > 0
        unsafe_setindex!(seq, other[j+=1], i)
    end
    return seq
end

function Base.setindex!{A}(seq::BioSequence{A}, other::BioSequence{A}, ::Colon)
    return setindex!(seq, other, 1:endof(seq))
end

function Base.setindex!(seq::BioSequence, x, i::Integer)
    checkbounds(seq, i)
    orphan!(seq)
    return unsafe_setindex!(seq, x, i)
end

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A}, x, locs::AbstractVector{T})
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    for i in locs
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function Base.setindex!{A}(seq::BioSequence{A}, x, locs::AbstractVector{Bool})
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    i = j = 0
    while (i = findnext(locs, i + 1)) > 0
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

Base.setindex!{A}(seq::BioSequence{A}, x, ::Colon) = setindex!(seq, x, 1:endof(seq))

# this is "unsafe" because of no bounds check and no orphan! call
@inline function unsafe_setindex!{A}(seq::BioSequence{A}, x, i::Integer)
    bin = enc64(seq, x)
    return encoded_setindex!(seq, bin, i)
end

@inline function encoded_setindex!{A}(seq::BioSequence{A}, bin::UInt64, i::Integer)
    j, r = bitindex(seq, i)
    @inbounds seq.data[j] = (bin << r) | (seq.data[j] & ~(mask(A) << r))
    return seq
end

function Base.resize!{A}(seq::BioSequence{A}, size::Integer)
    if size < 0
        throw(ArgumentError("size must be non-negative"))
    end
    orphan!(seq, size)
    resize!(seq.data, seq_data_len(A, size + seq.part.start - 1))
    seq.part = seq.part.start:seq.part.start+size-1
    return seq
end

Base.empty!(seq::BioSequence) = resize!(seq, 0)

function Base.push!{A}(seq::BioSequence{A}, x)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    encoded_setindex!(seq, bin, endof(seq))
    return seq
end

function Base.unshift!{A}(seq::BioSequence{A}, x)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    copy!(seq, 2, seq, 1, length(seq) - 1)
    encoded_setindex!(seq, bin, 1)
    return seq
end

function Base.pop!(seq::BioSequence)
    if isempty(seq)
        throw(ArgumentError("sequence must be non-empty"))
    end
    x = seq[end]
    deleteat!(seq, endof(seq))
    return x
end

function Base.shift!(seq::BioSequence)
    if isempty(seq)
        throw(ArgumentError("sequence must be non-empty"))
    end
    x = seq[1]
    deleteat!(seq, 1)
    return x
end

function Base.insert!{A}(seq::BioSequence{A}, i::Integer, x)
    checkbounds(seq, i)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    copy!(seq, i + 1, seq, i, endof(seq) - i)
    encoded_setindex!(seq, bin, i)
    return seq
end

function Base.deleteat!{A,T<:Integer}(seq::BioSequence{A}, range::UnitRange{T})
    checkbounds(seq, range)
    copy!(seq, range.start, seq, range.stop + 1, length(seq) - range.stop)
    resize!(seq, length(seq) - length(range))
    return seq
end

function Base.deleteat!(seq::BioSequence, i::Integer)
    checkbounds(seq, i)
    copy!(seq, i, seq, i + 1, length(seq) - i)
    resize!(seq, length(seq) - 1)
    return seq
end

function Base.append!{A}(seq::BioSequence{A}, other::BioSequence{A})
    resize!(seq, length(seq) + length(other))
    copy!(seq, endof(seq) - length(other) + 1, other, 1)
    return seq
end

function Base.copy!{A}(seq::BioSequence{A}, doff::Integer,
                       src::Vector{UInt8},  soff::Integer, len::Integer)
    datalen = seq_data_len(A, len)
    return seq
end

function Base.copy!{A}(dst::BioSequence{A}, src::BioSequence{A})
    return copy!(dst, 1, src, 1)
end

function Base.copy!{A}(dst::BioSequence{A}, doff::Integer,
                       src::BioSequence{A}, soff::Integer)
    return copy!(dst, doff, src, soff, length(src) - soff + 1)
end

function Base.copy!{A}(dst::BioSequence{A}, doff::Integer,
                       src::BioSequence{A}, soff::Integer, len::Integer)
    checkbounds(dst, doff:doff+len-1)
    checkbounds(src, soff:soff+len-1)

    if dst.shared || (dst === src && doff > soff)
        orphan!(dst, length(dst), true)
    end

    id = bitindex(dst, doff)
    is = bitindex(src, soff)
    rest = len * bitsof(A)

    while rest > 0
        # move `k` bits from `src` to `dst`
        x = dst.data[index(id)]
        y = src.data[index(is)]
        if offset(id) < offset(is)
            y >>= offset(is) - offset(id)
            k = min(64 - offset(is), rest)
        else
            y <<= offset(id) - offset(is)
            k = min(64 - offset(id), rest)
        end
        m = mask(k) << offset(id)
        dst.data[index(id)] = y & m | x & ~m

        id += k
        is += k
        rest -= k
    end

    return dst
end

function Base.map!(f::Function, seq::BioSequence)
    orphan!(seq)
    for i in 1:endof(seq)
        unsafe_setindex!(seq, f(inbounds_getindex(seq, i)), i)
    end
    return seq
end

function Base.filter!{A}(f::Function, seq::BioSequence{A})
    orphan!(seq)

    len = 0
    next = bitindex(seq, 1)
    j = index(next)
    datum::UInt64 = 0
    for i in 1:endof(seq)
        x = inbounds_getindex(seq, i)
        if f(x)
            datum |= enc64(seq, x) << offset(next)
            len += 1
            next += bitsof(A)
            if index(next) != j
                seq.data[j] = datum
                datum = 0
                j = index(next)
            end
        end
    end
    if offset(next) > 0
        seq.data[j] = datum
    end
    resize!(seq, len)

    return seq
end

function Base.similar{A}(seq::BioSequence{A}, len::Integer=length(seq))
    return BioSequence{A}(len)
end

# actually, users don't need to create a copy of a sequence.
function Base.copy{A}(seq::BioSequence{A})
    # NOTE: no need to set `seq.shared = true` here
    # since `newseq` will be `orphan!`ed soon.
    newseq = BioSequence{A}(seq.data, 1:endof(seq), true)
    orphan!(newseq, length(seq), true)  # force orphan!
    @assert newseq.data !== seq.data
    return newseq
end

# Replace a BioSequence's data with a copy, copying only what's needed.
# The user should never need to call this, as it has no outward effect on the
# sequence.
function orphan!{A}(seq::BioSequence{A}, size::Integer=length(seq), force::Bool=false)
    if !seq.shared && !force
        return seq
    end

    j, r = bitindex(seq, 1)
    data = Vector{UInt64}(seq_data_len(A, size))

    if !isempty(seq)
        x = seq.data[j] >> r
        l = min(endof(data), index(bitindex(seq, endof(seq))) - j + 1)
        @inbounds @simd for i in 1:l-1
            y = seq.data[j+i]
            data[i] = x | y << (64 - r)
            x = y >> r
        end
        data[l] = x
    end

    seq.data = data
    seq.part = 1:length(seq)
    seq.shared = false
    return seq
end

# Transformations
# ---------------

function Base.reverse!(seq::BioSequence)
    orphan!(seq)
    for i in 1:div(endof(seq), 2)
        x = inbounds_getindex(seq, i)
        unsafe_setindex!(seq, inbounds_getindex(seq, endof(seq) - i + 1), i)
        unsafe_setindex!(seq, x, endof(seq) - i + 1)
    end
    return seq
end

Base.reverse(seq::BioSequence) = reverse!(copy(seq))

@generated function Base.reverse{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    n = bitsof(A)
    if n == 2
        nucrev = :nucrev2
    elseif n == 4
        nucrev = :nucrev4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        data = Vector{UInt64}(seq_data_len(A, length(seq)))
        i = 1
        next = bitindex(seq, endof(seq))
        stop = bitindex(seq, 0)
        r = rem(offset(next) + $n, 64)
        if r == 0
            @inbounds while next - stop > 0
                x = seq.data[index(next)]
                data[i] = $nucrev(x)
                i += 1
                next -= 64
            end
        else
            @inbounds while next - stop > 64
                j = index(next)
                x = (seq.data[j] << (64 - r)) | (seq.data[j-1] >> r)
                data[i] = $nucrev(x)
                i += 1
                next -= 64
            end
            if next - stop > 0
                data[i] = $nucrev(seq.data[index(next)] << (64 - r))
            end
        end
        return BioSequence{A}(data, 1:length(seq), false)
    end
end

@inline function nucrev2(x::UInt64)
     x = (x & 0x3333333333333333) <<  2 | (x & 0xCCCCCCCCCCCCCCCC) >>  2
     x = (x & 0x0F0F0F0F0F0F0F0F) <<  4 | (x & 0xF0F0F0F0F0F0F0F0) >>  4
     x = (x & 0x00FF00FF00FF00FF) <<  8 | (x & 0xFF00FF00FF00FF00) >>  8
     x = (x & 0x0000FFFF0000FFFF) << 16 | (x & 0xFFFF0000FFFF0000) >> 16
     x = (x & 0x00000000FFFFFFFF) << 32 | (x & 0xFFFFFFFF00000000) >> 32
     return x
end

@inline function nucrev4(x::UInt64)
     x = (x & 0x0F0F0F0F0F0F0F0F) <<  4 | (x & 0xF0F0F0F0F0F0F0F0) >>  4
     x = (x & 0x00FF00FF00FF00FF) <<  8 | (x & 0xFF00FF00FF00FF00) >>  8
     x = (x & 0x0000FFFF0000FFFF) << 16 | (x & 0xFFFF0000FFFF0000) >> 16
     x = (x & 0x00000000FFFFFFFF) << 32 | (x & 0xFFFFFFFF00000000) >> 32
     return x
end

"""
    complement!(seq)

Make a complement sequence of `seq` in place.
"""
function complement!{A<:Union{DNAAlphabet{2},RNAAlphabet{2}}}(seq::BioSequence{A})
    orphan!(seq)
    next = bitindex(seq, 1)
    stop = bitindex(seq, endof(seq) + 1)
    @inbounds while next < stop
        seq.data[index(next)] = ~seq.data[index(next)]
        next += 64
    end
    return seq
end

function complement!{A<:Union{DNAAlphabet{4},RNAAlphabet{4}}}(seq::BioSequence{A})
    orphan!(seq)
    next = bitindex(seq, 1)
    stop = bitindex(seq, endof(seq) + 1)
    @inbounds while next < stop
        x = seq.data[index(next)]
        seq.data[index(next)] = (
            ((x & 0x1111111111111111) << 3) | ((x & 0x8888888888888888) >> 3) |
            ((x & 0x2222222222222222) << 1) | ((x & 0x4444444444444444) >> 1))
        next += 64
    end
    return seq
end

"""
    complement(seq)

Make a complement sequence of `seq`.
"""
function complement{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    return complement!(copy(seq))
end

"""
    reverse_complement!(seq)

Make a reversed complement sequence of `seq` in place.

Ambiguous nucleotides are left as-is.
"""
function reverse_complement!{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    return complement!(reverse!(seq))
end

"""
    reverse_complement(seq)

Make a reversed complement sequence of `seq`.

Ambiguous nucleotides are left as-is.
"""
function reverse_complement{A<:Union{DNAAlphabet,RNAAlphabet}}(seq::BioSequence{A})
    return complement!(reverse(seq))
end

# Ambiguous nucleotides iterator
# ------------------------------

immutable AmbiguousNucleotideIterator{A<:Union{DNAAlphabet,RNAAlphabet}}
    seq::BioSequence{A}
end

ambiguous_positions(seq::BioSequence) = AmbiguousNucleotideIterator(seq)

Base.start(it::AmbiguousNucleotideIterator) = find_next_ambiguous(it.seq, 1)
Base.done(it::AmbiguousNucleotideIterator, nextpos) = nextpos == 0
function Base.next(it::AmbiguousNucleotideIterator, nextpos)
    return nextpos, find_next_ambiguous(it.seq, nextpos + 1)
end

Base.iteratorsize(::AmbiguousNucleotideIterator) = Base.SizeUnknown()

function find_next_ambiguous{A<:Union{DNAAlphabet{2},RNAAlphabet{2}}}(
        seq::BioSequence{A}, i::Integer)
    # no ambiguity
    return 0
end

function find_next_ambiguous{A<:Union{DNAAlphabet{4},RNAAlphabet{4}}}(
        seq::BioSequence{A}, from::Integer)
    for i in max(from, 1):endof(seq)
        nt = inbounds_getindex(seq, i)
        if isambiguous(nt)
            return i
        end
    end
    # no ambiguity
    return 0
end

# Mismatch counting
# -----------------

"""
    mismatches(seq1::BioSequence, seq2::BioSequence[, compatible=false])

Return the number of mismatches between `seq1` and `seq2`.

If `seq1` and `seq2` are of differing lengths, only the first `min(length(seq1),
length(seq2))` nucleotides are compared.  When `compatible` is `true`, sequence
symbols are comapred using `iscompatible`; otherwise using `==`.
"""
function mismatches{A<:Alphabet}(
        seq1::BioSequence{A},
        seq2::BioSequence{A},
        compatible::Bool=false)
    if ((bitsof(A) == 2 || bitsof(A) == 4) && !compatible) ||
        A == DNAAlphabet{2} ||
        A == RNAAlphabet{2}
        return bitparallel_mismatches(seq1, seq2)
    end

    mis = 0
    if compatible
        for (x, y) in zip(seq1, seq2)
            if !iscompatible(x, y)
                mis += 1
            end
        end
    else
        for (x, y) in zip(seq1, seq2)
            if x != y
                mis += 1
            end
        end
    end
    return mis
end

@generated function bitparallel_mismatches{A}(a::BioSequence{A}, b::BioSequence{A})
    n = bitsof(A)
    if n == 2
        bitpar_mismatches = :bitpar_mismatches2
    elseif n == 4
        bitpar_mismatches = :bitpar_mismatches4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        if length(a) > length(b)
            return bitparallel_mismatches(b, a)
        end
        @assert length(a) ≤ length(b)

        nexta = bitindex(a, 1)
        nextb = bitindex(b, 1)
        stopa = bitindex(a, endof(a) + 1)
        mismatches = 0

        # align reading position of `a.data` so that `offset(nexta) == 0`
        if nexta < stopa && offset(nexta) != 0
            x = a.data[index(nexta)] >> offset(nexta)
            y = b.data[index(nextb)] >> offset(nextb)
            if offset(nextb) > offset(nexta)
                y |= b.data[index(nextb)+1] << (64 - offset(nextb))
            end
            k = 64 - offset(nexta)
            m = mask(k)
            mismatches += $bitpar_mismatches(x & m, y & m)
            nexta += k
            nextb += k
        end
        @assert offset(nexta) == 0

        if offset(nextb) == 0  # data are aligned with each other
            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                mismatches += $bitpar_mismatches(x, y)
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                m = mask(stopa - nexta)
                mismatches += $bitpar_mismatches(x & m, y & m)
            end
        elseif nexta < stopa
            y = b.data[index(nextb)]
            nextb += 64

            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                z = b.data[index(nextb)]
                y = y >> offset(nextb) | z << (64 - offset(nextb))
                mismatches += $bitpar_mismatches(x, y)
                y = z
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                x = a.data[index(nexta)]
                y = y >> offset(nextb)
                if 64 - offset(nextb) < stopa - nexta
                    y |= a.data[index(nextb)] << (64 - offset(nextb))
                end
                m = mask(stopa - nexta)
                mismatches += $bitpar_mismatches(x & m, y & m)
            end
        end

        return mismatches
    end
end

# bit-parallel mismatch count algorithm for 2 and 4-bit encoding
@inline function bitpar_mismatches2(x::UInt64, y::UInt64)
    xyxor = x $ y
    mismatches = UInt64(0)
    mismatches |=  xyxor & 0x5555555555555555
    mismatches |= (xyxor & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(mismatches)
end

@inline function bitpar_mismatches4(x::UInt64, y::UInt64)
    xyxor = x $ y
    mismatches = UInt64(0)
    mismatches |=  xyxor & 0x1111111111111111
    mismatches |= (xyxor & 0x2222222222222222) >> 1
    mismatches |= (xyxor & 0x4444444444444444) >> 2
    mismatches |= (xyxor & 0x8888888888888888) >> 3
    return count_ones(mismatches)
end

# Shuffle
# -------

function Base.shuffle(seq::BioSequence)
    return shuffle!(copy(seq))
end

function Base.shuffle!(seq::BioSequence)
    orphan!(seq)
    # Fisher-Yates shuffle
    for i in 1:endof(seq)-1
        j = rand(i:endof(seq))
        seq[i], seq[j] = seq[j], seq[i]
    end
    return seq
end

# Encoding
# --------

function encode_copy!{A}(dst::BioSequence{A},
                         src::Union{AbstractVector,AbstractString})
    return encode_copy!(dst, 1, src, 1)
end

function encode_copy!{A}(dst::BioSequence{A},
                         doff::Integer,
                         src::Union{AbstractVector,AbstractString},
                         soff::Integer)
    return encode_copy!(dst, doff, src, soff, length(src) - soff + 1)
end

function encode_copy!{A}(dst::BioSequence{A},
                         doff::Integer,
                         src::Union{AbstractVector,AbstractString},
                         soff::Integer,
                         len::Integer)
    if soff != 1 && !isascii(src)
        throw(ArgumentError("source offset ≠ 1 is not supported for non-ASCII string"))
    end

    checkbounds(dst, doff:doff+len-1)
    if length(src) < soff + len - 1
        throw(ArgumentError("source string does not contain $len elements from $soff"))
    end

    orphan!(dst)
    next = bitindex(dst, doff)
    stop = bitindex(dst, doff + len)
    i = soff
    while next < stop
        x = UInt64(0)
        j = index(next)
        while index(next) == j && next < stop
            char, i = Base.next(src, i)
            x |= enc64(dst, convert(Char, char)) << offset(next)
            next += bitsof(A)
        end
        dst.data[j] = x
    end
    return dst
end

function encode_copy!{A<:Union{DNAAlphabet{4},RNAAlphabet{4}}}(
        dst::BioSequence{A}, doff::Integer,
        src::AbstractVector{UInt8}, soff::Integer, len::Integer)
    checkbounds(dst, doff:doff+len-1)
    if length(src) < soff + len - 1
        throw(ArgumentError("source string does not contain $len elements from $soff"))
    end

    orphan!(dst)
    charmap = A <: DNAAlphabet ? char_to_dna : char_to_rna
    i = soff
    next = bitindex(dst, doff)
    stop = bitindex(dst, doff + len)

    # head
    if offset(next) != 0
        for d in 0:div(64 - offset(next), 4)-1
            dst[doff+d] = charmap[src[i+d]+1]
        end
        i += div(64 - offset(next), 4)
        next += 64 - offset(next)
    end

    # body
    D = 16
    while next < (stop - offset(stop))
        x::UInt64 = 0
        check = 0x00
        @inbounds for d in 0:D-1
            y = reinterpret(UInt8, charmap[src[i+d]+1])
            x |= UInt64(y) << 4d
            check |= y
        end
        if check & 0x80 != 0
            # invalid byte(s) is detected
            for d in 0:D-1
                if !isvalid(charmap[src[i+d]+1])
                    error("cannot encode $(src[i+d])")
                end
            end
        end
        dst.data[index(next)] = x
        i += D
        next += 64
    end

    # tail
    for d in 0:div(stop - next, 4)-1
        dst[doff+i-soff+d] = charmap[src[i+d]+1]
    end

    return dst
end

function enc64{A}(::BioSequence{A}, x)
    return UInt64(encode(A, convert(eltype(A), x)))
end

# Sequences to Matrix
# -------------------

"""
    seqmatrix{A<:Alphabet}(vseq::Vector{BioSequence{A}}, major::Symbol)

Construct a matrix of nucleotides or amino acids from a vector of `BioSequence`s.

If parameter `major` is set to `:site`, the matrix is created such that one
nucleotide from each sequence is placed in each column i.e. the matrix is laid
out in site-major order.
This means that iteration over one position of many sequences is efficient,
as julia arrays are laid out in column major order.

If the parameter `major` is set to `:seq`, the matrix is created such that each
sequence is placed in one column i.e. the matrix is laid out in sequence-major
order.
This means that iteration across each sequence in turn is efficient,
as julia arrays are laid out in column major order.

# Examples
```julia
julia> seqs = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
4-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 3nt DNA Sequence:
AAA
 3nt DNA Sequence:
TTT
 3nt DNA Sequence:
CCC
 3nt DNA Sequence:
GGG

julia> seqmatrix(seqs, :site)
4x3 Array{Bio.Seq.DNANucleotide,2}:
 DNA_A  DNA_A  DNA_A
 DNA_T  DNA_T  DNA_T
 DNA_C  DNA_C  DNA_C
 DNA_G  DNA_G  DNA_G

 julia> seqmatrix(seqs, :seq)
 3x4 Array{Bio.Seq.DNANucleotide,2}:
  DNA_A  DNA_T  DNA_C  DNA_G
  DNA_A  DNA_T  DNA_C  DNA_G
  DNA_A  DNA_T  DNA_C  DNA_G
```
"""
function seqmatrix{A<:Alphabet}(vseq::AbstractVector{BioSequence{A}}, major::Symbol)
    nseqs = length(vseq)
    @assert nseqs > 0 throw(ArgumentError("Vector of BioSequence{$A} is empty."))
    nsites = length(vseq[1])
    @inbounds for i in 2:nseqs
        length(vseq[i]) == nsites || throw(ArgumentError("Sequences in vseq must be of same length"))
    end
    if major == :site
        mat = Matrix{eltype(A)}(nseqs, nsites)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[seq, site] = vseq[seq][site]
        end
        return mat
    elseif major == :seq
        mat = Matrix{eltype(A)}(nsites, nseqs)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[site, seq] = vseq[seq][site]
        end
        return mat
    else
        throw(ArgumentError("major must be :site or :seq"))
    end
end

"""
    seqmatrix{A<:Alphabet,T}(::Type{T}, vseq::Vector{BioSequence{A}}, major::Symbol)

Construct a matrix of `T` from a vector of `BioSequence`s.

If parameter `major` is set to `:site`, the matrix is created such that one
nucleotide from each sequence is placed in each column i.e. the matrix is laid
out in site-major order.
This means that iteration over one position of many sequences is efficient,
as julia arrays are laid out in column major order.

If the parameter `major` is set to `:seq`, the matrix is created such that each
sequence is placed in one column i.e. the matrix is laid out in sequence-major
order.
This means that iteration across each sequence in turn is efficient,
as julia arrays are laid out in column major order.

# Examples
```julia
julia> seqs = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
4-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 3nt DNA Sequence:
AAA
 3nt DNA Sequence:
TTT
 3nt DNA Sequence:
CCC
 3nt DNA Sequence:
GGG

julia> seqmatrix(seqs, :site, UInt8)
4×3 Array{UInt8,2}:
 0x01  0x01  0x01
 0x08  0x08  0x08
 0x02  0x02  0x02
 0x04  0x04  0x04

julia> seqmatrix(seqs, :seq, UInt8)
3×4 Array{UInt8,2}:
 0x01  0x08  0x02  0x04
 0x01  0x08  0x02  0x04
 0x01  0x08  0x02  0x04
```
"""
function seqmatrix{T,A<:Alphabet}(::Type{T}, vseq::AbstractVector{BioSequence{A}}, major::Symbol)
    nseqs = length(vseq)
    @assert nseqs > 0 throw(ArgumentError("Vector of BioSequence{$A} is empty."))
    nsites = length(vseq[1])
    @inbounds for i in 2:nseqs
        length(vseq[i]) == nsites || throw(ArgumentError("Sequences in vseq must be of same length."))
    end
    if major == :site
        mat = Matrix{T}(nseqs, nsites)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[seq, site] = convert(T, vseq[seq][site])
        end
        return mat
    elseif major == :seq
        mat = Matrix{T}(nsites, nseqs)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[site, seq] = convert(T, vseq[seq][site])
        end
        return mat
    else
        throw(ArgumentError("major must be :site or :seq"))
    end
end

# Consensus
# ---------

"""
    majorityvote{A<:NucleotideAlphabet}(seqs::AbstractVector{BioSequence{A}})

Construct a sequence that is a consensus of a vector of sequences.

The consensus is established by a simple majority vote rule, where amiguous
nucleotides cast an equal vote for each of their possible states.
For each site a winner(s) out of A, T(U), C, or G is determined, in the cases
of ties the ambiguity symbol that unifies all the winners is returned.
E.g if A and T tie, then W is inserted in the consensus. If all A, T, C, and G
tie at a site, then N is inserted in the consensus. Note this means that if a
nucletide e.g. 'C' and a gap '-' draw, the nucleotide will always win over the
gap, even though they tied.

# Examples
```julia
julia> seqs = [dna"CTCGATCGATCC", dna"CTCGAAAAATCA", dna"ATCGAAAAATCG", dna"ATCGGGGGATCG"]

4-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 CTCGATCGATCC
 CTCGAAAAATCA
 ATCGAAAAATCG
 ATCGGGGGATCG

julia> majorityvote(seqs)
12nt DNA Sequence:
MTCGAAARATCG
```
"""
function majorityvote{A<:NucleotideAlphabet}(seqs::AbstractVector{BioSequence{A}})
    mat = seqmatrix(UInt8, seqs, :site)
    nsites = size(mat, 2)
    nseqs = size(mat, 1)
    result = BioSequence{A}(nsites)
    votes = Array{Int}(16)
    @inbounds for site in 1:nsites
        fill!(votes, 0)
        for seq in 1:nseqs
            nuc = mat[seq, site]
            votes[1] += nuc == 0x00
            votes[2] += (nuc & 0x01) != 0x00
            votes[3] += (nuc & 0x02) != 0x00
            votes[5] += (nuc & 0x04) != 0x00
            votes[9] += (nuc & 0x08) != 0x00
        end
        m = maximum(votes)
        merged = 0x00
        for i in 0x01:0x10
            merged |= ifelse(votes[i] == m, i - 0x01, 0x00)
        end
        result[site] = reinterpret(eltype(A), merged)
    end
    return result
end


text:
#= ForwardDiff is removed as a dependency from this version
import MechanicalSketch
import MechanicalSketch: empty_figure, PALETTE, O, HE, WI, EM, finish, ∙, settext
import MechanicalSketch: @import_expand, set_scale_sketch, color_with_lumin
import MechanicalSketch: generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: clamped_velocity_matrix, matrix_to_function
import MechanicalSketch: lenient_max
import MechanicalSketch: place_image, @layer
import MechanicalSketch: lic_matrix_current
import MechanicalSketch: ColorSchemes, color_with_lumin, convolution_matrix, clamped_velocity_matrix
import MechanicalSketch: BinLegend, BinLegendVector, streamlines_matrix, draw_legend
import MechanicalSketch: leonardo, fontsize
import ColorSchemes:     Paired_6, Greys_9, isoluminant_cgo_80_c38_n256
import Base.show

let
if !@isdefined m²
    @import_expand(~m, s, °)
end

include("test_functions_33.jl")

empty_figure(filename = joinpath(@__DIR__, "test_33.png");
    backgroundcolor = color_with_lumin(PALETTE[8], 10))

# Scaling and placement
physwidth = 10.0m
physheight = 4.0m
framedheight = physheight * 1.2
totheight = framedheight * 3
set_scale_sketch(totheight, HE)
totwidth = totheight * WI / HE
Δx = totwidth / 4
Δy = framedheight
Oadj = O + (0, EM /3)

# Reused velocity field from earlier tests
max_velocity = 0.5m/s
velocity_matrix = clamped_velocity_matrix(ϕ_33; physwidth, physheight, cutoff = max_velocity);

v_xy = matrix_to_function(velocity_matrix)

streamlinepixels = streamlines_matrix(v_xy; physwidth, physheight)

cent1 = Oadj + (-Δx, Δy)
legend1 = BinLegend(;maxlegend = 1.0, minlegend = 0.0, noofbins = 2,
                colorscheme = leonardo,
                nan_color = color_with_lumin(PALETTE[1], 80), name = Symbol("Value{Float64}"))
bw(x::Bool) = legend1(Float64(x))
ulp, _ = place_image(cent1, bw.(streamlinepixels))
@layer begin
    fontsize(25)
    draw_legend(ulp + (EM, EM), legend1)
end
str = "    1: Streamline pixels"
settext(str, ulp, markup = true)

complex_convolution_matrix = convolution_matrix(velocity_matrix, streamlinepixels)
legend2 = BinLegendVector(;operand_example = first(complex_convolution_matrix),
        max_magn_legend = lenient_max(complex_convolution_matrix),
        noof_ang_bins = 40, noof_magn_bins = 5,
        name = Symbol("Convolution matrix"))
cent2 = Oadj + (Δx, Δy)
ulp, _ = place_image(cent2, legend2.(complex_convolution_matrix))
@layer begin
    fontsize(25)
    draw_legend(ulp + (EM, EM), legend2)
end
str = "    2: Phase and amplitude for the animation"
settext(str, ulp, markup = true)

cent3 = Oadj + (-Δx, 0.0m)
curmat = lic_matrix_current(complex_convolution_matrix, 0, zeroval = -0.0)
legend3 = BinLegend(;maxlegend = 1.0, minlegend = -1.0, noofbins = 256,
                       colorscheme = reverse(Greys_9),
                       name = Symbol(" "))

ulp, _ = place_image(cent3, legend3.(curmat))
@layer begin
    fontsize(25)
    draw_legend(ulp + (EM, EM), legend3; max_vert_height = 3.0m, background_opacity = 0.4)
end
str = "    3: A single frame taken from 2"
settext(str, ulp, markup = true)

cent4 = Oadj + (Δx, 0.0m)
speedmatrix = hypot.(velocity_matrix)
binwidth = 0.1m/s
lum = 50
cols = isoluminant_cgo_80_c38_n256.colors .|> co -> color_with_lumin(co, lum)
legend4 = BinLegend(maxlegend = max_velocity, binwidth = binwidth,
          colorscheme = cols,
          name = :Speed)
ulp, _ = place_image(cent4, legend4.(speedmatrix));
@layer begin
    fontsize(25)
    draw_legend(ulp + (EM, EM), legend4)
end
str = "    4: isoluminant_cgo_80_c38_n256 at lumin = $lum"
settext(str, ulp, markup = true)

cent5 = Oadj + (-Δx, -Δy)
mixmat = map(legend4.(speedmatrix), curmat) do col, lu
    color_with_lumin(col, 50 + 50 * lu)
end;
ulp, _ = place_image(cent5, mixmat);
str = "    5: Chroma and hue from 4, lumin as in 3."
settext(str, ulp, markup = true)

cent6 = Oadj + (Δx, -Δy)
legend6 = BinLegend(;maxlegend = max_velocity, binwidth,
        colorscheme = Paired_6, name = :Speed)
mixmat = map(legend6.(speedmatrix), curmat) do col, lu
    color_with_lumin(col, 50 + lu * 50)
end;
ulp, _ = place_image(cent6, mixmat)
@layer begin
    fontsize(25)
    draw_legend(ulp + (EM, EM), legend6)
end
str = "    6: Colorscheme: Paired_6"
settext(str, ulp, markup = true)

set_scale_sketch(m)
finish()

end # let
=#


text:
Base.start(ll::LinkedList) = ll.head
Base.done(ll::LinkedList{T}, st::AbstractNode{T}) where T = st isa EmptyNode
Base.next(ll::LinkedList{T}, st::AbstractNode{T}) where T = st.data, st.next

lst = LinkedList{Int}()
push!(lst, 1)
push!(lst, 2)
push!(lst, 3)

for n in lst
    print(n, " ")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using TestDB
using Test

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

max_stars_repo_name: greenspaceexplorer/TestDB

max_issues_repo_name: greenspaceexplorer/TestDB

max_forks_repo_name: greenspaceexplorer/TestDB


text:
#------------------------------------------------------------------------------
"""
    printTeXTable(fh::IO,x,colNames=[],rowNames=[];width=12,prec=2)

Create LaTeX table from numeric matrix x
"""
function printTeXTable(fh::IO,x,colNames=[],rowNames=[];width=12,prec=2)

  (m,n) = (size(x,1),size(x,2))

  isempty(rowNames) && (rowNames = [string("r",i) for i = 1:m])   #create row names "r1"
  isempty(colNames) && (colNames = [string("c",i) for i = 1:n])   #create column names "c1"

  str = """
  \\begin{table}
    \\begin{tabular}{l"""
  str = string(str,"r"^n,"}\n","    ")

  str = string(str," & ")                              #empty cell above rowNames[1]
  for i = 1:n-1                                        #column names
    str = string(str,colNames[i]," & ")
  end
  str = string(str,colNames[n]," \\\\ \\hline \n")

  for i = 1:m                                                          #loop over rows
    str = string(str,"   ",rowNames[i]," & ")                          #row name
    for j = 1:n-1
      str = string(str,fmtNumPs(x[i,j],width,prec,"right"), " & ")     #x[i,1:end-1]
    end
    str = string(str,fmtNumPs(x[i,n],width,prec,"right")," \\\\ \n")   #x[i,end]
  end

  endstr = """
    \\hline
    \\end{tabular}
  \\end{table}"""
  str = string(str,endstr)

  print(fh,str,"\n")                                                    #print
  return nothing

end
                                      #when fh is not supplied: printing to screen
printTeXTable(x,colNames=[],rowNames=[];width=12,prec=2) =
printTeXTable(stdout::IO,x,colNames,rowNames,width=width,prec=prec)
#------------------------------------------------------------------------------

max_stars_repo_name: snowdj/JuliaTutorial

max_issues_repo_name: snowdj/JuliaTutorial

max_forks_repo_name: snowdj/JuliaTutorial


text:
"""
	mapproject(cmd0::String="", arg1=nothing, kwargs...)

Forward and inverse map transformations, datum conversions and geodesy.

Full option list at [`mapproject`]($(GMTdoc)mapproject.html)

Parameters
----------

- $(GMT.opt_R)
- $(GMT.opt_J)

- **A** | **origin** :: [Type => Str]    ``Arg = b|B|f|F|o|O[lon0/lat0][+v]``

    Calculate azimuth along track or to the optional fixed point set with lon0/lat0.
    ($(GMTdoc)mapproject.html#a)
- **C** | **center** :: [Type => Str | List | []]     ``Arg = [dx/dy]``

    Set center of projected coordinates to be at map projection center [Default is lower left corner].
    ($(GMTdoc)mapproject.html#c)
- **D** | **override_units** :: [Type => Str]    ``Arg = c|i|p``

    Temporarily override PROJ_LENGTH_UNIT and use c (cm), i (inch), or p (points) instead.
    ($(GMTdoc)mapproject.html#d)
- **E** | **geod2ecef** | **ecef** :: [Type => Str | []]    ``Arg = [datum]``

    Convert from geodetic (lon, lat, height) to Earth Centered Earth Fixed (ECEF) (x,y,z) coordinates.
    ($(GMTdoc)mapproject.html#e)
- **F** | **one2one** :: [Type => Str | []]    ``Arg = [unit]``

    Force 1:1 scaling, i.e., output (or input, see I) data are in actual projected meters.
    ($(GMTdoc)mapproject.html#f)
- **G** | **track_distances** :: [Type => Str | List]    ``Arg = [lon0/lat0][+a][+i][+u[+|-]unit][+v]``

    Calculate distances along track or to the optional fixed point set with G="lon0/lat0".
    ($(GMTdoc)mapproject.html#g)
- **I** | **inverse** :: [Type => Bool]

    Do the Inverse transformation, i.e., get (longitude,latitude) from (x,y) data.
    ($(GMTdoc)mapproject.html#i)
- **L** | **dist2line** :: [Type => Str | NamedTuple]   ``Arg = line.xy[+u[+|-]unit][+p] | (line=Matrix, unit=x, fractional_pt=_,cartesian=true, projected=true)``

    Determine the shortest distance from the input data points to the line(s) given in the
    ASCII multisegment file line.xy.
    ($(GMTdoc)mapproject.html#l)
- **N** | **geod2aux** :: [Type => Str | []]       ``Arg = [a|c|g|m]``

    Convert from geodetic latitudes to one of four different auxiliary latitudes (longitudes are unaffected).
    ($(GMTdoc)mapproject.html#n)
- **Q** | **list** :: [Type => Str | []]           ``Arg = [d|e]``

    List all projection parameters. To only list datums, use Q=:d, to only list ellipsoids, use Q=:e.
    ($(GMTdoc)mapproject.html#q)
- **S** | **supress** :: [Type => Bool]

    Suppress points that fall outside the region.
    ($(GMTdoc)mapproject.html#s)
- **T** | **change_datum** :: [Type => Str]    ``Arg = [h]from[/to]``

    Coordinate conversions between datums from and to using the standard Molodensky transformation.
    ($(GMTdoc)mapproject.html#t)
- **W** | **map_size** :: [Type => Str | []]    ``Arg = [w|h]``

    Prints map width and height on standard output. No input files are read.
    ($(GMTdoc)mapproject.html#w)
- **Z** | **travel_times** :: [Type => Str | Number]    ``Arg = [speed][+a][+i][+f][+tepoch]``

    Calculate travel times along track as specified with -G.
    ($(GMTdoc)mapproject.html#z)

- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_p)
- $(GMT.opt_s)
- $(GMT.opt_swap_xy)
"""
function mapproject(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("mapproject", cmd0, arg1, arg2)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :V_params :b :d :e :f :g :h :i :j :o :p :s :yx])
	cmd  = parse_these_opts(cmd, d, [[:C :center], [:E :geod2ecef :ecef], [:I :inverse], [:S :supress], #[:L :dist2line],
	                                 [:T :change_datum], [:W :map_size], [:Z :travel_times]])
	cmd  = add_opt_1char(cmd, d, [[:D :override_units], [:F :one2one], [:Q :list], [:N :geod2aux]])

	cmd = add_opt(d, cmd, 'A', [:A :azim],
	              (fixed_pt=("", arg2str), back=("b", nothing, 1), back_geocentric=("B", nothing, 1), forward=("f", nothing, 1), forward_geocentric=("F", nothing, 1), orientation=("o", nothing, 1), orientation_geocentric=("O", nothing, 1), unit="+u1", var_pt="_+v"))
	cmd = add_opt(d, cmd, 'G', [:G :track_distances],
	              (fixed_pt=("", arg2str, 1), accumulated="_+a", incremental="_+i", unit="+u1", var_pt="_+v"))

    cmd, args, n, = add_opt(d, cmd, 'L', [:L :dist2line], :line, Array{Any,1}([arg1, arg2]),
                            (unit="+u1", cartesian="_+uc", project="_+uC", fractional_pt="_+p"))
	if (n > 0)
		arg1, arg2 = args[:]
	end

	if (!occursin("-G", cmd))
		map = occursin(" -W", cmd) ? true : false
		cmd, = parse_J(d, cmd, " ", map)		# Do not append a default fig size
	end

	if (occursin(" -W", cmd))				# No input data in this case
		(dbg_print_cmd(d, cmd) !== nothing) && return cmd
		gmt("mapproject " * cmd)
	else
	    #cmd, got_fname, arg1, arg2 = find_data(d, cmd0, cmd, arg1, arg2)
		R = common_grd(d, cmd0, cmd, "mapproject ", arg1, arg2)
        (isa(R, Vector{GMTdataset}) && contains(cmd, " -I")) && (R[1].proj4 = prj4WGS84)
        R
	end
end

# ---------------------------------------------------------------------------------------------------
mapproject(arg1, arg2=nothing, cmd0::String=""; kw...) = mapproject(cmd0, arg1, arg2; kw...)

#mapproject(, G=(fixed_pt=[1 2], unit=:n,accumulated=1,incremental=1), azim=(fixed_pt=(3,4),forward=1), Vd=2) = "mapproject  -Af3/4 -G1/2+un+a+i"


text:
"""
Functions that return all objects required for a particular kind of simulation. E.g.,

    ```julia
    function rotor_sweep_template() <: Template
        ...
        return system, parameters, actions, freestream_function, environment_function, objective_function, step_range
    end
    ```
"""
abstract type Template end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
module Ledger

# package code goes here
include("engine.jl"); export LedgerEngine, run;

end # module


text:
immutable TinySparseMat{T} <: AbstractLinearMap{T}
    f::Function
    fc::Function
    M::Int
    N::Int
    _nnz::Int
    dims::Tuple{Int, Int}
    rowind::TinyIntVec
    colind::TinyIntVec
    nzind::TinyIntVec
    nzval::Vector{T}
    rowind_wblock::Vector{UInt32}
    colind_wblock::Vector{UInt32}
    nzind_wblock::Vector{UInt32}
    _issymmetric::Bool
    _ishermitian::Bool
    _isposdef::Bool
end

# Properties
Base.size(A::TinySparseMat) = A.dims
Base.nnz(A::TinySparseMat) = A._nnz
Base.issymmetric(A::TinySparseMat) = A._issymmetric
Base.ishermitian(A::TinySparseMat) = A._ishermitian
Base.isposdef(A::TinySparseMat) = A._isposdef

# Multiplication with vector
function Base.A_mul_B!(y::AbstractVector, A::TinySparseMat, x::AbstractVector)
    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch())
    A.f(y,A,x)
    return y
end

function *(A::TinySparseMat, x::AbstractVector)
    length(x) == A.N || throw(DimensionMismatch())
    y = similar(x, promote_type(eltype(A), eltype(x)), A.M)
    A.f(y,A,x)
    return y
end

function Base.At_mul_B!(y::AbstractVector, A::TinySparseMat, x::AbstractVector)
    issymmetric(A) && return Base.A_mul_B!(y, A, x)
    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())
    if !isreal(A)
        x = conj(x)
    end
    A.fc(y,A,x)
    if !isreal(A)
        conj!(y)
    end
    return y
end
function Base.At_mul_B(A::TinySparseMat, x::AbstractVector)
    issymmetric(A) && return A*x
    length(x) == A.M || throw(DimensionMismatch())
    if !isreal(A)
        x = conj(x)
    end
    y = similar(x, promote_type(eltype(A), eltype(x)), A.N)
    A.fc(y,A,x)
    if !isreal(A)
        conj!(y)
    end
    return y
end

function Base.Ac_mul_B!(y::AbstractVector, A::TinySparseMat, x::AbstractVector)
    ishermitian(A) && return Base.A_mul_B!(y,A,x)
    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())
    return A.fc(y,A,x)
end

function Base.Ac_mul_B(A::TinySparseMat, x::AbstractVector)
    ishermitian(A) && return A*x
    length(x) == A.M || throw(DimensionMismatch())
    y = similar(x, promote_type(eltype(A), eltype(x)), A.N)
    A.fc(y,A,x)
    return y
end


text:
"""
DistributedFactorGraphs.jl provides a flexible factor graph API for use in the Caesar.jl ecosystem.

The package supplies:

- A standardized API for interacting with factor graphs
- Implementations of the API for in-memory and database-driven operation
- Visualization extensions to validate the underlying graph

"""
module DistributedFactorGraphs

##==============================================================================
## imports
##==============================================================================

using Base
using DocStringExtensions
using Requires
using Dates
using Distributions
using Reexport
using JSON
using Unmarshal
using JSON2 # JSON2 requires all properties to be in correct sequence, can't guarantee that from DB.
using LinearAlgebra
using SparseArrays
using UUIDs

##==============================================================================
# Exports
##==============================================================================
const DFG = DistributedFactorGraphs
export DFG
##------------------------------------------------------------------------------
## DFG
##------------------------------------------------------------------------------
export AbstractDFG
export AbstractParams, NoSolverParams

# accessors & crud
export getUserId, getRobotId, getSessionId
export getDFGInfo
export getDescription, setDescription!,
       getSolverParams, setSolverParams!,
       getUserData, setUserData!,
       getRobotData, setRobotData!,
       getSessionData, setSessionData!,
       getAddHistory

# TODO Not sure these are needed or should work everywhere, implement in cloud?
export updateUserData!, updateRobotData!, updateSessionData!, deleteUserData!, deleteRobotData!, deleteSessionData!
export emptyUserData!, emptyRobotData!, emptySessionData!

# Graph Types: exported from modules or @reexport
export InMemoryDFGTypes, DefaultDFG
# LightDFG
# GraphsDFG

# CloudGraphsDFG
# FileDFG

# AbstractDFG Interface
export exists,
       addVariable!, addFactor!,
       getVariable, getFactor,
       updateVariable!, updateFactor!,
       deleteVariable!, deleteFactor!,
       listVariables, listFactors,
       getVariables, getFactors,
       isVariable, isFactor

export isConnected

export getBiadjacencyMatrix

#summary structure
export DFGSummary

export getSummary, getSummaryGraph

# Abstract Nodes
export DFGNode, AbstractDFGVariable, AbstractDFGFactor

# Variables
export DFGVariable, DFGVariableSummary, SkeletonDFGVariable

# Factors
export DFGFactor, DFGFactorSummary, SkeletonDFGFactor

# Common
export getSolvable, setSolvable!, isSolvable
export getInternalId
export getVariableLabelNumber

# accessors
export getLabel, getTimestamp, setTimestamp, setTimestamp!, getTags, setTags!

# Node Data
export isSolveInProgress, getSolveInProgress

# CRUD & SET
export listTags, mergeTags!, removeTags!, emptyTags!

#this isn't acttually implemented. TODO remove or implement
export addTags!

##------------------------------------------------------------------------------
# Variable
##------------------------------------------------------------------------------
# Abstract Variable Data
export InferenceVariable

# accessors
export getSolverDataDict, setSolverData!
export getSofttype, getSofttypename

export getSolverData

export getVariableType

export getSmallData, setSmallData!

# CRUD & SET
export getVariableSolverData,
       addVariableSolverData!,
       updateVariableSolverData!,
       deleteVariableSolverData!,
       listVariableSolverData,
       mergeVariableSolverData!

# PPE
##------------------------------------------------------------------------------
# PPE TYPES
export AbstractPointParametricEst
export MeanMaxPPE
export getMaxPPE, getMeanPPE, getSuggestedPPE, getLastUpdatedTimestamp

# accessors

export getPPEDict, getVariablePPEDict, getVariablePPE

# CRUD & SET
export getPPE,
       getVariablePPE,
       addPPE!,
       updatePPE!,
       deletePPE!,
       listPPEs,
       mergePPEs! #TODO look at rename to just mergePPE to match with cloud?

# Variable Node Data
##------------------------------------------------------------------------------
export VariableNodeData, PackedVariableNodeData

export packVariableNodeData, unpackVariableNodeData

export getSolvedCount, isSolved, setSolvedCount!, isInitialized

export getNeighborhood, getNeighbors, _getDuplicatedEmptyDFG
export findFactorsBetweenNaive
export copyGraph!, deepcopyGraph, deepcopyGraph!, buildSubgraph, mergeGraph!
# Big Data
##------------------------------------------------------------------------------
export addBigDataEntry!, getBigDataEntry, updateBigDataEntry!, deleteBigDataEntry!, getBigDataEntries, getBigDataKeys, hasDataEntry, hasBigDataEntry
# convenience wrappers
export addDataEntry!, getDataEntryElement
# aliases
export addData!, fetchData, fetchDataEntryElement

##------------------------------------------------------------------------------
# Factors
##------------------------------------------------------------------------------
# Factor Data
export GenericFunctionNodeData, PackedFunctionNodeData, FunctionNodeData
export InferenceType, PackedInferenceType, FunctorInferenceType, FactorOperationalMemory
export FunctorSingleton, FunctorPairwise, FunctorPairwiseMinimize

# accessors
export getVariableOrder
export getFactorType, getFactorFunction

# Node Data
export mergeVariableData!, mergeGraphVariableData!

##------------------------------------------------------------------------------
## Other utility functions
##------------------------------------------------------------------------------

# Sort
export natural_lt, sortDFG

# Validation
export isValidLabel

## List
export ls, lsf, ls2
export lsTypes, lsfTypes, lsTypesDict, lsfTypesDict
export lsWho, lsfWho
export isPrior, lsfPriors
export hasTags, hasTagsNeighbors

## Finding
export findClosestTimestamp, findVariableNearTimestamp

# Serialization
export packVariable, unpackVariable, packFactor, unpackFactor
export rebuildFactorMetadata!

# File import and export
export saveDFG, loadDFG!
export toDot, toDotFile

# Comparisons
export
    compare,
    compareField,
    compareFields,
    compareAll,
    compareAllSpecial,
    compareVariable,
    compareFactor,
    compareAllVariables,
    compareSimilarVariables,
    compareSubsetFactorGraph,
    compareSimilarFactors,
    compareFactorGraphs

## Deprecated exports should be listed in Deprecated.jl if possible, otherwise here
#TODO remove export in DFG v0.8.0
# export ConvolutionObject

## CustomPrinting.jl
export printFactor, printVariable, printNode

##==============================================================================
## Files Includes
##==============================================================================

# Entities

include("entities/AbstractDFG.jl")

include("entities/DFGFactor.jl")

include("entities/DFGVariable.jl")

include("entities/AbstractDFGSummary.jl")

include("services/AbstractDFG.jl")

# In Memory Types
# include("GraphsDFG/GraphsDFG.jl")
# @reexport using .GraphsDFGs
include("LightDFG/LightDFG.jl")
@reexport using .LightDFGs
#supported in Memory fg types
# const InMemoryDFGTypes = Union{GraphsDFG, LightDFG}
const InMemoryDFGTypes = LightDFG
const DefaultDFG = LightDFG

# Common includes
include("services/CommonAccessors.jl")
include("services/Serialization.jl")
include("services/DFGVariable.jl")
include("services/DFGFactor.jl")
include("Deprecated.jl")
include("services/CompareUtils.jl")

# Include the FilesDFG API.
include("FileDFG/FileDFG.jl")

# Custom show and printing for variable factor etc.
include("services/CustomPrinting.jl")

# To be moved as necessary.
include("Common.jl")

# Big data extensions
include("BigData/BigData.jl")

function __init__()
    @require GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" begin
        @info "Including Plots"
        include("DFGPlots/DFGPlots.jl")
        @reexport using .DFGPlots
    end

    if get(ENV, "DFG_USE_CGDFG", "") == "true"
        @info "Including CloudGraphsDFG"
        Base.include(DistributedFactorGraphs, joinpath(@__DIR__, "CloudGraphsDFG/CloudGraphsDFG.jl"))
        @eval DistributedFactorGraphs @reexport using .CloudGraphsDFGs
    end
end

end