text:
nrows(ifd::IFD) = Int(ifd[IMAGELENGTH].data)::Int
ncols(ifd::IFD) = Int(ifd[IMAGEWIDTH].data)::Int
nsamples(ifd::IFD) = Int(ifd[SAMPLESPERPIXEL].data)::Int

"""
    interpretation(ifd)

For a given IFD, determine the proper colorimetric interpretation of the data.
It returns subtypes of `Colorant` depending on the values of the tiff tags and
whether there are extrasamples it doesn't know how to deal with.
"""
function interpretation(ifd::IFD)
    interp = PhotometricInterpretations(ifd[PHOTOMETRIC].data)
    extras = EXTRASAMPLE_UNSPECIFIED
    if EXTRASAMPLES in ifd
        try
            extras = ExtraSamples(ifd[EXTRASAMPLES].data)
        catch
            extras = EXTRASAMPLE_ASSOCALPHA
        end
    end
    interpretation(interp, extras, nsamples(ifd))
end

# dummy color type for palette colored images to dispatch on
struct Palette{T} <: Colorant{T, 1}
    i::T
end
Base.reinterpret(::Type{Palette{T}}, arr::A) where {T, N, S, A <: AbstractArray{S, N}} = arr

interpretation(p::PhotometricInterpretations) = interpretation(Val(p))
interpretation(::Val{PHOTOMETRIC_RGB}) = RGB
interpretation(::Val{PHOTOMETRIC_MINISBLACK}) = Gray
interpretation(::Val{PHOTOMETRIC_PALETTE}) = Palette
interpretation(::Val{PHOTOMETRIC_YCBCR}) = YCbCr
interpretation(::Val{PHOTOMETRIC_CIELAB}) = Lab

function interpretation(p::PhotometricInterpretations, extrasamples::ExtraSamples, samplesperpixel::Int)
    interp = interpretation(p)::Type{<:Colorant}
    len = length(interp)::Int
    if len + 1 == samplesperpixel
        return interpretation(p, extrasamples, Val(samplesperpixel))
    elseif len == samplesperpixel
        return interp, false
    elseif len < samplesperpixel
        return interp, true
    else
        error("TIFF file says it contains $interp values, but only has $samplesperpixel samples per pixel instead of the minimum required $len")
    end
end
_pad(::Type{RGB}) = RGBX
_pad(::Type{T}) where {T} = T

interpretation(p::PhotometricInterpretations, extrasamples::ExtraSamples, nsamples::Val) = interpretation(p, Val(extrasamples), nsamples)
interpretation(p::PhotometricInterpretations, ::Val{EXTRASAMPLE_UNSPECIFIED}, @nospecialize(::Val)) = interpretation(p), true
interpretation(p::PhotometricInterpretations, ::Val{EXTRASAMPLE_UNSPECIFIED}, ::Val{4}) = _pad(interpretation(p)), false
interpretation(p::PhotometricInterpretations, ::Val{EXTRASAMPLE_ASSOCALPHA}, @nospecialize(::Val)) = coloralpha(interpretation(p)), false
interpretation(p::PhotometricInterpretations, ::Val{EXTRASAMPLE_UNASSALPHA}, nsamples::Val) = interpretation(p, Val(EXTRASAMPLE_ASSOCALPHA), nsamples)

_mappedtype(::Type{T}) where {T} = T
_mappedtype(::Type{T}) where {T <: Unsigned} = Normed{T, sizeof(T) * 8}
_mappedtype(::Type{T}) where {T <: Signed}   = Fixed{T, sizeof(T) * 8 - 1}

function rawtype(ifd::IFD)
    samplesperpixel = nsamples(ifd)
    bitsperpixel = ifd[BITSPERSAMPLE].data
    sampleformats = fill(UInt16(0x01), samplesperpixel)
    if SAMPLEFORMAT in ifd
        sampleformats = ifd[SAMPLEFORMAT].data
    end
    rawtype_mapping[SampleFormats(first(sampleformats)), first(bitsperpixel)]
end

"""
    $(SIGNATURES)

Allocate a cache for this IFD with correct type and size.
"""
function getcache(ifd::IFD)
    T = rawtype(ifd)
    if T === Bool
        return BitArray(undef, ncols(ifd), nrows(ifd))
    end
    colortype, extras = interpretation(ifd)
    return Array{colortype{_mappedtype(T)}}(undef, ncols(ifd), nrows(ifd))
end


text:
using SparseArrays
include("./src/preprocessing.jl")
r, c = SparseReadData("../auto.mtx")

sr = sparse(r, c, ones(length(r),1), max(c...),max(c...))
num_tr = sum((sr*sr).*sr)

ext: jl


text:
function test_edges_to_syn_dicts()
	tbl = Array{Any,2}(1,3);
	tbl[1,:] = [1, [1,2], [1,1,1]];
	pre, post = edges_to_syn_dicts(tbl)
	@test length(pre) == 1
	@test length(post) == 1
	@test pre[1] == 1
	@test post[1] == 2
end

function test_map_synapses()
	tbl_A = Array{Any,2}(2,3);
	tbl_B = Array{Any,2}(2,3);
	tbl_A[1,:] = [1, [1,2], [1,1,1]]
	tbl_A[2,:] = [5, [2,3], [1,2,1]]
	tbl_B[1,:] = [10, [20,21], [1,1,1]]
	tbl_B[2,:] = [1, [1,2], [1,1,5]]
	A_to_B = map_synapses(tbl_A, tbl_B)
	@test size(A_to_B) == (3,2)
	@test A_to_B[1,1] == 1
	@test A_to_B[1,2] == 10
	@test A_to_B[2,1] == 5
	@test A_to_B[2,2] == 0
	@test A_to_B[3,1] == 0
	@test A_to_B[3,2] == 1
end

function test_get_indexed_seg_IDs()
	tbl = Array{Any,2}(2,3);
	tbl[1,:] = [1, [1,2], [1,1,1]]
	tbl[2,:] = [5, [2,5], [1,2,1]]
	seg_to_index, index_to_seg = get_indexed_seg_IDs(tbl)
	@test length(seg_to_index) == 3
	@test seg_to_index[1] == 1
	@test seg_to_index[2] == 2
	@test seg_to_index[5] == 3
	@test !haskey(seg_to_index, 3)

	@test index_to_seg[1] == 1
	@test index_to_seg[2] == 2
	@test index_to_seg[3] == 5
	seg_to_index, index_to_seg = get_indexed_seg_IDs(tbl, 1)
	@test index_to_seg[2] == 1
	@test index_to_seg[3] == 2
	@test index_to_seg[4] == 5
end

function test_build_count_table()
	tbl_A = Array{Any,2}(2,3);
	tbl_B = Array{Any,2}(2,3);
	tbl_A[1,:] = [1, [1,2], [1,1,1]]
	tbl_A[2,:] = [5, [2,3], [1,2,1]]
	tbl_B[1,:] = [10, [20,21], [1,1,1]]
	tbl_B[2,:] = [1, [1,2], [1,1,5]]
	count_table, A_to_inds, B_to_inds = build_count_table(tbl_A, tbl_B)
	@test size(count_table) == (3+1,4+1)
	@test count_table[1,1] == 0
	@test count_table[1,2] == 1
	@test count_table[1,3] == 1
	@test count_table[1,4] == 0
	@test count_table[1,5] == 0
	@test count_table[2,1] == 0
	@test count_table[2,2] == 0
	@test count_table[2,3] == 0
	@test count_table[2,4] == 1
	@test count_table[2,5] == 0
	@test count_table[3,1] == 1
	@test count_table[3,2] == 0
	@test count_table[3,3] == 0
	@test count_table[3,4] == 0
	@test count_table[3,5] == 1
	@test count_table[4,1] == 1
	@test count_table[4,2] == 0
	@test count_table[4,3] == 0
	@test count_table[4,4] == 0
	@test count_table[4,5] == 0

end

function create_dummy_count_table()
	tbl_A = Array{Any,2}(3,3);
	tbl_B = Array{Any,2}(3,3);
	tbl_A[1,:] = [1, [1,2], [1,1,1]]
	tbl_A[2,:] = [5, [2,3], [1,2,1]]
	tbl_A[3,:] = [7, [1,3], [1,3,3]]
	tbl_B[1,:] = [10, [20,21], [1,1,1]]
	tbl_B[2,:] = [1, [1,2], [1,1,5]]
	tbl_B[3,:] = [4, [20,2], [1,3,3]]
	return build_count_table(tbl_A, tbl_B)
end

function test_merge_columns()
	count_table, _, _ = create_dummy_count_table()
	@test_throws AssertionError merge_columns(count_table, 1, 2)
	@test_throws AssertionError merge_columns(count_table, 6, 2)
	@test_throws AssertionError merge_columns(count_table, 2, 1)
	@test_throws AssertionError merge_columns(count_table, 2, 6)
	count_table = merge_columns(count_table, 3, 5)
	@test size(count_table) == (4,4)
	@test count_table[1,3] == 1
	@test count_table[2,3] == 0
	@test count_table[3,3] == 1
	@test count_table[4,3] == 1
end

function test_split_column()
	count_table, _, _ = create_dummy_count_table()
	@test_throws AssertionError split_column(count_table, 1)
	@test_throws AssertionError split_column(count_table, 6)
	count_table = split_column(count_table, 4)
	@test size(count_table) == (4,6)
	@test count_table[1,4] == 0
	@test count_table[2,4] == 1
	@test count_table[3,4] == 0
	@test count_table[4,4] == 0
	@test count_table[1,5] == 0
	@test count_table[2,5] == 1
	@test count_table[3,5] == 0
	@test count_table[4,5] == 0
	@test count_table[2,1] == 0
end

function test_remove_synapse()
	count_table, _, _ = create_dummy_count_table()
	@test_throws AssertionError remove_synapse(count_table, (1,2), (2,3))
	@test_throws AssertionError remove_synapse(count_table, (6,2), (2,3))
	@test_throws AssertionError remove_synapse(count_table, (2,1), (2,3))
	@test_throws AssertionError remove_synapse(count_table, (2,6), (2,3))
	count_table = remove_synapse(count_table, (2,4), (3,5))
	@test size(count_table) == (4,5)
	@test count_table[1,1] == 0
	@test count_table[2,1] == 1
	@test count_table[3,1] == 2
	@test count_table[4,1] == 1
	@test count_table[1,4] == 0
	@test count_table[2,4] == 1
	@test count_table[3,4] == 0
	@test count_table[4,4] == 0
	@test count_table[1,5] == 0
	@test count_table[2,5] == 0
	@test count_table[3,5] == 0
	@test count_table[4,5] == 0
end

function test_add_synapse()
	count_table, _, _ = create_dummy_count_table()
	@test_throws AssertionError add_synapse(count_table, (1,2))
	@test_throws AssertionError add_synapse(count_table, (6,2))
	@test_throws AssertionError add_synapse(count_table, (2,1))
	@test_throws AssertionError add_synapse(count_table, (2,6))
	count_table = add_synapse(count_table, (2,4))
	@test size(count_table) == (4,5)
	@test count_table[1,1] == 0
	@test count_table[1,2] == 2
	@test count_table[1,3] == 1
	@test count_table[1,4] == 1
	@test count_table[1,5] == 0
end

test_edges_to_syn_dicts()
test_map_synapses()
test_get_indexed_seg_IDs()
test_build_count_table()
test_merge_columns()
test_split_column()
test_remove_synapse()
test_add_synapse()


text:
struct SDRClassifier
  predictionSize::Int
  targetSize::Int
  α::Float32
  W::Matrix{Float32}
  history_pred::Matrix{Float32}
end
SDRClassifier(predictionSize,targetSize; α,buffer_length)=
    SDRClassifier(predictionSize,targetSize,α,
        zeros(targetSize,predictionSize), zeros(targetSize,buffer_length))

"""
    predict!(classifier::SDRClassifier, Π,target::Int; enable_learning=true)=

Predict probability that the output represents each bucket of an arithmetic encoder.
0<=predict!<=1
"""
predict!(classifier::SDRClassifier, Π,target::Int; enable_learning=true)=
    predict!(classifier, Π,bitarray(target,classifier.targetSize), enable_learning=enable_learning)
function predict!(classifier::SDRClassifier, Π,target::BitArray{1}; enable_learning=true)
  A(Π)= classifier.W*Π
  saturate(x)= isnan(x) ? 1 : x
  sumsat(x)= ( s= sum(x); isinf(s) ? maximum(x) : s )
  nonlinearity(A)= saturate.(ℯ.^A ./ sumsat(ℯ.^A))
  adapt()= -classifier.α .* (classifier.history_pred[:,end] .- target)*Π[Π]'
  circshift!(history,prediction)= begin
    history[:,2:end].= history[:,1:end-1]
    history[:,1].= prediction
  end
  prediction= Π|>A|>nonlinearity
  enable_learning &&
      ( classifier.W[:,Π]+= adapt() )
  circshift!(classifier.history_pred,prediction)
  return prediction
end
"""
    reverse_simpleArithmetic(bucketProbDist, algorithm,params)

Reverses the simple arithmetic encoder, given the same parameters.
Inputs a probability distribution across all buckets and collapses it to a single arithmetic
value, representing the most likely estimation in the timeseries' domain (like a defuzzifier)

# Arguments
- bucketProbDist [nbucket]: discrete probability of each bucket [0-1]
- algorithm: {'mean','mode','highmean'} 'highmean' is the mean of the highest-estimated values
"""
function reverse_simpleArithmetic(bucketProbDist::Vector, algorithm,params)
  bucketL= params.range / (2*params.buckets)
  bucketCenter= (params.in_min + bucketL .+ 2 .* (0:params.buckets-1) .* bucketL)#|> collect
  # normalize
  bucketProbDist= (bucketProbDist .+ eps(1/length(bucketCenter))) ./ sum(bucketProbDist)
  # Only the likeliest prediction, if it stands out
  if algorithm == "mode"
    if maximum(bucketProbDist) - minimum(bucketProbDist) > 0.2
      mostLikely= bucketCenter[findmax(bucketProbDist)[2]]
    else
      mostLikely= bucketCenter'bucketProbDist
    end
  # Mean of the entire distribution
  elseif algorithm == "mean"
    mostLikely= bucketCenter'bucketProbDist
  # Mean of the most likely predictions (like "mean", but filtering)
  elseif algorithm == "highmean"
    highprob= bucketProbDist .>= quantile(bucketProbDist,0.85)
    bucketProbDist[highprob].= bucketProbDist[highprob] ./ sum(bucketProbDist[highprob])
    mostLikely= bucketCenter[highprob]'bucketProbDist[highprob]
  else
    error("[reverse_simpleArithmetic] Algorithm $algorithm doesn't exist")
  end
  return mostLikely
end


text:
# column
unit_cost_c = [6 4 3 2 1]
# row
amount_r = [2, 3, 4, 6, 12]

# row * column adds up to a sum
unit_cost_c * amount_r

using LinearAlgebra
# define as row
unit_cost = [6, 4, 3, 2, 1]
# diagnol matrix times a vector gives a vector
Diagonal(unit_cost) * amount_r

#this obviously gives same answer as row 7 since transpose makes it a column
transpose(unit_cost) * amount_r

# matrix with different orders in eacch column
amounts =  [6 1 10;
            4 1 10;
            3 1 10;
            2 1 10;
            1 1 10]

# get the summed cost for each column order
transpose(unit_cost) * amounts

# first row calculates sum based on unit costs (as before)
# second row calculates average number of units (divide by 5 and sum)
operations =   [6.0 4.0 3.0 2.0 1.0;
                0.2 0.2 0.2 0.2 0.2]
operations * amounts

# identity operation
Diagonal([1,1,1]) * [8,12,3]
# I does it better
I*[8,12,3]
I*[8,12,3,5]

# reverse operation - since this is part of base not sure why Engheim doesnt mention simple operator
R =  [0 0 1;
      0 1 0;
      1 0 0]

R * [8, 12, 3]
reverse([8, 12, 3])
# can also reverse a row
reverse([8 12 3])
# you can use multiplication on the matrix
2I*[8, 12, 3]
# more examples
z = [1 2 3;
     4 5 6;
     7 8 9]

I*z
2*z
reverse(z)
2*reverse(z)
# be careful with inverse
#inv(z) # commented out see below
# this gives you a singular exception. If you take the determinant you see it is ) which explains why
det(z)
# a simple change solves the problem
zf = [1 2 3;
      4 5 6;
      7 8 10]
det(zf)
inv(zf)
# rounding error for floating point doesnt create perfect identity
If = zf * inv(zf)
# but its close enough
I ≈ If  # yields true

# Back to Engheim
A =   [1  1 -2;
       1 -1 -1;
       1  1  3]

B = [3, 0, 12]

X = inv(A)*B

A*X
B
A*X ≈ B  # yields true

# concatenation
C = [1 3 5]
D = [7 9 11]
E = [1, 3, 5]
F = [7, 9, 11]
G = [1 3 5; 2 4 6]
H = [7 9 11; 8 10 12]

# horizontal - dims indiccate which dimention
hcat(C, D)
cat(C, D, dims=2)
hcat(E, F)
cat(E, F, dims=2)
hcat(G, H)
cat(G, H, dims=2)
hcat(G, 2G, 3G, 4G)

# vertical
vcat(C, D)
cat(C, D, dims=1)
vcat(E, F)
cat(E, F, dims=1)
vcat(G, H)
cat(G, H, dims=1)
vcat(C, D, C, D)

# creating matrices
zeros(Int8, 3, 4)
ones(Int8, 3, 4)
rand(Int8, 3, 4)
zeros(Float64, 3, 4)
ones(Float64, 3, 4)
rand(Float64, 3, 4)
fill(12, 3, 4)  #Int64
fill(12.0, 3, 4) #Float64
fill(1//2, 3, 4) #Rational{Int64}
fill(0x4, 3, 4) #Int8

# in the example on page 249 for vector addition this is what you are really seeing:
u = [4,2]  # from the origin [0,0]
v = [-1,1] # from the starting point [4,2]
u + v
# yields:
# 2-element Vector{Int64}:
# 3
# 3
# scalar addition needs to be done per element
u .+ 1
# yields:
# 2-element Vector{Int64}:
# 5
# 3
# assignment of vectors
w = v
w == v # true
w[1] = 1 # change w
w == v # true
v[2] = 2 # change v
w == v # true
# element wise assignment has different result in that w is NOT pointing to unit
# w .= u is equivalent to w[1] = 4 and w[2] = 2 which obviously has no impact on u
w .= u
w == v # true
w == u # true temporarily
v == u # true temporarily
# but
w[1] = 1 # change w
w == v # true
w == u # false
v == u # false
# subtraction
v = [-1, 1]
# next two give same result. This also shows multiplying vector by scaler works as expected
u + (-v)
u - v
# vector rotation
x = [3, 4]
xm = 3*[1, 0] + 4*[0, 1]
x == xm # true
# rotate x 90 degrees
xmr = 3*[0, 1] + 4*[-1, 0]
# do it as matrix multiplication - NB: Engheim got the matrix wrong
# [0,1] and [-1, 0] are columns not rows so the matrix looks like A below
A = [ 0 -1;
      1 0]
# a vector is a column not a row
# multiply each row in the vector x times a column in A and add the two vectors together (as in xmr)
# 3 is first row times [0, 1] + 4 is second row * [-1,0] giving vector (column) [-4, 3]
A*x == xmr # true

#rotation(degrees)
#Create a rotation matrix for rotating a 2D vector `deg` degrees.

function rotation(deg::Real)
  rad = deg2rad(deg)
  cosθ = cos(rad)
  sinθ = sin(rad)

  [cosθ -sinθ;
   sinθ  cosθ]
end

# create a Matrix for rotating a vector 90 degrees
M = rotation(90)

# create a saaple [x,y] vector and rotate 90 degrees
v = [3, 4]
M*v
v2 = [3, 3]
M*v2
v3 = [3, 0]
M*v3
v4 = [0, 3]
M*v4

# # Dot Product
# # important you add "using LinearAlgebra" before running the following
# a = [2, 3]
# b = [3, 1]

# # to get the dot product character type "a \cdot[tab] b" (i.e. thit the tab key)
# a ⋅ b
# transpose(a)*b
# dot(a, b)

# # Work equations
# r = [3, 1]
# # the norm is the sqrt of r ⋅ r which gives the length of the vector r
# # u is a unit vector with same direction as r (divide r by its lenght
# # to get a vector lenght 1 with same direction as r)
# u = r / norm(r)
# # in the graph in the book Engheim shifts teh vector from the 0origin
# F = [2, 3]
# # the work
# F ⋅ u
# # from this equality “𝐅 ⋅ 𝐮 = |𝐅| cos θ " we can calulate the angle θ
# θ = acos(F ⋅ u/norm(F))
# # check the equality
# F ⋅ u == norm(F) * cos(θ)
# # substitute back the defintion of u
# F ⋅ (r/norm(r)) == norm(F) * cos(θ)
# # multiply norm(r) to both sides of the equation (there are rounding errors)
# F ⋅ r ≈ norm(F) * norm(r) * cos(θ)
# # Cross product
# u = [6, 0, 0]
# v = [3, 4, 0]
# # ⨯ is \times
# w = u ⨯ v

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
mutable struct TreeSwaptionEngine{S <: ShortRateModel, T <: ShortRateTree} <: LatticeShortRateModelEngine{S}
  model::S
  timeSteps::Int
  common::LatticeShortRateModelEngineCommon{T}
  latticeBuilt::Bool
  # ts::Y

  # function call{S, I}(::Type{TreeSwaptionEngine}, m::S, tsteps::I)
  #   t = new{S, I, YieldTermStructure}(m, tsteps)
  #   add_observer!(m, t)
  #
  #   return t
  # end
  #
  # call{S, I}(::Type{TreeSwaptionEngine}, m::S, tsteps::I, l::LatticeShortRateModelEngineCommon) = new{S, I, YieldTermStructure}(m, tsteps, l)
  #
  # call{S, I, Y}(::Type{TreeSwaptionEngine}, m::S, tsteps::I, l::LatticeShortRateModelEngineCommon, ts::Y) = new{S, I, T, Y}(m, tsteps, l, ts)
  function TreeSwaptionEngine{S, T}(model::S,
                              timeSteps::Int,
                              common::LatticeShortRateModelEngineCommon{T},
                              latticeGen::Bool = true) where {S, T}
    ts = new{S, T}(model, timeSteps, common, latticeGen)
    add_observer!(model, ts)

    return ts
  end

  # function TreeSwaptionEngine{S}(model::S, timeSteps::Int)
  #   # create empty tree
  #   tg = TimeGrid([1.0], 1)
  #   lattice = tree(model, tg)
  #   common = LatticeShortRateModelEngineCommon{typeof(lattice)}(tg, lattice)
  #   ts = new{S, typeof(lattice)}(model, timeSteps, common, false)
  #   add_observer!(model, ts)
  #
  #   return ts
  # end
end

function TreeSwaptionEngine(model::S, tg::TimeGrid) where {S <: ShortRateModel}
  lattice = tree(model, tg)
  return TreeSwaptionEngine{S, typeof(lattice)}(model, 0, LatticeShortRateModelEngineCommon(tg, lattice))
end

function TreeSwaptionEngine(model::S, timeSteps::Int) where {S <: ShortRateModel}
  # create empty tree
  tg = TimeGrid([1.0], 1)
  lattice = tree(model, tg)
  common = LatticeShortRateModelEngineCommon{typeof(lattice)}(tg, lattice)
  return TreeSwaptionEngine{S, typeof(lattice)}(model, timeSteps, common, false)
end

# methods
function _calculate!(pe::TreeSwaptionEngine, swaption::Swaption)
  tsmodel = pe.model

  refDate = reference_date(tsmodel.ts)
  dc = tsmodel.ts.dc

  dSwaption = DiscretizedSwaption(swaption, refDate, dc, pe.common.lattice.treeLattice)

  if ~pe.latticeBuilt
    times = mandatory_times(dSwaption)
    tg = TimeGrid(times, pe.timeSteps)
    pe.common.tg = tg
    update!(pe)
    pe.latticeBuilt = true
  end
  lattice = pe.common.lattice

  stoppingTimes = zeros(length(swaption.exercise.dates))
  # @simd for i in eachindex(stoppingTimes)
  #   @inbounds stoppingTimes[i] = year_fraction(dc, refDate, swaption.exercise.dates[i])
  # end
  map!(x -> year_fraction(dc, refDate, x), stoppingTimes, swaption.exercise.dates)

  initialize!(dSwaption, lattice.treeLattice, stoppingTimes[end])

  nextExerciseIdx = findnext(greater_than_or_equal_to, stoppingTimes, 1, 0.0)

  nextExercise = nextExerciseIdx != 0 ? stoppingTimes[nextExerciseIdx] : stoppingTimes[end]

  rollback!(dSwaption, nextExercise)
  swaption.results.value = present_value(dSwaption)
end


text:
function instantiate_plot(θ, p_θ, E, γ)

    phasespace_position = Node([SVector(0.0, 0.0, 0.0)])
    onany(θ, p_θ, E) do θ_val, p_θ_val, E_val
        phasespace_position[] = [SVector(θ_val, p_θ_val, E_val)]
    end

    physical_position = Node([SVector(0.0, 0.0), SVector(0.0, 0.0)])
    on(θ) do θ_val
        physical_position[] = [SVector(0.0, 0.0), SVector(r(θ_val)...)]
    end

    θ_ax   = range(-1.5π, 1.5π, length = 100)
    p_θ_ax = range(-1.5p_θ₀(2), 1.5p_θ₀(2), length = 100)

    s1 = Scene()
    surface!(s1, θ_ax, p_θ_ax, energy, alpha=0.1, transparency = true)
    surf = s1[end]
    contour3d!(s1, θ_ax, p_θ_ax, energy, levels = [E.val], color=:orange, linewidth=5)
    cont = s1[end]
    map!(E -> [E], cont.attributes.levels, E)
    scatter!(s1, phasespace_position, color=:white, markersize=0.2)
    scat = s1[end]

    s1[Axis].attributes.names.axisnames[] = ("θ", "pθ", "E")
    surf[:colormap] = :RdYlBu
    scale!(s1, 1, 1/2, 1/20)

    s2 = Scene(resolution=(500, 500))
    scatter!(s2, physical_position, limits=FRect(-1.0, -1.0, 2.0, 2.0))
    lines!(s2, physical_position, limits=FRect(-1.0, -1.0, 2.0, 2.0))

    s3, slider_obs = AbstractPlotting.textslider(0.01:0.01:3, "Energía")
    Observables.connect!(slider_obs, γ)

    scene = vbox(s1, hbox(s2,s3))
    display(scene)

    on(scene.events.window_area) do val
        update_cam!(scene.children[1], FRect3D(Vec3f0(0), Vec3f0(1)))
    end

    return scene
end

max_stars_repo_name: tlqmj/phasespace

max_issues_repo_name: tlqmj/phasespace

max_forks_repo_name: tlqmj/phasespace


text:
module TestUnordered
include("../preamble.jl")
using Transducers: append_unordered!, channel_unordered, reduce_commutative

@testset "reduce_commutative" begin
    @testset "Map(inc)(+)" begin
        @test reduce_commutative(+, Map(inc), 1:100) == sum((1:100) .+ 1)
    end
end

@testset "append_unordered!" begin
    input = Channel(Map(identity), 1:100)
    output = Channel{Int}(Inf)
    append_unordered!(output, Map(inc), input)
    close(output)
    ys = sort!(collect(output))
    @test ys == (1:100) .+ 1
end

@testset "append_unordered!(output, (x + 1 for x in input))" begin
    input = Channel(Map(identity), 1:100)
    output = Channel{Int}(Inf)
    append_unordered!(output, (x + 1 for x in input))
    close(output)
    ys = sort!(collect(output))
    @test ys == (1:100) .+ 1
end

@testset "channel_unordered" begin
    @testset "Map(inc)" begin
        input = Channel(Map(identity), 1:100)
        ys = sort!(collect(channel_unordered(Map(inc), input)))
        @test ys == (1:100) .+ 1
    end
    @testset "x + 1 for x in input" begin
        input = Channel(Map(identity), 1:100)
        ys = sort!(collect(channel_unordered(x + 1 for x in input)))
        @test ys == (1:100) .+ 1
    end
    @testset "eduction(x + 1 for x in input)" begin
        input = Channel(Map(identity), 1:100)
        ys = sort!(collect(channel_unordered(eduction(x + 1 for x in input))))
        @test ys == (1:100) .+ 1
    end
    @testset "ReduceIf(isodd): early termination" begin
        input_data = fill(2, 100)
        input_data[10] = 1
        input = Channel(Map(identity), input_data)
        ys = collect(channel_unordered(ReduceIf(isodd), input))
        @test count(x -> x == 1, ys) == 1
        @test sort!(unique(ys)) == [1, 2]
        @test length(ys) < length(input_data)
    end
    @testset "error handling (no dangling tasks)" begin
        input_length = 100
        input = Channel(Map(identity), 1:input_length)
        trace = Channel(input_length)
        ntasks = 10
        xf = Map() do x
            put!(trace, x)

            # Try to make sure all tasks hits above line.  This may
            # not be required since `put!` would invoke task switch.
            # But let's play on the safe side:
            sleep(0.1)

            x
        end |> Map() do x
            error("Throwing error with x = $x")
        end
        @test_throws Exception take!(channel_unordered(xf, input; ntasks=ntasks))
        close(trace)
        consumed = sort!(collect(trace))
        @test consumed == 1:ntasks
    end
    @testset "error handling (terminate other tasks)" begin
        input_length = 100
        input = Channel(Map(identity), 1:input_length)
        trace = Channel(input_length)
        ntasks = 10
        xf = Map() do x
            put!(trace, x)
            sleep(0.01)
            x
        end |> Map() do x
            if x == 1
                sleep(0.02)
                error("Throwing error with x = $x")
            end
            x
        end
        @test_throws Exception collect(channel_unordered(xf, input; ntasks=ntasks))
        close(trace)
        consumed = sort!(collect(trace))
        @test length(consumed) > ntasks
        @test length(consumed) < input_length
    end
end

end  # module


text:
function maxpool1d(ip, kernel, strides=1, pads=0)
    l = size(ip)[2]
    ip1 = ip[:, 1]
    temp = []
    for x=1:length(ip1)-kernel+1
        push!(temp, findmax(ip1[x:x+kernel])[1])
    end
    return temp
end


text:
import Base: .+, .-, .*, ./, .^, max, min, .==, .>, .>=, .<, .<=, +, -, *, /, \

cuda01 = [
("add",".+","s+xi"),
("sub",".-","s-xi"),
("mul",".*","s*xi"),
("div","./","s/xi"),
("pow",".^","pow(s,xi)"),
("max","max","(xi>s?xi:s)"),
("min","min","(xi<s?xi:s)"),
("eq",".==","xi==s"),
("gt",".>","xi>s"),
("ge",".>=","xi>=s"),
("lt",".<","xi<s"),
("le",".<=","xi<=s"),
# "hypot",
# "rhypot",
# "atan2",
# "frexp",
# "ldexp",
# "scalbn",
# "scalbln",
# "jn",
# "yn",
# "fmod",
# "remainder",
# "mod",
# "fdim",
]

# ambiguity fixes:
max{T<:Real,S<:Real}(a::KnetArray{T},s::S)=max(T(s),a)
max{T<:Real,S<:Real}(s::S,a::KnetArray{T})=max(T(s),a)
min{T<:Real,S<:Real}(a::KnetArray{T},s::S)=min(T(s),a)
min{T<:Real,S<:Real}(s::S,a::KnetArray{T})=min(T(s),a)
max{T<:Real,S<:Number}(a::KnetArray{T},s::S)=max(T(s),a)
max{T<:Real,S<:Number}(s::S,a::KnetArray{T})=max(T(s),a)
min{T<:Real,S<:Number}(a::KnetArray{T},s::S)=min(T(s),a)
min{T<:Real,S<:Number}(s::S,a::KnetArray{T})=min(T(s),a)
(+)(a::KnetArray{Bool},s::Bool)=(.+)(s,a)
(+)(s::Bool,a::KnetArray{Bool})=(.+)(s,a)
(-)(a::KnetArray{Bool},s::Bool)=(.+)(-s,a)
(-)(s::Bool,a::KnetArray{Bool})=(.-)(s,a)
.^(x::Base.Irrational{:e}, a::KnetArray)=.^(float(x),a)

# For array,scalar we can get some for free:
# Only type corrected number,array need implementing for basic arithmetic:
(.+){T}(a::KnetArray{T},s::Number)=(.+)(T(s),a)
(.+){T}(s::Number,a::KnetArray{T})=(.+)(T(s),a)
(.-){T}(a::KnetArray{T},s::Number)=(.+)(T(-s),a)
(.-){T}(s::Number,a::KnetArray{T})=(.-)(T(s),a)
(.*){T}(a::KnetArray{T},s::Number)=(.*)(T(s),a)
(.*){T}(s::Number,a::KnetArray{T})=(.*)(T(s),a)
(./){T}(a::KnetArray{T},s::Number)=(.*)(T(1/s),a)
(./){T}(s::Number,a::KnetArray{T})=(./)(T(s),a)
#(.^){T}(a::KnetArray{T},s::Number) # cannot convert to an s,a operation
(.^){T}(s::Number,a::KnetArray{T})=(.^)(T(s),a)
max{T}(a::KnetArray{T},s::Number)=max(T(s),a)
max{T}(s::Number,a::KnetArray{T})=max(T(s),a)
min{T}(a::KnetArray{T},s::Number)=min(T(s),a)
min{T}(s::Number,a::KnetArray{T})=min(T(s),a)

.=={T}(a::KnetArray{T},s::Number)=(T(s).==a)
.=={T}(s::Number,a::KnetArray{T})=(T(s).==a)
.>{T}(a::KnetArray{T},s::Number)=(T(s).<a)
.>{T}(s::Number,a::KnetArray{T})=(T(s).>a)
.>={T}(a::KnetArray{T},s::Number)=(T(s).<=a)
.>={T}(s::Number,a::KnetArray{T})=(T(s).>=a)
.<{T}(a::KnetArray{T},s::Number)=(T(s).>a)
.<{T}(s::Number,a::KnetArray{T})=(T(s).<a)
.<={T}(a::KnetArray{T},s::Number)=(T(s).>=a)
.<={T}(s::Number,a::KnetArray{T})=(T(s).<=a)

# familiar aliases for broadcasting operations of array & scalar (#7226):
(+){T}(a::KnetArray{T},s::Number)=(.+)(T(s),a)
(+){T}(s::Number,a::KnetArray{T})=(.+)(T(s),a)
(-){T}(a::KnetArray{T},s::Number)=(.+)(T(-s),a)
(-){T}(s::Number,a::KnetArray{T})=(.-)(T(s),a)
(*){T}(a::KnetArray{T},s::Number)=(.*)(T(s),a)
(*){T}(s::Number,a::KnetArray{T})=(.*)(T(s),a)
(/){T}(a::KnetArray{T},s::Number)=(.*)(T(1/s),a)
(\){T}(s::Number,a::KnetArray{T})=(.*)(T(1/s),a)
#(/){T}(s::Number,a::KnetArray{T})=(.*)(T(1/s),a) # not defined in base
#(^){T}(a::KnetArray{T},s::Number)=(.^)(a,T(s)) # linalg
#(^){T}(s::Number,a::KnetArray{T})=(.^)(T(s),a) # linalg

function cuda01def(f, j=f, o...)
    J=Symbol(j)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$(S)_01"
        @eval begin
            function $J(s::$T,x::KnetArray{$T})
                y = similar(x)
                ccall(($F,$libknet8),Void,(Cint,$T,Ptr{$T},Ptr{$T}),length(y),s,x,y)
                return y
            end
        end
    end
end

#if isdefined(:libknet8)
    for f in cuda01
        isa(f,Tuple) || (f=(f,))
        cuda01def(f...)
    end
#end


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

using Documenter, LAR

makedocs(
	format = :html,
	sitename = "LAR.jl",
	assets = ["assets/lar.css", "assets/logo.png"],
	pages = [
		"Home" => "index.md",
	]
)


text:
# Adapted from "Example: Neural network" by Kjell Magne Fauske:
# https://texample.net/tikz/examples/neural-network/
# which is licensed under CC BY 2.5:
# https://creativecommons.org/licenses/by/2.5/

TikzPicture(
    chain::Chain;
    inputs::Vector{String}=Vector{String}(),
    outputs::Vector{String}=Vector{String}(),
    outputname="Output",
    kwargs...,
) = TikzPicture(Network(chain, inputs, outputs, outputname); kwargs...)

function TikzPicture(network::Network; xscale=1.3, yscale=1, options="", kwargs...)
    inputs = escaped.(network.inputs)
    outputs = escaped.(network.outputs)
    widths = [layer.inputs for layer in network.layers]
    names = escaped.([layer.name for layer in network.layers])

    nndepth = length(network.layers)
    nnwidth = maximum(widths)

    x(i) = xscale * i
    y(i, width) = yscale * (0.5 + 0.5width - i)

    xstart = 0
    xend = x(nndepth - 1)
    ystart = -y(1, nnwidth)
    yend = y(1, nnwidth)

    yannot = yend + 1

    xymin = (xstart - 2, ystart - 1)
    xymax = (xend + 2, yend + 1.75)

    ypositions(width) = y.(1:width, width)
    hiddenlist(width) = tikzformat(enumerate(ypositions(width)), tuplefmt)
    connectlist(layer) = tikzformat(layer.connect, tuplefmt)

    inputlist = enumerate(zip(ypositions(widths[1]), inputs))
    hiddenlists = enumerate(zip(x.(1:nndepth - 2), hiddenlist.(widths[2:end - 1])))
    outputlist = enumerate(zip(ypositions(widths[end]), outputs))
    connectlists = enumerate(connectlist.(network.layers))
    annotlist = zip(x.(0:length(network.layers) - 1), names)

    tikzinputs = tikzformat(inputlist, inputfmt)
    tikzhidden = tikzformat(hiddenlists, layerfmt)
    tikzoutputs = tikzformat(outputlist, outputfmt)
    tikzconnect = tikzformat(connectlists, connectfmt)
    tikzannot = tikzformat(annotlist, tuplefmt)

    data = """
        \\fill[white, use as bounding box] $xymin rectangle $xymax;

        \\foreach \\j/\\y/\\input in {$tikzinputs}
            \\node[input neuron,pin=left:\\input] (1-\\j) at (0,\\y) {};

        \\foreach \\i/\\x/\\tikzhidden in {$tikzhidden}
            \\foreach \\j/\\y in \\tikzhidden
                \\node[hidden neuron] (\\i-\\j) at (\\x,\\y) {};

        \\foreach \\j/\\y/\\output in {$tikzoutputs}
            \\node[output neuron,pin=right:\\output] ($nndepth-\\j) at ($xend,\\y) {};

        \\foreach \\i/\\k/\\connectlist in {$tikzconnect}
            \\foreach \\j/\\l in \\connectlist
                \\path (\\i-\\j) edge (\\k-\\l);

        \\foreach \\x/\\txt in {$tikzannot}
            \\node[annot] at (\\x,$yannot) {\\txt};
    """
    options = """
        neuron/.style={circle,draw,minimum size=17pt},
        input neuron/.style={neuron,fill=green!50},
        hidden neuron/.style={neuron,fill=blue!50},
        output neuron/.style={neuron,fill=red!50},
        annot/.style={text width=4em, align=center},
        $(isempty(options) ? '%' : options)
    """
    TikzPicture(data; options=options, kwargs...)
end

escaped(str) = replace(str, '#' => "\\#")

tikzformat(A, fmt) = join(fmt.(A), ',')

tuplefmt((a, b)) = "$a/$b"

inputfmt((j, (y, input))) = "$j/$y/\\pgftext[right] {$input}"

layerfmt((i, (x, tikzhidden))) = "$(i + 1)/$x/{$tikzhidden}"

outputfmt((j, (y, output))) = "$j/$y/\\pgftext[left] {$output}"

connectfmt((i, layerconnect)) = "$i/$(i + 1)/{$layerconnect}"


text:
module TypeTrees
##
# Generate a text graphic of Julia modules type tree
##

# The node type holds the type of the cuurent node and a dict of subtypes
type TTNode
    strname::String
    typ::Type
    subtypes::Dict{String, TTNode}

    TTNode(sname::String, t::Type) = new(sname, t, Dict{String, TTNode}())
end

# Add a node to a dict if not added
function add_ttnode(subtypes::Dict{String, TTNode}, sname::String, tnode::TTNode)
    ret = get(subtypes, sname, nothing)
    (nothing == ret) && (ret = subtypes[sname] = tnode)
    ret
end

function add_ttnode(subtypes::Dict{String, TTNode}, sname::String, t::Type)
    ret = get(subtypes, sname, nothing)
    (nothing == ret) && (subtypes[sname] = ret = TTNode(sname, t))
    ret
end

# Get a string name for the type
typ_name(t::UnionType) = string(t)
typ_name(t::TypeConstructor) = string(t)
typ_name(t) = string(t.name)

# Store a type and its type hierarchy chain
# Recurse till we reach the top level type
function store_type(sname::String, t::UnionType)
    suptype = UnionType
    tnode = TTNode(sname, t)

    # store unions under UnionType type
    subtypes = store_type(typ_name(suptype), suptype)
    add_ttnode(subtypes, sname, tnode)

    # unions are also in a sense related to the types of their components
    for suptype = t.types
        subtypes = store_type(typ_name(suptype), suptype)
        add_ttnode(subtypes, sname, tnode)
    end

    return tnode.subtypes
end

function store_type(sname::String, t::TypeConstructor)
    suptype = t.body
    subtypes = store_type(typ_name(suptype), suptype)
    tnode = add_ttnode(subtypes, sname, t)
    return tnode.subtypes
end

function store_type(sname::String, t::DataType)
    suptype = super(t)
    subtypes = (suptype != t) ? store_type(typ_name(suptype), suptype) : types_tree
    tnode = add_ttnode(subtypes, sname, t)
    return tnode.subtypes
end

function store_type(sname::String, t::Tuple)
    tnode = add_ttnode(types_tree, sname, t)
    return tnode.subtypes
end

function store_type(sname::String, t)
    suptype = super(t)
    subtypes = (suptype != t) ? store_type(typ_name(suptype), suptype) : types_tree
    tnode = add_ttnode(subtypes, sname, t)
    return tnode.subtypes
end

# examine all symbols in module and store those that are types
function store_all_from(m::Module)
    for expr = names(m,true)
        try
            t = eval(m,expr)
            isa(t, Type) && store_type(string(expr), t)
        #catch ex
        #    println("Error adding ", string(expr), m, " (", ex, ")")
        end
    end
end

type_props(typ) = ""
type_props(typ::DataType) = string("<<",
                                 typ.abstract    ? " abstract"    : " concrete",
                                 typ.mutable     ? " mutable"     : " immutable",
                                 typ.pointerfree ? " pointerfree" : "",
                                 " size:", typ.size,
                                 " >>")

function print_tree(subtypes::Dict{String, TTNode}, pfx::String="")
    for n in sort!([keys(subtypes)...])
        v = subtypes[n]
        if(n == string(v.typ))
            println(pfx, "+- ", n, " ", type_props(v.typ))
        else
            println(pfx, "+- ", n, " = ",  v.typ, " ", type_props(v.typ))
        end
        print_tree(v.subtypes, pfx * ".  ")
    end
end

# TODO: optionally take module names in command line
# TODO: sort output
# TODO: option to list subtrees of type tree, or other symbol types
const types_tree = Dict{String, TTNode}()

for m in (Base, Core, Main)
    store_all_from(m)
end

# print_tree(types_tree)

end # module

max_stars_repo_name: alanedelman/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
# TODO: this should be more robust, ideally.
using Logging
import WebIO, HTTP, WebSockets

const pages = Dict{String,Any}()
const server = Ref{Any}()
const port = Ref{Int}(9000)

function isrunning(server)
    isassigned(server) && !istaskdone(server[].serve_task)
end

# Be ready for the deprecation of tohtml :P
# once the no unsafe-script lands
if isdefined(WebIO, :tohtml)
    const tohtml = WebIO.tohtml
else
    tohtml(io, app) = show(io, MIME"text/html"(), app)
end

function routepages(req)
    target = req.target[2:end]

    haskey(pages, target) || return missing

    io = IOBuffer()

    print(io, """
        <!doctype html>
        <html>
        <head>
        <meta charset="UTF-8">
        <script>window._webIOWebSocketURL = 'ws://localhost:$(port[])/webio-juno-websocket';</script>
        <script src=$(repr(WebIO.bundle_key))></script>
        <style>
            body {
                width: 100vw;
                height: 100vh;
                margin: 0;
            }
            .webio-scope, webio-mountpoint {
                width: 100%;
                height: 100%;
            }
        </style>
        </head>
        <body>
    """)
    tohtml(io, pages[target])
    print(io, """
        </body>
        </html>
    """)

    return HTTP.Response(
        200,
        ["Content-Type" => "text/html"],
        body = String(take!(io))
    )
end

function Base.show(io::IO, ::MIME"application/prs.juno.plotpane+html", n::Union{WebIO.Node, WebIO.Scope, WebIO.AbstractWidget})
    global pages, server
    id = rand(UInt128)
    pages[string(id)] = WebIO.render(n)

    if !isrunning(server)
        server[] = setup_server()
    end
    print(io, "<meta http-equiv=\"refresh\" content=\"0; url=http://localhost:$(port[])/$(id)\"/>")
end

function setup_server()
    port[] = rand(8000:9000)

    server = WebIO.WebIOServer(
        routepages,
        http_port = port[],
        singleton = false,
        websocket_route = "/webio-juno-websocket"
    )

    return server
end


text:
#-------------------------------------------------------------# Cross Validation
"""
`StatLearnCV(o::StatLearn, xtest, ytest)`

Automatically tune the regularization parameter λ for `o` by minimizing loss on
test data `xtest`, `ytest`.

```julia
sl = StatLearn(size(x, 2), LassoPenalty(.1))
o = StatLearnCV(sl, xtest, ytest)
fit!(o, x, y)
```
"""
type StatLearnCV{T<:Real, S<:Real, W<:Weight} <: OnlineStat{XYInput}
    o::StatLearn
    burnin::Int
    xtest::AMat{T}
    ytest::AVec{S}
    weight::W
end
function StatLearnCV(o::StatLearn, xtest::AMat, ytest::AVec, burnin = 1000; wgt::LearningRate = LearningRate())
    @assert length(o.β) == size(xtest, 2) "number of predictors doesn't match test data"
    StatLearnCV(o, burnin, xtest, ytest, wgt)
end

function _fit!(o::StatLearnCV, x::AVec, y::Real, γ::Float64)
    if nobs(o) < o.burnin
        _fit!(o.o, x, y, γ)
    else
        tuneλ!(o.o, x, y, γ, o.xtest, o.ytest)
    end
    o
end
# function _fit!(o::StatLearnCV, x::AMat, y::AVec, γ::Float64)
#     if nobs(o) < o.burnin
#         _fit!(o.o, x, y, γ)
#     else
#         tuneλ!(o.o, x, y, γ, o.xtest, o.ytest)
#     end
#     o
# end

function tuneλ!{A<:Algorithm, M<:ModelDefinition}(
        o::StatLearn{A, M, NoPenalty}, x, y, γ, xtest, ytest
    )
    fit!(o, x, y)
end
function tuneλ!(o::StatLearn, x, y, γ, xtest, ytest)
    o_l = copy(o)
    o_h = copy(o)
    o_l.penalty.λ = max(0.0, o_l.penalty.λ - γ)
    o_h.penalty.λ += γ

    fit!(o, x, y)
    fit!(o_l, x, y)
    fit!(o_h, x, y)

    loss_low = loss(o_l, xtest, ytest)
    loss_mid = loss(o, xtest, ytest)
    loss_hi = loss(o_h, xtest, ytest)
    best = min(loss_low, loss_mid, loss_hi)

    if best == loss_low # o_l is winner
        o.β0 = o_l.β0
        o.β = o_l.β
        o.algorithm = o_l.algorithm
        o.penalty.λ = o_l.penalty.λ
    elseif best == loss_hi # o_h is winner
        o.β0 = o_h.β0
        o.β = o_h.β
        o.algorithm = o_h.algorithm
        o.penalty.λ = o_h.penalty.λ
    end
end
nobs(o::StatLearnCV) = nobs(o.o)
nups(o::StatLearnCV) = nups(o.o)
value(o::StatLearnCV) = value(o.o)
coef(o::StatLearnCV) = coef(o.o)
predict(o::StatLearnCV, x) = predict(o.o, x)
loss(o::StatLearnCV, x, y) = loss(o.o, x, y)
loss(o::StatLearnCV) = loss(o, o.xtest, o.ytest)

function Base.show(io::IO, o::StatLearnCV)
    printheader(io, "StatLearnCV")
    print_item(io, "burnin", o.burnin)
    print_item(io, "loss", loss(o))
    show(io, o.o)
end


text:
using Test, OffsetArrays

@testset "Tridiag" begin
    d  = 2 .+ rand(5)
    dl = rand(4)
    du = rand(4)
    M = Tridiagonal(dl, d, du)
    F = lu(M)
    src = rand(5,5,5)
    for dim = 1:3
        dest1 = mapslices(x->ldiv!(F, x), copy(src), dims=dim)
        dest2 = similar(src)
        AxisAlgorithms.A_ldiv_B_md!(dest2, F, src, dim)
        @test dest1 ≈ dest2
    end

    src = OffsetArray(src, -2:2, 1:5, 0:4)
    dest1 = mapslices(x->ldiv!(F, x), copy(src), dims=2)
    dest2 = similar(src)
    AxisAlgorithms.A_ldiv_B_md!(dest2, F, src, 2)
    @test dest1 ≈ dest2
end

max_stars_repo_name: timholy/FilterMD

max_issues_repo_name: timholy/FilterMD

max_forks_repo_name: timholy/FilterMD


text:
using Graphs, SimpleWeightedGraphs, StatsBase

"""
"""
function get_adj_dict(G, A)
    adjacency_dict = Dict()
    for i in range(1, nv(G))
        adjacency_dict[i] = findall(>(0), A[i,:])
    end
    return adjacency_dict
end

"""
    Compute all k-edge walks from start_node
"""
function all_walks(adjacency_dict, start_node, k)
    curr_walks = Vector{Tuple}()
    neighbors = adjacency_dict[start_node]
    for ne in neighbors
        edge = Tuple((start_node, ne))
        push!(curr_walks, Tuple(edge))
    end
    for i in range(2, k)
        new_walks = Vector{Tuple}()
        for walk in curr_walks
	    neighbors = adjacency_dict[walk[end]]
            if length(neighbors) > 0
                for ne in neighbors
                    push!(new_walks, (walk..., ne))
                end
            end
        end
        curr_walks = new_walks
    end
    return curr_walks
end

"""
    Compute walks_per_node k-edge random walks through G
    starting from start_node.
"""
function random_walks(G, start_node::Integer, k::Integer, walks_per_node::Integer, verbose::Bool=false)
    curr_walks = Vector{Tuple}()
    max_tries = walks_per_node*4
    for i in range(1, walks_per_node)
	tries = 0
        walk = Tuple(randomwalk(G, start_node, k+1))
        while length(walk) < k+1
            walk = randomwalk(G, start_node, k+1)
	    tries += 1
	    if tries == max_tries && verbose
		println("Couldn't find a walk in $max_tries tries. \
			Returning nothing.")
		return nothing
	    end
        end
	push!(curr_walks, Tuple(walk))
    end
    return curr_walks
end

"""
    Get the out-degree weighted probability of each walk
"""
function get_walk_probabilities!(walks, walk_weight_lookup, out_degrees, weight_walks)
    if weight_walks
	    walk_weights = Vector{Float64}()
	    for w in walks
		if !haskey(walk_weight_lookup, w)
		    prod = out_degrees[w[1]]
		    for i in range(2, length(w)-1)
			prod *= out_degrees[i]
		    end
		    walk_weight_lookup[w] = prod
		end
		push!(walk_weights, walk_weight_lookup[w])
	    end
	    walk_weights /= sum(walk_weights)
    else
        walk_probabilities = Vector{Float64}(range(1, length(walks)))
        fill!(walk_probabilities, 1. / length(walks))
    end
    return walk_weights
    end

function sample_start_node(nodes, node_probabilities, bias_nodes)
    if bias_nodes
        start_node = sample(nodes, node_probabilities)[1]
    else
        start_node = rand(nodes)[1]
    end
    return start_node
end

"""
"""
function filter_nodes(G, k)
    # Compute path counts from adjacency matrix
    A = adjacency_matrix(G)
    A_bin = adjacency_matrix(G)
    A_bin[findall(>(0), A)] .= 1
    path_counts = sum(A^k, dims=2)
    # We are only interested in nodes who have k-edge walks
    nodes = [u for u in range(1, nv(G)) if path_counts[u] > 0]
    println("Total nodes: $(nv(G)) nodes with k-edge walks $(length(nodes))")
    return nodes, A, A_bin, path_counts
end

function compute_node_probs(nodes, path_counts)
    node_weights = Vector()
    for node in nodes
	push!(node_weights, path_counts[node])
    end
    node_probabilities = ProbabilityWeights(node_weights / sum(node_weights))
    return node_probabilities
end

"""
Helper function for uniform walk sampler.
Mutates walks_by_node.
"""
function update_rw_sample!(G, start_node, k, walks_per_node, walks_by_node, path_counts)
    if !haskey(walks_by_node, start_node)
        # Compute the first set of walks from start_node
	curr_walks = random_walks(G, start_node, k, walks_per_node)
	while curr_walks == nothing
	    curr_walks = random_walks(G, start_node, k, walks_per_node)
	end
    	walks_by_node[start_node] = curr_walks
    else
        # add more walks
	if length(walks_by_node[start_node]) < path_counts[start_node]
	    curr_walks = random_walks(G, start_node, k, walks_per_node)
            while curr_walks == nothing
                curr_walks = random_walks(G, start_node, k, walks_per_node)
            end
            # Only add new walks
            curr_walks = [w for w in curr_walks if !(w in walks_by_node[start_node])]
            if length(curr_walks) > 0
	    	append!(walks_by_node[start_node], curr_walks)
            end
        end
    end
    curr_walks = walks_by_node[start_node]
    return curr_walks
end

"""
    Sample from the complete DeBruijn graph defined by G.
"""
function uniform_walk_sample(G::SimpleWeightedDiGraph, k::Integer, num_walks::Integer,
        walks_per_node::Integer, weight_walks::Bool, bias_nodes::Bool, verbose::Bool=false)
    # Filter nodes based on whether they
    # originate any k-edge walks
    nodes, A, A_bin, path_counts = filter_nodes(G, k)
    adjacency_dict = get_adj_dict(G, A)
    # dims=2 is row sum (out deg/walks)
    out_degrees = sum(A_bin, dims=2)
    # Compute weights for nodes
    node_probabilities = compute_node_probs(nodes, path_counts)

    walks = Vector{Tuple}(undef, num_walks)
    walks_by_node = Dict()
    walk_weight_lookup = Dict()
    for w in range(1, num_walks)
	start_node = sample_start_node(nodes, node_probabilities, bias_nodes)
        if walks_per_node > 0
	    curr_walks = update_rw_sample!(G, start_node, k, walks_per_node, walks_by_node, path_counts)
        else
            curr_walks = all_walks(adjacency_dict, start_node, k)
            walks_by_node[start_node] = curr_walks
        end

        if length(curr_walks) > 0
            # Sample one walk from walks_by_node
	    walk_probabilities = get_walk_probabilities!(curr_walks, walk_weight_lookup,
                                                     out_degrees, weight_walks)

            walk = sample(curr_walks, ProbabilityWeights(walk_probabilities))
            walks[w] = walk
	end
    end
    if verbose
	println("Computed $(length(walks)) length-$k walks out of $num_walks requested.")
    end
    return walks
end


text:
create table categories (
  id integer not null,
  catname varchar(40) not null,
  primary key(id)
);

create table quotes (
  id      integer not null,
  cid     integer not null,
  author  varchar(100),
  quoname varchar(250) not null,
  primary key(id)
);

max_stars_repo_path: Module 3/Chapter09/apt.jl

max_issues_repo_path: Module 3/Chapter09/apt.jl

max_forks_repo_path: Module 3/Chapter09/apt.jl


text:
testOnAllContexts("PSO") do ctx
    println("\tsolves a simple problem")
    testSimpleProblem(
        comparer ->
        begin
            pso = PSO(ctx, 60, comparer, dna -> dna)
            @test pso.dnaSize == 60
            @test pso.populationSize == 100
            @test pso.vMax == 0.1f0
            @test pso.w == 0.8f0
            @test pso.c1 == 2f0
            @test pso.c2 == 2f0
            @test pso.mutationChance == 0.02f0
            @test pso.mutationStrength == 0.05f0
            pso
        end)

    println("\tsolves a problem with synthesizer")
    testSimpleProblem(
        ctx,
        (comparer, synth) ->
        begin
            pso = PSO(
                ctx,
                dnaSize!(synth),
                comparer,
                dna -> decode(synth, dna),
                vMax = 0.9f0,
                w = 0.4f0,
                mutationChance = 0.01f0,
                mutationStrength = 0.02f0)
            @test pso.dnaSize == dnaSize!(synth)
            @test pso.populationSize == 100
            @test pso.vMax == 0.9f0
            @test pso.w == 0.4f0
            @test pso.c1 == 2f0
            @test pso.c2 == 2f0
            @test pso.mutationChance == 0.01f0
            @test pso.mutationStrength == 0.02f0
            pso
        end)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
function get_scen_vars(pt::Point, scenario::String)
    a = Iterators.filter(x->contains(x.name, scenario), keys(pt))
    return Point(Dict([k=>pt[k] for k in a]))
end

function get_volt_vars(pt::Point)
    a = Iterators.filter(x->contains(x.name, "VOLT"), keys(pt))
    return Point(Dict([k=>pt[k] for k in a]))
end

function get_bin_vars(pt::Point)
    a = Iterators.filter(x->contains(x.name, "Bin"), keys(pt))
    return Point(Dict([k=>pt[k] for k in a]))
end

function get_prod_vars(pt::Point)
    a = Iterators.filter(x->contains(x.name, "Sgen"), keys(pt))
    return Point(Dict([k=>pt[k] for k in a]))
end

function get_delta_var(pt::Point, scenario::String)
    a = Iterators.filter(x->contains(x.name, get_delta_varname(scenario)), keys(pt))
    return Point(Dict([k=>pt[k] for k in a]))
end

function get_splitted_Cpt(pt, scenario)
    vars = get_scen_vars(pt, scenario)
    volt_vars = real2cplx(get_volt_vars(vars))
    bin_vars = get_bin_vars(vars)
    prod_vars = real2cplx(get_prod_vars(vars))

    return volt_vars, bin_vars, prod_vars
end

###################
## Plot functions
###################

function plot_Volt_vars(point, scenarios)
    r1 = 0.9
    plot(r1*cos.(-π:0.01:π), r1*sin.(-π:0.01:π), label=:Vmin)
    r2 = 1.1
    plot!(r2*cos.(-π:0.01:π), r2*sin.(-π:0.01:π), label=:Vmax)
    for sc in scenarios
        println(sc)
        volt_vars, _, _ = get_splitted_Cpt(point, sc)
        ptsV = Array{Complex}(collect(values(volt_vars)))
        scatter!(real(ptsV), imag(ptsV), lab=sc)
    end
    gui()
end

function plot_Sgen_vars(pt_global, scenarios)
    fig = plot()
    pt = real2cplx(get_prod_vars(pt_global))
    Sgens = Dict()
    # Get prod by scenarios
    for scenario in scenarios
        _, _, bc_prod_vars = get_splitted_Cpt(pt_knitro, "BaseCase")
        Sgens[scenario] = bc_prod_vars
    end

    # One color per bus/generator
    init_scen = sort(collect(keys(Sgens)))[1]
    i=1
    cols = Dict()
    for (varsc, valsc) in Sgens[init_scen]
        scenario, busid, genname, varname = split(varsc.name, "_")
        cols[(busid, genname)] = i
        i+=1
    end

    # Loop on generators
    for (varsc, valsc) in Sgens[init_scen]
        base_scenario, busid, genname, varname = split(varsc.name, "_")
        println("$base_scenario, $busid, $genname, $varname")

        gen = OPFpbs[base_scenario].ds.bus[bus_name(parse(busid))][genname]
        Smin, Smax = gen.power_min, gen.power_max
        domain = [ Smin, real(Smin)+im*imag(Smax), Smax, real(Smax)+im*imag(Smin), Smin]
        plot!(real(domain), imag(domain), color=cols[(busid, genname)], label="Cstr_$(busid)_$genname")

        for scenario in setdiff(scenarios, init_scen)
            gen_cur = OPFpbs[scenario].ds.bus[bus_name(parse(busid))][genname]
            Smin_cur, Smax_cur = gen_cur.power_min, gen_cur.power_max
            (Smin == Smin_cur) && (Smax == Smax_cur) || warn("bus $busid, gen $genname: between $base_scenario and $scenario, different bounds($Smin, $Smax vs. $Smin_cur, $Smax_cur)")

            var = Variable(variable_name("Sgen", bus_name(parse(busid)), String(genname), scenario), Complex)
            scatter!([pt[var]], color=cols[(busid, genname)], label=var.name)
        end
    end
    plot!(xlabel="P gen", ylabel = "Q gen")
    return fig
end

function plot_ViVj_vars(pt, scenarios)
    r1 = 0.9^2
    plt = plot(r1*cos.(-π:0.01:π), r1*sin.(-π:0.01:π), label=:Vmin2)
    r2 = 1.1^2
    plot!(r2*cos.(-π:0.01:π), r2*sin.(-π:0.01:π), label=:Vmax2)

    cols = Dict()
    i = 1
    for scen in scenarios
        cols[scen] = i
        i+=1
    end

    for scen in scenarios
        volt_vars, _, _ = get_splitted_Cpt(pt_knitro, scen)
        V = Vector{Complex}([val for (var, val) in volt_vars])
        ViVj = collect(Set(V*V'))
        scatter!(ViVj, label=scen, color=cols[scen])
    end
    return plt
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Trains
using Test

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


text:
using Base.Test
using BioSymbols
using BioSequences
using GenomicFeatures
using SuffixArrays
using Gadfly
using Automa

using BufferedStreams
using Libz

using PyCall, PyPlot

using ScikitLearn
using ScikitLearn.Utils: meshgrid

@sk_import ensemble: IsolationForest

@pyimport matplotlib.font_manager as fm
@pyimport scipy.stats as stats

import Automa
import Automa.RegExp: @re_str
import Compat: take!

importall BioSymbols
importall GenomicFeatures

include("../src/pairs.jl")
include("../src/energy.jl")
include("../src/rnaduplex.jl")
include("../src/intervals.jl")
include("../src/traverse.jl")
include("../src/mafreader.jl")
include("../src/gtrmodel.jl")
include("../src/newick.jl")
include("../src/suffix.jl")
include("../src/evoduplex.jl")
include("../src/io.jl")
include("../src/regions.jl")
include("../src/train.jl")

const pair_set   = [AU_PAIR, UA_PAIR, CG_PAIR, GC_PAIR, GU_PAIR, UG_PAIR]
const mismat_set = [AA_MISMATCH, AG_MISMATCH, AC_MISMATCH, CA_MISMATCH, CC_MISMATCH,
                    CU_MISMATCH, GA_MISMATCH, GG_MISMATCH, UC_MISMATCH, UU_MISMATCH]
const bulge_set  = [AB_BULGE, CB_BULGE, GB_BULGE, UB_BULGE, BA_BULGE, BC_BULGE, BG_BULGE, BU_BULGE]

@testset "Nucleotide Pairs" begin

   @test convert(RNAPair, RNA_A, RNA_U) == AU_PAIR
   @test convert(RNAPair, RNA_U, RNA_A) == UA_PAIR
   @test convert(RNAPair, RNA_C, RNA_G) == CG_PAIR
   @test convert(RNAPair, RNA_G, RNA_C) == GC_PAIR
   @test convert(RNAPair, RNA_G, RNA_U) == GU_PAIR
   @test convert(RNAPair, RNA_U, RNA_G) == UG_PAIR

   @test convert(RNAMismatch, RNA_A, RNA_A) == AA_MISMATCH
   @test convert(RNAMismatch, RNA_A, RNA_G) == AG_MISMATCH
   @test convert(RNAMismatch, RNA_A, RNA_C) == AC_MISMATCH
   @test convert(RNAMismatch, RNA_C, RNA_A) == CA_MISMATCH
   @test convert(RNAMismatch, RNA_C, RNA_C) == CC_MISMATCH
   @test convert(RNAMismatch, RNA_C, RNA_U) == CU_MISMATCH
   @test convert(RNAMismatch, RNA_G, RNA_A) == GA_MISMATCH
   @test convert(RNAMismatch, RNA_G, RNA_G) == GG_MISMATCH
   @test convert(RNAMismatch, RNA_U, RNA_C) == UC_MISMATCH
   @test convert(RNAMismatch, RNA_U, RNA_U) == UU_MISMATCH

   @test convert(RNABulge, RNA_A, RNA_Gap) == AB_BULGE
   @test convert(RNABulge, RNA_C, RNA_Gap) == CB_BULGE
   @test convert(RNABulge, RNA_G, RNA_Gap) == GB_BULGE
   @test convert(RNABulge, RNA_U, RNA_Gap) == UB_BULGE
   @test convert(RNABulge, RNA_Gap, RNA_A) == BA_BULGE
   @test convert(RNABulge, RNA_Gap, RNA_C) == BC_BULGE
   @test convert(RNABulge, RNA_Gap, RNA_G) == BG_BULGE
   @test convert(RNABulge, RNA_Gap, RNA_U) == BU_BULGE

   unique = Dict{Int,Bool}()
   ind_sum    = 0
   for val in [pair_set; mismat_set]
      i = index(val)
      @test !haskey(unique, i)
      unique[i] = true
      ind_sum += i
   end
   @test ind_sum == sum(collect(1:16))

   @test split(AU_PAIR) == (1,4)
   @test split(CG_PAIR) == (2,3)
   @test split(GU_PAIR) == (3,4)

   @test AU_PAIR == flip(UA_PAIR)
   @test UA_PAIR == flip(AU_PAIR)
   @test CG_PAIR == flip(GC_PAIR)
   @test GC_PAIR == flip(CG_PAIR)
   @test UG_PAIR == flip(GU_PAIR)
   @test GU_PAIR == flip(UG_PAIR)

   @test is_purine_pyrimidine( AU_PAIR ) == true
   @test is_purine_pyrimidine( GC_PAIR ) == true
   @test is_purine_pyrimidine( GU_PAIR ) == true
   @test is_purine_pyrimidine( UA_PAIR ) == false
   @test is_purine_pyrimidine( CG_PAIR ) == false
   @test is_purine_pyrimidine( UG_PAIR ) == false

   @test AA_MISMATCH == flip(AA_MISMATCH)
   @test AG_MISMATCH == flip(GA_MISMATCH)
   @test AC_MISMATCH == flip(CA_MISMATCH)
   @test CC_MISMATCH == flip(CC_MISMATCH)
   @test CU_MISMATCH == flip(UC_MISMATCH)
   @test GG_MISMATCH == flip(GG_MISMATCH)
   @test UU_MISMATCH == flip(UU_MISMATCH)

   @test AB_BULGE == flip(BA_BULGE)
   @test CB_BULGE == flip(BC_BULGE)
   @test GB_BULGE == flip(BG_BULGE)
   @test UB_BULGE == flip(BU_BULGE)
   @test BA_BULGE == flip(AB_BULGE)
   @test BC_BULGE == flip(CB_BULGE)
   @test BG_BULGE == flip(GB_BULGE)
   @test BU_BULGE == flip(UB_BULGE)

   @test isfiveprime( AB_BULGE ) == true
   @test isfiveprime( CB_BULGE ) == true
   @test isfiveprime( GB_BULGE ) == true
   @test isfiveprime( UB_BULGE ) == true
   @test isfiveprime( BA_BULGE ) == false
   @test isfiveprime( BC_BULGE ) == false
   @test isfiveprime( BG_BULGE ) == false
   @test isfiveprime( BU_BULGE ) == false

   for b in bulge_set
      @test !ispair(b)
      @test isbulge(b)
      @test !ismismatch(b)
   end
   @test nbulges(bulge_set) == length(bulge_set)
   @test nmismatches(bulge_set) == 0
   @test npairs(bulge_set) == 0

   for m in mismat_set
      @test !ispair(m)
      @test !isbulge(m)
      @test ismismatch(m)
   end
   @test nmismatches(mismat_set) == length(mismat_set)
   @test nbulges(mismat_set) == 0
   @test npairs(mismat_set) == 0

   for p in pair_set
      @test ispair(p)
      @test !isbulge(p)
      @test !ismismatch(p)
   end
   @test npairs(pair_set) == length(pair_set)
   @test nbulges(pair_set) == 0
   @test nmismatches(pair_set) == 0

   @test split(AB_BULGE) == 1
   @test split(CB_BULGE) == 2
   @test split(GB_BULGE) == 3
   @test split(UB_BULGE) == 4
   @test split(BA_BULGE) == 1
   @test split(BC_BULGE) == 2
   @test split(BG_BULGE) == 3
   @test split(BU_BULGE) == 4

   mypath = [AU_PAIR, CG_PAIR, AB_BULGE]
   @test five_three(mypath, 1:3) == (1,4)
   @test five_three(mypath, 1:3, last=true) == (1,3)
   mypath = [BC_BULGE, CG_PAIR, BA_BULGE]
   @test five_three(mypath, 1:3) == (2,2)
   @test five_three(mypath, 1:3, last=true) == (2,1)

end

@testset "Energy Parameters" begin
   for i in pair_set, j in pair_set
      if i == GU_PAIR && j == UG_PAIR || i == UG_PAIR && j == GU_PAIR
         @test TURNER_1998_STACK[index(i), index(j)] > 0.0
      else
         @test TURNER_1998_STACK[index(i), index(j)] < 0.0
      end

      @test length(TURNER_2004_INTERNAL_TWO[ index(i), index(j) ]) == 16
      @test sum(TURNER_2004_INTERNAL_TWO[ index(i), index(j) ]) > 0.0

      for k in 1:4
         @test length(TURNER_2004_INTERNAL_THREE[ index(i), index(j), k ]) == 16
         @test sum(TURNER_2004_INTERNAL_THREE[ index(i), index(j), k ]) > 0.0
      end

      @test length(TURNER_2004_INTERNAL_FOUR[ index(i), index(j) ]) == 256
      @test sum(TURNER_2004_INTERNAL_FOUR[ index(i), index(j) ]) > 0.0
   end
end

@testset "Duplex Energy Functions" begin

   duplex = RNADuplex()

   # Test simple 1x1 lookup
   duplex.path = [AU_PAIR, AA_MISMATCH, AU_PAIR]
   @test internal_1x1_energy( duplex, 1:3 ) == 1.9
   duplex.path = [UG_PAIR, UU_MISMATCH, UG_PAIR]
   @test internal_1x1_energy( duplex, 1:3 ) == 1.6
   @test motif_energy( duplex, 3 )          == 1.6

   # Test lookup with 1x2 table.
   duplex.path = [AU_PAIR, AA_MISMATCH, BA_BULGE, AU_PAIR]
   @test internal_1x2_energy( duplex, 1:4 ) == 3.7
   # Reverse bulge/mismatch order
   duplex.path = [AU_PAIR, BA_BULGE, AA_MISMATCH, AU_PAIR]
   @test internal_1x2_energy( duplex, 1:4 ) == 3.7
   # Invert bulge side
   duplex.path = [UA_PAIR, AB_BULGE, AA_MISMATCH, UA_PAIR]
   @test internal_1x2_energy( duplex, 1:4 ) == 3.7
   # Invert and reverse
   duplex.path = [UA_PAIR, AA_MISMATCH, AB_BULGE, UA_PAIR]
   @test internal_1x2_energy( duplex, 1:4 ) == 3.7
   # Test range
   duplex.path = [UG_PAIR, UU_MISMATCH, BU_BULGE, UG_PAIR]
   @test internal_1x2_energy( duplex, 1:4 ) == 3.0
   @test motif_energy( duplex, 4 )          == 3.0

   # Test lookup with 2x2 table.
   duplex.path = [AU_PAIR, AA_MISMATCH, AA_MISMATCH, AU_PAIR]
   @test internal_2x2_energy( duplex, 1:4 ) == 2.8
   duplex.path = [UG_PAIR, UU_MISMATCH, UU_MISMATCH, UG_PAIR]
   @test internal_2x2_energy( duplex, 1:4 ) == 3.1
   @test motif_energy( duplex, 4 )          == 3.1

   # Test helix symmetry
   duplex.path = [AU_PAIR, UA_PAIR, CG_PAIR, GC_PAIR] # no_symmetry
   @test helix_symmetry( duplex ) == 0.0
   duplex.path = [AU_PAIR, CG_PAIR, GC_PAIR, UA_PAIR] # even_symmetry
   @test helix_symmetry( duplex ) == TURNER_1998_HELICAL_SYMMETRY
   duplex.path = [AU_PAIR, CG_PAIR, GU_PAIR, GC_PAIR, UA_PAIR] # odd_symmetry
   @test helix_symmetry( duplex ) == TURNER_1998_HELICAL_SYMMETRY

end

@testset "Duplex Calculations" begin

   # Duplex unit-test examples from http://rna.urmc.rochester.edu/NNDB

   #= Non-self complimentary duplex

     5'GCACG 3'
     3'CGUGC 5'  =#

    duplex = RNADuplex()
    push!( duplex, [GC_PAIR, CG_PAIR, AU_PAIR, CG_PAIR, GC_PAIR] )
    @test length( duplex.path ) == 5
    @test length( duplex.energy ) == length( duplex.path ) + 1
    @test energy( duplex ) == -6.01

   #= Self complimentary duplex

     5'AGCGCU 3'
     3'UCGCGA 5'  =#

    duplex = RNADuplex()
    push!( duplex, [AU_PAIR, GC_PAIR, CG_PAIR, GC_PAIR, CG_PAIR, UA_PAIR] )
    @test length( duplex.path ) == 6
    @test length( duplex.energy ) == length( duplex.path ) + 1
    @test energy( duplex ) == -7.98

   #= 2x2 internal loop

         GA
     5'CA  CG
     3'GU  GC
         AG  =#

    duplex = RNADuplex()
    push!( duplex, [CG_PAIR,AU_PAIR,GA_MISMATCH,AG_MISMATCH,CG_PAIR,GC_PAIR] )
    @test length( duplex.path ) == 6
    @test length( duplex.energy ) == length( duplex.path ) + 1
    @test energy( duplex ) == -1.51

   #= 1x5 internal loop

          G
     5'CA     CG
     3'GU     GC
         GAAAG   =#

     duplex = RNADuplex()
     push!( duplex, [CG_PAIR, AU_PAIR, GG_MISMATCH, BA_BULGE, BA_BULGE, BA_BULGE,
                     BG_BULGE, CG_PAIR, GC_PAIR] )
     @test length( duplex.path ) == 9
     @test length( duplex.energy ) == length( duplex.path ) + 1
     @test energy( duplex ) == 4.69

   #= 2x3 internal loop with stabilizing mismatches

         GA
     5'CA   CG
     3'GU   GC
         GAG  =#

     duplex = RNADuplex()
     push!( duplex, [CG_PAIR, AU_PAIR, GG_MISMATCH, AA_MISMATCH, BG_BULGE, CG_PAIR,
                     GC_PAIR] )
     @test length( duplex.path ) == 7
     @test length( duplex.energy ) == length( duplex.path ) + 1
     @test energy( duplex ) == 0.89

   #= Single exception bulge with multiple states

          C
     5'GCC G
     3'CGG C  =#

     duplex = RNADuplex()
     push!( duplex, [GC_PAIR, CG_PAIR, CG_PAIR, CB_BULGE, GC_PAIR] )
     @test length( duplex.path ) == 5
     @test length( duplex.energy ) == length( duplex.path ) + 1
     @test energy( duplex ) == -2.7867

   #= Multiple nucleotide bulge

         ACA
     5'GA   G
     3'CU   C  =#

     duplex = RNADuplex()
     push!( duplex, [GC_PAIR, AU_PAIR, AB_BULGE, CB_BULGE, AB_BULGE, GC_PAIR] )
     @test length( duplex.path ) == 6
     @test length( duplex.energy ) == length( duplex.path ) + 1
     @test energy( duplex ) == 5.34

end

@testset "Duplex Interval Collection" begin

   # Test pushing, overwriting, and ordering
   dup = RNADuplex()
   push!(dup, [CG_PAIR,AU_PAIR,AB_BULGE,GC_PAIR])
   c = DuplexCollection{String}()
   i = DuplexInterval(Interval("c",1,6,'+',"genea"),
                      Interval("c",50,55,'+',"genea"), dup)
   push!( c, i )
   @test length(c.names["c"]) == 1
   @test length(c.names["c"][one(UInt64)]) == 1

   dup2 = deepcopy( dup )
   push!( dup2, CG_PAIR )
   @test energy(dup2) < energy(dup)
   i2 = DuplexInterval(Interval("c",2,7,'+',"genea"),
                       Interval("c",20,25,'+',"genea"), dup2)
   push!( c, i2 )
   @test length(c.names["c"]) == 1
   @test length(c.names["c"][one(UInt64)]) == 2

   i3 = DuplexInterval(Interval("c",1,7,'+',"genea"),
                       Interval("c",49,55,'+',"genea"), dup2)
   push!( c, i3 )
   @test length(c.names["c"]) == 1
   @test length(c.names["c"][one(UInt64)]) == 2

   @test isordered(c)

   # Test addition-of-larger overwriting duplex
   c = DuplexCollection{String}()
   dup = RNADuplex()
   push!(dup, [CG_PAIR,GC_PAIR,CG_PAIR])
   i = DuplexInterval(Interval("c",1,5,'+',"genea"),
                      Interval("c",100,105,'+',"genea"), dup)
   dup2 = deepcopy(dup)
   push!(dup2, [GC_PAIR,CG_PAIR])
   @test energy(dup2) < energy(dup)
   i2 = DuplexInterval(Interval("c",6,10,'+',"genea"),
                       Interval("c",106,110,'+',"genea"), dup2)
   push!( c, i  )
   push!( c, i2 )
   @test length(c.names["c"]) == 1
   @test length(c.names["c"][one(UInt64)]) == 2

   worse = RNADuplex()
   push!(worse, [GU_PAIR,AB_BULGE,GB_BULGE,UA_PAIR,UG_PAIR,CB_BULGE,UA_PAIR])
   @test energy(worse) > energy(dup)
   i3 = DuplexInterval(Interval("c",1,10,'+',"genea"),
                       Interval("c",100,110,'+',"genea"), worse)

   push!( c, i3 ) # less favorable energy should not get pushed.
   @test length(c.names["c"]) == 1
   @test length(c.names["c"][one(UInt64)]) == 2
   firstint = c.names["c"][one(UInt64)][1]
   @test firstint.first.first == 1 && firstint.first.last == 5

   better = RNADuplex()
   push!(better, [CG_PAIR,GC_PAIR,CG_PAIR,CG_PAIR,GC_PAIR,CG_PAIR])
   @test energy(better) < energy(dup2)
   i3 = DuplexInterval(Interval("c",1,10,'+',"genea"),
                       Interval("c",100,110,'+',"genea"), better)
   push!( c, i3 ) # better energy should overwrite overlapping
   @test length(c.names["c"]) == 1
   firstint = c.names["c"][one(UInt64)][1]
   @test firstint.first.first == 1 && firstint.first.last == 10
   @test length(c.names["c"][one(UInt64)]) == 1

   # Test addition of larger overwriting duplex, and subsequent deletion because of smaller higher energy duplex
   c = DuplexCollection{String}()
   push!( c, i  )
   push!( c, i2 )
   middle = RNADuplex()
   push!(middle, [CG_PAIR,GC_PAIR,CG_PAIR, GC_PAIR])
   @test energy(dup2) < energy(middle) < energy(dup)
   i3 = DuplexInterval(Interval("c",1,10,'+',"genea"),
                       Interval("c",100,110,'+',"genea"), middle)
   push!( c, i3 )
   @test length(c.names["c"]) == 1
   firstint = c.names["c"][one(UInt64)][1]
   @test firstint.first.first == 6 && firstint.first.last == 10
   @test length(c.names["c"][one(UInt64)]) == 1

end

@testset "Duplex Intervals Stitching" begin

   # test duplex strings function

   c = DuplexCollection{String}()
   dup = RNADuplex()
   push!(dup, [CG_PAIR,GC_PAIR,CG_PAIR,CG_PAIR,GC_PAIR,GC_PAIR])
   i = DuplexInterval(Interval("c",  1,  6,'+',"genea"),
                      Interval("c",105,110,'+',"genea"), dup)
   dup2 = RNADuplex()
   push!(dup2, [GC_PAIR,GC_PAIR,UG_PAIR,AU_PAIR,CG_PAIR,GC_PAIR])
   i2 = DuplexInterval(Interval("c",  5, 10,'+',"genea"),
                       Interval("c",101,106,'+',"genea"), dup2)

   @test !isnull( stitch(i, i2, 3, 3) )
   @test !isnull( stitch(i2, i, 3, 3) )

   # test stitching with bulges
   c = DuplexCollection{String}()
   dup = RNADuplex()
   push!(dup, [CG_PAIR,GC_PAIR,CG_PAIR,CG_PAIR,GC_PAIR,AB_BULGE,GC_PAIR])
   i = DuplexInterval(Interval("c",  1,  7,'+',"genea"),
                      Interval("c",105,110,'+',"genea"), dup)
   dup2 = RNADuplex()
   push!(dup2, [GC_PAIR,AB_BULGE,GC_PAIR,UG_PAIR,AU_PAIR,CG_PAIR,GC_PAIR])
   i2 = DuplexInterval(Interval("c",  5, 11,'+',"genea"),
                       Interval("c",101,106,'+',"genea"), dup2)

   @test !isnull( stitch(i, i2, 3, 3) )
   @test !isnull( stitch(i2, i, 3, 3) )

   # test negative stitching cases

end

@testset "Load Bed & BedGraph" begin
   g = loadbedgraph( "head.bedGraph.gz" )
   @test length(g) == 999
   i = loadbed( "small.bed" )
   g = loadbedgraph( "head.bedGraph.gz", regions=i, regionbool=true )
   @test length(g) == 10
end

@testset "MAF Parser" begin
   mafheader = """
##maf version=1 scoring=tba.v8
 # tba.v8 (((human chimp) baboon) (mouse rat))
 # multiz.v7
 # maf_project.v5 _tba_right.maf3 mouse _tba_C
 # single_cov2.v4 single_cov2 /dev/stdin
"""
   maflines = """
a score=23262.0
s hg16.chr7    27578828 38 + 158545518 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s panTro1.chr6 28741140 38 + 161576975 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s baboon.chr1    116834 38 +   4622798 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
s mm4.chr6     53215344 38 + 151104725 -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
s rn3.chr4     81344243 40 + 187371129 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG

a score=5062.0
s hg16.chr7    27699739 6 + 158545518 TAAAGA
s panTro1.chr6 28862317 6 + 161576975 TAAAGA
s baboon.chr1    241163 6 +   4622798 TAAAGA
s mm4.chr6     53303881 6 + 151104725 TAAAGA
s rn3.chr4     81444246 6 + 187371129 taagga

a score=6636.0
s hg16.chr7    27707221 13 + 158545518 gcagctgaaaaca
s panTro1.chr6 28869787 13 + 161576975 gcagctgaaaaca
s baboon.chr1    249182 13 +   4622798 gcagctgaaaaca
s mm4.chr6     53310102 13 + 151104725 ACAGCTGAAAATA
"""

   maffile = mafheader * "\n" * maflines

   # Single buffer:

   headreader = MAFReader(BufferedInputStream(IOBuffer(maffile)))
   noheadreader = MAFReader(BufferedInputStream(IOBuffer(maflines)))

   discard_header!(headreader)
   @test headreader.stream.buffer[1:10] == noheadreader.stream.buffer[1:10]

   headreader = MAFReader(BufferedInputStream(IOBuffer(maffile)))

   while !done(headreader)
      headrec   = MAFRecord()
      noheadrec = MAFRecord()
      read!(headreader, headrec)
      read!(noheadreader, noheadrec)
      @test length(noheadrec) == length(headrec)
      @test length(noheadrec) > 1
      for i in 1:length(headrec)
         @test headrec[i] == noheadrec[i]
      end
   end
   @test done(headreader)
   @test done(noheadreader)

   # Multi buffer

   longreader = MAFReader(BufferedInputStream(IOBuffer(maffile)))
   shortreader = MAFReader(BufferedInputStream(IOBuffer(maffile), 512))

   while !done(longreader)
      longrec  = MAFRecord()
      shortrec = MAFRecord()
      read!(longreader, longrec)
      read!(shortreader, shortrec)
      @test length(shortrec) == length(longrec)
      @test length(shortrec) > 1
      for i in 1:length(shortrec)
         @test shortrec[i] == longrec[i]
      end
   end
   @test done(shortreader)
   @test done(longreader)

   # Bad Records
   maffile = """
##maf version=1 scoring=roast.v3.3
a score=781.000000
s hg19.chr22                16092278 1 +  51304566 C
s panTro4.chr22             14470458 1 +  49737984 C
s nomLeu3.chrUn_GL397501_2   1216368 1 -   1485088 C
s chlSab1.chrUn_KE147573       66869 1 -     85383 C
s monDom5.chr4              96780965 1 - 435153693 T
s danRer7.chr3                323838 1 -  63268876 C
s petMar2.GL501154             43910 1 -     71231 t

a score=-50707.000000
s hg19.chr22                     16092279 32 +  51304566 AGTCCCCAGGTGGT------CATGACACCTCAATTGGA
s panTro4.chr22                  14470459 32 +  49737984 AGTCCCCAGGTGGT------CATGACACCTCAATTGGA
s nomLeu3.chrUn_GL397501_2        1216369 32 -   1485088 AGTCCCCAGGTGGT------CATGACACTTCAACTGGA
s chlSab1.chrUn_KE147573            66870 32 -     85383 AGTCCCCAGGTGAT------CATGACACCTTAACTGGA
s monDom5.chr4                   96780966 24 - 435153693 TGACTCCAGCCAGA------GACGGAGCCC--------
s falChe1.KB397375                  11884 29 +     36355 GACCCTTGTGAAAA------CATGAGGGCTGTGAA---
s falPer1.KB390849                  40484 13 +     49439 AGCCCTTGTCTGA-------------------------
s ficAlb2.chrUn_KE165385           131290 31 -    141484 -cctccccgagtat------cCTGCGACCCCAGATCGA
s geoFor1.JH741012                   5677 32 -     10893 GGGCCCTTTGGGAA------GAGGACACTGCGACTGGA
s taeGut2.chr26_EQ833276_random      8281 32 -     85973 AGCCCCTCAGCTGT------CCCCACACCTCTCCTGCA
s pseHum1.KB221494                  45158 24 +     77128 gaccccacaatgac-------------cctcggttga-
s galGal4.chrUn_JH375878              275 32 -      4229 AGCCCAGTggggga------ccccggagccaaagggga
s pelSin1.JH208549                  94224 38 +    588315 GGCCCCAGGGCCAAATCCGTCCCTGGGTCCATATTCAA
s danRer7.chr3                     323839 32 -  63268876 ATTATCCAACTAGA------CTTGGCTGTTTAGAAGAC
s petMar2.GL501154                  43911 30 -     71231 --agttcacgtgag------ctctacaactcaggtgaa
"""

   noshort = MAFReader(BufferedInputStream(IOBuffer(maffile)))
   noshort.minlen = 5
   noshortrec = MAFRecord()
   read!(noshort, noshortrec)
   @test done(noshort)
   @test length(noshortrec) == 15

   noweak = MAFReader(BufferedInputStream(IOBuffer(maffile)))
   noweak.minscore = 0.0
   noweakrec = MAFRecord()
   read!(noweak, noweakrec)
   @test !done(noweak)
   @test length(noweakrec) == 7
   res = read!(noweak, MAFRecord())
   @test done(noweak)
   @test length(res) == 0

   #Trailing newline
   withnew = MAFReader(BufferedInputStream(IOBuffer(maffile * "\n")))
   withnew.minlen = 5
   withnewrec = MAFRecord()
   read!(withnew, withnewrec)
   @test done(withnew)
   @test length(withnewrec) == 15

   #Test MAFRecord gap deletion
   for i in 1:length(withnewrec)
      @test length(withnewrec[i].sequence) == length("AGTCCCCAGGTGGT------CATGACACCTCAATTGGA")
   end
   deletegaps!(withnewrec)
   for i in 1:length(withnewrec)
      @test length(withnewrec[i].sequence) == length("AGTCCCCAGGTGGTCATGACACCTCAATTGGA")
   end
end

@testset "Parsing MAF File & Stitching" begin
   # Stitching adjacent MAF records
   maflines = """
   a score=23262.0
   s hg16.chr7    27571001 38 + 158545518 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
   s panTro1.chr6 28741140 38 + 161576975 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
   s baboon.chr1    116834 38 +   4622798 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
   s rn3.chr4     81344243 40 + 187371129 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG

   a score=5062.0
   s hg16.chr7    27571039 6 + 158545518 TAAAGA
   s panTro1.chr6 28862317 6 + 161576975 TAAAGA
   s baboon.chr1    241163 6 +   4622798 TAAAGA
   s rn3.chr4     81444246 6 + 187371129 taagga

   a score=6636.0
   s hg16.chr7    27571045 13 + 158545518 gcagctgaaaaca
   s panTro1.chr6 28869787 13 + 161576975 gcagctgaaaaca
   s baboon.chr1    249182 13 +   4622798 gcagctgaaaaca
   s mm4.chr6     53310102 13 + 151104725 ACAGCTGAAAATA
   """

   tree = parsenewick("(((((hg16:0.006591,panTro1:0.006639):0.002184,baboon:0.009411):0.009942,
                            mm4:0.018342):0.014256,rn3:0.036199):0.021496,unused:0.04);")

   reader = MAFReader(BufferedInputStream(IOBuffer(maflines)))

   col = readmaf!( reader, tree.index )
   for i in col
      for j in i.metadata.species
         @test length(j.sequence) == length(i.metadata.species[1].sequence)
      end
   end

   println(col)

   #=MAFRecord(MAFSpecies[
   MAFSpecies("hg16","chr7",27571001,true,   AAAGGGAATGTTAACCAAATGAATTGTCTCTTACGGTGTAAAGAGCAGCTGAAAACA),
   MAFSpecies("panTro1","chr6",28741140,true,AAAGGGAATGTTAACCAAATGAATTGTCTCTTACGGTGTAAAGAGCAGCTGAAAACA),
   MAFSpecies("baboon","chr1",116834,true,   AAAGGGAATGTTAACCAAATGAGTTGTCTCTTATGGTGTAAAGAGCAGCTGAAAACA),
   MAFSpecies("mm4","chr6",53215344,true,    -AAGGGAATGTTAAGCAAACGAATTGTCTCTCAGTGTGTAAAGAACAGCTGAAAATA),
   MAFSpecies("rn3","chr4",81344243,true,    -AAGGGGATGCTAAGCCAATGAGTTGTCTCTCAATGTGTAAGGA-------------)
   =#
end

@testset "Newick tree building and Felsenstein's algorithm" begin
   tree = parsenewick("(((((hg19:0.006591,panTro2:0.006639):0.002184,gorGor1:0.009411):0.009942,
                            ponAbe2:0.018342):0.014256,rheMac2:0.036199):0.021496,papHam1:0.04);")
   names = collect(tree.root)
   expnames = ["hg19", "panTro2", "gorGor1", "ponAbe2", "rheMac2", "papHam1"]
   @test length(names) == length(expnames) == length(tree.order)
   for i in 1:length(names)
      @test names[i] == expnames[i] == tree.order[i]
   end
   # Test rate matrix -> probability
   Q = [-3.0 1.0 1.0 1.0
        1.0 -3.0 1.0 1.0
        1.0 1.0 -3.0 1.0
        1.0 1.0 1.0 -3.0]
   set_prob_mat!( tree, Q )
   hg19_node = tree.root.left.value.left.value.left.value.left.value.left.value
   @test hg19_node.label == "hg19"
   @test hg19_node.prob == expm(Q*0.006591)

   # Test likelihood function
   tree   = parsenewick("((hg19:0.006591,panTro2:0.006639):0.002184,gorGor1:0.5);")
   set_prob_mat!( tree, GTR_SINGLE_Q )

   #=
   hg19 length=0.006591
   A->G 0.0125902
   C->G 0.00378768
   G->G 0.985462
   T->G 0.00314851
   =#

   hg19 = expm(GTR_SINGLE_Q*0.006591) * [0,0,1,0]
   @test length(hg19) == 4

   #=
   panTro length=0.006639
   A->G 0.0126803
   C->G 0.00381511
   G->G 0.985357
   T->G 0.00317139
   =#

   panTro = expm(GTR_SINGLE_Q*0.006639) * [0,0,1,0]
   @test length(panTro) == 4

   #=
   hgpanAnc length=0.002184
   A->A 0.000158486
   A->C 2.07831e-7
   A->G 4.64153e-7
   A->T 1.8804e-7

   C->A 2.28759e-8
   C->C 1.43597e-5
   C->G 1.51679e-8
   C->T 5.20155e-8

   G->A 0.00409759
   G->C 0.00122328
   G->G 0.966317
   G->T 0.0010147

   T->A 1.47604e-8
   T->C 3.7086e-8
   T->G 8.97429e-9
   T->T 9.92702e-6

   *->A 0.0042561136363
   *->C 0.001237884617
   *->G 0.9663174882951899
   *->T 0.0010248670755000002
   =#

   hgpanAnc = expm(0.002184*GTR_SINGLE_Q) * (hg19 .* panTro)
   @test length(hgpanAnc) == 4

   @test round(0.0042561136363, 6) == round(hgpanAnc[1], 6)
   @test round(0.001237884617, 6) == round(hgpanAnc[2], 6)
   @test round(0.9663174882951899, 6) == round(hgpanAnc[3], 6)
   @test round(0.0010248670755000002, 6) == round(hgpanAnc[4], 6)

   @test round(likelihood( tree, DNA[DNA_G, DNA_G, DNA_A] ),4) == 0.0539

   smtree   = parsenewick("((hg19:0.006591,panTro2:0.006639):0.002184,gorGor1:0.15);")
   pairtree = deepcopy(smtree)
   set_prob_mat!( smtree,   GTR_SINGLE_Q )
   set_prob_mat!( pairtree, GTR_PAIRED_Q )

   single_p = likelihood( smtree, DNA[DNA_G, DNA_G, DNA_G] ) * likelihood( smtree, DNA[DNA_C, DNA_C, DNA_G] )
   paired_p = likelihood( pairtree, DNA[DNA_G, DNA_G, DNA_G], DNA[DNA_C, DNA_C, DNA_G] )
   @test single_p > paired_p

   single_p = likelihood( smtree, DNA[DNA_G, DNA_G, DNA_G] ) * likelihood( smtree, DNA[DNA_C, DNA_C, DNA_T] )
   paired_p = likelihood( pairtree, DNA[DNA_G, DNA_G, DNA_G], DNA[DNA_C, DNA_C, DNA_T] )
   @test paired_p > single_p

end

@testset "RNADuplexArray Building and Traversal" begin

   # test depth and sai composition
   seq = dna"AAATGATGCCGCAGGGGGGGGGGTGCGGCAATCATTT"
   rsuf = RNASuffixArray{DNAAlphabet{2},UInt8,UInt8}( seq, 10 )
   for i in 1:length(rsuf.sai)-1
      @test <=( rsuf.depth, i, rsuf.depth, i+1 )
   end

   push!( rsuf, reverse(seq) )
   for i in 1:length(rsuf.sai)-1
      @test <=( rsuf.depth, i, rsuf.depth, i+1 )
   end

   # test basic duplex array building and traversal
   rda = RNADuplexArray{DNAAlphabet{2},UInt8,UInt8}( seq, 25 )
   @test length(collect(traverse( rda, 10:50, bulge_max=0 ), minenergy=-15.0)) == 0
   @test length(collect(traverse( rda, 10:50, bulge_max=1 ), minenergy=-15.0)) == 1
   @test length(traverse( rda, 10:50, bulge_max=1 )) == 1
   res = shift!(collect(traverse( rda, 10:50, bulge_max=1 )))
   @test first(res.first) == 1 && last(res.first) == 13
   @test first(res.last) == 24 && last(res.last) == length(seq)

   # test stitching of smaller duplex collection
   rda = RNADuplexArray{DNAAlphabet{2},UInt8,UInt8}( seq, 10 )
   res = collect(traverse(rda, bulge_max=1, minfold=-6.0))
   @test length(res) == 5
   i = 1
   for j in res
      @test j.first.first == i
      i += 1
   end
   sti = stitch(traverse(rda, bulge_max=1, minfold=-6.0))
   @test length(collect(sti)) == 1

   # TODO
   # test fwd and rev sequences for intermolecular constructor
   rda  = RNADuplexArray{DNAAlphabet{2},UInt8,UInt8}( seq, seq, 25 )
   res  = collect(traverse( rda, bulge_max=1 ), minenergy=-15.0)

   # test push! additional sequences
   prda = RNADuplexArray{DNAAlphabet{2},UInt8,UInt8}( seq, 25 )
   push!( prda, seq )
   pres = collect(traverse( prda, bulge_max=1 ), minenergy=-15.0)

   # test vector constructor

   # test larger scale example hnRNP D and genome offsets

   # test recursive stitching
   # test iterval collection and filtering

end

@testset "RNADuplexArray EvoDuplexes: building, traversal, stitching, scoring" begin
      maffile = """
##maf version=1 scoring=roast.v3.3
a score=-50707.000000
s hg19.chr22                     16092279 32 +  51304566 AAATGATGCCGCAGGGG------GGGGGGTGCGGCAATCATTT
s panTro4.chr22                  14470459 32 +  49737984 AAATGATGCCGCAGGGG------GGGGGGTGCGGCAATCATTT
"""
   mafread = MAFReader(BufferedInputStream(IOBuffer(maffile * "\n")))
   mafrec = MAFRecord()
   read!(mafread, mafrec)
   @test done(mafread)
   @test length(mafrec) == 2

   deletegaps!(mafrec)
   tree = parsenewick("(hg19:0.006591,panTro4:0.006639);")

   # correct match of tree species and maf file
   phylosuf = RNASuffixArray{DNAAlphabet{4}, UInt16, UInt8}( mafrec, tree, 10 )
   @test phylosuf.depth == phylosuf.species[1]
   @test length(phylosuf.species) == 1

   # missing species end of tree from maf
   tree = parsenewick("((hg19:0.006591,panTro4:0.006639):0.002184,gorGor1:0.009411);")
   phylosuf = RNASuffixArray{DNAAlphabet{4}, UInt16, UInt8}( mafrec, tree, 10 )
   @test phylosuf.depth == phylosuf.species[1]
   @test phylosuf.depth != phylosuf.species[2]
   @test length(phylosuf.species) == 2

   # missing species middle of tree from maf
   mafrec[2].name = "gorGor1"
   tree = parsenewick("((hg19:0.006591,panTro4:0.006639):0.002184,gorGor1:0.009411);")
   phylosuf = RNASuffixArray{DNAAlphabet{4}, UInt16, UInt8}( mafrec, tree, 10 )
   @test phylosuf.depth != phylosuf.species[1]
   @test phylosuf.depth == phylosuf.species[2]
   @test length(phylosuf.species) == 2

   # compare to single seq
   seq = dna"AAATGATGCCGCAGGGGGGGGGGTGCGGCAATCATTT"
   rda = RNADuplexArray{DNAAlphabet{2},UInt8,UInt8}( seq, 10 )
   fwd_dep = rda.fwd.depth
   @test fwd_dep == phylosuf.depth

   # test push! MAF record
   push!( phylosuf, mafrec, tree )
   @test length(phylosuf.depth[1]) == length(phylosuf.species[2][1])
   for i in 1:length(phylosuf.sai)-1
      @test <=( phylosuf.depth, i, phylosuf.depth, i+1 )
      @test <=( phylosuf.species[2], i, phylosuf.species[2], i+1 )
   end

   # test build RNADuplexArray from MAF file...
   smtree   = parsenewick("((hg19:0.006591,panTro4:0.006639):0.002184,gorGor1:0.15);")
   pairtree = deepcopy(smtree)
   set_prob_mat!( smtree, GTR_SINGLE_Q )
   set_prob_mat!( pairtree, GTR_PAIRED_Q )

   # test stitching of EvoDuplexes with toy case
   rda = RNADuplexArray{DNAAlphabet{2},UInt8,UInt8}( mafrec, smtree, 10 )
   res = collect(traverse( rda, bulge_max=1, single=smtree, minfold=-6.0 ))
   for i in 1:length(res)-1
      sti = stitch( res[i], res[i+1], 3, 3 )
      println(sti)
      @test !isnull(sti)
      left,right = strings(sti.value.duplex)
      left_nts   = collect(DNASequence(RNASequence(left)))
      right_nts  = collect(DNASequence(RNASequence(right)))
      @test sti.value.duplex.alignment[1,1:length(left_nts)] == left_nts
      @test sti.value.duplex.alignment[1,length(left_nts)+1:end] == reverse(right_nts)
      @test length(sti.value.duplex.bracket) == length(left_nts) + length(right_nts)
      @test sti.value.duplex.first == length(left_nts)
   end

   # try stitching of EvoDuplexes with complex real-world test case
   maffile = """
##maf version=1 scoring=roast.v3.3
a score=-50707.000000
s hg19.chr22                     16092279 32 +  51304566 ACTCTATCTTAAGCTTTTCTGCTTTTTAATTATCCTGAAGTAAAGATCTTTGCTGATCTTCTGACTTTAGTGAACCTATTAATGTGCTGCAGGCCCCAGTCAAAACTGGAACCAGGGATATAGTAACTATTGGAATCAAG--------------AGACCAGTTCTTGGAGTTATATCCTTTCTTAGGTGACTAGGCCTGCTGCACAATAATAGGTTAATTAAAGTCAGAAGAAGGTCAGCAAAGATGGATTGGGTGAGATTGGGGCCCTTTTCTTAGAAGGGCAGAGATACTAAGCACTG
s panTro4.chr22                  14470459 32 +  49737984 ACTCTATCTTAAGCTTTTCTGCTTTTTAATTATCCTGAAGTAAAGATCTTTGCTGATCTTCTGACTTTAGTGAACCTATTAATGTGCTGCAGGCCCCAGTCAAAACTGGAACCAGGGATATAGTAACTATTGGAATCAAG--------------AGACCAGTTCTTGGAGTTATATCCTTTCTTAGGTGACTAGGCCTGCTGCACAATAATAGGTTAATTAAAGTCAGAAGAAGGTCAGCAAAGATGGATTGGGTGAGATTGGGGCCCTTTTCTTAGAAGGGCAGAGATACTAAGCACTG
"""
   mafread = MAFReader(BufferedInputStream(IOBuffer(maffile * "\n")))
   mafrec = MAFRecord()
   read!(mafread, mafrec)
   @test done(mafread)
   @test length(mafrec) == 2

   rda = RNADuplexArray{DNAAlphabet{2},UInt16,UInt16}( mafrec, smtree, 50 )
   res = collect(traverse( rda, bulge_max=3, mismatch_max=3, single=smtree, minfold=-6.0 ))
   for i in 1:length(res)-1
      sti = stitch( res[i], res[i+1], 3, 3 )
      isnull(sti) && continue
      println(res[i])
      println(res[i+1])
      println(sti)
      left,right = strings(sti.value.duplex)
      left_nts   = collect(DNASequence(RNASequence(left)))
      right_nts  = collect(DNASequence(RNASequence(right)))
      @test sti.value.duplex.alignment[1,1:length(left_nts)] == left_nts
      @test sti.value.duplex.alignment[1,length(left_nts)+1:end] == reverse(right_nts)
      @test length(sti.value.duplex.bracket) == length(left_nts) + length(right_nts)
      @test sti.value.duplex.first == length(left_nts)
   end

   # score using EvoFold phylo-likelihood model
   # str,unstr = score!(res[1].duplex, smtree, pairtree)
   # @test str - unstr > 0
end

@testset "Distance-specific IsolationForests: train, predict" begin

   data = open(readdlm, "hnrnpd.jlt")
   df = DistanceForest()

   train!( df, data )
   p = predict( df, data )

   @test length(p)/size(data, 1) < 0.02
   @test data[p[1], 1] == 859

end


text:
#helper functions
@inline function termination_X(k, set, t, X_store_d, blocks, diff_blocks, block_options, termination)
    terminate, skip, reach_set_intersected = termination(k, set, t)
    X_store = getX_store(reach_set_intersected, X_store_d, block_options, blocks, diff_blocks)
    return terminate, skip, X_store
end

@inline function getX_store(X_store, X_store_d, block_options, blocks, diff_blocks)
    rs_oa = Approximations.overapproximate(X_store, CartesianProductArray, block_options)
    return combine_cpas(rs_oa, X_store_d, blocks, diff_blocks)
end

@inline function deco_post_sparse(b, blocks, Whatk_blocks, partition,
                                  ϕpowerk, Xhatk, Xhat0, output_function, overapproximate)
    for i in 1:b
        bi = partition[blocks[i]]
        Xhatk_bi = ZeroSet(length(bi))
        for (j, bj) in enumerate(partition)
            block = ϕpowerk[bi, bj]
            if !iszero(block)
                Xhatk_bi = Xhatk_bi + block * Xhat0[j]
            end
        end
        Xhatk_bi_lazy = (U == nothing ? Xhatk_bi : Xhatk_bi + Whatk[i])
        Xhatk[i] = (output_function == nothing) ?
            overapproximate(blocks[i], Xhatk_bi_lazy) :
            Xhatk_bi_lazy
    end

    array = CartesianProductArray(copy(Xhatk))

    X_store = (output_function == nothing) ?
                  array_d :
                  box_approximation(output_function(array))

    return X_store
end

function deco_post_dense(b, blocks, Whatk_blocks, partition, ϕpowerk, arr::Vector{LazySet{NUM}},
                                  arr_length, U, Xhat0, Xhatk, output_function, overapproximate) where {NUM}
    for i in 1:b
        bi = partition[blocks[i]]
        for (j, bj) in enumerate(partition)
            arr[j] = ϕpowerk[bi, bj] * Xhat0[j]
        end
        if U != nothing
            arr[arr_length] = Whatk_blocks[i]
        end

        Xhatk[i] = (output_function == nothing) ?
                 overapproximate(blocks[i], MinkowskiSumArray(arr)) :
                 MinkowskiSumArray(copy(arr))

    end
    array = CartesianProductArray(copy(Xhatk))

    X_store = (output_function == nothing) ?
                array :
                box_approximation(output_function(array))

    return X_store
end

# sparse
function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
                       Xhat0::Vector{<:LazySet{NUM}},
                       U::Union{ConstantInput, Nothing},
                       overapproximate::Function,
                       overapproximate_inputs::Function,
                       n::Int,
                       N::Union{Int, Nothing},
                       output_function::Union{Function, Nothing},
                       blocks::AbstractVector{Int},
                       partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
                       dimensions::Vector{Int},
                       δ::NUM,
                       guards_proj::Vector{<:LazySet{NUM}},
                       block_options,
                       vars::Vector{Int},
                       termination::Function,
                       progress_meter::Union{Progress, Nothing},
                       res::Vector{<:SparseReachSet}
                       )::Tuple{Int, Bool} where {NUM}
    ProgressMeter.update!(progress_meter, 1)
    array = CartesianProductArray(Xhat0[blocks])
    X_store = (output_function == nothing) ?
        array :
        box_approximation(output_function(array))
    t0 = zero(δ)
    t1 = δ

    diff_blocks = setdiff(1:length(partition), blocks)
    all_dimensions = 1:n

    b = length(blocks)
    bd = length(diff_blocks)

    array_d = CartesianProductArray(Xhat0[diff_blocks])
    X_store_d = (output_function == nothing) ?
            array_d :
            box_approximation(output_function(array_d))

    terminate, skip, X_store = termination_X(1, X_store, t0, X_store_d, blocks,
                                             diff_blocks, block_options, termination)
    store!(res, 1, X_store, t0, t1, all_dimensions, N)
    if terminate
        return 1, skip
    end

    b = length(blocks)
    Xhatk = Vector{LazySet{NUM}}(undef, b)
    Xhatk_d = Vector{LazySet{NUM}}(undef, bd)
    ϕpowerk = copy(ϕ)

    if U != nothing
        Whatk_blocks = Vector{LazySet{NUM}}(undef, b)
        Whatk_diff_blocks = Vector{LazySet{NUM}}(undef, bd)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk_blocks[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
        @inbounds for i in 1:bd
            bi = partition[diff_blocks[i]]
            Whatk_diff_blocks[i] = overapproximate_inputs(1, diff_blocks[i], proj(bi, n) * inputs)
        end
    end

    k = 2
    @inbounds while true
        ProgressMeter.update!(progress_meter, k)
        X_store = deco_post_sparse(b, blocks, Whatk_blocks, partition, ϕpowerk,
                                   Xhatk, Xhat0, output_function, overapproximate)

        t0 = t1
        t1 += δ

        terminate, skip, reach_set_intersected = termination(k, X_store, t0)
        if reach_set_intersected isa EmptySet
            store!(res, k, X_store, t0, t1, dimensions, N)
            break
        end
        if isdisjoint(reach_set_intersected, UnionSetArray(guards_proj))
            # keep sparse reach set
            store!(res, k, X_store, t0, t1, dimensions, N)
        else
            # compute full-dimensional reach set
            X_store_d = deco_post_sparse(bd, diff_blocks, Whatk_diff_blocks, partition,
                                         ϕpowerk, Xhatk_d, Xhat0, output_function, overapproximate)
            X_store = getX_store(reach_set_intersected, X_store_d, block_options, blocks, diff_blocks)
            store!(res, k, X_store, t0, t1, all_dimensions, N)
        end

        if terminate
            break
        end

        if U != nothing
            for i in 1:b
                bi = partition[blocks[i]]
                Whatk_blocks[i] = overapproximate_inputs(k, blocks[i],
                    Whatk_blocks[i] + row(ϕpowerk, bi) * inputs)
            end
            for i in 1:bd
                bi = partition[diff_blocks[i]]
                Whatk_diff_blocks[i] = overapproximate_inputs(k, diff_blocks[i],
                    Whatk_diff_blocks[i] + row(ϕpowerk, bi) * inputs)
            end
        end

        ϕpowerk = ϕpowerk * ϕ
        k += 1
    end

    return k, skip
end

# dense
function reach_blocks!(ϕ::AbstractMatrix{NUM},
                       Xhat0::Vector{<:LazySet{NUM}},
                       U::Union{ConstantInput, Nothing},
                       overapproximate::Function,
                       overapproximate_inputs::Function,
                       n::Int,
                       N::Union{Int, Nothing},
                       output_function::Union{Function, Nothing},
                       blocks::AbstractVector{Int},
                       partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
                       dimensions::Vector{Int},
                       δ::NUM,
                       guards_proj::Vector{<:LazySet{NUM}},
                       block_options,
                       vars::Vector{Int},
                       termination::Function,
                       progress_meter::Union{Progress, Nothing},
                       res::Vector{<:SparseReachSet}
                       )::Tuple{Int, Bool} where {NUM}
    ProgressMeter.update!(progress_meter, 1)
    array = CartesianProductArray(Xhat0[blocks])
    X_store = (output_function == nothing) ?
        array :
        box_approximation(output_function(array))
    t0 = zero(δ)
    t1 = δ

    diff_blocks = setdiff(1:length(partition), blocks)
    all_dimensions = 1:n

    b = length(blocks)
    bd = length(diff_blocks)

    array_d = CartesianProductArray(Xhat0[diff_blocks])
    X_store_d = (output_function == nothing) ?
        array_d :
        box_approximation(output_function(array_d))

    terminate, skip, X_store = termination_X(1, X_store, t0, X_store_d, blocks,
                                             diff_blocks, block_options, termination)
    store!(res, 1, X_store, t0, t1, all_dimensions, N)
    if terminate
        return 1, skip
    end

    Xhatk = Vector{LazySet{NUM}}(undef, b)
    Xhatk_d = Vector{LazySet{NUM}}(undef, bd)
    ϕpowerk = copy(ϕ)
    ϕpowerk_cache = similar(ϕ)

    if U != nothing
        Whatk_blocks = Vector{LazySet{NUM}}(undef, b)
        Whatk_diff_blocks = Vector{LazySet{NUM}}(undef, bd)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk_blocks[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
        @inbounds for i in 1:bd
            bi = partition[diff_blocks[i]]
            Whatk_diff_blocks[i] = overapproximate_inputs(1, diff_blocks[i], proj(bi, n) * inputs)
        end
    end

    arr_length = (U == nothing) ? length(partition) : length(partition) + 1
    arr = Vector{LazySet{NUM}}(undef, arr_length)
    k = 2
    @inbounds while true
        ProgressMeter.update!(progress_meter, k)
        X_store = deco_post_dense(b, blocks, Whatk_blocks, partition, ϕpowerk,
                                  arr, arr_length, U, Xhat0, Xhatk, output_function, overapproximate)
        t0 = t1
        t1 += δ
        terminate, skip, reach_set_intersected = termination(k, X_store, t0)
        if reach_set_intersected isa EmptySet
            store!(res, k, X_store, t0, t1, dimensions, N)
            break
        end

        if isdisjoint(reach_set_intersected, UnionSetArray(guards_proj))
            # keep sparse reach set
            store!(res, k, X_store, t0, t1, dimensions, N)
        else
            # compute full-dimensional reach set
            X_store_d = deco_post_dense(bd, diff_blocks, Whatk_diff_blocks, partition, ϕpowerk, arr,
                                        arr_length, U, Xhat0, Xhatk_d, output_function, overapproximate)

            X_store = getX_store(reach_set_intersected, X_store_d, block_options, blocks, diff_blocks)
            store!(res, k, X_store, t0, t1, all_dimensions, N)
        end

        if terminate
            break
        end

        if U != nothing
            for i in 1:b
                bi = partition[blocks[i]]
                Whatk_blocks[i] = overapproximate_inputs(k, blocks[i],
                    Whatk_blocks[i] + row(ϕpowerk, bi) * inputs)
            end
            for i in 1:bd
                bi = partition[diff_blocks[i]]
                Whatk_diff_blocks[i] = overapproximate_inputs(k, diff_blocks[i],
                    Whatk_diff_blocks[i] + row(ϕpowerk, bi) * inputs)
            end
        end

        mul!(ϕpowerk_cache, ϕpowerk, ϕ)
        copyto!(ϕpowerk, ϕpowerk_cache)

        k += 1
    end

    return k, skip
end


text:
# Test implementation for builtin types

@testset "Array" begin
    tree = Any[1,Any[2,3]]

    T = Vector{Any}  # This is printed as "Array{Any,1}" in older versions of Julia
    @test repr_tree(tree) == "$T\n├─ 1\n└─ $T\n   ├─ 2\n   └─ 3\n"

    @test collect(Leaves(tree)) == [1,2,3]
    @test collect(Leaves(tree)) isa Vector{Int}
    @test collect(PostOrderDFS(tree)) == Any[1,2,3,Any[2,3],Any[1,Any[2,3]]]
    @test collect(StatelessBFS(tree)) == Any[Any[1,Any[2,3]],1,Any[2,3],2,3]

    tree2 = Any[Any[1,2],Any[3,4]]
    @test collect(PreOrderDFS(tree2)) == Any[tree2,Any[1,2],1,2,Any[3,4],3,4]
end