text:
# This is an experiment in adding indexing to vectors. The idea is
# that if a DataFrame column is indexed, the following will have fast
# lookups:
#     df[:(2 .< col1 .< 7), :]
#     df[:( datecol .>= "2011-01-01" ), "col3"]
#     df[:( col .== "red" ), :]
#
# Keeping indexing with columns has advantages:
#     - Column ordering is less likely to be messed up inadvertantly.
#     - Multiple columns can be indexed; there's no main key grouping.
#     - More flexible logic combinations are possible.
#     - You can mix and match keyed and non-keyed comparisons.
#     - The DataFrame structure and indexing are not disturbed.
#     - It's probably less coding and maintenance.
# The disadvantages are:
#     - It's yet another vector type.
#     - It's somewhat slower than having one overall set of DataFrame
#       keys in a sorted DataFrame. It should still be pretty fast (no
#       speed checks, yet).
#     - You can't do data.table/pandas shortcuts like df["A"] for
#       df[:( keycol .== "A" ), :]. (But, df["A"] is less descriptive if
#       you don't know what df's keys are.)

# An IndexedVector is a pointer to the original column along with an
# indexing vector equal to `order(orig)`. A comparison operation like
# `idv .> 3` returns an Indexer type. The Indexer type includes a
# pointer to the IndexedVector along with a vector of Range1's.
# DataVector's and DataFrame's can be indexed with Indexers. It's fast
# because you're using a slice of the already indexed vector.

# Indexer's can be combined with `|` and `&`. In the case where the
# IndexedVector is the same, the Indexer is reduced or expanded as
# appropriate. This includes: `0 .< idv .< 3` or `idv .== 1 | idv .==
# 4`. If Indexers have different IndexedVectors (like `idv1 .== 1 |
# idv2 .== 1`), then the result is converted to a BitVector.

# We handle NA's by excluding them from the index. `order` puts NA's
# at the front. We use the following function to exclude those. This
# does make the indexing a little trickier as the length of the index
# can be less than the length of the DataArray.

# Note that the following leaves out NA's.
indexorder(x) = sortperm(x)
function indexorder(v::AbstractDataVector)
    Nna = sum(isna(v))
    sortperm(v)[Nna + 1 : end]
end

type IndexedVector{T,S<:AbstractVector} <: AbstractVector{T}
    x::S
    idx::Vector{Int}
end
IndexedVector{T}(x::AbstractVector{T}) = IndexedVector{T,typeof(x)}(x, indexorder(x))

Base.getindex{I<:Real}(v::IndexedVector,i::AbstractVector{I}) = IndexedVector(v.x[i])
Base.getindex{I<:Real}(v::IndexedVector,i::I) = v.x[i]
Base.getindex(v::IndexedVector,i::Real) = v.x[i]
Base.setindex!(v::IndexedVector, val::Any, i::Real) = IndexedVector(setindex!(v.x, val, i))
Base.setindex!(v::IndexedVector, val::Any, inds::AbstractVector) = IndexedVector(setindex!(v.x, val, inds))
Base.reverse(v::IndexedVector) = IndexedVector(reverse(v.x), reverse(v.idx))
Base.similar(v::IndexedVector, T, dims::Dims) = similar(v.x, T, dims)

vecbind_type(x::IndexedVector) = vecbind_type(x.x)

# to make assign in a DataFrame work:
upgrade_vector(v::IndexedVector) = v

Base.sortperm{S,A<:AbstractDataVector}(x::IndexedVector{S,A}) = [findin(x.x, [NA]), x.idx]
Base.sortperm(x::IndexedVector) = x.idx
Base.sortperm(x::IndexedVector, ::Base.Sort.ReverseOrdering) = reverse(x.idx)
Base.sort(x::IndexedVector) = x.x[sortperm(x)]
Base.sort(x::IndexedVector, ::Base.Sort.ReverseOrdering) = x.x[reverse(sortperm(x))]
Base.Perm{O<:Base.Sort.Ordering}(o::O, v::IndexedVector) = FastPerm(o, v)

type Indexer
    r::Vector{Range1}
    iv::IndexedVector
    function Indexer(r::Vector{Range1}, iv::IndexedVector)
        for i in 1:length(r)
            ri = r[i]
            if length(ri) < 1 || ri[1] < 1 || ri[end] > length(iv)
                delete!(r, i)
            end
        end
        new(r, iv)
    end
end

function Base.intersect(v1::Vector{Range1}, v2::Vector{Range1})
    # Assumes that the Range1's in each vector are sorted and don't overlap.
    # (Actually, it doesn't, but it should to get more speed.)
    res = Range1[]
    # This does more work than it needs to.
    for r1 in v1
        for r2 in v2
            isoverlap = !((r1[1] > r2[end]) | (r2[1] > r1[end]))
            if isoverlap
                push!(res, max(r1[1], r2[1]):min(r1[end], r2[end]))
            end
        end
    end
    res
end

function Base.union(v1::Vector{Range1}, v2::Vector{Range1})
    # Assumes that the Range1's in each vector are sorted and don't overlap.
    # (Actually, it doesn't, but it should to get more speed.)
    # TODO Check for zero length.
    compare(v1, v2) = (v1[end][end] > v2[end][end] ? v2 : v1,
                       v1[end][end] > v2[end][end] ? v1 : v2)
    res = Range1[]
    v1 = copy(v1) # Destructively operate on these
    v2 = copy(v2)
    while length(v1) > 0 && length(v2) > 0
        # right part, working right to left
        (left, right) = compare(v1, v2)
        while length(right) > 0 && right[end][1] > left[end][end]
            push!(res, pop!(right))
        end
        if length(right) == 0 break; end
        # overlap
        r_end = max(right[end][end], left[end][end])
        overlap = false
        while length(left) > 0 && length(right) > 0 && right[end][end] > left[end][1]
            r_start = min(right[end][1], left[end][1])
            overlap = true
            (left, right) = compare(v1, v2)
            if right[end][1] >= r_start
                pop!(right)
            end
        end
        if overlap
            if length(right) == 0 && length(left) > 0 && r_start <= left[end][1]
                r_start = left[end][1]
                pop!(left)
            end
            push!(res, r_start : r_end)
        end
    end
    # Rest of v1 or v2 (no overlaps here)
    while length(v1) > 0
        push!(res, pop!(v1))
    end
    while length(v2) > 0
        push!(res, pop!(v2))
    end
    reverse(res)
end

function Base.(:!)(x::Indexer)   # Negate the Indexer
    res = Range1[1 : x.r[1][1] - 1]
    for i in 1:length(x.r) - 1
        push!(res, x.r[i][end] + 1 : x.r[i + 1][1] - 1)
    end
    push!(res, x.r[end][end] + 1 : length(x.iv.idx))
    Indexer(res, x.iv)
end

function Base.(:&)(x1::Indexer, x2::Indexer)
    if is(x1.iv, x2.iv)
        Indexer(intersect(x1.r, x2.r), x1.iv)
    else
        bool(x1) & bool(x2)
    end
end
Base.(:&)(x1::BitVector, x2::Indexer) = x1 & bool(x2)
Base.(:&)(x1::Indexer, x2::BitVector) = x2 & bool(x1)

function Base.(:|)(x1::Indexer, x2::Indexer)
    if is(x1.iv, x2.iv)
        Indexer(union(x1.r, x2.r), x1.iv)
    else
        bool(x1) | bool(x2)
    end
end
Base.(:|)(x1::BitVector, x2::Indexer) = x1 | bool(x2)
Base.(:|)(x1::Indexer, x2::BitVector) = x2 | bool(x1)

function Base.bool(ix::Indexer)
    res = falses(length(ix.iv))
    for i in ix.iv.idx[[ix.r...]]
        res[i] = true
    end
    res
end

# `getindex` -- each Range1 of the Indexer (i.r...) is applied to the indexing vector (i.iv.idx)

Base.getindex(x::IndexedVector, i::Indexer) = x[i.iv.idx[[i.r...]]]
Base.getindex(x::AbstractVector, i::Indexer) = x[i.iv.idx[[i.r...]]]
Base.getindex(x::AbstractDataVector, i::Indexer) = x[i.iv.idx[[i.r...]]]

# df[MultiRowIndex, SingleColumnIndex] => (Sub)?AbstractDataVector
function Base.getindex(df::DataFrame, row_inds::Indexer, col_ind::ColumnIndex)
    selected_column = df.colindex[col_ind]
    return df.columns[selected_column][row_inds.iv.idx[[row_inds.r...]]]
end

# df[MultiRowIndex, MultiColumnIndex] => (Sub)?DataFrame
function Base.getindex{T <: ColumnIndex}(df::DataFrame, row_inds::Indexer, col_inds::AbstractVector{T})
    selected_columns = df.colindex[col_inds]
    new_columns = {dv[row_inds.iv.idx[[row_inds.r...]]] for dv in df.columns[selected_columns]}
    return DataFrame(new_columns, Index(df.colindex.names[selected_columns]))
end

typealias ComparisonTypes Union(Number, String)   # defined mainly to avoid warnings

# element-wise (in)equality operators
# these may need range checks
# Should these results be sorted? Could be a counting sort.
Base.(:.==){T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : search_sorted_last(a.x, v, a.idx)], a)
Base.(:.==){T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : search_sorted_last(a.x, v, a.idx)], a)
Base.(:.>=){T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : length(a.idx)], a)
Base.(:.<=){T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[1 : search_sorted_last(a.x, v, a.idx)], a)
Base.(:.>=){T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[1 : search_sorted_last(a.x, v, a.idx)], a)
Base.(:.<=){T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : length(a.idx)], a)
Base.(:.>){T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[search_sorted_first_gt(a.x, v, a.idx) : length(a.idx)], a)
Base.(:.<){T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[1 : search_sorted_last_lt(a.x, v, a.idx)], a)
Base.(:.<){T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[search_sorted_first_gt(a.x, v, a.idx) : length(a.idx)], a)
Base.(:.>){T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[1 : search_sorted_last_lt(a.x, v, a.idx)], a)

function search_sorted_first_gt{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    res = search_sorted_last(a, x, idx)
    if res == 0 return 1 end
    if res == length(a) && a[idx[res]] != x return(length(a)+1) end
    a[idx[res]] == x ? res + 1 : res
end
function search_sorted_last_lt{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    res = search_sorted_first(a, x, idx)
    if res > length(idx) return length(idx) end
    if res == 1 && a[idx[res]] != x return(0) end
    a[idx[res]] == x ? res - 1 : res
end

Base.findin(a::IndexedVector,r::Range1) = findin(a, [r])
Base.findin(a::IndexedVector,v::Real) = findin(a, [v])
function Base.findin(a::IndexedVector, b::AbstractVector)
    ## Returns an Indexer with the elements in "a" that appear in "b"
    res = a .== b[1]
    for i in 2:length(b)
        res = res | (a .== b[i])
    end
    res
end

Base.size(a::IndexedVector) = size(a.x)
Base.length(a::IndexedVector) = length(a.x)
Base.ndims(a::IndexedVector) = 1
Base.eltype(a::IndexedVector) = eltype(a.x)

## print(io, a::IndexedVector) = print(io, a.x)
function Base.show(io::IO, a::IndexedVector)
    print(io, "IndexedVector: ")
    show(io, a.x)
end
function Base.repl_show(io::IO, a::IndexedVector)
    print(io, "IndexedVector: ")
    repl_show(io, a.x)
end

function search_sorted_last{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    ## Index of the last value of vector a that is less than or equal to x.
    ## Returns 0 if x is less than all values of a.
    ## idx is an indexing vector equal in length to a that sorts a
    ## @assert length(a) == length(idx)
    lo = 0
    hi = length(idx) + 1
    while lo < hi-1
        i = (lo+hi)>>>1
        if isless(x,a[idx[i]])
            hi = i
        else
            lo = i
        end
    end
    lo
end

function search_sorted_first{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    ## Index of the first value of vector a that is greater than or equal to x.
    ## Returns length(a) + 1 if x is greater than all values in a.
    ## idx is an indexing vector equal in length to a that sorts a
    ## @assert length(a) == length(idx)
    lo = 0
    hi = length(idx) + 1
    while lo < hi-1
        i = (lo+hi)>>>1
        if isless(a[idx[i]],x)
            lo = i
        else
            hi = i
        end
    end
    hi
end

# the following is needed for show(df)
maxShowLength(v::IndexedVector) = length(v) > 0 ? maximum([length(_string(x)) for x = v.x]) : 0

# Methods to speed up grouping and merging
function DataArrays.PooledDataArray{R}(d::IndexedVector, ::Type{R})
    refs = zeros(R, size(d))
    oneval = one(R)
    local idx::Int
    ## local lastval::T
    local poolidx::R
    pool = Array(eltype(d), 0)
    # skip over NAs
    nna = length(d) - length(d.x)
    if nna == length(d)
        return PooledDataArray(DataArrays.RefArray(refs), pool)
    end
    lastval = d.x[d.idx[nna+1]]
    push!(pool, d.x[d.idx[nna+1]])
    poolidx = oneval
    for i = 1 : length(d.idx)
        idx = d.idx[i]
        val = d.x[idx]
        if val != lastval
            push!(pool, val)
            poolidx += oneval
            lastval = val
        end
        refs[idx] = poolidx
    end
    return PooledDataArray(DataArrays.RefArray(refs), pool)
end
DataArrays.PooledDataArray(d::IndexedVector) = PooledDataArray(d, DEFAULT_POOLED_REF_TYPE)

DataArrays.DataArray(d::IndexedVector) = DataArray(x.x)

function DataArrays.PooledDataVecs{S}(v1::IndexedVector{S},
                           v2::IndexedVector{S})
    return PooledDataVecs(PooledDataArray(v1),
                          PooledDataArray(v2))
end


text:
"""
    gpdfitnew(x)

Estimate the paramaters for the Generalized Pareto Distribution (GPD). Returns empirical Bayes estimate for the parameters of the two-parameter generalized Parato distribution given the data.

# Arguments
* `x::AbstractArray`: One dimensional data array.

# Returns
* `k::Float64, sigma::Float64`: Estimated parameter values.

# Notes
* This function returns a negative of Zhang and Stephens's k, because it is more common parameterization.
"""
function gpdfitnew(x::AbstractArray)
    if ndims(x) != 1 || length(x) <= 1
        throw(DimensionMismatch())
    end
    sort!(x)

    n = length(x)
    m = 80 + trunc(Int, sqrt(n))
    prior = 3

    bs = collect(1:m) .- 0.5
    bs .= (1 .- sqrt.(m ./ bs)) ./ prior ./ x[trunc(Int64, n/4 + 0.5)] .+ 1 / x[end]

    ks = mean(log1p.(-bs .* x'), dims=2)[:,1]
    L = n * (log.(bs ./ -ks) .- ks .- 1)
    w = similar(L)
    w .= (1 ./ sum(exp.(L .- L'), dims=1)')[:,1]

    # remove negligible weights
    dii = w .>= (10 * eps(Float64))
    if ~all(dii)
        w = w[dii]
        bs = bs[dii]
    end
    # normalise w
    w ./= sum(w)
    # posterior mean for b
    b = sum(bs .* w)
    # Estimate for k, note that we return a negative of Zhang and
    # Stephens's k, because it is more common parameterisation.
    k = mean(log1p.(-b .* x))
    # estimate for sigma
    sigma = -k / b

    return k, sigma

end

ext: jl


text:
using Oceananigans.Operators
using Oceananigans.Solvers: FFTBasedPoissonSolver, FourierTridiagonalPoissonSolver, solve!
using Oceananigans.Architectures: device_event
using Oceananigans.Distributed: DistributedFFTBasedPoissonSolver

using PencilArrays: Permutation

#####
##### Calculate the right-hand-side of the non-hydrostatic pressure Poisson equation.
#####

const ZXYPermutation = Permutation{(3, 1, 2), 3}
const ZYXPermutation = Permutation{(3, 2, 1), 3}

@kernel function calculate_pressure_source_term_fft_based_solver!(rhs, grid, Δt, U★)
    i, j, k = @index(Global, NTuple)
    @inbounds rhs[i, j, k] = divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
end

@kernel function calculate_permuted_pressure_source_term_fft_based_solver!(rhs, grid, Δt, U★, ::ZXYPermutation)
    i, j, k = @index(Global, NTuple)
    @inbounds rhs[k, i, j] = divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
end

@kernel function calculate_permuted_pressure_source_term_fft_based_solver!(rhs, grid, Δt, U★, ::ZYXPermutation)
    i, j, k = @index(Global, NTuple)
    @inbounds rhs[k, j, i] = divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
end

@kernel function calculate_pressure_source_term_fourier_tridiagonal_solver!(rhs, grid, Δt, U★)
    i, j, k = @index(Global, NTuple)
    @inbounds rhs[i, j, k] = Δzᶜᶜᶜ(i, j, k, grid) * divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
end

#####
##### Solve for pressure
#####

function solve_for_pressure!(pressure, solver::DistributedFFTBasedPoissonSolver, Δt, U★)
    rhs = parent(first(solver.storage))
    arch = architecture(solver)
    grid = solver.local_grid

    rhs_event = launch!(arch, grid, :xyz, calculate_permuted_pressure_source_term_fft_based_solver!,
                        rhs, grid, Δt, U★, solver.input_permutation, dependencies = device_event(arch))

    wait(device(arch), rhs_event)

    # Solve pressure Poisson equation for pressure, given rhs
    solve!(pressure, solver)

    return pressure
end

function solve_for_pressure!(pressure, solver::FFTBasedPoissonSolver, Δt, U★)

    # Calculate right hand side:
    rhs = solver.storage
    arch = architecture(solver)
    grid = solver.grid

    rhs_event = launch!(arch, grid, :xyz, calculate_pressure_source_term_fft_based_solver!,
                        rhs, grid, Δt, U★, dependencies = device_event(arch))

    wait(device(arch), rhs_event)

    # Solve pressure Poisson given for pressure, given rhs
    solve!(pressure, solver, rhs)

    return nothing
end

function solve_for_pressure!(pressure, solver::FourierTridiagonalPoissonSolver, Δt, U★)

    # Calculate right hand side:
    rhs = solver.source_term
    arch = architecture(solver)
    grid = solver.grid

    rhs_event = launch!(arch, grid, :xyz, calculate_pressure_source_term_fourier_tridiagonal_solver!,
                        rhs, grid, Δt, U★, dependencies = device_event(arch))

    wait(device(arch), rhs_event)

    # Pressure Poisson rhs, scaled by Δzᶜᶜᶜ, is stored in solver.source_term:
    solve!(pressure, solver)

    return nothing
end


text:
using Dates

function tiktak_test()
    @testset "TikTakRule" begin
        nPoints = 17;
        nParams = 3;
        ttr = TikTakRule(; nPoints);
        h = mdl.make_test_history(nPoints, nParams);
        solnV = first(h.guessesV) .+ 0.1;
        fVal = 9.8;
        for j = 1 : nPoints
            startV = mdl.new_guess(ttr, h, j);
            # println(startV);
            @test size(startV) == size(solnV);
            p = Point(startV, solnV .+ (1/j), fVal + 1/j, fVal - 1/j,
                :success, Second(120 + j), 10 * j);
            mdl.add_to_history!(h, j, p);
        end
    end
end

@testset "UpdateRule" begin
    tiktak_test();
end

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


text:
using RLEVectors
using DataFrames

macro timeit(ex)
# like @time, but returning the timing rather than the computed value
  return quote
    #gc_disable()
    local val = $ex # compile
    local t0 = time()
    for i in 1:1e4 val = $ex end
    local t1 = time()
    #gc_enable()
    t1-t0
  end
end

foo = IntegerRle(Vector{Int32}(collect(1:1000)), Vector{Int32}(collect(5:5:5000)))
length(foo)
timings = DataFrame()
timings[:language] = "julia"
timings[:language_version] = VERSION
timings[:date] = chomp(readstring(`date "+%Y-%m-%d"`))
timings[:indexing] = @timeit foo[100]
timings[:range_indexing] = @timeit foo[801:900]
timings[:setting] = @timeit foo[800] = 5
timings[:range_setting] = @timeit foo[801:900] = 1:100
timings[:scalar_add] = @timeit foo + 4
timings[:length] = @timeit length(foo)
timings[:nrun] = @timeit nrun(foo)
timings[:mean] = @timeit mean(foo)
timings[:max] = @timeit maximum(foo)
timings[:width] = @timeit rwidth(foo)
timings[:last] = @timeit rlast(foo)
timings[:first] = @timeit rfirst(foo)
timings[:add_two_rles] = @timeit foo + foo
timings[:disjoin] = @timeit disjoin(foo.runends,foo.runends)
timings[:which] = @timeit findin(foo,[800,200,357])
timings[:scalar_less] = @timeit foo .< 3
timings[:median] = @timeit median(foo)
timings[:which_max] = @timeit findmax(foo)

bdf = DataFrames.readtable("/Users/phaverty/.julia/v0.4/RLEVectors/benchmark/rle.timings.csv",header=true);

for n in names(timings)
  if !(n in names(bdf))
    bdf[n] = NaN
  end
end

bdf = vcat(bdf,timings)

writetable( "/Users/phaverty/.julia/v0.4/RLEVectors/benchmark/rle.timings.csv",
             bdf, separator=',',header=true)

jdf = timings
rdf = bdf[ bdf[:language] .== "R",:];
rdf = rdf[end,:]

for n in names(bdf)[4:end]
       println(n)
       println( rdf[1,n] / jdf[1,n])
end
r_over_julia = zeros(ncol(bdf)-3)

for (i in 1:length(r_over_julia))
  r_over_julia[i] = log2(rdf[1,i+3] / jdf[1,i+3])
end

### Plotting
using Gadfly

## Performance Relative to R
bench_plot = plot(x=names(bdf)[4:end],y=r_over_julia, Geom.bar, Guide.ylabel("Elapsed Time: log2(R/julia)"),
     Guide.xticks(orientation=:vertical),Scale.color_continuous(minvalue=-15,maxvalue=15),color=r_over_julia,
     Guide.title("Relative Performance of R and Julia Rle Vectors"),Geom.hline(color="black"),yintercept=[0],Guide.xlabel(""))

date = jdf[1,:date]
relative_perf_file = "/Users/phaverty/.julia/v0.4/RLEVectors/benchmark/plots/benchmark_rle_vectors.$(date).png"
draw(PNG(relative_perf_file,8inch,5inch),bench_plot )
current_relative_perf_file = "/Users/phaverty/.julia/v0.4/RLEVectors/benchmark/plots/benchmark_rle_vectors.png"
cp(relative_perf_file, current_relative_perf_file, remove_destination=true)

## Performance over time
jdf = bdf[ bdf[:,:language] .== "julia", 3:end ]
melted_bdf = melt(jdf, :date)
timeline_plot = plot(melted_bdf, x="date", y="value", color="variable", Guide.xlabel("Date"), Geom.line,
                     Scale.y_log10, Guide.ylabel("log2 elapsed seconds (1e4 runs)"))
timeline_file = "/Users/phaverty/.julia/v0.4/RLEVectors/benchmark/plots/benchmark_rle_vectors.$(date).timeline.png"
draw(PNG(timeline_file,10inch,6inch),timeline_plot )
current_timeline_file = "/Users/phaverty/.julia/v0.4/RLEVectors/benchmark/plots/benchmark_rle_vectors.timeline.png"
cp(timeline_file, current_timeline_file, remove_destination=true)

## Profiling
using ProfileView
foo + foo; Profile.clear(); @profile for i in 1:1e6 foo + foo end; ProfileView.view()
foo .< 3; Profile.clear(); @profile for i in 1:1e4 foo .< 3 end; ProfileView.view()
foo .+ 3; Profile.clear(); @profile for i in 1:1e4 foo .+ 3 end; ProfileView.view()
sum(foo); Profile.clear(); @profile for i in 1:1e4 sum(foo) end; ProfileView.view()
findin(foo,[800,300,357]); Profile.clear(); @profile for i in 1:1e5 findin(foo,[800,300,357]) end; ProfileView.view()
median(foo); Profile.clear(); @profile for i in 1:1e4 median(foo) end; ProfileView.view()
rfirst(foo); Profile.clear(); @profile for i in 1:1e6 rfirst(foo) end; ProfileView.view()
rwidth(foo); Profile.clear(); @profile for i in 1:1e6 rwidth(foo) end; ProfileView.view()
foo[800] = 5; Profile.clear(); @profile for i in 1:1e6 foo[800] = 5 end; ProfileView.view()
foo[802] = 5; Profile.clear(); @profile for i in 1:1e6 foo[802] = 5 end; ProfileView.view()
foo[801:900] = 1:100; Profile.clear(); @profile for i in 1:1e5 foo[801:900] = 1:100 end; ProfileView.view()

x = collect(2:6:5000)
y = collect(6:4:5000)

foo = disjoin(x, y)
goo = disjoin2(x, y)
goo = goo[1:1667]
foo == goo

@time for i in 1:1e6 disjoin(x, y) end
@time for i in 1:1e6 disjoin2(x, y) end
Profile.clear(); @profile for i in 1:1e6 disjoin2(x, y) end; ProfileView.view()


text:
@testset "compute_int_J.jl" begin

end


text:
module MyCalculus

export EulerMethod

"""
    EulerMethod(f, x₀, y₀, step= 0.5, n=100.0)

Function to approximate solution to an ordinary differential equation using the Euler's Method.
dy/dx = f(x, y) where y(x₀) = y₀ and x₀ = x₀.

# Arguments
    f: function to approximate solution to, of the form f(x, y).
        This is the right hand side of the differential equation.
    x₀: initial x value condition
    y₀: initial y value condition
    step: step size for the Euler's Method
    n: number of steps to take

# Returns
    x: array of x values
    y: array of y values

"""
function EulerMethod(f::Function, x₀::Real, y₀::Real, step::Real, n::Real=100.0)

    y =[y₀] #initialize y array with intiial condition
    x = [x for x in range(x₀, step=step, length=n)] # array of x values starting with initial condition
    m = [] # this array will hold slope (right hand side of differential equation)

    push!(m, f(x[1], y[1])) # push initial slope onto array

    for k in 2:length(x) # iterate through x values starting at 2nd element

        push!(y, y[k-1] + m[k-1]*step) # push new y value to array using Euler's Method
        push!(m, f(x[k], y[k])) # push slope to array
    end

    return x, y

end

end


text:
export Differential, expand_derivatives, @derivatives

struct Differential <: Function
    x::Expression
end

Base.show(io::IO, D::Differential) = print(io, "(D'~", D.x, ")")
Base.convert(::Type{Expr}, D::Differential) = D

(D::Differential)(x::Operation) = Operation(D, Expression[x])
function (D::Differential)(x::Variable)
    D.x === x             && return Constant(1)
    has_dependent(x, D.x) || return Constant(0)
    return Operation(D, Expression[x])
end
(::Differential)(::Any) = Constant(0)
Base.:(==)(D1::Differential, D2::Differential) = isequal(D1.x, D2.x)

function expand_derivatives(O::Operation)
    @. O.args = expand_derivatives(O.args)

    if O.op isa Differential
        D = O.op
        o = O.args[1]
        isa(o, Operation) || return O
        return simplify_constants(sum(i->derivative(o,i)*expand_derivatives(D(o.args[i])),1:length(o.args)))
    end

    return O
end
expand_derivatives(x) = x

# Don't specialize on the function here
derivative(O::Operation, idx) = derivative(O.op, (O.args...,), Val(idx))

# Pre-defined derivatives
import DiffRules, SpecialFunctions, NaNMath
for (modu, fun, arity) ∈ DiffRules.diffrules()
    for i ∈ 1:arity
        @eval function derivative(::typeof($modu.$fun), args::NTuple{$arity,Any}, ::Val{$i})
            M, f = $(modu, fun)
            partials = DiffRules.diffrule(M, f, args...)
            dx = @static $arity == 1 ? partials : partials[$i]
            convert(Expression, dx)
        end
    end
end

function count_order(x)
    @assert !(x isa Symbol) "The variable $x must have an order of differentiation that is greater or equal to 1!"
    n = 1
    while !(x.args[1] isa Symbol)
        n = n+1
        x = x.args[1]
    end
    n, x.args[1]
end

_repeat_apply(f, n) = n == 1 ? f : f ∘ _repeat_apply(f, n-1)
function _differential_macro(x)
    ex = Expr(:block)
    lhss = Symbol[]
    x = flatten_expr!(x)
    for di in x
        @assert di isa Expr && di.args[1] == :~ "@derivatives expects a form that looks like `@derivatives D''~t E'~t`"
        lhs = di.args[2]
        rhs = di.args[3]
        order, lhs = count_order(lhs)
        push!(lhss, lhs)
        expr = :($lhs = $_repeat_apply(Differential($rhs), $order))
        push!(ex.args,  expr)
    end
    push!(ex.args, Expr(:tuple, lhss...))
    ex
end

macro derivatives(x...)
    esc(_differential_macro(x))
end

function calculate_jacobian(eqs,vars)
    Expression[Differential(vars[j])(eqs[i]) for i in 1:length(eqs), j in 1:length(vars)]
end


text:
import HTTP.Handlers: HandlerFunction, handle

function handle(::HandlerFunction{typeof(handle_rpcs)}, http::HTTP.Stream)
    HTTP.WebSockets.upgrade(handle_rpcs, http)
end

function version(req::HTTP.Request)
    host, port = INSTANCE[]
    resp = Dict(
        "Protocol-Version" => "1.3",
        "webSocketDebuggerUrl" => "ws://$host:$port/ws",
        )
    return HTTP.Response(200, JSON2.write(resp))
end

function main(req::HTTP.Request)
    host, port = INSTANCE[]
    resp = [
        Dict(
             "webSocketDebuggerUrl" => "ws://$host:$port/ws"
            ),
            ]
    return HTTP.Response(200, JSON2.write(resp))
end


text:
# The original license on the t-tests from HypothesisTests.jl
# is as follows. We're including it here because we're mostly copying the
# documentation of the functions.
#
# t.jl
# Various forms of t-tests
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

import HypothesisTests.OneSampleTTest
import HypothesisTests.EqualVarianceTTest
import HypothesisTests.UnequalVarianceTTest

"""
    EqualVarianceTTest(d1::AbstractUncertainValue, d2::AbstractUncertainValue,
        n::Int = 1000; μ0::Real = 0) -> EqualVarianceTTest

Consider two samples `s1` and `s2`, each consisting of `n` random draws from the
distributions furnishing `d1` and `d2`, respectively.

This function performs a two-sample t-test of the null hypothesis that `s1` and
`s2` come from distributions with equal means and variances against the
alternative hypothesis that the distributions have different means but equal
variances.
"""
function EqualVarianceTTest(d1::AbstractUncertainValue,
        d2::AbstractUncertainValue,
        n::Int = 1000;
        μ0::Real = 0)

    x = resample(d1, n)
    y = resample(d2, n)
    EqualVarianceTTest(x, y, μ0)
end

"""
    EqualVarianceTTestPooled(d1::UncertainDataset, d2::UncertainDataset,
        n::Int = 1000; μ0::Real = 0) -> EqualVarianceTTest

Consider two samples `s1[i]` and `s2[i]`, each consisting of `n` random draws
from the distributions furnishing the uncertain values `d1[i]` and `d2[i]`,
respectively. Gather all `s1[i]` in a pooled sample `S1`, and all `s2[i]` in
a pooled sample `S2`.

Perform a two-sample t-test of the null hypothesis that `S1` and `S2` come from
distributions with equal means and variances against the alternative hypothesis
that the distributions have different means but equal variances.
"""
function EqualVarianceTTestPooled(d1::UncertainDataset,
        d2::UncertainDataset,
        n::Int = 1000;
        μ0::Real = 0)

    x = vcat(resample(d1, n)...,)
    y = vcat(resample(d2, n)...,)

    EqualVarianceTTest(x, y, μ0)
end

"""
    EqualVarianceTTestElementWise(d1::UncertainDataset, d2::UncertainDataset,
        n::Int = 1000; μ0::Real = 0) -> Vector{EqualVarianceTTest}

Consider two samples `s1[i]` and `s2[i]`, each consisting of `n` random draws
from the distributions furnishing the uncertain values `d1[i]` and `d2[i]`,
respectively. This function performs an elementwise `EqualVarianceTTest`
on the pairs `(s1[i], s2[i])`. Specifically:

Performs an pairwise two-sample t-test of the null hypothesis that `s1[i]` and
`s2[i]` come from distributions with equal means and variances against the
alternative hypothesis that the distributions have different means but equal
variances.
"""
function EqualVarianceTTestElementWise(d1::UncertainDataset,
        d2::UncertainDataset,
        n::Int = 1000;
        μ0::Real = 0)

    N = length(d1)
    [EqualVarianceTTest(resample(d1[i], n), resample(d2[i], n), μ0) for i = 1:N]
end

"""
    UnequalVarianceTTest(d1::AbstractUncertainValue, d2::AbstractUncertainValue,
        n::Int = 1000; μ0::Real = 0) -> UnequalVarianceTTest

Consider two samples `s1` and `s2`, each consisting of `n` random draws from the
distributions furnishing `d1` and `d2`, respectively.

Perform an unequal variance two-sample t-test of the null hypothesis that `s1`
and `s2` come from distributions with equal means against the alternative
hypothesis that the distributions have different means.
"""
function UnequalVarianceTTest(d1::AbstractUncertainValue,
        d2::AbstractUncertainValue,
        n::Int = 1000;
        μ0::Real = 0)

    x = resample(d1, n)
    y = resample(d2, n)
    UnequalVarianceTTest(x, y, μ0)
end

"""
    UnequalVarianceTTestPooled(d1::UncertainDataset, d2::UncertainDataset,
        n::Int = 1000; μ0::Real = 0) -> UnequalVarianceTTest

Consider two samples `s1[i]` and `s2[i]`, each consisting of `n` random draws
from the distributions furnishing the uncertain values `d1[i]` and `d2[i]`,
respectively. Gather all `s1[i]` in a pooled sample `S1`, and all `s2[i]` in
a pooled sample `S2`.

This function performs an unequal variance two-sample t-test of the null
hypothesis that `S1` and `S2` come from distributions with equal means against
the alternative hypothesis that the distributions have different means.

"""
function UnequalVarianceTTestPooled(d1::UncertainDataset,
        d2::UncertainDataset,
        n::Int = 1000;
        μ0::Real = 0)

    x = vcat(resample(d1, n)...,)
    y = vcat(resample(d2, n)...,)

    UnequalVarianceTTest(x, y, μ0)
end

"""
    UnequalVarianceTTestElementWise(d1::UncertainDataset, d2::UncertainDataset,
        n::Int = 1000; μ0::Real = 0) -> Vector{UnequalVarianceTTest}

Consider two samples `s1[i]` and `s2[i]`, each consisting of `n` random draws
from the distributions furnishing the uncertain values `d1[i]` and `d2[i]`,
respectively. This function performs an elementwise `EqualVarianceTTest`
on the pairs `(s1[i], s2[i])`. Specifically:

Performs an pairwise unequal variance two-sample t-test of the null hypothesis
that `s1[i]` and `s2[i]` come from distributions with equal means against the
alternative hypothesis that the distributions have different means.

This test is sometimes known as Welch's t-test. It differs from the equal
variance t-test in that it computes the number of degrees of freedom of the
test using the Welch-Satterthwaite equation:

```math
    ν_{χ'} ≈ \\frac{\\left(\\sum_{i=1}^n k_i s_i^2\\right)^2}{\\sum_{i=1}^n
        \\frac{(k_i s_i^2)^2}{ν_i}}
```
"""
function UnequalVarianceTTestElementWise(d1::UncertainDataset,
        d2::UncertainDataset,
        n::Int = 1000;
        μ0::Real = 0)

    N = length(d1)
    [UnequalVarianceTTest(resample(d1[i], n), resample(d2[i], n), μ0) for i = 1:N]
end

"""
    OneSampleTTest(d::AbstractUncertainValue, n::Int = 1000;
        μ0::Real = 0) -> OneSampleTTest

Perform a one sample t-test of the null hypothesis that the uncertain value has
a distribution with mean `μ0` against the alternative hypothesis that its
distribution does not have mean `μ0`. `n` indicates the number of draws during
resampling.
"""
function OneSampleTTest(d::AbstractUncertainValue, n::Int = 1000; μ0::Real = 0)
    x = resample(d, n)
    OneSampleTTest(x, μ0)
end

"""
    OneSampleTTest(d1::AbstractUncertainValue,
                    d2::AbstractUncertainValue,
                    n::Int = 1000; μ0::Real=0) -> OneSampleTTest

Perform a paired sample t-test of the null hypothesis that the differences
between pairs of uncertain values in `d1` and `d2` come from a distribution
with mean `μ0` against the alternative hypothesis that the distribution does
not have mean `μ0`.
"""
function OneSampleTTest(d1::AbstractUncertainValue,
                        d2::AbstractUncertainValue,
                        n::Int = 1000; μ0::Real = 0)
    x = resample(d1, n)
    y = resample(d2, n)
    OneSampleTTest(x, y, μ0)
end

"""
    OneSampleTTestPooled(d::UncertainDataset, n::Int = 1000;
        μ0::Real=0) -> OneSampleTTest

First, sample `n` draws of each uncertain value in `d1`, then pooling the draws
together. Then, perform a one sample t-test of the null hypothesis that the
uncertain values have a pooled distribution with mean `μ0` against the
alternative hypothesis that its pooled distribution does not have mean `μ0`.
`n` indicates the number of draws during resampling.
"""
function OneSampleTTestPooled(d::UncertainDataset, n::Int = 1000; μ0::Real = 0)
    x = vcat([resample(uncertainval, n) for uncertainval in d]...,)
    OneSampleTTest(x, μ0)
end

"""
    OneSampleTTestElementWise(d::UncertainDataset, n::Int = 1000;
        μ0::Real = 0) -> Vector{OneSampleTTest}

First, sample `n` draws of each uncertain value in `d`, then pooling the draws
together. Then, perform a one sample t-test of the null hypothesis that the
uncertain values have a pooled distribution with mean `μ0` against the
alternative hypothesis that its pooled distribution
does not have mean `μ0`. `n` indicates the number of draws during resampling.
"""
function OneSampleTTestElementWise(d::UncertainDataset, n::Int = 1000; μ0::Real = 0)
    [OneSampleTTest(resample(uncertainval, n), μ0) for uncertainval in d]
end

"""
    OneSampleTTestPooled(d1::UncertainDataset,
        d2::UncertainDataset,
        n::Int = 1000; μ0::Real = 0) -> OneSampleTTest

First, sample `n` draws of each uncertain value in each dataset, pooling the
draws from the elements of `d1` and the draws from the elements of `d2` separately.
Then, perform a paired sample t-test of the null hypothesis that the differences
between pairs of uncertain values in `d1` and `d2` come from a distribution
with mean `μ0` against the alternative hypothesis that the distribution does
not have mean `μ0`.
"""
function OneSampleTTestPooled(d1::UncertainDataset,
                                d2::UncertainDataset,
                                n::Int = 1000; μ0::Real = 0)
    x = vcat(resample(d1, n)...,)
    y = vcat(resample(d2, n)...,)
    OneSampleTTest(x, y, μ0)
end

"""
    OneSampleTTestElementWise(d1::UncertainDataset,
        d2::UncertainDataset,
        n::Int = 1000; μ0::Real = 0) -> Vector{OneSampleTTest}

Perform a one sample t-test of the null hypothesis that the uncertain value has
a distribution with mean `μ0` against the alternative hypothesis that its
distribution does not have mean `μ0` for uncertain value in `d`.

`n` indicates the number of draws during resampling.
"""
function OneSampleTTestElementWise(d1::UncertainDataset,
                                    d2::UncertainDataset,
                                    n::Int = 1000; μ0::Real = 0)
    N = length(d1)
    [OneSampleTTest(resample(d1[i], n), resample(d2[i], n), μ0) for i in 1:N]
end

export
EqualVarianceTTest,
EqualVarianceTTestPooled,
EqualVarianceTTestElementWise,

UnequalVarianceTTest,
UnequalVarianceTTestPooled,
UnequalVarianceTTestElementWise,

OneSampleTTest,
OneSampleTTestPooled,
OneSampleTTestElementWise


text:
module AcuteML

import EzXML.Node

include("utilities.jl")
include("templating.jl")

# main macro
export @aml
# types
export Node, Document
# literals
export @xd_str, @hd_str, @sc_str, @a_str
################################################################
"""
  @aml typedef

# Type defnition
Use `@aml` macro to define a Julia type, and then the package automatically creates a xml or html associated with the defined type.

### Document Defnition
* Use `xd""` or `hd""` to define a XML or HTML document:
```julia
@aml mutable struct Doc xd""
# add fields (elements) here
end
```

### Nodes (Elements) Defnition
* Specify the html/xml struct name as a string after the struct name after a space
```julia
@aml mutable struct Person "person"
# add fields (elements) here
end
```
* If the html/xml name is the same as struct name, you can use `"~"` instead
```julia
@aml mutable struct person "~"
# add fields (elements) here
end
```

### Fields Names
* Sepecify the html/xml field name as a string in front of the field after `,`
```julia
field, "study-field"
```
* If the html/xml name is the same as variable name, you can use `"~"` instead
```julia
age::UInt, "~"
```

### Attributes
* If the value is going to be an attribute put `a` before its name
```julia
id::Int64, a"~"
```

### Default Value
* You can specify the default value for an argument by using `= defVal` syntax
```julia
GPA::Float64 = 4.5, "~"
```

### Value Types
You can use Julia types or  defined types for values.

* If you don't specify the type of a variable, it is considered to be string for aml manipulations:
```julia
field, "study-field"
```
However, for a high performance code specify String type (`field::String, "study-field"`)

* For already `@aml` defined types, name should be the same as the defined type root name
```julia
university::University, "university"
```

### Value Checking
You can define any restriction for values using functions.

* To define any restrictions for the values of one field, define a function that checks a criteria for the field value and returns Bool, and put its name after a `,` after the field name:
```julia
GPA::Float64, "~", GPAcheck
```

* To define any restrictions for multiple values of a struct, define a function that gets all the variables and checks a criteria and returns Bool, and put its name after a `,` after the struct name:
```julia
@aml mutable struct Person "person", courseCheck
# ...
end
```

Refer to https://aminya.github.io/AcuteML.jl/dev/valueChecking/ for some of these functions examples.

### Optional Fields
* If a field is optional, don't forget to define its type as `UN{}` (Union with Nothing), and set the default value as `nothing`.
```julia
residence::UN{String}=nothing, "residence-stay" # optional with nothing as default value
```
```julia
funds::UN{String}, "financial-funds"   # optional, but you should pass nothing manually in construction
```

### Empty Elements (Self-Closing) Defnition
* Use `sc"name"` to define a self-closing (empty) element (e.g. `<rest />`)
```julia
@aml struct rest sc"~"
end
```

-------------------------------------------------------

# Example - Struct Definition

First, we define the structs using `@aml` to store the data in:

```julia
using AcuteML

# Types definition

# Person Type
@aml mutable struct Person "person", courseCheck
    age::UInt64, "~"
    field, "study-field"
    GPA::Float64 = 4.5, "~", GPAcheck
    courses::Vector{String}, "taken-courses"
    id::Int64, a"~"
end

@aml mutable struct University "university"
    name, a"university-name"
    people::Vector{Person}, "person"
end

@aml mutable struct Doc xd""
    university::University, "~"
end

```

```julia
# Value Checking Functions
GPAcheck(x) = x <= 4.5 && x >= 0

function courseCheck(age, field, GPA, courses, id)

    if field == "Mechanical Engineering"
        relevant = ["Artificial Intelligence", "Robotics", "Machine Design"]
    elseif field == "Computer Engineering"
        relevant = ["Julia", "Algorithms"]
    else
        error("study field is not known")
    end

    return any(in.(courses, Ref(relevant)))
end
```
-------------------------------------------------------

# Example - Constructor

After we defined the structs, we can create instances of them by passing our data to the fields:

```julia

P1 = Person(age=24, field="Mechanical Engineering", courses=["Artificial Intelligence", "Robotics"], id = 1)
P2 = Person(age=18, field="Computer Engineering", GPA=4, courses=["Julia"], id = 2)

P2.GPA=4.2 # mutability support

U = University(name="Julia University", people=[P1, P2])

D = Doc(university = U)
```

```julia
# An example that doesn't meet the critertia function for GPA because GPA is more than 4.5
P3 = Person(age=99, field="Macro Wizard", GPA=10, courses=["Julia Magic"], id = 3)
julia>
GPA doesn't meet criteria function
```

```html
julia> print(P1.aml)
<person id="1">
  <age>24</age>
  <study-field>Mechanical Engineering</study-field>
  <GPA>4.5</GPA>
  <taken-courses>Artificial Intelligence</taken-courses>
  <taken-courses>Robotics</taken-courses>
</person>

julia> print(U.aml)
<university university-name="Julia University">
  <person id="1">
    <age>24</age>
    <study-field>Mechanical Engineering</study-field>
    <GPA>4.5</GPA>
    <taken-courses>Artificial Intelligence</taken-courses>
    <taken-courses>Robotics</taken-courses>
  </person>
  <person id="2">
    <age>18</age>
    <study-field>Computer Engineering</study-field>
    <GPA>4.2</GPA>
    <taken-courses>Julia</taken-courses>
  </person>
</university>

julia> print(D.aml)
<?xml version="1.0" encoding="UTF-8" standalone="yes"?><!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd">
<university university-name="Julia University">
  <person id="1">
    <age>24</age>
    <study-field>Mechanical Engineering</study-field>
    <GPA>4.5</GPA>
    <taken-courses>Artificial Intelligence</taken-courses>
    <taken-courses>Robotics</taken-courses>
  </person>
  <person id="2">
    <age>18</age>
    <study-field>Computer Engineering</study-field>
    <GPA>4.2</GPA>
    <taken-courses>Julia</taken-courses>
  </person>
</university>
```
-------------------------------------------------------

# Example - Extractor

After we defined the structs, we can automatically extract and store the data in their fields:

```julia
using AcuteML

xml = parsexml(\"\"\"
<?xml version="1.0" encoding="UTF-8" standalone="yes"?><!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd">
<university university-name="Julia University">
  <person id="1">
    <age>24</age>
    <study-field>Mechanical Engineering</study-field>
    <GPA>4.5</GPA>
    <taken-courses>Artificial Intelligence</taken-courses>
    <taken-courses>Robotics</taken-courses>
  </person>
  <person id="2">
    <age>18</age>
    <study-field>Computer Engineering</study-field>
    <GPA>4</GPA>
    <taken-courses>Julia</taken-courses>
  </person>
</university>
\"\"\")

# extract Doc
D = Doc(xml) # StructName(xml) like Doc(xml) extracts the data and stores them in proper format

# Now you can access all of the data by calling the fieldnames

# extract University
U = D.university

julia>U.name
"Julia University"

# extract Person

P1 = U.people[1]

julia>P1.age
24

julia>P1.field
Mechanical Engineering

julia>P1.GPA
4.5

julia>P1.courses
["Artificial Intelligence", "Robotics"]

julia>P1.id
1
```
"""
macro aml(expr)
    expr = macroexpand(__module__, expr) # to expand @static

    #  check if aml is used before struct
    expr isa Expr && expr.head == :struct || error("Invalid usage of @aml")

    mutability = expr.args[1]
    T = expr.args[2] # Type name +(curly braces)

    aml = Symbol(:aml)

    # expr.args[3] # arguments
    # argParams.args # empty
    expr.args[3], argParams, argDefVal, argTypes, argVars, argNames, argFun, amlTypes, amlName, docOrElmType, amlFun = _aml(expr)

    # defining outter constructors
    # Only define a constructor if the type has fields, otherwise we'll get a stack
    # overflow on construction
    if !isempty(argVars)

        ################################################################
        # Arguments methods

        # Non-aml arguments are ignored
        idAmlArgs = (!).(ismissing.(argNames)) # missing argName means not aml

        amlNames = argNames[idAmlArgs]
        amlFuns = argFun[idAmlArgs]
        amlVars = argVars[idAmlArgs]
        amlParams = argParams[idAmlArgs]
        amlDefVal = argDefVal[idAmlArgs]
        amlTypes = amlTypes[idAmlArgs]
        argTypes  = argTypes[idAmlArgs]

        numAml = length(amlVars)

        amlconst=Vector{Expr}(undef,numAml)
        amlext=Vector{Expr}(undef,numAml)
        amlmutability=Vector{Expr}(undef,numAml)

        amlVarsCall = Vector{Expr}(undef,numAml)

        ##########################
        # Each argument of the struct
        for i=1:numAml
            argTypesI = argTypes[i]
            amlVarsI = amlVars[i]
            amlNamesI = amlNames[i]
            amlTypesI = amlTypes[i]
            amlFunsI=amlFuns[i]
            amlSymI=QuoteNode(amlVarsI)
            ##########################
            # call Expr - For mutability
            amlVarsCall[i] = :(str.$amlVarsI)

            # Vector
            if (isa(argTypesI, Expr) && argTypesI.args[1] == :Vector) || (!isa(argTypesI, Union{Symbol, Expr}) && argTypesI <: Array)

                # Function missing
                if ismissing(amlFunsI)

                    amlconst[i]=:(addelementVect!(aml, $amlNamesI, $amlVarsI, $amlTypesI))

                    amlext[i]=:($amlVarsI = findallcontent($(esc(argTypesI)), $amlNamesI, aml, $amlTypesI))

                    if mutability
                        amlmutability[i] = quote
                            if name == $amlSymI
                                updateallcontent!(value, $amlNamesI, str.aml, $amlTypesI)
                            end
                        end
                    end

                # Function provided
                else
                    amlconst[i]=quote
                        if ($(esc(amlFunsI)))($amlVarsI)
                            addelementVect!(aml, $amlNamesI, $amlVarsI, $amlTypesI)
                        else
                            error("$($amlNamesI) doesn't meet criteria function")
                        end
                    end

                    amlext[i]=quote

                        $amlVarsI = findallcontent($(esc(argTypesI)), $amlNamesI, aml, $amlTypesI)

                        if !(($(esc(amlFunsI)))($amlVarsI))
                            error("$($amlNamesI) doesn't meet criteria function")
                        end
                    end

                    if mutability
                        amlmutability[i] = quote
                            if name == $amlSymI
                                if ($(esc(amlFunsI)))($(amlVarsCall[i]))
                                    updateallcontent!(value, $amlNamesI, str.aml, $amlTypesI)
                                else
                                    error("$($amlNamesI) doesn't meet criteria function")
                                end
                            end
                        end
                    end

                end

            ##########################
            # Non Vector
            elseif isa(argTypesI, Symbol) || (isa(argTypesI, Expr) && argTypesI.args[1] == :Union ) || (isa(argTypesI, Expr) && argTypesI.args[1] == :UN) || !(argTypesI <: Array)

                # Function missing
                if ismissing(amlFunsI)

                    amlconst[i]=:(addelementOne!(aml, $amlNamesI, $amlVarsI, $amlTypesI))
                    amlext[i]=:($amlVarsI = findfirstcontent($(esc(argTypesI)), $amlNamesI, aml, $amlTypesI))

                    if mutability

                        amlmutability[i] = quote
                            if name == $amlSymI
                                updatefirstcontent!(value, $amlNamesI, str.aml, $amlTypesI)
                            end
                        end
                    end

                # Function provided
                else
                    amlconst[i]=quote
                        if ($(esc(amlFunsI)))($amlVarsI)
                            addelementOne!(aml, $amlNamesI, $amlVarsI, $amlTypesI)
                        else
                            error("$($amlNamesI) doesn't meet criteria function")
                        end
                    end

                    amlext[i]=quote

                        $amlVarsI = findfirstcontent($(esc(argTypesI)), $amlNamesI, aml, $amlTypesI)

                        if !(($(esc(amlFunsI)))($amlVarsI))
                            error("$($amlNamesI) doesn't meet criteria function")
                        end
                    end
                    if mutability
                        amlmutability[i] = quote
                            if name == $amlSymI
                                if ($(esc(amlFunsI)))($(amlVarsCall[i]))
                                    updatefirstcontent!(value, $amlNamesI, str.aml, $amlTypesI)
                                else
                                    error("$($amlNamesI) doesn't meet criteria function")
                                end
                            end
                        end
                    end

                end

            end

        end # endfor

        ################################################################
        # aml Function
        if !ismissing(amlFun[1])
            F=amlFun[1]
            amlFunChecker = quote
                if !( ($(esc(F)))($(argVars...)) )
                    error("struct criteria function ($($(esc(F)))) isn't meet")
                end
            end
        else
            amlFunChecker = nothing
        end

        if mutability
            if !ismissing(amlFun[1])
                F=amlFun[1]
                amlFunCheckerMutability = quote
                    if !( ($(esc(F)))($(amlVarsCall...)) )
                        error("struct criteria function ($($(esc(F)))) isn't meet")
                    end
                end
            else
                amlFunCheckerMutability = nothing
            end
        end

        ################################################################
        # Type name is a single name (symbol)
        if T isa Symbol

            docOrElmconst = :( aml = docOrElmInit($docOrElmType, $amlName) )

            typeDefinition =:($expr)

            amlConstructor = quote
                function ($(esc(T)))(; $(argParams...))
                    $amlFunChecker
                    $docOrElmconst
                    $(amlconst...)
                    return ($(esc(T)))($(argVars...), aml)
                end
            end

            amlExtractor = quote
                function ($(esc(T)))(aml::Union{Document, Node})
                    $(amlext...)
                    $amlFunChecker
                    return ($(esc(T)))($(argVars...), aml)
                end

            end

            nothingMethod = :( ($(esc(T)))(::Nothing) = nothing )
            # convertNothingMethod = :(Base.convert(::Type{($(esc(T)))}, ::Nothing) = nothing) # for passing nothing to function without using Union{Nothing, T} in the definition
            selfMethod = :( ($(esc(T)))(in::$(esc(T))) = $(esc(T))(in.aml) )

            if mutability
                mutabilityExp = quote
                     function Base.setproperty!(str::($(esc(T))),name::Symbol, value)
                         setfield!(str,name,value)
                         $amlFunCheckerMutability
                         $(amlmutability...)
                     end
                 end
            else
                mutabilityExp = nothing
            end

            out = quote
                Base.@__doc__($(esc(typeDefinition)))
                $amlConstructor
                $amlExtractor
                $nothingMethod
                # $convertNothingMethod
                $selfMethod
                $mutabilityExp
            end

        ################################################################
        # Parametric type structs
        elseif T isa Expr && T.head == :curly
            # if T == S{A<:AA,B<:BB}, define two methods
            #   S(...) = ...
            #   S{A,B}(...) where {A<:AA,B<:BB} = ...
            S = T.args[1]
            P = T.args[2:end]
            Q = [U isa Expr && U.head == :<: ? U.args[1] : U for U in P]
            SQ = :($S{$(Q...)})

            docOrElmconst = :( aml = docOrElmInit($docOrElmType, $amlName) )

            typeDefinition =:($expr)

            amlConstructor = quote
                function ($(esc(S)))(; $(argParams...))
                    $amlFunChecker
                    $docOrElmconst
                    $(amlconst...)
                    return ($(esc(S)))($(argVars...), aml)
                end
            end

            amlConstructorCurly = quote
                function ($(esc(SQ)))(; $(argParams...)) where {$(esc.(P)...)}
                    $amlFunChecker
                    $docOrElmconst
                    $(amlconst...)
                    return ($(esc(SQ)))($(argVars...), aml)
                end
            end

            amlExtractor = quote
                function ($(esc(S)))(aml::Union{Document, Node})
                    $(amlext...)
                    $amlFunChecker
                    return ($(esc(S)))($(argVars...), aml)
                end

            end

            amlExtractorCurly = quote
                function ($(esc(SQ)))(aml::Union{Document, Node}) where {$(esc.(P)...)}
                    $(amlext...)
                    $amlFunChecker
                    return ($(esc(SQ)))($(argVars...), aml)
                end

            end

            nothingMethod = :( ($(esc(S)))(::Nothing) = nothing )
            # convertNothingMethod = :(Base.convert(::Type{($(esc(S)))}, ::Nothing) = nothing) # for passing nothing to function without using Union{Nothing, ...} in the definition
            selfMethod = :( ($(esc(S)))(in::$(esc(S))) = $(esc(S))(in.aml) )

            if mutability
                mutabilityExp = quote
                     function Base.setproperty!(str::($(esc(T))),name::Symbol, value)
                         setfield!(str,name,value)
                         $amlFunCheckerMutability
                         $(amlmutability...)
                     end
                 end
             else
                 mutabilityExp = nothing
            end

             out = quote
                 Base.@__doc__($(esc(typeDefinition)))
                 $amlConstructor
                 $amlConstructorCurly
                 $amlExtractor
                 $amlExtractorCurly
                 $nothingMethod
                 # $convertNothingMethod
                 $selfMethod
                 $mutabilityExp
             end
        ################################################################
        else
            error("Invalid usage of @aml")
        end
    else
        out = nothing
    end

    return out
end

################################################################
# @aml Parser Function
function _aml(expr)

    # reminder:
    # var is a symbol
    # var::T or anything more complex is an expression

    argExpr = expr.args[3] # arguments of the type
    T = expr.args[2] # Type name +(curly braces)

    argParams = Union{Expr,Symbol}[] # Expr(:parameters)[]
    argVars = Union{Expr,Symbol}[]
    argDefVal = Any[]
    argTypes = Union{Missing,Type, Symbol, Expr}[]
    argNames = Union{Missing,String}[]
    argFun = Union{Missing, Symbol, Function}[]
    amlTypes = Int64[]
    amlName = "my type"
    docOrElmType = 0
    amlFun = Array{Union{Missing, Symbol, Function},0}(undef)

    for i in eachindex(argExpr.args) # iterating over arguments of each type argument
        ei = argExpr.args[i] # type argument element i

        ########################
        # Line number skipper
        if typeof(ei) == LineNumberNode
            continue
        end

        ########################
        # Single struct name - "aml name"
        if isa(ei, String)

            amlFun[1]=missing # function

            # Self-name checker
            if ei == "~"
                if T isa Symbol
                    amlName = string(T)
                elseif T isa Expr && T.head == :curly
                    amlName = string(T.args[1]) # S
                end
            else
                amlName = ei  # Type aml name
            end

            argExpr.args[i]= nothing # removing "aml name" from expr args
            docOrElmType = 0

        ################################################################
        # Literal Struct name - xd/hd"aml name"
        elseif isa(ei, Tuple)
            ################################################################
            # Struct aml
            ########################
            # Literal only  xd/hd"aml name"
            if isa(ei, Tuple{Int64,String})

                amlFun[1]=missing # function

                # Self-name checker
                if ei[2] == "~"
                    if T isa Symbol
                        amlName = string(T)
                    elseif T isa Expr && T.head == :curly
                        amlName = string(T.args[1]) # S
                    end
                else
                    amlName = ei[2] # Type aml name
                end

                docOrElmType = ei[1]

                argExpr.args[i]= nothing # removing "aml name" from expr args
            end

        elseif ei.head == :tuple
                ########################
                # Struct Function - "aml name", F
            if isa(ei.args[1], String) && isa(ei.args[2], Union{Symbol,Function}) # "aml name", F

                amlFun[1]=ei.args[2] # function

                # Self-name checker
                if ei.args[1] == "~"
                    if T isa Symbol
                        amlName = string(T)
                    elseif T isa Expr && T.head == :curly
                        amlName = string(T.args[1]) # S
                    end
                else
                    amlName = ei.args[1] # Type aml name
                end

                docOrElmType = 0
                argExpr.args[i]= nothing # removing "aml name" from expr args

                ########################
                # Literal and Struct Function - xd/hd"aml name", F
            elseif isa(ei.args[1], Tuple)  && isa(ei.args[2], Union{Symbol,Function})

                amlFun[1]=ei.args[2] # function

                # Self-name checker
                if ei.args[1][2] == "~"
                    if T isa Symbol
                        amlName = string(T)
                    elseif T isa Expr && T.head == :curly
                        amlName = string(T.args[1]) # S
                    end
                else
                    amlName = ei.args[1][2] # Type aml name
                end

                docOrElmType = ei.args[1][1]
                argExpr.args[i]= nothing # removing "aml name" from expr args
        ################################################################
        # Arguments
            ########################
            # No Def Value
            elseif ei.args[1] isa Union{Symbol,Expr} # var/var::T, "name"

                # Def Value
                push!(argDefVal, missing)

                # Type Checker
                lhs = ei.args[1]
                if lhs isa Symbol #  var, "name"

                    var = ei.args[1]

                    push!(argTypes, String) # consider String as the type
                    push!(argParams, var)
                    push!(argVars, var)

                    argExpr.args[i]=var  # removing "name",...

                elseif lhs isa Expr && lhs.head == :(::) && lhs.args[1] isa Symbol # var::T, "name"

                    var = lhs.args[1]
                    varType = lhs.args[2] # Type

                    push!(argTypes, varType)
                    push!(argParams, var)
                    push!(argVars, var)

                    argExpr.args[i]=lhs  # removing "name",...

                end

                # Literal Checker
                if length(ei.args[2]) == 2 # literal

                    elmType = ei.args[2][1]
                    push!(amlTypes, elmType) # literal type

                    ni = ei.args[2][2]

                    # Self-name checker
                    if ni == "~"
                        push!(argNames,string(var))
                    else
                        push!(argNames,ni)
                    end

                else
                    push!(amlTypes, 0) # non-literal

                    ni = ei.args[2]

                    # Self-name checker
                    if ni == "~"
                        push!(argNames,string(var))
                    else
                        push!(argNames,ni)
                    end
                end

                # Function Checker
                if length(ei.args) == 3 && isa(ei.args[3], Union{Function, Symbol}) #  var/var::T, "name", f

                    fun = ei.args[3]   # function
                    push!(argFun, fun)

                else # function name isn't given
                    push!(argFun, missing)
                end

            end  # end Tuple sub possibilities
        ################################################################
        # Def Value
        elseif ei.head == :(=) # def value provided

            # aml name Checker
            if ei.args[2].head == :tuple # var/var::T = defVal, "name"

                # Def Value
                defVal = ei.args[2].args[1]

                push!(argDefVal, defVal)

                lhs = ei.args[1]

                argExpr.args[i]=lhs # remove =defVal for type definition

                # Type Checker
                if lhs isa Symbol #  var = defVal, "name"

                    var = ei.args[1]

                    push!(argTypes, String) # consider String as the type
                    push!(argParams, Expr(:kw, var, defVal))
                    push!(argVars, var)

                    argExpr.args[i]=var  # removing "name",...

                elseif lhs isa Expr && lhs.head == :(::) && lhs.args[1] isa Symbol # var::T = defVal, "name"

                    var = lhs.args[1]
                    varType = lhs.args[2] # Type

                    push!(argTypes, varType)
                    push!(argParams, Expr(:kw, var, defVal)) # TODO also put type expression
                    push!(argVars, var)

                    argExpr.args[i]=lhs  # removing "name",...

                end

                # Literal Checker
                if length(ei.args[2].args[2]) == 2 # literal

                    elmType = ei.args[2].args[2][1]
                    push!(amlTypes, elmType) # literal type

                    ni = ei.args[2].args[2][2]

                    # Self-name checker
                    if ni == "~"
                        push!(argNames,string(var))
                    else
                        push!(argNames,ni)
                    end

                else
                    push!(amlTypes, 0) # non-literal

                    ni = ei.args[2].args[2]

                    # Self-name checker
                    if ni == "~"
                        push!(argNames,string(var))
                    else
                        push!(argNames,ni)
                    end

                end

                # Function Checker
                if length(ei.args[2].args) == 3 && isa(ei.args[2].args[3], Union{Function, Symbol}) #  var/var::T  = defVal, "name", f

                    fun = ei.args[2].args[3]  # function
                    push!(argFun, fun)

                else # function name isn't given
                    push!(argFun, missing)
                end

            ########################
            #  No aml Name - But defVal
            else # var/var::T = defVal # ignored for creating aml

                # Type Checker
                lhs = ei.args[1]
                if lhs isa Symbol #  var = defVal

                    defVal = ei.args[2]

                    push!(argDefVal, defVal)
                    push!(argNames,missing) # ignored for creating aml
                    push!(argFun, missing) # ignored for creating aml

                    var = ei.args[1]

                    push!(argTypes, Any)
                    push!(argParams, Expr(:kw, var, defVal))
                    push!(argVars, var)

                    argExpr.args[i]=var # remove =defVal for type definition

                elseif lhs isa Expr && lhs.head == :(::) && lhs.args[1] isa Symbol # var::T = defVal

                    defVal = ei.args[2]

                    push!(argDefVal, defVal)
                    push!(argNames,missing) # ignored for creating aml
                    push!(argFun, missing) # ignored for creating aml

                    var = lhs.args[1]
                    varType = lhs.args[2] # Type

                    push!(argTypes, varType)
                    push!(argParams, Expr(:kw, var, defVal)) # TODO also put type expression
                    push!(argVars, var)

                    argExpr.args[i]=lhs # remove =defVal for type definition

                else
                    # something else, e.g. inline inner constructor
                    #   F(...) = ...
                    continue
                end

            end

        ################################################################
        # No aml name - No defVal
        else  # var/var::T  # ignored for creating aml

            # Type Checker
            if ei isa Symbol #  var
                push!(argNames, missing) # argument ignored for aml
                push!(argFun, missing) # ignored for creating aml

                push!(argTypes, String)

                var = ei

                push!(argParams, var)
                push!(argVars, var)

            elseif ei.head == :(::) && ei.args[1] isa Symbol # var::T
                push!(argNames, missing) # argument ignored for aml
                push!(argFun, missing) # ignored for creating aml

                var = ei.args[1]
                varType = ei.args[2] # Type

                push!(argTypes, varType)
                push!(argParams, var)
                push!(argVars, var)

            elseif ei.head == :block  # anything else should be evaluated again
                # can arise with use of @static inside type decl
                argExpr, argParams, argDefVal, argTypes, argVars, argNames, argFun, amlTypes, amlName, docOrElmType, amlFun = _aml(expr)
            else
                continue
            end

        end # end ifs
    end # endfor

    ########################
    # self closing tags checker
    if  docOrElmType == 10
        # add a field with nothing type
        push!(argNames, "content") # argument ignored for aml
        push!(argTypes, Nothing)
        push!(argFun,missing)
        push!(amlTypes,10)
        push!(argParams, Expr(:kw, :content, nothing))
        push!(argVars, :content)
        push!(argDefVal, nothing)
        push!(argExpr.args,:(content::Nothing))
        # argParams, argDefVal, argTypes, argVars, argNames, amlTypes, amlName, docOrElmType = _aml(argExpr)
    end

    ########################
    # aml::Node adder
    push!(argExpr.args,:(aml::Union{Document,Node}))

    return argExpr, argParams, argDefVal, argTypes, argVars, argNames, argFun, amlTypes, amlName, docOrElmType, amlFun
end

################################################################
# Literal Macros:
# html document
macro hd_str(s)
    docOrElmType = -1
    return docOrElmType, s
end

# xml document
macro xd_str(s)
    docOrElmType = -2
    return docOrElmType, s
end

# self-closing
macro sc_str(s)
    docOrElmType= 10
    return docOrElmType, s
end

# attribute
macro a_str(s)
    elmType = 2
    return elmType, s
end

# namespace
macro ns_str(s)
    elmType = 18
    return elmType, s
end
################################################################
# I/O
import EzXML: parsexml, parsehtml, readxml, readhtml
# from EzXML
funs = [:parsexml, :parsehtml, :readxml, :readhtml]
for fun in funs

    @eval begin

        # doc
        @doc (@doc $(fun)) $(fun)

        # exporting
        export $(fun)
    end
end

end


text:
using Circuitscape
using Test
import Circuitscape: compute_single, compute_cholmod, compute_parallel,
                     runtests

# Unit tests for internals
@testset "Unit tests" begin
    include("internal.jl")
end

for f in (compute, compute_single, compute_cholmod, compute_parallel)
    runtests(f)
end


text:
import MixedModels: AbstractReMat

"""
    ExpandedReMat

Represent a random effects matrix that can't be represented in the compact
format of the standard ReMat.  The z field here is the full Z' matrix, and
linear algebra methods delegate to the (transposed?) matrix.

"""

immutable ExpandedReMat{T<:AbstractFloat,V,R} <: AbstractReMat
    f::Union{NullableCategoricalVector{V,R},CategoricalVector{V,R},PooledDataVector{V,R}}
    z::AbstractMatrix{T}
    fnm::Symbol
    cnms::Vector
end

ExpandedReMat(A::AbstractReMat) = ExpandedReMat(A.f, full(A), A.fnm, A.cnms)

Base.size(A::ExpandedReMat) = size(A.z)
Base.size(A::ExpandedReMat, i::Integer) = size(A.z, i)

Base.full(r::ExpandedReMat) = isa(r.z, AbstractSparseMatrix) ? full(r.z) : r.z
Base.sparse(r::ExpandedReMat) = isa(r.z, AbstractSparseMatrix) ? r.z : sparse(r.z)


text:
using CVCore
using CVVideoIO
using Base.Test

@testset "video capture" begin
    cap = VideoCapture(0)
    @test isOpened(cap)
    release(cap)
    @test !isOpened(cap)
end


text:
@testset "Betweenness" begin
    # self loops
    s2 = DiGraph(3)
    add_edge!(s2, 1, 2); add_edge!(s2, 2, 3); add_edge!(s2, 3, 3)
    s1 = Graph(s2)
    g3 = PathGraph(5)

    gint = loadgraph(joinpath(testdir, "testdata", "graph-50-500.jgz"), "graph-50-500")

    c = vec(readcsv(joinpath(testdir, "testdata", "graph-50-500-bc.txt")))
    for g in testdigraphs(gint)
        z  = @inferred(betweenness_centrality(g))
        zp = @inferred(parallel_betweenness_centrality(g))
        @test all(isapprox.(z, zp))
        @test map(Float32, z)  == map(Float32, c)

        y  = @inferred(betweenness_centrality(g, endpoints=true, normalize=false))
        yp = parallel_betweenness_centrality(g, endpoints=true, normalize=false)
        @test all(isapprox(y, yp))
        @test round.(y[1:3], 4) ==
            round.([122.10760591498584, 159.0072453120582, 176.39547945994505], 4)

        x  = @inferred(betweenness_centrality(g, 3))
        xp = parallel_betweenness_centrality(g, 3)

        x2  = @inferred(betweenness_centrality(g, collect(1:20)))
        xp2 = parallel_betweenness_centrality(g, collect(1:20))

        @test length(x) == 50
        @test length(xp) == 50
        @test length(x2) == 50
        @test length(xp2) == 50
    end

    @test @inferred(betweenness_centrality(s1)) == [0, 1, 0]
    @test parallel_betweenness_centrality(s1)   == [0, 1, 0]

    @test @inferred(betweenness_centrality(s2)) == [0, 0.5, 0]
    @test parallel_betweenness_centrality(s2)   == [0, 0.5, 0]

    g = Graph(2)
    add_edge!(g, 1, 2)
    z  = @inferred(betweenness_centrality(g; normalize=true))
    zp = parallel_betweenness_centrality(g; normalize=true)
    @test all(isapprox(z, zp))
    @test z[1] == z[2] == 0.0
    z2  = @inferred(betweenness_centrality(g, vertices(g)))
    zp2 = parallel_betweenness_centrality(g, vertices(g))
    @test all(isapprox(z2, zp2))
    z3  = @inferred(betweenness_centrality(g, [vertices(g);]))
    zp3 = parallel_betweenness_centrality(g, [vertices(g);])
    @test all(isapprox(z3, zp3))

    @test z == z2 == z3

    z  = @inferred(betweenness_centrality(g3; normalize=false))
    zp = parallel_betweenness_centrality(g3; normalize=false)
    @test all(isapprox(z, zp))
    @test z[1] == z[5] == 0.0

    # Weighted Graph tests
    g=Graph(6)
    add_edge!(g,1,2)
    add_edge!(g,2,3)
    add_edge!(g,3,4)
    add_edge!(g,2,5)
    add_edge!(g,5,6)
    add_edge!(g,5,4)

    distmx = [ 0.0  2.0  0.0  0.0  0.0  0.0;
               2.0  0.0  4.2  0.0  1.2  0.0;
               0.0  4.2  0.0  5.5  0.0  0.0;
               0.0  0.0  5.5  0.0  0.9  0.0;
               0.0  1.2  0.0  0.9  0.0  0.6;
               0.0  0.0  0.0  0.0  0.6  0.0;]

    @test isapprox(betweenness_centrality(g,vertices(g), distmx; normalize=false), [0.0,6.0,0.0,0.0,6.0,0.0])
    @test isapprox(betweenness_centrality(g,vertices(g), distmx; normalize=false, endpoints=true), [5.0,11.0,5.0,5.0,11.0,5.0])
    @test isapprox(betweenness_centrality(g,vertices(g), distmx; normalize=true), [0.0, 0.6000000000000001, 0.0, 0.0, 0.6000000000000001, 0.0])
    @test isapprox(betweenness_centrality(g,vertices(g), distmx; normalize=true, endpoints=true), [0.5,1.1,0.5,0.5,1.1,0.5])

    adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
    a2 = DiGraph(adjmx2)
    for g in testdigraphs(a2)
      distmx2 = [Inf 2.0 Inf; 3.2 Inf 4.2; 5.5 6.1 Inf]
      c2 = [0.24390243902439027,0.27027027027027023,0.1724137931034483]
      @test isapprox(betweenness_centrality(g,vertices(g), distmx2; normalize=false), [0.0,1.0,0.0])
      @test isapprox(betweenness_centrality(g,vertices(g), distmx2; normalize=false, endpoints=true), [4.0,5.0,4.0])
      @test isapprox(betweenness_centrality(g,vertices(g), distmx2; normalize=true), [0.0,0.5,0.0])
      @test isapprox(betweenness_centrality(g,vertices(g), distmx2; normalize=true, endpoints=true), [2.0,2.5,2.0])
    end
end


text:
#=
# Maximum Variance Unfolding (MVU)
This example uses `Chordal.jl` to perfomr the MVU dimensionality reduction
technique on a toy dataset.
=#

#-
using Chordal
using JuMP, Hypatia
using SparseArrays, LinearAlgebra, Random
using Plots

#=
Some helper functions are below. This section can be skipped but is included
for completeness.
=#
#-
## Generate swissroll dataset
function swissroll(n)
    ## NOTE: these could be randomized too
    theta = 1.5 .* pi .+ 2.5 .* pi .* LinRange(0, 1, n)
    z = [(i % 2)-1/2 for i in 1:n]
    return [theta' .* [cos.(theta'); sin.(theta')]; z'], theta
end

function show_data(x, theta; c=(70,70))
    Mx = maximum(abs.(x[1,:]))
    My = maximum(abs.(x[2,:]))
    Mz = maximum(abs.(x[3,:]))
    MM = max(Mx, My, Mz)
    Mx = max(Mx*1.1, MM*.9)
    My = max(My*1.1, MM*.3)
    Mz = max(Mz*1.1, MM*.3)
    return scatter(x[1,:], x[2,:], x[3,:],
        legend=false,
        zcolor=theta,
        xlims=(-Mx,Mx),
        ylims=(-My, My),
        zlims=(-Mz,Mz),
        markersize=7,
        camera=c,
    )
end

## Constructs Euclidean Distance Matrix given vectors x = [x_1 ... x_n]
function edm(x)
    G = x' * x
    dg = diag(G)
    return -2*G + dg*ones(size(G,1))' + ones(size(G,1))*dg'
end

## Given edge weights E, Eij = d(xi, xj), finds k nearest neighbors
function knn(E, k)
    nn = [Vector{Int}(undef,0) for _ in 1:size(E, 1)]
    for i in 1:size(E, 1)
        nn[i] = partialsortperm(@view(E[:,i]), 1:k+1)
    end
    return nn
end

## from nearest neighbors, constructs the MVU edge list
function knn_edges(nn, clique=true)
    k = length(nn) - 1
    E = CartesianIndex{2}[]

    for i in 1:k+1
        for j in i+1:k+1
            i == j && continue
            !clique && i > 1 && continue

            push!(E, CartesianIndex(nn[i], nn[j]))
        end
    end
    return E
end

## Returns edge list for k nearest neighbors given vectors X
function edges(x, k, clique=true)
    nn = knn(edm(x), k)
    return reduce(vcat, [knn_edges(nn[i], clique) for i in 1:length(nn)])
end

## Constructs V matrix (basis)
function get_V(n)
    return spdiagm(n, n-1, 0 => ones(n-1), -1 => -ones(n-1))
end

## Constructs constrain matrices (as in Vandenberghe and Andersen pg 401)
function get_constraint_matrices(x, k)
    n = size(x, 2)
    D = edm(x)

    E = edges(x, k)
    Ei = [CartesianIndex(e[1], e[1]) for e in E]
    Ej = [CartesianIndex(e[2], e[2]) for e in E]

    b = [D[e] for e in E]
    A = [sparse(
            [e[1], e[2], e[1], e[2]],
            [e[1], e[2], e[2], e[1]],
            [1, 1, -1, -1], n, n
         ) for e in E]
    V = get_V(n)
    VA = [V'*Ai*V for Ai in A]

    return V, VA, b
end

## Returns a SparseMatrixCSC with values of Zref
##   Zref is a SparseAxisArray, which functions like a dictionary
##   - Requires optimizer to have been run
function reconstruct_from_sparse_varref(Zref, n)
    Zcv = value.(Zref)
    Z_uncomp = spzeros(Float64, n, n)
    for ind in eachindex(Zcv)
        Z_uncomp[ind[1], ind[2]] = Zcv[ind]
    end
    return Z_uncomp
end

#=
## Problem Setup
First, we construct a "swiss roll" test dataset, which is 3d but has intrinsic
dimension 2
=#
#-
## Construct Dataset
n, k = 80, 2
x, theta = swissroll(n)
show_data(x, theta)

#=
Next, we setup the MVU SDP
=#

#-
## Construct problem data
V, VA, b = get_constraint_matrices(x, k)
C = V'*V;

#=
## Chordal Decomposition and Solve
=#
## Chordal setup
_, _, _, cliques = get_selectors(sparsity_pattern([VA..., C]); verbose=false, ret_cliques=true, perm=nothing)
sp = Chordal.get_sparsity_pattern_from_cliques(cliques)
sp_inds = zip(findnz(sp)[1:2]...)

## -- Construct Model --
model = Model(Hypatia.Optimizer)
Z = @variable(model, Z[i=1:n-1, j=1:n-1; (i,j) in sp_inds])

## Objective
@objective(model, Max, 1/length(sp_inds) * Chordal.tr(C, Z))

## Distrance-preserving equality constraints
for i in 1:length(b)
    @constraint(model, Chordal.tr(VA[i], Z) == b[i])
end

## Chordal decomposition constraints (replace big PSD constraints with many small ones)
for cl in cliques
    len_cl = length(cl)
    Zp = Matrix{VariableRef}(undef, len_cl, len_cl)
    for i in 1:len_cl, j in 1:len_cl
        Zp[i,j] = Z[cl[i], cl[j]]
    end
    @constraint(model, Zp in PSDCone())
end

## Optimize
optimize!(model)
solution_summary(model)
#=
## Solution Reconstruction
=#
## Reconstruct solution with min rank PSD completion
Z_ = Chordal.reconstruct_from_sparse_varref(Z, n-1)
Z_ = 0.5*(Z_ + Z_') # Deals with numerical error
Z_factor = Chordal.minrank_completion(Z_)
W, sv, _ = svd(V*Z_factor)
yc = Diagonal(sv[1:3])*W[:, 1:3]'
show_data(yc, theta; c =(40,70))


text:
# calculate P(z|t)
# P(z|t) = P(z,t) / P(t)

module CondZtBar
    using DataFrames, CSV, Plots, FreqTables
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # pivot
        freqs = freqtable(df_new, :lidar, :hour)
        probs = freqs / length(df_org.lidar)

        # sum
        keys = names(probs, 2)
        probs_t = [sum(probs[begin:end, Name(key)]) for key in keys]

        # P(z|t)
        cond_z_t = probs / probs_t[1] # P(t=0)
        z = names(probs, 1)
        println(cond_z_t)

        # bar plot
        bar(cond_z_t[begin:end, Name(6)], xticks=z, alpha=0.3, color=:blue, label="Hour:6")
        bar!(cond_z_t[begin:end, Name(14)], xticks=z, alpha=0.3, color=:red, label="Hour:14")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/cond_z_t_bar/cond_z_t.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# File to check and compare qdldl routine against built-in cholesky routine
include("qdldl.jl")

using QDLDL, Base.Test

nDense = 1000
nSparse = 1000
numProblems = nDense + nSparse

rng = Random.MersenneTwister(131123)

function generateSparsePosDefMatrix(rng,n::Int64,density::Float64)
      X = sprand(rng,n,n,density)
      X = 0.5(X+X')
      # make random matrix pos def
      X = X+(2*n+1)*speye(n);
      return X
  end

function generatePosDefMatrix(rng,n::Int64)
      X = rand(rng,n,n)
      Q, R = qr(X)
      eigs = rand(rng,n).*(5 .-0.1) .+ 0.1
      X = Q*Matrix(Diagonal(eigs))*Q'
      X = 0.5*(X+X')
      return X
  end

solveTimesQDLDL = zeros(numProblems)
solveTimesCHOL = zeros(numProblems)

@testset "QDLDL routine for random problems" begin

  for iii=1:numProblems

    dim = rand(rng,50:1:200);

    if iii <=nSparse
      density=rand(rng,0.1:0.1:0.4)
      A = generateSparsePosDefMatrix(rng,dim,density)
    else
      A = sparse(generatePosDefMatrix(rng,dim))
    end
    xtrue = rand(rng,dim)
    btrue = A*xtrue

    # solve with QDLDL
    tic()
    F = QDLDL.qdldl(A)
    x_qdldl = solve(F,btrue)
    solveTimesQDLDL[iii] = toq()

    # solve with cholesky
    tic()
    W = cholfact(A)
    x_chol = W\btrue
    solveTimesCHOL[iii] = toq()

    # check results
    @test norm(xtrue-x_qdldl,Inf) < 1e-9
    @test norm(x_chol-x_qdldl,Inf) < 1e-9

  end
end

# print mean solve times for both routines
println("\n\n$(nSparse) SPARSE PROBLEMS\n"*"-"^30*"\nMean solve time QLDL routine: \t$(mean(solveTimesQDLDL[1:nSparse])*1000)ms,\nMean solve Cholesky routine: \t$(mean(solveTimesCHOL[1:nSparse])*1000)ms.\n"*"-"^30)
println("\n$(nDense) DENSE PROBLEMS\n"*"-"^30*"\nMean solve time QLDL routine: \t$(mean(solveTimesQDLDL[nSparse+1:end])*1000)ms,\nMean solve Cholesky routine: \t$(mean(solveTimesCHOL[nSparse+1:end])*1000)ms.\n"*"-"^30)


text:
using SafeTestsets
using Test

@safetestset vas_sample = "VAS samples with direct" begin
using Random: MersenneTwister
using Fleck: VectorAdditionModel, MarkovDirect, vas_initial, send, vas_input
    rng = MersenneTwister(2930472)
    vas = VectorAdditionModel(transitions, rates)
    sampler = MarkovDirect()

    input_process = vas_initial(vas, [1, 1, 0])
    for i in 1:10
        y = send(vas, input_process)
        sample = send(sampler, y)
        input_process = vas_input(vas, sample)
    end
end


text:
using Statistics
N = 10_000         # number of nodes
m = 30             # number of active links per active node
mu = 0.1            # probability I -> S
tmax = 100         # maximum duration of the simulation
Ni = Int(0.1 * N)     # initial number of infects
ActMin = 0.001      # minimum activation probability
ActMax = 1          # maximum activation probability
gamma = -2.22       # exponent of the power law

AgentAct = rand(N, 1);
AgentAct = ((ActMax.^(gamma + 1) - ActMin.^(gamma + 1)) .* AgentAct .+
            ActMin.^(gamma + 1)).^(1.0 / (gamma + 1))
AvgDegree = 2 * m * mean(AgentAct)

NTrials = 5         # number of trials for each value of lambda

tolerance = 1e-15   # termination tolerance
TWindow = 1000      # time window for evaluating steady state
TMinConvergence = 500

# health state of agents S=0 I=1
global AgentState = zeros(N)
AgentState[1:Ni] .= 1.0

# diffferent lambda
# lambda=0.1:0.1:1;
# IAvgLambda=zeros(length(lambda))
# beta=zeros(length(lambda))

lambda = 0.2
termination = 1.0
t = 1;
ITrend = zeros(0)

 for t = 1:tmax
    AdjMatrix = zeros(N, N)
    ActiveIndicesBool = [rand() < agentAct for agentAct in AgentAct]
    indx = [k[1] for k in findall(ActiveIndicesBool)]
    global AgentState=AgentState
    for i in indx
        pick_from = setdiff(1:N, i)
        ContactedAgents = rand(pick_from, m)
        AdjMatrix[i, ContactedAgents] .= 1
        AdjMatrix[ContactedAgents, i] .= 1 # connections are bilateral!
    end

    global AgentStateNext = AgentState

    InfectedIndicesBool = [agentState == 1 for agentState in AgentState]
    InfectedIndices = [k[1] for k in findall(InfectedIndicesBool)]

    for j in InfectedIndices
        ConnectedSusceptibles = Bool.(AdjMatrix[j, :]) .& (.~Bool.(AgentState)) # find individuals connected to an infected that are susceptible
        ConnSuscIndices = findall(ConnectedSusceptibles) # find their indices
        AgentStateNext[ConnSuscIndices] = [rand() < lambda for v in ConnSuscIndices] # infection process S->I on the susceptible and connected;
    end

    InfectedIndices = [k[1] for k in findall(Bool.(AgentStateNext)) ]# find Infected again and
    AgentStateNext[InfectedIndices] = [~(rand() < mu) for v in InfectedIndices] # process I->S (it is independent from the network connections)
    AgentState = AgentStateNext

    push!(ITrend, Int(sum(AgentState)))
    if t > (TMinConvergence + TWindow)
        termination = (abs(mean(ITrend[t - TWindow:t - 1]) - mean(ITrend[t - TWindow + 1:t])) < tolerance);
        if termination
            break
        end
    end

end


text:
mutable struct Reactions
  IDs
  names
  S
  reversible
  grRules
  subsystems
  assGenes
  inputColours
  weights
  leadingGenes
  outputColours
  outputFrequency
  combinedOutput

  function Reactions(IDs,
	  				 names,
	  				 S,
					 reversible,
					 grRules,
					 subsystems,
					 assGenes,
					 inputColours,
					 weights,
					 leadingGenes,
					 outputColours,
					 outputFrequency,
					 combinedOutput)
    new(IDs,
		names,
		S,
		reversible,
		grRules,
		subsystems,
		assGenes,
		inputColours,
		weights,
  		leadingGenes,
		outputColours,
		outputFrequency,
		combinedOutput)
  end

  function Reactions(IDs,
	  				 names,
	  				 S,
					 reversible,
					 grRules,
					 subsystems,
					 assGenes)
    Reactions(IDs,
		names,
		S,
		reversible,
		grRules,
		subsystems,
		assGenes,
		[],
		[],
  		[],
		[],
		[],
		[])
  end
end

function nReactions(reactions::Reactions)
	return length(reactions.IDs)
end

function getAssGenes(rules)
	assGenes = repeat([[]], length(rules))
	for i in 1:length(rules)
	  assGenes[i] = unique(parse.(Int64, collect(m.match for m in eachmatch(r"\d+", rules[i]))))
	end

	return assGenes
end

"""
Determine colours and weights for reactions
"""
function colour!(reactions::Reactions, genes::Genes, pThresh, alpha, precision)
	nReac = nReactions(reactions)
	colours = zeros(nReac)
	weights = fill(weight(0, pThresh, alpha, precision), nReac)
	leadingGenes = repeat([zero(Int64)], nReac)

	for reacInd in 1:nReac
		currAssGenes = reactions.assGenes[reacInd]
		if !isempty(currAssGenes)
			currColours = genes.colours[currAssGenes]
			currWeights = genes.weights[currAssGenes]
			if (1 in currColours) & (-1 in currColours)
				weights[reacInd] = weight(0.5, pThresh, alpha, precision)
			else
				colours[reacInd] = sign(sum(currColours))
				(weights[reacInd], leadGeneInd) = findmax(currWeights)
				if weights[reacInd] > 0
					leadingGenes[reacInd] = currAssGenes[leadGeneInd]
				end
			end
		end
	end

	reactions.inputColours = colours
	reactions.weights = weights
	reactions.leadingGenes = leadingGenes
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
function PeriodicOrbit(sys::DiffSys; ps=(2,0.0), eps=10^(-10), transient=10^4, steps = 2^12)
    n_dim = sys.n_dim
    x0 = rand(n_dim)
    dt = 0.01
    next_step(x, dt) = RungeKutta(sys.func, dt, x)
    for _ in 1:dt:transient
        x0 = next_step(x0, dt)
    end
    x_prev, x_now = x0, next_step(x0, dt)
    idx, section = ps
    while true
        if (x_prev[idx]>section) & (x_now[idx]<section)
            if dt < eps
                break
            else
                dt *= 0.5
                x_now = next_step(x_prev, dt)
                continue
            end
        end
        x_prev = x_now
        x_now = next_step(x_prev, dt)
    end
    x_dest = x_now
    @info "Poincare section for the system crossing x[$idx]=$section from positive to negative\nperiodic section" x_dest
    t, dt = 0.0, 0.01
    x_prev = x_now
    x_now = next_step(x_prev, dt)
    while true
        if (x_prev[idx]>section) & (x_now[idx]<section)
            if dt < eps
                t += dt
                break
            else
                dt *= 0.5
                x_now = next_step(x_prev, dt)
                continue
            end
        end
        x_prev = x_now
        x_now = next_step(x_prev, dt)
        t += dt
    end
    period = t
    dt = period / steps
    ts = [0:dt:period;]
    x = x_dest
    orbits = zeros(steps+1, n_dim)
    orbits[1,:] = x
    for i in 1:steps
        x = next_step(x, dt)
        orbits[i+1,:] = x
    end
    distance = norm(orbits[begin,:]-orbits[end,:])
    @info "Periodicity checked: distance between start and end of orbit = $distance"
    return period, ts, orbits
end

PeriodicOrbit(sys::SpecialDiffSys; ps=(2,0.0), eps=10^(-10), transient=10^4, steps = 2^12) = PeriodicOrbit(DiffSys(sys), ps=ps, eps=eps, transient=transient, steps = steps)

function PhaseSensitivity(sys::SpecialDiffSys, period::Float64, periodic_orbits::AbstractArray; eps=10^(-12))
    func(x) = vec_field(sys, x)
    jac(x) = jacobian(sys, x)
    ω = 2π/period
    steps, n_dim = size(periodic_orbits)
    dt = period/(steps-1)
    periodic_orbits_reverse = reverse(periodic_orbits, dims=1)
    adjoint_func(z, x) = -jac(x)'*z
    Z_now = zeros(steps, n_dim)
    z = randn(n_dim)
    # contour integral
    for i in 1:steps
        Z_now[i,:] = z
        z -= dt*adjoint_func(z, periodic_orbits_reverse[i,:])
    end
    # normalization
    for i in 1:steps
        @views coeff = ω / dot(Z_now[i,:], func(periodic_orbits_reverse[i, :]))
        Z_now[i,:] *= coeff
    end
    Z_prev = copy(Z_now)
    while true
        Z_now[1,:] = Z_prev[end,:]
        # contour integral
        for i in 1:steps
            Z_now[i,:] = z
            z -= dt*adjoint_func(z, periodic_orbits_reverse[i,:])
        end
        # normalization
        for i in 1:steps
            @views coeff = ω / dot(Z_now[i,:], func(periodic_orbits_reverse[i, :]))
            Z_now[i,:] *= coeff
        end
        # checking convergence
        L∞ = maximum(Z_now-Z_prev)
        if L∞ < eps
            @info "L^∞  norm between phase sensitivity function = $L∞ \n well converged"
            break
        else
            # @info "L^∞ norm between phase sensitivity function = $L∞ \n one more contour integral"
        end
        Z_prev = copy(Z_now)
    end
    return reverse(Z_now, dims=1)
end

function PhaseSensitivity(sys::SpecialDiffSys; steps=2^12)
    period, _, periodic_orbits = PeriodicOrbit(sys, steps=steps)
    @info "Period of periodic orbit" period
    return PhaseSensitivity(sys, period, periodic_orbits)
end

function PhaseCoupling(phasesensitivity::AbstractArray, coupling::Function, periodic_orbit1::AbstractArray, periodic_orbit2::AbstractArray)
    steps = size(phasesensitivity, 1)
    step1, step2 = size(periodic_orbit1, 1), size(periodic_orbit1, 1)
    @assert step1==steps "Time step of phasesensitivity function ($steps) differs from first periodic orbit ($step1)"
    @assert step2==steps "Time step of phasesensitivity function ($steps) differs from first periodic orbit ($step2)"
    dθ = 2π/(steps-1)
    shift_phase(orbits, j) = begin
        tmp = circshift(orbits[begin:end-1,:], -j)
        return vcat(tmp, tmp[1,:]')
    end
    Γ_arr = zeros(steps)
    trapezoidal_weight = ones(steps); trapezoidal_weight[[begin,end]] .= 0.5
    for i in 1:steps
        shifted_periodic_orbit2 = shift_phase(periodic_orbit2, i)
        @views Gjks = [coupling(periodic_orbit1[l,:], shifted_periodic_orbit2[l,:]) for l in 1:steps]
        @views integrands = [dot(phasesensitivity[l,:], Gjks[l]) for l in 1:steps]
        Γ_arr[i] = dot(trapezoidal_weight, integrands)*dθ/2π
    end
    return Γ_arr
end

function PhaseCoupling(sys::CoupledDiffSys)
    periodic_orbit_data = [PeriodicOrbit(s) for s in sys.SelfDiffSyss]
    periodic_orbits = [d[3] for d in periodic_orbit_data]
    phasesensitivitys = [PhaseSensitivity(sys.OriginalSpecialDiffSyss[i], periodic_orbit_data[i][1], periodic_orbit_data[i][3]) for i in 1:sys.n_body]
    phasecoupling_dict = Dict()
    for (key, coupling) in sys.CouplingFuncDict
        sender, reciever = key
        Γ_arr = PhaseCoupling(phasesensitivitys[sender], coupling, periodic_orbits[sender], periodic_orbits[reciever])
        phasecoupling_dict[key] = Γ_arr
    end
    return phasecoupling_dict
end