text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# BGZF* bgzf_dopen(int fd, const char *mode);
@defun bgzf_dopen(fd::Cint, mode::Cstring)::Ptr{BGZF}

# BGZF* bgzf_open(const char* path, const char *mode);
@defun bgzf_open(path::Cstring, mode::Cstring)::Ptr{BGZF}

# BGZF* bgzf_hopen(struct hFILE *fp, const char *mode);
@defun bgzf_hopen(fp::Ptr{hFILE}, mode::Cstring)::Ptr{BGZF}

# int bgzf_close(BGZF *fp);
@defun bgzf_close(fp::Ptr{BGZF})::Cint

# ssize_t bgzf_read(BGZF *fp, void *data, size_t length) HTS_RESULT_USED;
@defun bgzf_read(fp::Ptr{BGZF}, data::Ptr{Cvoid}, length::Csize_t)::Cssize_t

# ssize_t bgzf_write(BGZF *fp, const void *data, size_t length) HTS_RESULT_USED;
@defun bgzf_write(fp::Ptr{BGZF}, data::Ptr{Cvoid}, length::Csize_t)::Cssize_t

# ssize_t bgzf_block_write(BGZF *fp, const void *data, size_t length);
@defun bgzf_block_write(fp::Ptr{BGZF}, data::Ptr{Cvoid}, length::Csize_t)::Cssize_t

# int bgzf_peek(BGZF *fp);
@defun bgzf_peek(fp::Ptr{BGZF})::Cint

# ssize_t bgzf_raw_read(BGZF *fp, void *data, size_t length) HTS_RESULT_USED;
@defun bgzf_raw_read(fp::Ptr{BGZF}, data::Ptr{Cvoid}, length::Csize_t)::Cssize_t

# ssize_t bgzf_raw_write(BGZF *fp, const void *data, size_t length) HTS_RESULT_USED;
@defun bgzf_raw_write(fp::Ptr{BGZF}, data::Ptr{Cvoid}, length::Csize_t)::Cssize_t

# int bgzf_flush(BGZF *fp) HTS_RESULT_USED;
@defun bgzf_flush(fp::Ptr{BGZF})::Cint

# int64_t bgzf_seek(BGZF *fp, int64_t pos, int whence) HTS_RESULT_USED;
@defun bgzf_seek(fp::Ptr{BGZF}, pos::Int64, whence::Cint)::Int64

# int bgzf_check_EOF(BGZF *fp);
@defun bgzf_check_EOF(fp::Ptr{BGZF})::Cint

# int bgzf_compression(BGZF *fp);
@defun bgzf_compression(fp::Ptr{BGZF})::Cint

# int bgzf_is_bgzf(const char *fn) HTS_DEPRECATED("Use bgzf_compression() or hts_detect_format() instead");
@defun bgzf_is_bgzf(fn::Cstring)::Cint

# void bgzf_set_cache_size(BGZF *fp, int size);
@defun bgzf_set_cache_size(fp::Ptr{BGZF}, size::Cint)::Cvoid

# int bgzf_flush_try(BGZF *fp, ssize_t size) HTS_RESULT_USED;
@defun bgzf_flush_try(fp::Ptr{BGZF}, size::Cssize_t)::Cint

# int bgzf_getc(BGZF *fp);
@defun bgzf_getc(fp::Ptr{BGZF})::Cint

# int bgzf_getline(BGZF *fp, int delim, struct kstring_t *str);
@defun bgzf_getline(fp::Ptr{BGZF}, delim::Cint, str::Ptr{kstring_t})::Cint

# int bgzf_read_block(BGZF *fp) HTS_RESULT_USED;
@defun bgzf_read_block(fp::Ptr{BGZF})::Cint

# int bgzf_thread_pool(BGZF *fp, struct hts_tpool *pool, int qsize);
@defun bgzf_thread_pool(fp::Ptr{BGZF}, pool::Ptr{hts_tpool}, qsize::Cint)::Cint

# int bgzf_mt(BGZF *fp, int n_threads, int n_sub_blks);
@defun bgzf_mt(fp::Ptr{BGZF}, n_threads::Cint, n_sub_blks::Cint)::Cint

# int bgzf_compress(void *dst, size_t *dlen, const void *src, size_t slen, int level);
@defun bgzf_compress(dst::Ptr{Cvoid}, dlen::Ptr{Csize_t}, src::Ptr{Cvoid}, slen::Csize_t, level::Cint)::Cint

# int bgzf_useek(BGZF *fp, off_t uoffset, int where) HTS_RESULT_USED;
@defun bgzf_useek(fp::Ptr{BGZF}, uoffset::Coff_t, where::Cint)::Cint

# off_t bgzf_utell(BGZF *fp);
@defun bgzf_utell(fp::Ptr{BGZF})::Coff_t

# int bgzf_index_build_init(BGZF *fp);
@defun bgzf_index_build_init(fp::Ptr{BGZF})::Cint

# int bgzf_index_load(BGZF *fp, const char *bname, const char *suffix) HTS_RESULT_USED;
@defun bgzf_index_load(fp::Ptr{BGZF}, bname::Cstring, suffix::Cstring)::Cint

# int bgzf_index_load_hfile(BGZF *fp, struct hFILE *idx, const char *name) HTS_RESULT_USED;
@defun bgzf_index_load_hfile(fp::Ptr{BGZF}, idx::Ptr{hFILE}, name::Cstring)::Cint

# int bgzf_index_dump(BGZF *fp, const char *bname, const char *suffix) HTS_RESULT_USED;
@defun bgzf_index_dump(fp::Ptr{BGZF}, bname::Cstring, suffix::Cstring)::Cint

# int bgzf_index_dump_hfile(BGZF *fp, struct hFILE *idx, const char *name) HTS_RESULT_USED;
@defun bgzf_index_dump_hfile(fp::Ptr{BGZF}, idx::Ptr{hFILE}, name::Cstring)::Cint


text:
immutable SharedDict{K,V} <: Associative{K,V}
    slots::SharedArray{UInt8,1}
    keys::SharedArray{K,1}
    vals::SharedArray{V,1}
    count::Int
    idxfloor::Int  # an in
    pids::AbstractVector{Int}
end

function share{T}(a::AbstractArray{T}; kwargs...)
    sh = SharedArray{T}(size(a); kwargs...)
    @simd for i=1:length(a)
        sh.s[i] = a[i]
    end
    return sh
end
share(A::SharedArray, pids::AbstractVector{Int}) = (pids==A.pids ? A : share(sdata(A), pids))

share(D::Dict) = share(D::Dict, workers())
function share(D::Union{Dict,SharedDict}, pids::AbstractVector{Int})
    isa(D, SharedDict) && (pids == D.pids) && (return D)

    slots = share(D.slots, pids=pids)
    keys = share(D.keys, pids=pids)
    vals = share(D.vals, pids=pids)

    SharedDict(slots, keys, vals, D.count, D.idxfloor, pids)
end

hashindex(key, sz) = ((hash(key)%Int) & (sz-1)) + 1

isslotempty(h::SharedDict, i::Int) = h.slots[i] == 0x0
isslotfilled(h::SharedDict, i::Int) = h.slots[i] == 0x1
isslotmissing(h::SharedDict, i::Int) = h.slots[i] == 0x2

# get the index where a key is stored, or -1 if not present
function ht_keyindex{K,V}(h::SharedDict{K,V}, key)
    sz = length(h.keys)
    iter = 0
    maxprobe = max(16, sz>>6)
    index = hashindex(key, sz)
    keys = h.keys

    while true
        if isslotempty(h,index)
            break
        end
        if !isslotmissing(h,index) && isequal(key,keys[index])
            return index
        end

        index = (index & (sz-1)) + 1
        iter+=1
        iter > maxprobe && break
    end

    return -1
end

function getindex{K,V}(h::SharedDict{K,V}, key)
    index = ht_keyindex(h, key)
    return (index<0) ? throw(KeyError(key)) : h.vals[index]::V
end

haskey(h::SharedDict, key) = (ht_keyindex(h, key) >= 0)
in{T<:SharedDict}(key, v::KeyIterator{T}) = (ht_keyindex(v.dict, key) >= 0)

function skip_deleted(h::SharedDict, i)
    L = length(h.slots)
    while i<=L && !isslotfilled(h,i)
        i += 1
    end
    return i
end

function start(t::SharedDict)
    i = skip_deleted(t, t.idxfloor)
    t.idxfloor = i
    return i
end
done(t::SharedDict, i) = i > length(t.vals)
next{K,V}(t::SharedDict{K,V}, i) = (Pair{K,V}(t.keys[i],t.vals[i]), skip_deleted(t,i+1))

isempty(t::SharedDict) = (t.count == 0)
length(t::SharedDict) = t.count

next{T<:SharedDict}(v::KeyIterator{T}, i) = (v.dict.keys[i], skip_deleted(v.dict,i+1))
next{T<:SharedDict}(v::ValueIterator{T}, i) = (v.dict.vals[i], skip_deleted(v.dict,i+1))

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    WeightedCrossValidation(weigthing, folding; lambda=1.0, loss=Dict())

An error estimation method which samples are weighted with
`weighting` method and split into folds with `folding` method.
Weights are raised to `lambda` power in `[0,1]`. Optionally,
specify `loss` function from `LossFunctions.jl` for some of
the variables.

## References

* Sugiyama et al. 2006. [Importance-weighted cross-validation for
  covariate shift](https://link.springer.com/chapter/10.1007/11861898_36)
* Sugiyama et al. 2007. [Covariate shift adaptation by importance weighted
  cross validation](http://www.jmlr.org/papers/volume8/sugiyama07a/sugiyama07a.pdf)
"""
struct WeightedCrossValidation{W<:WeightingMethod,
                               F<:FoldingMethod,
                               T<:Real} <: ErrorEstimationMethod
  weighting::W
  folding::F
  lambda::T
  loss::Dict{Symbol,SupervisedLoss}

  function WeightedCrossValidation{W,F,T}(weighting, folding,
                                          lambda, loss) where {W,F,T}
    @assert 0 ≤ lambda ≤ 1 "lambda must lie in [0,1]"
    new(weighting, folding, lambda, loss)
  end
end

WeightedCrossValidation(weighting::W, folding::F;
                        lambda::T=one(T), loss=Dict()) where {W,F,T} =
  WeightedCrossValidation{W,F,T}(weighting, folding, lambda, loss)

function error(solver, problem, method::WeightedCrossValidation)
  # retrieve problem info
  sdata = _foldable(problem)
  ovars = _outputvars(problem)
  loss  = method.loss
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[var][1])
    end
  end

  # weight all samples
  ws = weight(sdata, method.weighting) .^ method.lambda

  # folds for cross-validation
  fs = folds(sdata, method.folding)

  # error for a fold
  function ε(f)
    # setup and solve sub-problem
    subproblem = _subproblem(problem, f)
    solution   = _solve(subproblem, solver)

    # holdout set and weights
    holdout = _holdout(problem, f)
    weights = view(ws, f[2])

    # loss for each variable
    losses = map(ovars) do var
      y = holdout[var]
      ŷ = solution[var]
      ℒ = value(loss[var], y, ŷ, AggMode.WeightedSum(weights)) / length(y)
      var => ℒ
    end

    Dict(losses)
  end

  # compute error for each fold in parallel
  εs = foldxt(vcat, Map(ε), fs)

  # combine error from different folds
  Dict(var => mean(get.(εs, var, 0)) for var in ovars)
end

# foldable data of the problem
_foldable(p::EstimationProblem) = data(p)
_foldable(p::LearningProblem) = sourcedata(p)

# output variables of the problem
_outputvars(p::EstimationProblem) = name.(variables(p))
_outputvars(p::LearningProblem) = outputvars(task(p))

# subproblem for a given fold
function _subproblem(p::EstimationProblem, f)
  sdat = view(data(p), f[1])
  sdom = view(domain(data(p)), f[2])
  vars = name.(variables(p))
  EstimationProblem(sdat, sdom, vars)
end
function _subproblem(p::LearningProblem, f)
  source = view(sourcedata(p), f[1])
  target = view(sourcedata(p), f[2])
  LearningProblem(source, target, task(p))
end

# solution for a given fold
function _solve(p::EstimationProblem, solver)
  vars = name.(variables(p))
  sol  = solve(p, solver)
  reduce(hcat, [sol[v] for v in vars])
end
_solve(p::LearningProblem, solver) = solve(p, solver)

# holdout set for a given fold
_holdout(p::EstimationProblem, f) = view(data(p), f[2])
_holdout(p::LearningProblem, f) = view(sourcedata(p), f[2])


text:
function SMD_settings(D_ul::Int, D_ll::Int, fnum::Int)
    settings = zeros(Int32, 6)

    ccall((:blb18_cop_settings, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Ptr{Int32},  Int32),
          D_ul,  D_ll, settings, fnum)

    # returns p, q, r, s, lenG,leng
    return tuple(settings...)
end

function SMD_solutions(D_ul::Int, D_ll::Int, fnum::Int)

    x = zeros(D_ul)
    y = zeros(D_ll)

    ccall((:blb18_cop_solutions, bilevelBenchmark),
              Cvoid,
            (Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
            D_ul, D_ll, x, y, fnum)

    return x, y
end

function SMD_ranges(D_ul::Int, D_ll::Int, fnum::Int)
    bounds_ul = zeros(2D_ul)
    bounds_ll = zeros(2D_ll)

    ccall((:blb18_cop_ranges, bilevelBenchmark),
              Cvoid,
            (Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
            D_ul, D_ll, bounds_ul, bounds_ll, fnum)

    return Array(reshape(bounds_ul, D_ul, :)'), Array(reshape(bounds_ll, D_ll, :)')
end

function SMD_leader(x::Array{Float64}, y::Array{Float64}, fnum::Int)
    D_ul = length(x)
    D_ll = length(y)
    F = [0.0]
    G = [0.0]
    lenG = 0

    if fnum > 8
        p, q, r, s, lenG,leng = SMD_settings(D_ul, D_ll, fnum)
        G = zeros(lenG)
    end

    ccall((:blb18_leader_cop, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
        1, D_ul, D_ll, x, y, F, G, fnum)

    if lenG == 0
      return F[1]
    end

    return F[1], G, [0.0]
end

function SMD_follower(x::Array{Float64}, y::Array{Float64}, fnum::Int)
    D_ul = length(x)
    D_ll = length(y)
    f = [0.0]
    g = [0.0]
    leng = 0

    if fnum > 8
        p, q, r, s, lenG,leng = SMD_settings(D_ul, D_ll, fnum)
        g = zeros(leng)
    end

    ccall((:blb18_follower_cop, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
        1, D_ul, D_ll, x, y, f, g, fnum)

    if leng == 0
        return f[1]
    end

    return f[1], g, [0.0]
end

function SMD_leader(x::Array{Float64}, y::Array{Float64}, fnum::Int,p::Int,q::Int,r::Int,s::Int = 0)
    F = [0.0]

    if fnum == 1
        ccall((:SMD1_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 2
        ccall((:SMD2_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 3
        ccall((:SMD3_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 4
        ccall((:SMD4_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 5
        ccall((:SMD5_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 6
        ccall((:SMD6_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, s, x, y, F)
    elseif fnum == 7
        ccall((:SMD7_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 8
        ccall((:SMD8_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    else
        error("Not a valid number")
    end

    return F[1]
end

function SMD_follower(x::Array{Float64},y::Array{Float64},fnum::Int,p::Int,q::Int,r::Int,s::Int = 0)
    D_ul = length(x)
    D_ll = length(y)
    F = [0.0]

    if fnum == 1
        ccall((:SMD1_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 2
        ccall((:SMD2_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 3
        ccall((:SMD3_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 4
        ccall((:SMD4_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 5
        ccall((:SMD5_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 6
        ccall((:SMD6_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, s, x, y, F)
    elseif fnum == 7
        ccall((:SMD7_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    elseif fnum == 8
        ccall((:SMD8_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          p, q, r, x, y, F)
    else
        error("Not a valid number")
    end

    return F[1]
end

function SMD_Ψ(x, D_ll, fnum)
    D_ul = length(x)
    y = zeros(D_ll)

    p, q, r, s, lenG,leng = SMD_settings(D_ul, D_ll, fnum)
    # @show p, q, r, s

    if fnum == 1
        y[q+1:end] = atan.(x[p+1:end])
    elseif fnum == 2
        y[q+1:end] = exp.(x[p+1:end])
    elseif fnum == 3
        y[q+1:end] = atan.(x[p+1:end].^2)
    elseif fnum == 4
        y[q+1:end] = exp.(abs.(x[p+1:end])) .- 1.0
    elseif fnum == 5
        y[1:q] .= 1.0
        y[q+1:end] = sqrt.(abs.(x[p+1:end]))
    elseif fnum == 6
        # q = 1; s = 2;
        y[q+s+1:end] = x[p+1:end]
    elseif fnum == 7
        y[q+1:end] = exp.(x[p+1:end])
    elseif fnum == 8
        xx = x[p+1:end]
        y[q+1:end] = sign.(xx) .* abs.(xx) .^(1/3)
    end

    y
end

function SMD_optimum(fnum)
    if 1 <= fnum <= 9
        return 0.0, 0.0
    elseif fnum == 10
        return 4.0, 3.0
    elseif fnum == 11
        return -1.0, 1.0
    elseif fnum == 12
        return 3.0, 4.0
    else
        error("SMD_optimum(fnum) with fnum = 1,2,...,12")
    end

end


text:
export gadplot_axes

split_gadargs(args...) = (),args
split_gadargs(ax::AxisId,args...) = (ax,),args
split_gadargs(ax1::AxisId,ax2::AxisId,args...) = (ax1,ax2),args
split_gadargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,args...) = (ax1,ax2,ax3),args
split_gadargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,ax4::AxisId,args...) =
  (ax1,ax2,ax3,ax4),args
function split_gadargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,ax4::AxisId,
                      ax5::AxisId,args...)
  error("Plotting data using 5 or more axes is not supported by Gadfly backend.")
end

function gadplot_axes(data,args...;kwds...)
  ax, args = split_gadargs(args...)
  df, axes = asplotable(data,ax...;kwds...)
  if eltype(df.value) <: Complex
    df.value = abs.(df.value)
    @warn("Ignoring phase of complex value")
  end

  gadplot_axes_(df,axes,names(df)[2:end]...;args=args)
end

available_backends[:gadfly] = gadplot_axes
set_backend!(:gadfly)

function gadplot_axes_(df,axes,x;args)
  plot(df,x=x,y=:value,Geom.line,args...)
end

function gadplot_axes_(df,axes,x,y;args)
  (xmin,xmax) = extrema(df[!,x])
  (ymin,ymax) = extrema(df[!,y])
  plot(df,x=x,y=y,color=:value,Geom.rectbin,
       Coord.cartesian(xmin=xmin,xmax=xmax),args...)
end

function gadplot_axes_(df,axes,x,y,z;args)
  xmin,xmax = extrema(df[!,x])
  ymin,ymax = extrema(df[!,y])
  plot(df,x=x,y=y,color=:value,xgroup=z,
       Geom.subplot_grid(Geom.rectbin,
                         Coord.cartesian(xmin=xmin,xmax=xmax,
                                         ymin=ymin,ymax=ymax)),args...)
end

function gadplot_axes_(df,axes,x,y,z,w;args)
  xmin,xmax = extrema(df[!,x])
  ymin,ymax = extrema(df[!,y])
  plot(df,x=x,y=y,color=:value,xgroup=z,ygroup=w,
       Geom.subplot_grid(Geom.rectbin,
                         Coord.cartesian(xmin=xmin,xmax=xmax,
                                         ymin=ymin,ymax=ymax)),args...)
end
function gadplot_axes_(df,axes,args...;kwds...)
  error("Plotting data with $(length(axes)) dims along
        $(length(args)) axes is not supported.")
end


text:
module test_examples
include("example_pointsets.jl")
include("example_domains.jl")

#
# Called by runtest.
#
function test()
    # Call all functions in this module whose names start with "example"
    for symbol in names(@__MODULE__,all=true)
        if findfirst("example_", String(symbol))==1:length("example_")
            println("Testing $(String(symbol))...")
            getfield(@__MODULE__,symbol)()
        end
    end
    true
end

# End of module
end


text:
abstract type CartesianEquilibrium <: AnalyticEquilibrium end
abstract type CartesianPerturbation <: AnalyticPerturbation end

CartesianEquPert = Union{CartesianEquilibrium, CartesianPerturbation}

function X(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    x[1]
end

function Y(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    x[2]
end

function Z(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    x[3]
end

function J(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    one(T)
end

function g₁₁(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    one(T)
end

function g₂₂(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    one(T)
end

function g₃₃(x::AbstractArray{T,1}, equ::CartesianEquPert) where {T <: Number}
    one(T)
end


text:
######### StanVariational Bernoulli example  ###########

using StanVariational, Test

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Use tmpdir across multiple runs to prevent re-compilation
#tmpdir = joinpath(@__DIR__, "tmp")

sm = VariationalModel("variational", bernoulli_model)

rc = stan_variational(sm; data=bernoulli_data)

if success(rc)

  (samples, cnames) = read_variational(sm)

end


text:
# EXAMPLE: gather
# here we have 3 ranks

using MPI

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)

if rank == 0
    data = [7.0,8.0,9.0,10.0,11.0,12.0]
else
    data = Float64[]
end

data = MPI.Scatter(data, 2, 0, comm) #each rank has 2 elements

data_gather = MPI.Gather(data, 2, 0, comm)

#print output
if rank == 0
    @show data_gather
end

MPI.Finalize()


text:
@testset "CoreModel type" begin
    cp = test_LP()
    @test cp isa CoreModel
end


text:
module PPTs

using StatsBase, LinearAlgebra
using ExprRules

export
        PPT,
        PPTNode,
        get_child,
        probabilities,
        probability,
        prune!

"""
    PPT

Probabilistic Prototype Tree.
# Arguments:
- `w_terminal::Float64`: probability of selecting a terminal on initialization
- `w_nonterm::Float64`: probability of selecting a non-terminal on initialization
"""
struct PPT
    w_terminal::Float64
    w_nonterm::Float64
end
PPT(w_terminal::Float64=0.6) = PPT(w_terminal, 1-w_terminal)

struct PPTNode
    ps::Dict{Symbol,Vector{Float64}}
    children::Vector{PPTNode}
end

"""
    PPTNode(pp::PPT, grammar::Grammar)

Node of a PPT.
"""
function PPTNode(pp::PPT, grammar::Grammar)
    ps = Dict(typ=>normalize!([isterminal(grammar, i) ?
                               pp.w_terminal : pp.w_nonterm for i in grammar[typ]], 1)
             for typ in nonterminals(grammar))
    PPTNode(ps, PPTNode[])
end

"""
    nchildren(node::PPTNode)

Returns the number of children of a node.
"""
ExprRules.nchildren(node::PPTNode) = length(node.children)

"""
    probabilities(node::PPTNode, typ::Symbol)

Returns the probability vector of a node.
"""
probabilities(node::PPTNode, typ::Symbol) = node.ps[typ]

"""
    get_child(pp::PPT, node::PPTNode, grammar::Grammar, i::Int)

Returns child node i of a ppt node.  Will construct and initialize it if it doesn't already exist.
"""
function get_child(pp::PPT, node::PPTNode, grammar::Grammar, i::Int)
    if i > length(node.children)
        push!(node.children, PPTNode(pp, grammar))
    end
    node.children[i]
end

"""
    rand(pp::PPT, node::PPTNode, grammar::Grammar, typ::Symbol)

Randomly sample an expression tree from ppt model using parameters pp, grammar 'grammar', and start symbol typ.
"""
function Base.rand(pp::PPT, pptnode::PPTNode, grammar::Grammar, typ::Symbol)
    rules = grammar[typ]
    rule_index = sample(rules, Weights(pptnode.ps[typ]))
    ctypes = child_types(grammar, rule_index)
    node = iseval(grammar, rule_index) ?
        RuleNode(rule_index, Core.eval(grammar, rule_index), Vector{RuleNode}(undef,length(ctypes))) :
        RuleNode(rule_index, Vector{RuleNode}(undef,length(ctypes)))

    for (i,typ) in enumerate(ctypes)
        node.children[i] = rand(pp, get_child(pp, pptnode, grammar, i), grammar, typ)
    end
    node
end

"""
    probability(pp::PPT, node::PPTNode, grammar::Grammar, expr::RuleNode)

Compute the probability of an expression tree expr given the model ppt using parameters pp and grammar 'grammar'.
"""
function probability(pp::PPT, node::PPTNode, grammar::Grammar, expr::RuleNode)
    typ = return_type(grammar, expr)
    i = something(findfirst(isequal(expr.ind),grammar[typ]), 0)
    retval = node.ps[typ][i]
    for (i,c) in enumerate(expr.children)
        retval *= probability(pp, get_child(pp, node, grammar, i), grammar, c)
    end
    retval
end

"""
    prune!(node::PPTNode, grammar::Grammar, p_threshold::Float64)

Prune the ppt of decisions with probability less than p_threshold.
"""
function prune!(node::PPTNode, grammar::Grammar, p_threshold::Float64)
    kmax, pmax = :None, 0.0
    for (k, p) in node.ps
        pmax′ = maximum(p)
        if pmax′ > pmax
            kmax, pmax = k, pmax′
        end
    end
    if pmax > p_threshold
        i = indmax(node.ps[kmax])
        if isterminal(grammar, i)
            empty!(node.children)
        else
            max_arity_for_rule = maximum(nchildren(grammar, r) for
                                         r in grammar[kmax])
            while length(node.children) > max_arity_for_rule
                pop!(node.children)
            end
        end
    end
    return node
end

end #module


text:
""" Data structures for graphs, based on C-sets.

Provides the category theorist's four basic kinds of graphs: (directed) graphs,
symmetric graphs, reflexive graphs, and symmetric reflexive graphs. Also defines
half-edge graphs.
"""
module BasicGraphs
export AbstractGraph, Graph, nv, ne, src, tgt, edges, vertices,
  has_edge, has_vertex, add_edge!, add_edges!, add_vertex!, add_vertices!,
  rem_edge!, rem_edges!, rem_vertex!, rem_vertices!,
  neighbors, inneighbors, outneighbors, all_neighbors,
  AbstractSymmetricGraph, SymmetricGraph, inv,
  AbstractReflexiveGraph, ReflexiveGraph, refl,
  AbstractSymmetricReflexiveGraph, SymmetricReflexiveGraph,
  AbstractHalfEdgeGraph, HalfEdgeGraph, vertex, half_edges,
  add_dangling_edge!, add_dangling_edges!

using Compat: only

import Base: inv
import LightGraphs: SimpleGraph, SimpleDiGraph,
  nv, ne, src, dst, edges, vertices, has_edge, has_vertex,
  add_edge!, add_vertex!, add_vertices!, rem_edge!, rem_vertex!, rem_vertices!,
  neighbors, inneighbors, outneighbors, all_neighbors

using ...Present, ...CSetDataStructures

# Graphs
########

@present TheoryGraph(FreeSchema) begin
  V::Ob
  E::Ob
  src::Hom(E,V)
  tgt::Hom(E,V)
end

const AbstractGraph = AbstractACSetType(TheoryGraph)
const Graph = CSetType(TheoryGraph, index=[:src,:tgt])

function (::Type{T})(nv::Int) where T <: AbstractGraph
  g = T(); add_vertices!(g, nv); g
end

nv(g::AbstractACSet) = nparts(g, :V)
ne(g::AbstractACSet) = nparts(g, :E)
ne(g::AbstractACSet, src::Int, tgt::Int) =
  count(subpart(g, e, :tgt) == tgt for e in incident(g, src, :src))

src(g::AbstractACSet, args...) = subpart(g, args..., :src)
tgt(g::AbstractACSet, args...) = subpart(g, args..., :tgt)
dst(g::AbstractACSet, args...) = tgt(g, args...) # LightGraphs compatibility

vertices(g::AbstractACSet) = parts(g, :V)
edges(g::AbstractACSet) = parts(g, :E)
edges(g::AbstractACSet, src::Int, tgt::Int) =
  (e for e in incident(g, src, :src) if subpart(g, e, :tgt) == tgt)

has_vertex(g::AbstractACSet, v) = has_part(g, :V, v)
has_edge(g::AbstractACSet, e) = has_part(g, :E, e)
has_edge(g::AbstractACSet, src::Int, tgt::Int) = tgt ∈ outneighbors(g, src)

add_vertex!(g::AbstractACSet; kw...) = add_part!(g, :V; kw...)
add_vertices!(g::AbstractACSet, n::Int; kw...) = add_parts!(g, :V, n; kw...)

add_edge!(g::AbstractGraph, src::Int, tgt::Int; kw...) =
  add_part!(g, :E; src=src, tgt=tgt, kw...)

function add_edges!(g::AbstractGraph, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(g, :E, n; src=srcs, tgt=tgts, kw...)
end

rem_vertex!(g::AbstractACSet, v::Int; kw...) = rem_vertices!(g, v:v; kw...)

function rem_vertices!(g::AbstractGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    es = reduce(vcat, [incident(g,vs,:src); incident(g,vs,:tgt)], init=Int[])
    rem_parts!(g, :E, unique!(sort!(es)))
  end
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractGraph, e::Int) = rem_part!(g, :E, e)
rem_edge!(g::AbstractACSet, src::Int, tgt::Int) =
  rem_edge!(g, first(edges(g, src, tgt)))

rem_edges!(g::AbstractGraph, es) = rem_parts!(g, :E, es)

neighbors(g::AbstractGraph, v::Int) = outneighbors(g, v)
inneighbors(g::AbstractGraph, v::Int) = subpart(g, incident(g, v, :tgt), :src)
outneighbors(g::AbstractGraph, v::Int) = subpart(g, incident(g, v, :src), :tgt)
all_neighbors(g::AbstractGraph, v::Int) =
  Iterators.flatten((inneighbors(g, v), outneighbors(g, v)))

# Symmetric graphs
##################

@present TheorySymmetricGraph <: TheoryGraph begin
  inv::Hom(E,E)

  compose(inv,inv) == id(E)
  compose(inv,src) == tgt
  compose(inv,tgt) == src
end

# Don't index `inv` because it is self-inverse and don't index `tgt`
# because `src` contains the same information due to symmetry of graph.
const AbstractSymmetricGraph = AbstractACSetType(TheorySymmetricGraph)
const SymmetricGraph = CSetType(TheorySymmetricGraph, index=[:src])

function (::Type{T})(nv::Int) where T <: AbstractSymmetricGraph
  g = T(); add_vertices!(g, nv); g
end

inv(g::AbstractACSet, args...) = subpart(g, args..., :inv)

add_edge!(g::AbstractSymmetricGraph, src::Int, tgt::Int; kw...) =
  add_edges!(g, src:src, tgt:tgt; kw...)

function add_edges!(g::AbstractSymmetricGraph, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  k = nparts(g, :E)
  add_parts!(g, :E, 2n; src=vcat(srcs,tgts), tgt=vcat(tgts,srcs),
             inv=vcat((k+n+1):(k+2n),(k+1):(k+n)), kw...)
end

function rem_vertices!(g::AbstractSymmetricGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    es = reduce(vcat, incident(g, vs, :src), init=Int[])
    rem_parts!(g, :E, unique!(sort!([es; inv(g, es)])))
  end
  # FIXME: Vertex removal is inefficient because `rem_parts!` still searches for
  # edges with given targets but `tgt` is not indexed.
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractSymmetricGraph, e::Int) = rem_edges!(g, e:e)
rem_edges!(g::AbstractSymmetricGraph, es) =
  rem_parts!(g, :E, unique!(sort!([es; inv(g, es)])))

neighbors(g::AbstractSymmetricGraph, v::Int) =
  subpart(g, incident(g, v, :src), :tgt)
inneighbors(g::AbstractSymmetricGraph, v::Int) = neighbors(g, v)
outneighbors(g::AbstractSymmetricGraph, v::Int) = neighbors(g, v)
all_neighbors(g::AbstractSymmetricGraph, v::Int) = neighbors(g, v)

# Reflexive graphs
##################

@present TheoryReflexiveGraph <: TheoryGraph begin
  refl::Hom(V,E)

  compose(refl, src) == id(V)
  compose(refl, tgt) == id(V)
end

const AbstractReflexiveGraph = AbstractACSetType(TheoryReflexiveGraph)
const ReflexiveGraph = CSetType(TheoryReflexiveGraph, index=[:src,:tgt])

function (::Type{T})(nv::Int) where T <: AbstractReflexiveGraph
  g = T(); add_vertices!(g, nv); g
end

refl(g::AbstractACSet, args...) = subpart(g, args..., :refl)

add_vertex!(g::AbstractReflexiveGraph; kw...) =
  only(add_vertices!(g, 1; kw...))

function add_vertices!(g::AbstractReflexiveGraph, n::Int; kw...)
  vs = add_parts!(g, :V, n; kw...)
  es = add_parts!(g, :E, n, src=vs, tgt=vs)
  set_subpart!(g, vs, :refl, es)
  vs
end

add_edge!(g::AbstractReflexiveGraph, src::Int, tgt::Int; kw...) =
  add_part!(g, :E; src=src, tgt=tgt, kw...)

function add_edges!(g::AbstractReflexiveGraph, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(g, :E, n; src=srcs, tgt=tgts, kw...)
end

function rem_vertices!(g::AbstractReflexiveGraph, vs; keep_edges::Bool=false)
  es = if keep_edges
    sort(refl(g, vs)) # Always delete reflexive edges.
  else
    es = reduce(vcat, [incident(g,vs,:src); incident(g,vs,:tgt)], init=Int[])
    unique!(sort!(es))
  end
  rem_parts!(g, :E, es)
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractReflexiveGraph, e::Int) = rem_part!(g, :E, e)
rem_edges!(g::AbstractReflexiveGraph, es) = rem_parts!(g, :E, es)

# Symmetric reflexive graphs
############################

@present TheorySymmetricReflexiveGraph <: TheorySymmetricGraph begin
  refl::Hom(V,E)

  compose(refl, src) == id(V)
  compose(refl, tgt) == id(V)
  compose(refl, inv) == refl # Reflexive loop fixed by involution.
end

const AbstractSymmetricReflexiveGraph =
  AbstractACSetType(TheorySymmetricReflexiveGraph)
const SymmetricReflexiveGraph =
  CSetType(TheorySymmetricReflexiveGraph, index=[:src])

function (::Type{T})(nv::Int) where T <: AbstractSymmetricReflexiveGraph
  g = T(); add_vertices!(g, nv); g
end

add_vertex!(g::AbstractSymmetricReflexiveGraph; kw...) =
  only(add_vertices!(g, 1; kw...))

function add_vertices!(g::AbstractSymmetricReflexiveGraph, n::Int; kw...)
  vs = add_parts!(g, :V, n; kw...)
  es = add_parts!(g, :E, n, src=vs, tgt=vs)
  set_subpart!(g, vs, :refl, es)
  set_subpart!(g, es, :inv, es)
  vs
end

add_edge!(g::AbstractSymmetricReflexiveGraph, src::Int, tgt::Int; kw...) =
  add_edges!(g, src:src, tgt:tgt; kw...)

function add_edges!(g::AbstractSymmetricReflexiveGraph,
                    srcs::AbstractVector{Int}, tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  k = nparts(g, :E)
  add_parts!(g, :E, 2n; src=vcat(srcs,tgts), tgt=vcat(tgts,srcs),
             inv=vcat((k+n+1):(k+2n),(k+1):(k+n)), kw...)
end

function rem_vertices!(g::AbstractSymmetricReflexiveGraph, vs;
                       keep_edges::Bool=false)
  es = if keep_edges
    sort(refl(g, vs)) # Always delete reflexive edges.
  else
    es = reduce(vcat, incident(g, vs, :src), init=Int[])
    unique!(sort!([es; inv(g, es)]))
  end
  rem_parts!(g, :E, es)
  # FIXME: Vertex removal is inefficient for same reason as `SymmetricGraph`.
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractSymmetricReflexiveGraph, e::Int) = rem_edges!(g, e:e)
rem_edges!(g::AbstractSymmetricReflexiveGraph, es) =
  rem_parts!(g, :E, unique!(sort!([es; inv(g, es)])))

# Half-edge graphs
##################

@present TheoryHalfEdgeGraph(FreeSchema) begin
  V::Ob
  H::Ob
  vertex::Hom(H,V)
  inv::Hom(H,H)

  compose(inv, inv) == id(H)
end

const AbstractHalfEdgeGraph = AbstractACSetType(TheoryHalfEdgeGraph)
const HalfEdgeGraph = CSetType(TheoryHalfEdgeGraph, index=[:vertex])

function (::Type{T})(nv::Int) where T <: AbstractHalfEdgeGraph
  g = T(); add_vertices!(g, nv); g
end

vertex(g::AbstractACSet, args...) = subpart(g, args..., :vertex)

half_edges(g::AbstractACSet) = parts(g, :H)
half_edges(g::AbstractACSet, v) = incident(g, v, :vertex)

function half_edge_pairs(g::AbstractACSet, src::Int, tgt::Int)
  hs = half_edges(g, src)
  hs′ = inv(g, hs)
  has_tgt = vertex(g, hs′) .== tgt
  (hs[has_tgt], hs′[has_tgt])
end

@inline add_edge!(g::AbstractHalfEdgeGraph, src::Int, tgt::Int; kw...) =
  add_half_edge_pair!(g, src, tgt; kw...)

@inline add_edges!(g::AbstractHalfEdgeGraph, srcs::AbstractVector{Int},
                   tgts::AbstractVector{Int}; kw...) =
  add_half_edge_pairs!(g, srcs, tgts; kw...)

function add_half_edge_pair!(g::AbstractACSet, src::Int, tgt::Int; kw...)
  k = nparts(g, :H)
  add_parts!(g, :H, 2; vertex=[src,tgt], inv=[k+2,k+1], kw...)
end

function add_half_edge_pairs!(g::AbstractACSet, srcs::AbstractVector{Int},
                              tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  k = nparts(g, :H)
  add_parts!(g, :H, 2n; vertex=vcat(srcs,tgts),
             inv=vcat((k+n+1):(k+2n),(k+1):(k+n)), kw...)
end

add_dangling_edge!(g::AbstractACSet, v::Int; kw...) =
  add_part!(g, :H; vertex=v, inv=nparts(g,:H)+1)

function add_dangling_edges!(g::AbstractACSet, vs::AbstractVector{Int}; kw...)
  n, k = length(vs), nparts(g, :H)
  add_parts!(g, :H, n; vertex=vs, inv=(k+1):(k+n), kw...)
end

function rem_vertices!(g::AbstractHalfEdgeGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    hs = reduce(vcat, incident(g, vs, :vertex), init=Int[])
    rem_parts!(g, :H, unique!(sort!([hs; inv(g, hs)])))
  end
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractHalfEdgeGraph, src::Int, tgt::Int) =
  rem_parts!(g, :H, sort!(unique(first.(half_edge_pairs(g, src, tgt)))))

rem_edge!(g::AbstractHalfEdgeGraph, h::Int) = rem_edges!(g, h:h)
rem_edges!(g::AbstractHalfEdgeGraph, hs) =
  rem_parts!(g, :H, unique!(sort!([hs; inv(g, hs)])))

# LightGraphs constructors
##########################

function (::Type{LG})(g::Union{AbstractGraph,AbstractSymmetricGraph}) where
    LG <: Union{SimpleGraph,SimpleDiGraph}
  lg = LG(nv(g))
  for e in edges(g)
    add_edge!(lg, src(g,e), tgt(g,e))
  end
  lg
end

function SimpleGraph(g::AbstractHalfEdgeGraph)
  lg = SimpleGraph(nv(g))
  for e in half_edges(g)
    e′ = inv(g,e)
    if e <= e′
      add_edge!(lg, vertex(g,e), vertex(g,e′))
    end
  end
  lg
end

end


text:
## State updates ##

function update!(interpreter::Interpreter,
                 domain::Domain, state::State, diff::Diff)
    error("Not implemented.")
end

function update!(interpreter::Interpreter,
                 domain::Domain, state::State, diff::ConditionalDiff)
    if isempty(diff.diffs) return state end
    combined = empty(diff.diffs[1])
    for (cs, d) in zip(diff.conds, diff.diffs)
        satisfy(domain, state, cs) && combine!(combined, d)
    end
    return update!(interpreter, domain, state, combined)
end

function update(interpreter::Interpreter,
                domain::Domain, state::State, diff::Diff)
    return update!(interpreter, domain, copy(state), diff)
end


text:
export sensekey, sensekeys

function sensekey(db::DB, ss::Synset, lem::Lemma)
    db.sensekeys[(ss.offset, lem.word)]
end

function sensekeys(db::DB, lem::Lemma)
    [db.sensekeys[(ss_offset, lem.word)] for ss_offset in lem.synset_offsets]
end


text:
module Jecco

using HDF5
using LinearAlgebra
using FFTW

include("types.jl")
include("fornberg.jl")
include("deriv_operators.jl")
include("grid.jl")
include("input.jl")
include("output.jl")
include("startup.jl")
include("spectral.jl")
include("filtering.jl")
include("utils.jl")

export AbstractPartition, FlattenedVector

export AbstractCoord, CartesianCoord, GaussLobattoCoord
export Chart, Atlas
export Cartesian, GaussLobatto

export CenteredDiff, EqualSizeStencilFD, ChebDeriv, FourierDeriv
export ChebInterpolator
export KO_Filter, Exp_Filter
export KO_Centered

export OpenPMDTimeSeries, get_field
export FieldTimeSeries, get_coords

include("KG_3_1/KG_3_1.jl")
include("AdS5_3_1/AdS5_3_1.jl")

end # module


text:
using Dates
using Downloads: download
using TimeZones: DEPS_DIR

const LATEST_FILE = joinpath(DEPS_DIR, "latest")
const LATEST_FORMAT = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS")
const LATEST_DELAY = Hour(1)  # In 1996 a correction to a release was made an hour later

function read_latest(io::IO)
    version = readline(io)
    retrieved_utc = DateTime(readline(io), LATEST_FORMAT)
    return version, retrieved_utc
end

function read_latest(filename::AbstractString)
    open(filename, "r") do io
        read_latest(io)
    end
end

function write_latest(io::IO, version::AbstractString, retrieved_utc::DateTime=now(Dates.UTC))
    write(io, version)
    write(io, "\n")
    write(io, Dates.format(retrieved_utc, LATEST_FORMAT))
end

const LATEST = let T = Tuple{AbstractString, DateTime}
    isfile(LATEST_FILE) ? Ref{T}(read_latest(LATEST_FILE)) : Ref{T}()
end

function set_latest_cached(version::AbstractString, retrieved_utc::DateTime=now(Dates.UTC))
    LATEST[] = version, retrieved_utc
    open(LATEST_FILE, "w") do io
        write_latest(io, version, retrieved_utc)
    end
end

function latest_cached(now_utc::DateTime=now(Dates.UTC))
    if isassigned(LATEST)
        latest_version, latest_retrieved_utc = LATEST[]

        if now_utc - latest_retrieved_utc < LATEST_DELAY
            return latest_version
        end
    end

    return nothing
end

"""
    tzdata_versions() -> Vector{String}

Retrieves all of the currently available tzdata versions from IANA. The version list is
ordered from earliest to latest.

# Examples
```julia
julia> last(tzdata_versions())  # Current latest available tzdata version
"2020a"
```
"""
function tzdata_versions()
    io = download("https://data.iana.org/time-zones/releases/", IOBuffer())
    html = String(take!(io))

    versions = [
        m[:version]
        for m in eachmatch(r"href=\"tzdata(?<version>(?:\d{2}){1,2}[a-z]?).tar.gz\"", html)
    ]

    # Correctly order releases which include 2-digit years (e.g. "96")
    sort!(versions, by=v -> length(v) < 5 ? "19$v" : v)

    return versions
end

"""
    tzdata_latest_version() -> String

Determine the latest version of tzdata available while limiting how often we check with
remote servers.
"""
function tzdata_latest_version()
    latest_version = latest_cached()

    # Retrieve the current latest version the cached latest has expired
    if latest_version === nothing
        latest_version = last(tzdata_versions())
        set_latest_cached(latest_version)
    end

    return latest_version
end


text:
# Add SC data to beginning of dataset assuming
# SC is equal to average SC prior to the release
function pad_data()

end

function prep_data_csv(filePath)
    df = CSV.read(filePath, DataFrame)
    mat = zeros(2, nrow(df))
    mat[1,:] = df.modelTime
    mat[2,:] = df.SpecificConductance
    return(mat)
end

function prep_data_db(flume, date, start_time, end_time)

end


text:
## Storage of constraints
#
# All `F`-in-`S` constraints are stored in a vector of `ConstraintEntry{F, S}`.
# The index in this vector of a constraint of index
# `ci::MOI.ConstraintIndex{F, S}` is given by `model.constrmap[ci.value]`. The
# advantage of this representation is that it does not require any dictionary
# hence it never needs to compute a hash.
#
# It may seem redundant to store the constraint index `ci` as well as the
# function and sets in the tuple but it is used to efficiently implement the
# getter for `MOI.ListOfConstraintIndices{F, S}`. It is also used to implement
# `MOI.delete`. Indeed, when a constraint is deleted, it is removed from the
# vector hence the index in the vector of all the functions that were stored
# after must be decreased by one. As the constraint index is stored in the
# vector, it readily gives the entries of `model.constrmap` that need to be
# updated.

struct VectorOfConstraints{
    F<:MOI.AbstractFunction,
    S<:MOI.AbstractSet
} <: MOI.ModelLike
    # FIXME: It is not ideal that we have `DataType` here, it might induce type
    #        instabilities. We should change `CleverDicts` so that we can just
    #        use `typeof(CleverDicts.index_to_key)` here.
    constraints::CleverDicts.CleverDict{
        MOI.ConstraintIndex{F,S},
        Tuple{F,S},
        typeof(CleverDicts.key_to_index),
        Base.Fix1{typeof(CleverDicts.index_to_key),DataType}
    }

    function VectorOfConstraints{F,S}() where {F,S}
        return new{F,S}(
            CleverDicts.CleverDict{MOI.ConstraintIndex{F,S},Tuple{F,S}}()
        )
    end
end

MOI.is_empty(v::VectorOfConstraints)  = isempty(v.constraints)
MOI.empty!(v::VectorOfConstraints) = empty!(v.constraints)

function MOI.add_constraint(
    v::VectorOfConstraints{F,S},
    func::F,
    set::S,
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    # We canonicalize the constraint so that solvers can avoid having to
    # canonicalize it most of the time (they can check if they need to with
    # `is_canonical`.
    # Note that the canonicalization is not guaranteed if for instance
    # `modify` is called and adds a new term.
    # See https://github.com/jump-dev/MathOptInterface.jl/pull/1118
    return CleverDicts.add_item(v.constraints, (canonical(func), copy(set)))
end

function MOI.is_valid(
    v::VectorOfConstraints{F,S},
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    return haskey(v.constraints, ci)
end

function MOI.delete(
    v::VectorOfConstraints{F,S},
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.throw_if_not_valid(v, ci)
    delete!(v.constraints, ci)
    return
end

function MOI.get(
    v::VectorOfConstraints{F,S},
    ::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.throw_if_not_valid(v, ci)
    return v.constraints[ci][1]
end

function MOI.get(
    v::VectorOfConstraints{F,S},
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.throw_if_not_valid(v, ci)
    return v.constraints[ci][2]
end

function MOI.set(
    v::VectorOfConstraints{F,S},
    ::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F,S},
    func::F,
) where {F,S}
    MOI.throw_if_not_valid(v, ci)
    v.constraints[ci] = (func, v.constraints[ci][2])
    return
end

function MOI.set(
    v::VectorOfConstraints{F,S},
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F,S},
    set::S,
) where {F,S}
    MOI.throw_if_not_valid(v, ci)
    v.constraints[ci] = (v.constraints[ci][1], set)
    return
end

function MOI.get(
    v::VectorOfConstraints{F,S},
    ::MOI.ListOfConstraints,
)::Vector{Tuple{DataType,DataType}} where {F,S}
    return isempty(v.constraints) ? [] : [(F, S)]
end

function MOI.get(
    v::VectorOfConstraints{F,S},
    ::MOI.NumberOfConstraints{F,S},
) where {F,S}
    return length(v.constraints)
end

function MOI.get(
    v::VectorOfConstraints{F,S},
    ::MOI.ListOfConstraintIndices{F,S},
) where {F,S}
    return keys(v.constraints)
end

function MOI.modify(
    v::VectorOfConstraints{F,S},
    ci::MOI.ConstraintIndex{F,S},
    change::MOI.AbstractFunctionModification,
) where {F,S}
    func, set = v.constraints[ci]
    v.constraints[ci] = (modify_function(func, change), set)
    return
end

# Deletion of variables in vector of variables

function _remove_variable(v::VectorOfConstraints, vi::MOI.VariableIndex)
    CleverDicts.map_values!(v.constraints) do func_set
        return remove_variable(func_set..., vi)
    end
    return
end
function _filter_variables(keep::Function, v::VectorOfConstraints)
    CleverDicts.map_values!(v.constraints) do func_set
        return filter_variables(keep, func_set...)
    end
    return
end

function throw_delete_variable_in_vov(vi::MOI.VariableIndex)
    message = string(
        "Cannot delete variable as it is constrained with other",
        " variables in a `MOI.VectorOfVariables`.",
    )
    return throw(MOI.DeleteNotAllowed(vi, message))
end

# Nothing to do as it's not `VectorOfVariables` constraints
_throw_if_cannot_delete(::VectorOfConstraints, vis, fast_in_vis) = nothing

function _throw_if_cannot_delete(
    v::VectorOfConstraints{MOI.VectorOfVariables,S},
    vis,
    fast_in_vis,
) where {S<:MOI.AbstractVectorSet}
    if MOI.supports_dimension_update(S)
        return
    end
    for fs in values(v.constraints)
        f = fs[1]::MOI.VectorOfVariables
        if length(f.variables) > 1 && f.variables != vis
            for vi in f.variables
                if vi in fast_in_vis
                    # If `supports_dimension_update(S)` then the variable
                    # will be removed in `_filter_variables`.
                    throw_delete_variable_in_vov(vi)
                end
            end
        end
    end
    return
end

function _delete_variables(
    ::Function,
    ::VectorOfConstraints,
    ::Vector{MOI.VariableIndex},
)
    return  # Nothing to do as it's not `VectorOfVariables` constraints
end

function _delete_variables(
    callback::Function,
    v::VectorOfConstraints{MOI.VectorOfVariables,S},
    vis::Vector{MOI.VariableIndex},
) where {S<:MOI.AbstractVectorSet}
    filter!(v.constraints) do p
        f = p.second[1]
        del = if length(f.variables) == 1
            first(f.variables) in vis
        else
            vis == f.variables
        end
        if del
            callback(p.first)
        end
        return !del
    end
    return
end

function _deleted_constraints(
    callback::Function,
    v::VectorOfConstraints,
    vi::MOI.VariableIndex,
)
    vis = [vi]
    _delete_variables(callback, v, vis)
    _remove_variable(v, vi)
    return
end

function _deleted_constraints(
    callback::Function,
    v::VectorOfConstraints,
    vis::Vector{MOI.VariableIndex},
)
    removed = Set(vis)
    _delete_variables(callback, v, vis)
    _filter_variables(vi -> !(vi in removed), v)
    return
end


text:
# https://github.com/open-mpi/ompi/blob/master/ompi/include/mpi.h.in

const MPI_Aint = Int
const MPI_Count = Int64
const MPI_Offset = Int64

const MPI_Status_size = 24
const MPI_Status_Source_offset = 0
const MPI_Status_Tag_offset = 4
const MPI_Status_Error_offset = 8

const MPI_Info = Ptr{Cvoid}
MPI_Info_f2c(c::Cint) = ccall(:MPI_Info_f2c,MPI_Info,(Cint,),c)
const MPI_INFO_NULL = Cint(0)

const MPI_Win = Ptr{Cvoid}
MPI_Win_f2c(c::Cint) = ccall(:MPI_Win_f2c,MPI_Win,(Cint,),c)
const MPI_WIN_NULL = Cint(0)

const MPI_Comm = Ptr{Cvoid}
MPI_Comm_f2c(c::Cint) = ccall(:MPI_Comm_f2c,MPI_Comm,(Cint,),c)
const MPI_COMM_NULL = Cint(2)
const MPI_COMM_SELF = Cint(1)
const MPI_COMM_WORLD = Cint(0)

const MPI_Request = Ptr{Cvoid}
MPI_Request_f2c(c::Cint) = ccall(:MPI_Request_f2c,MPI_Request,(Cint,),c)
const MPI_REQUEST_NULL = Cint(0)

const MPI_File = Ptr{Cvoid}
MPI_File_f2c(c::Cint) = ccall(:MPI_File_f2c,MPI_File,(Cint,),c)
const MPI_FILE_NULL = Cint(0)

const MPI_Op = Ptr{Cvoid}
MPI_Op_f2c(c::Cint) = ccall(:MPI_Op_f2c,MPI_Op,(Cint,),c)
const MPI_OP_NULL = Cint(0)
const MPI_BAND = Cint(6)
const MPI_BOR = Cint(8)
const MPI_BXOR = Cint(10)
const MPI_LAND = Cint(5)
const MPI_LOR = Cint(7)
const MPI_LXOR = Cint(9)
const MPI_MAX = Cint(1)
const MPI_MIN = Cint(2)
const MPI_SUM = Cint(3)
const MPI_PROD = Cint(4)
const MPI_MAXLOC = Cint(11)
const MPI_MINLOC = Cint(12)
const MPI_REPLACE = Cint(13)
const MPI_NO_OP = Cint(14)

const MPI_Datatype = Ptr{Cvoid}
MPI_Datatype_f2c(c::Cint) = ccall(:MPI_Type_f2c,MPI_Datatype,(Cint,),c)
const MPI_DATATYPE_NULL = Cint(0)
const MPI_PACKED = Cint(2)
const MPI_CHAR = Cint(34)
const MPI_SIGNED_CHAR = Cint(36)
const MPI_UNSIGNED_CHAR = Cint(35)
const MPI_WCHAR = Cint(33)
const MPI_BYTE = Cint(1)
const MPI_SHORT = Cint(37)
const MPI_UNSIGNED_SHORT = Cint(38)
const MPI_INT = Cint(39)
const MPI_UNSIGNED = Cint(40)
const MPI_LONG = Cint(41)
const MPI_UNSIGNED_LONG = Cint(42)
const MPI_FLOAT = Cint(45)
const MPI_DOUBLE = Cint(46)
const MPI_INT8_T = Cint(58)
const MPI_INT16_T = Cint(60)
const MPI_INT32_T = Cint(62)
const MPI_INT64_T = Cint(64)
const MPI_UINT8_T = Cint(59)
const MPI_UINT16_T = Cint(61)
const MPI_UINT32_T = Cint(63)
const MPI_UINT64_T = Cint(65)
const MPI_C_FLOAT_COMPLEX = Cint(69)
const MPI_C_DOUBLE_COMPLEX = Cint(70)

const MPI_THREAD_SINGLE = Cint(0)
const MPI_THREAD_FUNNELED = Cint(1)
const MPI_THREAD_SERIALIZED = Cint(2)
const MPI_THREAD_MULTIPLE = Cint(3)
const MPI_PROC_NULL = Cint(-2)
const MPI_ANY_SOURCE = Cint(-1)
const MPI_ANY_TAG = Cint(-1)
const MPI_SUCCESS = Cint(0)
const MPI_UNDEFINED = Cint(-32766)
const MPI_LOCK_EXCLUSIVE = Cint(1)
const MPI_LOCK_SHARED = Cint(2)
const MPI_MAX_INFO_KEY = Cint(36)
const MPI_MAX_INFO_VAL = Cint(256)
const MPI_TAG_UB = Cint(0)
const MPI_COMM_TYPE_SHARED = Cint(0)
const MPI_ORDER_C = Cint(0)
const MPI_ORDER_FORTRAN = Cint(1)
const MPI_UNIVERSE_SIZE = Cint(6)
const MPI_MAX_ERROR_STRING = Cint(256)
const MPI_SUCCESS = Cint(0)
const MPI_MODE_RDONLY = Cint(2)
const MPI_MODE_RDWR = Cint(8)
const MPI_MODE_WRONLY = Cint(4)
const MPI_MODE_CREATE = Cint(1)
const MPI_MODE_EXCL = Cint(64)
const MPI_MODE_DELETE_ON_CLOSE = Cint(16)
const MPI_MODE_UNIQUE_OPEN = Cint(32)
const MPI_MODE_SEQUENTIAL = Cint(256)
const MPI_MODE_APPEND = Cint(128)

const MPI_SEEK_SET = Cint(600)
const MPI_SEEK_CUR = Cint(602)
const MPI_SEEK_END = Cint(604)

const MPI_IDENT = Cint(0)
const MPI_CONGRUENT = Cint(1)
const MPI_SIMILAR = Cint(2)
const MPI_UNEQUAL = Cint(3)

const MPI_BOTTOM = reinterpret(SentinelPtr, 0)
const MPI_IN_PLACE = reinterpret(SentinelPtr, 1)
const MPI_STATUS_IGNORE = reinterpret(SentinelPtr, 0)
const MPI_STATUSES_IGNORE = reinterpret(SentinelPtr, 0)


text:
# Formulas for representing and working with linear-model-type expressions
# Original by Harlan D. Harris.  Later modifications by John Myles White
# and Douglas M. Bates.

## Formulas are written as expressions and parsed by the Julia parser.
## For example :(y ~ a + b + log(c))
## In Julia the & operator is used for an interaction.  What would be written
## in R as y ~ a + b + a:b is written :(y ~ a + b + a&b) in Julia.
## The equivalent R expression, y ~ a*b, is the same in Julia

## The lhs of a one-sided formula is 'nothing'
## The rhs of a formula can be 1

type Formula
    lhs::@compat(Union{Symbol, Expr, Void})
    rhs::@compat(Union{Symbol, Expr, Integer})
end

macro ~(lhs, rhs)
    ex = Expr(:call,
              :Formula,
              Base.Meta.quot(lhs),
              Base.Meta.quot(rhs))
    return ex
end

#
# TODO: implement contrast types in Formula/Terms
#
type Terms
    terms::Vector
    eterms::Vector        # evaluation terms
    factors::Matrix{Bool} # maps terms to evaluation terms
    ## An eterms x terms matrix which is true for terms that need to be "promoted"
    ## to full rank in constructing a model matrx
    is_non_redundant::Matrix{Bool}
# order can probably be dropped.  It is vec(sum(factors, 1))
    order::Vector{Int}    # orders of rhs terms
    response::Bool        # indicator of a response, which is eterms[1] if present
    intercept::Bool       # is there an intercept column in the model matrix?
end

type ModelFrame
    df::AbstractDataFrame
    terms::Terms
    msng::BitArray
    ## mapping from df keys to contrasts matrices
    contrasts::Dict{Symbol, ContrastsMatrix}
end

typealias AbstractFloatMatrix{T<:AbstractFloat} AbstractMatrix{T}

type ModelMatrix{T <: AbstractFloatMatrix}
    m::T
    assign::Vector{Int}
end

Base.size(mm::ModelMatrix) = size(mm.m)
Base.size(mm::ModelMatrix, dim...) = size(mm.m, dim...)

function Base.show(io::IO, f::Formula)
    print(io,
          string("Formula: ",
                 f.lhs == nothing ? "" : f.lhs, " ~ ", f.rhs))
end

## Return, as a vector of symbols, the names of all the variables in
## an expression or a formula
function allvars(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    cc=Symbol[]
    for i in ex.args[2:end] cc=append!(cc,allvars(i)) end
    cc
end
allvars(f::Formula) = unique(vcat(allvars(f.rhs), allvars(f.lhs)))
allvars(sym::Symbol) = [sym]
allvars(v::Any) = Array(Symbol, 0)

# special operators in formulas
const specials = Set([:+, :-, :*, :/, :&, :|, :^])

function dospecials(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in specials) return ex end
    excp = copy(ex)
    excp.args = vcat(a1,map(dospecials, ex.args[2:end]))
    if a1 != :* return excp end
    aa = excp.args
    a2 = aa[2]
    a3 = aa[3]
    if length(aa) > 3
        excp.args = vcat(a1, aa[3:end])
        a3 = dospecials(excp)
    end
    ## this order of expansion gives the R-style ordering of interaction
    ## terms (after sorting in increasing interaction order) for higher-
    ## order interaction terms (e.g. x1 * x2 * x3 should expand to x1 +
    ## x2 + x3 + x1&x2 + x1&x3 + x2&x3 + x1&x2&x3)
    :($a2 + $a2 & $a3 + $a3)
end
dospecials(a::Any) = a

## Distribution of & over +
const distributive = @compat Dict(:& => :+)

distribute(ex::Expr) = distribute!(copy(ex))
distribute(a::Any) = a
## apply distributive property in-place
function distribute!(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    [distribute!(a) for a in ex.args[2:end]]
    ## check that top-level can be distributed
    a1 = ex.args[1]
    if a1 in keys(distributive)

        ## which op is being DISTRIBUTED (e.g. &, *)?
        distributed_op = a1
        ## which op is doing the distributing (e.g. +)?
        distributing_op = distributive[a1]

        ## detect distributing sub-expression (first arg is, e.g. :+)
        is_distributing_subex(e) =
            typeof(e)==Expr && e.head == :call && e.args[1] == distributing_op

        ## find first distributing subex
        first_distributing_subex = findfirst(is_distributing_subex, ex.args)
        if first_distributing_subex != 0
            ## remove distributing subexpression from args
            subex = splice!(ex.args, first_distributing_subex)

            newargs = Any[distributing_op]
            ## generate one new sub-expression, which calls the distributed operation
            ## (e.g. &) on each of the distributing sub-expression's arguments, plus
            ## the non-distributed arguments of the original expression.
            for a in subex.args[2:end]
                new_subex = copy(ex)
                push!(new_subex.args, a)
                ## need to recurse here, in case there are any other
                ## distributing operations in the sub expression
                distribute!(new_subex)
                push!(newargs, new_subex)
            end
            ex.args = newargs
        end
    end
    ex
end
distribute!(a::Any) = a

const associative = Set([:+,:*,:&])       # associative special operators

## If the expression is a call to the function s return its arguments
## Otherwise return the expression
function ex_or_args(ex::Expr,s::Symbol)
    if ex.head != :call error("Non-call expression encountered") end
    if ex.args[1] == s
        ## recurse in case there are more :calls of s below
        return vcat(map(x -> ex_or_args(x, s), ex.args[2:end])...)
    else
        ## not a :call to s, return condensed version of ex
        return condense(ex)
    end
end
ex_or_args(a,s::Symbol) = a

## Condense calls like :(+(a,+(b,c))) to :(+(a,b,c))
function condense(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in associative) return ex end
    excp = copy(ex)
    excp.args = vcat(a1, map(x->ex_or_args(x,a1), ex.args[2:end])...)
    excp
end
condense(a::Any) = a

## always return an ARRAY of terms
getterms(ex::Expr) = (ex.head == :call && ex.args[1] == :+) ? ex.args[2:end] : Expr[ex]
getterms(a::Any) = Any[a]

ord(ex::Expr) = (ex.head == :call && ex.args[1] == :&) ? length(ex.args)-1 : 1
ord(a::Any) = 1

const nonevaluation = Set([:&,:|])        # operators constructed from other evaluations
## evaluation terms - the (filtered) arguments for :& and :|, otherwise the term itself
function evt(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    if !(ex.args[1] in nonevaluation) return ex end
    filter(x->!isa(x,Number), vcat(map(getterms, ex.args[2:end])...))
end
evt(a) = Any[a]

function Terms(f::Formula)
    rhs = condense(distribute(dospecials(f.rhs)))
    tt = unique(getterms(rhs))
    tt = tt[!(tt .== 1)]             # drop any explicit 1's
    noint = (tt .== 0) | (tt .== -1) # should also handle :(-(expr,1))
    tt = tt[!noint]
    oo = Int[ord(t) for t in tt]     # orders of interaction terms
    if !issorted(oo)                 # sort terms by increasing order
        pp = sortperm(oo)
        tt = tt[pp]
        oo = oo[pp]
    end
    etrms = map(evt, tt)
    haslhs = f.lhs != nothing
    if haslhs
        unshift!(etrms, Any[f.lhs])
        unshift!(oo, 1)
    end
    ev = unique(vcat(etrms...))
    sets = [Set(x) for x in etrms]
    facs = Bool[t in s for t in ev, s in sets]
    non_redundants = fill(false, size(facs)) # initialize to false
    Terms(tt, ev, facs, non_redundants, oo, haslhs, !any(noint))
end

## Default NA handler.  Others can be added as keyword arguments
function na_omit(df::DataFrame)
    cc = complete_cases(df)
    df[cc,:], cc
end

## Trim the pool field of da to only those levels that occur in the refs
function dropunusedlevels!(da::PooledDataArray)
    rr = da.refs
    uu = unique(rr)
    length(uu) == length(da.pool) && return da
    T = eltype(rr)
    su = sort!(uu)
    dict = Dict(zip(su, one(T):convert(T, length(uu))))
    da.refs = map(x -> dict[x], rr)
    da.pool = da.pool[uu]
    da
end
dropunusedlevels!(x) = x

is_categorical(::PooledDataArray) = true
is_categorical(::Any) = false

## Check for non-redundancy of columns.  For instance, if x is a factor with two
## levels, it should be expanded into two columns in y~0+x but only one column
## in y~1+x.  The default is the rank-reduced form (contrasts for n levels only
## produce n-1 columns).  In general, an evaluation term x within a term
## x&y&... needs to be "promoted" to full rank if y&... hasn't already been
## included (either explicitly in the Terms or implicitly by promoting another
## term like z in z&y&...).
##
## This modifies the Terms, setting `trms.is_non_redundant = true` for all non-
## redundant evaluation terms.
function check_non_redundancy!(trms::Terms, df::AbstractDataFrame)

    (n_eterms, n_terms) = size(trms.factors)

    encountered_columns = Vector{eltype(trms.factors)}[]

    if trms.intercept
        push!(encountered_columns, zeros(eltype(trms.factors), n_eterms))
    end

    for i_term in 1:n_terms
        for i_eterm in 1:n_eterms
            ## only need to check eterms that are included and can be promoted
            ## (e.g., categorical variables that expand to multiple mm columns)
            if Bool(trms.factors[i_eterm, i_term]) && is_categorical(df[trms.eterms[i_eterm]])
                dropped = trms.factors[:,i_term]
                dropped[i_eterm] = 0

                if findfirst(encountered_columns, dropped) == 0
                    trms.is_non_redundant[i_eterm, i_term] = true
                    push!(encountered_columns, dropped)
                end

            end
        end
        ## once we've checked all the eterms in this term, add it to the list
        ## of encountered terms/columns
        push!(encountered_columns, Compat.view(trms.factors, :, i_term))
    end

    return trms.is_non_redundant
end

const DEFAULT_CONTRASTS = DummyCoding

## Set up contrasts:
## Combine actual DF columns and contrast types if necessary to compute the
## actual contrasts matrices, levels, and term names (using DummyCoding
## as the default)
function evalcontrasts(df::AbstractDataFrame, contrasts::Dict = Dict())
    evaledContrasts = Dict()
    for (term, col) in eachcol(df)
        is_categorical(col) || continue
        evaledContrasts[term] = ContrastsMatrix(haskey(contrasts, term) ?
                                                contrasts[term] :
                                                DEFAULT_CONTRASTS(),
                                                col)
    end
    return evaledContrasts
end

function ModelFrame(trms::Terms, d::AbstractDataFrame;
                    contrasts::Dict = Dict())
    df, msng = na_omit(DataFrame(map(x -> d[x], trms.eterms)))
    names!(df, convert(Vector{Symbol}, map(string, trms.eterms)))
    for c in eachcol(df) dropunusedlevels!(c[2]) end

    evaledContrasts = evalcontrasts(df, contrasts)

    ## Check for non-redundant terms, modifying terms in place
    check_non_redundancy!(trms, df)

    ModelFrame(df, trms, msng, evaledContrasts)
end

ModelFrame(df::AbstractDataFrame, term::Terms, msng::BitArray) = ModelFrame(df, term, msng, evalcontrasts(df))
ModelFrame(f::Formula, d::AbstractDataFrame; kwargs...) = ModelFrame(Terms(f), d; kwargs...)
ModelFrame(ex::Expr, d::AbstractDataFrame; kwargs...) = ModelFrame(Formula(ex), d; kwargs...)

## modify contrasts in place
function setcontrasts!(mf::ModelFrame, new_contrasts::Dict)
    new_contrasts = Dict([ Pair(col, ContrastsMatrix(contr, mf.df[col]))
                      for (col, contr) in filter((k,v)->haskey(mf.df, k), new_contrasts) ])

    mf.contrasts = merge(mf.contrasts, new_contrasts)
    return mf
end
setcontrasts!(mf::ModelFrame; kwargs...) = setcontrasts!(mf, Dict(kwargs))

"""
    StatsBase.model_response(mf::ModelFrame)
Extract the response column, if present.  `DataVector` or
`PooledDataVector` columns are converted to `Array`s
"""
function StatsBase.model_response(mf::ModelFrame)
    if mf.terms.response
        convert(Array, mf.df[mf.terms.eterms[1]])
    else
        error("Model formula one-sided")
    end
end

## construct model matrix columns from model frame + name (checks for contrasts)
function modelmat_cols{T<:AbstractFloatMatrix}(::Type{T}, name::Symbol, mf::ModelFrame; non_redundant::Bool = false)
    if haskey(mf.contrasts, name)
        modelmat_cols(T, mf.df[name],
                      non_redundant ?
                      ContrastsMatrix{FullDummyCoding}(mf.contrasts[name]) :
                      mf.contrasts[name])
    else
        modelmat_cols(T, mf.df[name])
    end
end

modelmat_cols{T<:AbstractFloatMatrix}(::Type{T}, v::DataVector) = convert(T, reshape(v.data, length(v), 1))
modelmat_cols{T<:AbstractFloatMatrix}(::Type{T}, v::Vector) = convert(T, reshape(v, length(v), 1))

"""
    modelmat_cols{T<:AbstractFloatMatrix}(::Type{T}, v::PooledDataVector, contrast::ContrastsMatrix)

Construct `ModelMatrix` columns of type `T` based on specified contrasts, ensuring that
levels align properly.
"""
function modelmat_cols{T<:AbstractFloatMatrix}(::Type{T}, v::PooledDataVector, contrast::ContrastsMatrix)
    ## make sure the levels of the contrast matrix and the categorical data
    ## are the same by constructing a re-indexing vector. Indexing into
    ## reindex with v.refs will give the corresponding row number of the
    ## contrast matrix
    reindex = [findfirst(contrast.levels, l) for l in levels(v)]
    contrastmatrix = convert(T, contrast.matrix)
    return indexrows(contrastmatrix, reindex[v.refs])
end

indexrows(m::SparseMatrixCSC, ind::Vector{Int}) = m'[:, ind]'
indexrows(m::AbstractMatrix, ind::Vector{Int}) = m[ind, :]

"""
    expandcols{T<:AbstractFloatMatrix}(trm::Vector{T})
Create pairwise products of columns from a vector of matrices
"""
function expandcols{T<:AbstractFloatMatrix}(trm::Vector{T})
    if length(trm) == 1
        trm[1]
    else
        a = trm[1]
        b = expandcols(trm[2 : end])
        reduce(hcat, [broadcast(*, a, Compat.view(b, :, j)) for j in 1 : size(b, 2)])
    end
end

"""
    droprandomeffects(trms::Terms)
Expressions of the form `(a|b)` are "random-effects" terms and are not
incorporated in the ModelMatrix.  This function checks for such terms and,
if any are present, drops them from the `Terms` object.
"""
function droprandomeffects(trms::Terms)
    retrms = Bool[Meta.isexpr(t, :call) && t.args[1] == :| for t in trms.terms]
    if !any(retrms)  # return trms unchanged
        trms
    elseif all(retrms) && !trms.response   # return an empty Terms object
        Terms(Any[],Any[],Array(Bool, (0,0)),Array(Bool, (0,0)), Int[], false, trms.intercept)
    else
        # the rows of `trms.factors` correspond to `eterms`, the columns to `terms`
        # After dropping random-effects terms we drop any eterms whose rows are all false
        ckeep = !retrms                 # columns to retain
        facs = trms.factors[:, ckeep]
        rkeep = vec(sum(facs, 2) .> 0)
        Terms(trms.terms[ckeep], trms.eterms[rkeep], facs[rkeep, :],
              trms.is_non_redundant[rkeep, ckeep],
              trms.order[ckeep], trms.response, trms.intercept)
    end
end

"""
    dropresponse!(trms::Terms)
Drop the response term, `trms.eterms[1]` and the first row and column
of `trms.factors` if `trms.response` is true.
"""
function dropresponse!(trms::Terms)
    if trms.response
        ckeep = 2:size(trms.factors, 2)
        rkeep = vec(any(trms.factors[:, ckeep], 2))
        Terms(trms.terms, trms.eterms[rkeep], trms.factors[rkeep, ckeep],
              trms.is_non_redundant[rkeep, ckeep], trms.order[ckeep], false, trms.intercept)
    else
        trms
    end
end

"""
    ModelMatrix{T<:AbstractFloatMatrix}(mf::ModelFrame)
Create a `ModelMatrix` of type `T` (default `Matrix{Float64}`) from the
`terms` and `df` members of `mf`.

This is basically a map-reduce where terms are mapped to columns by `cols`
and reduced by `hcat`.  During the collection of the columns the `assign`
vector is created.  `assign` maps columns of the model matrix to terms in
the model frame.  It can also be considered as mapping coefficients to
terms and, hence, to names.

If there is an intercept in the model, that column occurs first and its
`assign` value is zero.

Mixed-effects models include "random-effects" terms which are ignored when
creating the model matrix.
"""
@compat function (::Type{ModelMatrix{T}}){T<:AbstractFloatMatrix}(mf::ModelFrame)
    dfrm = mf.df
    terms = droprandomeffects(dropresponse!(mf.terms))

    blocks = T[]
    assign = Int[]
    if terms.intercept
        push!(blocks, ones(size(dfrm, 1), 1))  # columns of 1's is first block
        push!(assign, 0)                       # this block corresponds to term zero
    end

    factors = terms.factors

    ## Map eval. term name + redundancy bool to cached model matrix columns
    eterm_cols = @compat Dict{Tuple{Symbol,Bool}, T}()
    ## Accumulator for each term's vector of eval. term columns.

    ## TODO: this method makes multiple copies of the data in the ModelFrame:
    ## first in term_cols (1-2x per evaluation term, depending on redundancy),
    ## second in constructing the matrix itself.

    ## turn each term into a vector of mm columns for its eval. terms, using
    ## "promoted" full-rank versions of categorical columns for non-redundant
    ## eval. terms:
    for (i_term, term) in enumerate(terms.terms)
        term_cols = T[]
        ## Pull out the eval terms, and the non-redundancy flags for this term
        ff = Compat.view(factors, :, i_term)
        eterms = Compat.view(terms.eterms, ff)
        non_redundants = Compat.view(terms.is_non_redundant, ff, i_term)
        ## Get cols for each eval term (either previously generated, or generating
        ## and storing as necessary)
        for (et, nr) in zip(eterms, non_redundants)
            if ! haskey(eterm_cols, (et, nr))
                eterm_cols[(et, nr)] = modelmat_cols(T, et, mf, non_redundant=nr)
            end
            push!(term_cols, eterm_cols[(et, nr)])
        end
        push!(blocks, expandcols(term_cols))
        append!(assign, fill(i_term, size(blocks[end], 2)))
    end

    if isempty(blocks)
        error("Could not construct model matrix. Resulting matrix has 0 columns.")
    end

    I = size(dfrm, 1)
    J = mapreduce(x -> size(x, 2), +, blocks)
    X = similar(blocks[1], I, J)
    i = 1
    for block in blocks
        len = size(block, 2)
        X[:, i:(i + len - 1)] = block
        i += len
    end
    ModelMatrix{T}(X, assign)
end
ModelMatrix(mf::ModelFrame) = ModelMatrix{Matrix{Float64}}(mf)

"""
    termnames(term::Symbol, col)
Returns a vector of strings with the names of the coefficients
associated with a term.  If the column corresponding to the term
is not a `PooledDataArray` a one-element vector is returned.
"""
termnames(term::Symbol, col) = [string(term)]
function termnames(term::Symbol, mf::ModelFrame; non_redundant::Bool = false)
    if haskey(mf.contrasts, term)
        termnames(term, mf.df[term],
                  non_redundant ?
                  ContrastsMatrix{FullDummyCoding}(mf.contrasts[term]) :
                  mf.contrasts[term])
    else
        termnames(term, mf.df[term])
    end
end

termnames(term::Symbol, col::Any, contrast::ContrastsMatrix) =
    ["$term: $name" for name in contrast.termnames]

function expandtermnames(term::Vector)
    if length(term) == 1
        return term[1]
    else
        return foldr((a,b) -> vec([string(lev1, " & ", lev2) for
                                   lev1 in a,
                                   lev2 in b]),
                     term)
    end
end

"""
    coefnames(mf::ModelFrame)
Returns a vector of coefficient names constructed from the Terms
member and the types of the evaluation columns.
"""
function coefnames(mf::ModelFrame)
    terms = droprandomeffects(dropresponse!(mf.terms))

    ## strategy mirrors ModelMatrx constructor:
    eterm_names = @compat Dict{Tuple{Symbol,Bool}, Vector{Compat.UTF8String}}()
    term_names = Vector{Compat.UTF8String}[]

    if terms.intercept
        push!(term_names, Compat.UTF8String["(Intercept)"])
    end

    factors = terms.factors

    for (i_term, term) in enumerate(terms.terms)

        ## names for columns for eval terms
        names = Vector{Compat.UTF8String}[]

        ff = Compat.view(factors, :, i_term)
        eterms = Compat.view(terms.eterms, ff)
        non_redundants = Compat.view(terms.is_non_redundant, ff, i_term)

        for (et, nr) in zip(eterms, non_redundants)
            if !haskey(eterm_names, (et, nr))
                eterm_names[(et, nr)] = termnames(et, mf, non_redundant=nr)
            end
            push!(names, eterm_names[(et, nr)])
        end
        push!(term_names, expandtermnames(names))
    end

    reduce(vcat, Vector{Compat.UTF8String}(), term_names)
end

function Formula(t::Terms)
    lhs = t.response ? t.eterms[1] : nothing
    rhs = Expr(:call,:+)
    if t.intercept
        push!(rhs.args,1)
    end
    append!(rhs.args,t.terms)
    Formula(lhs,rhs)
end


text:
# This script provides implementations of
# multi-step methods of different orders
# for solving a system of differential equations
# and compares them to the exact solutions and
# other methods while solving a 1-body problem

println('\n', " "^4, "> Loading the packages...")

using LinearAlgebra
using Printf

# Define the value of ϰ
const ϰ = -1

# Integrate equations of motion using the
# Euler's method, return the values of position
# and velocity on the last step
function euler(
    r₀₁::Vector{F},
    r₀₂::Vector{F},
    r₀₃::Vector{F},
    v₀₁::Vector{F},
    v₀₂::Vector{F},
    v₀₃::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F},Vector{F},Vector{F},Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀₁)
    # Prepare the output vectors
    r₁ = copy(r₀₁)
    r₂ = copy(r₀₂)
    r₃ = copy(r₀₃)
    v₁ = copy(v₀₁)
    v₂ = copy(v₀₂)
    v₃ = copy(v₀₃)
    # Compute the solutions (in-place)
    for _ in 1:n
        ρ₁₂ = norm(r₁ - r₂)^3
        ρ₁₃ = norm(r₁ - r₃)^3
        ρ₂₃ = norm(r₂ - r₃)^3
        for k in 1:N
            a₁₂ = ϰ * (r₁[k] - r₂[k]) / ρ₁₂
            a₁₃ = ϰ * (r₁[k] - r₃[k]) / ρ₁₃
            a₂₃ = ϰ * (r₂[k] - r₃[k]) / ρ₂₃
            a₂₁ = -a₁₂
            a₃₁ = -a₁₃
            a₃₂ = -a₂₃
            r₁[k] += h * v₁[k]
            r₂[k] += h * v₂[k]
            r₃[k] += h * v₃[k]
            v₁[k] += h * (a₁₂ + a₁₃)
            v₂[k] += h * (a₂₁ + a₂₃)
            v₃[k] += h * (a₃₁ + a₃₂)
        end
    end
    return r₁, r₂, r₃, v₁, v₂, v₃
end

# Integrate the three-body problem using the
# two-step Adams–Bashforth's method, return the
# values of position and velocity on the last step
function ab2(
    r₀₁::Vector{F},
    r₀₂::Vector{F},
    r₀₃::Vector{F},
    v₀₁::Vector{F},
    v₀₂::Vector{F},
    v₀₃::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F},Vector{F},Vector{F},Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀₁)
    # Prepare the output vectors
    r₁ = copy(r₀₁)
    r₂ = copy(r₀₂)
    r₃ = copy(r₀₃)
    v₁ = copy(v₀₁)
    v₂ = copy(v₀₂)
    v₃ = copy(v₀₃)
    # Prepare buffers for previous values
    r₁_ₖ₋₁ = copy(r₀₁)
    r₂_ₖ₋₁ = copy(r₀₂)
    r₃_ₖ₋₁ = copy(r₀₃)
    v₁_ₖ₋₁ = copy(v₀₁)
    v₂_ₖ₋₁ = copy(v₀₂)
    v₃_ₖ₋₁ = copy(v₀₃)
    # Compute the second value of the solution
    # by using the one-step Euler's method
    r₁, r₂, r₃, v₁, v₂, v₃ = euler(r₀₁, r₀₂, r₀₃, v₀₁, v₀₂, v₀₃, h, UInt(1))
    # Define a couple of independent coefficients
    k₁ = 3 / 2 * h
    k₂ = -1 / 2 * h
    # Compute the rest in two steps
    for _ in 2:n
        ρ₁₂_ₖ = norm(r₁ - r₂)^3
        ρ₁₃_ₖ = norm(r₁ - r₃)^3
        ρ₂₃_ₖ = norm(r₂ - r₃)^3
        ρ₁₂_ₖ₋₁ = norm(r₁_ₖ₋₁ - r₂_ₖ₋₁)^3
        ρ₁₃_ₖ₋₁ = norm(r₁_ₖ₋₁ - r₃_ₖ₋₁)^3
        ρ₂₃_ₖ₋₁ = norm(r₂_ₖ₋₁ - r₃_ₖ₋₁)^3
        for k in 1:N
            a₁₂_ₖ = k₁ * ϰ * (r₁[k] - r₂[k]) / ρ₁₂_ₖ
            a₁₃_ₖ = k₁ * ϰ * (r₁[k] - r₃[k]) / ρ₁₃_ₖ
            a₂₃_ₖ = k₁ * ϰ * (r₂[k] - r₃[k]) / ρ₂₃_ₖ
            a₁₂_ₖ₋₁ = k₂ * ϰ * (r₁_ₖ₋₁[k] - r₂_ₖ₋₁[k]) / ρ₁₂_ₖ₋₁
            a₁₃_ₖ₋₁ = k₂ * ϰ * (r₁_ₖ₋₁[k] - r₃_ₖ₋₁[k]) / ρ₁₃_ₖ₋₁
            a₂₃_ₖ₋₁ = k₂ * ϰ * (r₂_ₖ₋₁[k] - r₃_ₖ₋₁[k]) / ρ₂₃_ₖ₋₁

            a₂₁_ₖ = -a₁₂_ₖ
            a₃₁_ₖ = -a₁₃_ₖ
            a₃₂_ₖ = -a₂₃_ₖ
            a₂₁_ₖ₋₁ = -a₁₂_ₖ₋₁
            a₃₁_ₖ₋₁ = -a₁₃_ₖ₋₁
            a₃₂_ₖ₋₁ = -a₂₃_ₖ₋₁

            r₁_ₖ₋₁[k] = r₁[k]
            r₂_ₖ₋₁[k] = r₂[k]
            r₃_ₖ₋₁[k] = r₃[k]

            r₁[k] += k₁ * v₁[k] + k₂ * v₁_ₖ₋₁[k]
            r₂[k] += k₁ * v₂[k] + k₂ * v₂_ₖ₋₁[k]
            r₃[k] += k₁ * v₃[k] + k₂ * v₃_ₖ₋₁[k]

            v₁_ₖ₋₁[k] = v₁[k]
            v₂_ₖ₋₁[k] = v₂[k]
            v₃_ₖ₋₁[k] = v₃[k]

            v₁[k] += a₁₂_ₖ + a₁₃_ₖ + a₁₂_ₖ₋₁ + a₁₃_ₖ₋₁
            v₂[k] += a₂₁_ₖ + a₂₃_ₖ + a₂₁_ₖ₋₁ + a₂₃_ₖ₋₁
            v₃[k] += a₃₁_ₖ + a₃₂_ₖ + a₃₁_ₖ₋₁ + a₃₂_ₖ₋₁
        end
    end
    return r₁, r₂, r₃, v₁, v₂, v₃
end

# Integrate the three-body problem using the
# three-step Adams–Bashforth's method, return the
# values of position and velocity on the last step
function ab3(
    r₀₁::Vector{F},
    r₀₂::Vector{F},
    r₀₃::Vector{F},
    v₀₁::Vector{F},
    v₀₂::Vector{F},
    v₀₃::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F},Vector{F},Vector{F},Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀₁)
    # Prepare the output vectors
    r₁ = copy(r₀₁)
    r₂ = copy(r₀₂)
    r₃ = copy(r₀₃)
    v₁ = copy(v₀₁)
    v₂ = copy(v₀₂)
    v₃ = copy(v₀₃)
    # Prepare buffers for previous values
    r₁_ₖ₋₁ = copy(r₀₁)
    r₂_ₖ₋₁ = copy(r₀₂)
    r₃_ₖ₋₁ = copy(r₀₃)
    v₁_ₖ₋₁ = copy(v₀₁)
    v₂_ₖ₋₁ = copy(v₀₂)
    v₃_ₖ₋₁ = copy(v₀₃)
    r₁_ₖ₋₂ = copy(r₀₁)
    r₂_ₖ₋₂ = copy(r₀₂)
    r₃_ₖ₋₂ = copy(r₀₃)
    v₁_ₖ₋₂ = copy(v₀₁)
    v₂_ₖ₋₂ = copy(v₀₂)
    v₃_ₖ₋₂ = copy(v₀₃)
    # Compute the second value of the solution
    # by using the one-step Euler's method
    r₁, r₂, r₃, v₁, v₂, v₃ = euler(r₀₁, r₀₂, r₀₃, v₀₁, v₀₂, v₀₃, h, UInt(1))
    # Compute the third value of the solution
    # by using the two-step Adams–Bashforth's method
    r₁_ₖ₋₁, r₂_ₖ₋₁, r₃_ₖ₋₁, v₁_ₖ₋₁, v₂_ₖ₋₁, v₃_ₖ₋₁ = r₁, r₂, r₃, v₁, v₂, v₃
    r₁, r₂, r₃, v₁, v₂, v₃ = ab2(r₀₁, r₀₂, r₀₃, v₀₁, v₀₂, v₀₃, h, UInt(2))
    # Define a couple of independent coefficients
    k₁ = 23 / 12 * h
    k₂ = -4 / 3 * h
    k₃ = 5 / 12 * h
    # Compute the rest in two steps
    for _ in 3:n
        ρ₁₂_ₖ = norm(r₁ - r₂)^3
        ρ₁₃_ₖ = norm(r₁ - r₃)^3
        ρ₂₃_ₖ = norm(r₂ - r₃)^3
        ρ₁₂_ₖ₋₁ = norm(r₁_ₖ₋₁ - r₂_ₖ₋₁)^3
        ρ₁₃_ₖ₋₁ = norm(r₁_ₖ₋₁ - r₃_ₖ₋₁)^3
        ρ₂₃_ₖ₋₁ = norm(r₂_ₖ₋₁ - r₃_ₖ₋₁)^3
        ρ₁₂_ₖ₋₂ = norm(r₁_ₖ₋₂ - r₂_ₖ₋₂)^3
        ρ₁₃_ₖ₋₂ = norm(r₁_ₖ₋₂ - r₃_ₖ₋₂)^3
        ρ₂₃_ₖ₋₂ = norm(r₂_ₖ₋₂ - r₃_ₖ₋₂)^3
        for k in 1:N
            a₁₂_ₖ = k₁ * ϰ * (r₁[k] - r₂[k]) / ρ₁₂_ₖ
            a₁₃_ₖ = k₁ * ϰ * (r₁[k] - r₃[k]) / ρ₁₃_ₖ
            a₂₃_ₖ = k₁ * ϰ * (r₂[k] - r₃[k]) / ρ₂₃_ₖ
            a₁₂_ₖ₋₁ = k₂ * ϰ * (r₁_ₖ₋₁[k] - r₂_ₖ₋₁[k]) / ρ₁₂_ₖ₋₁
            a₁₃_ₖ₋₁ = k₂ * ϰ * (r₁_ₖ₋₁[k] - r₃_ₖ₋₁[k]) / ρ₁₃_ₖ₋₁
            a₂₃_ₖ₋₁ = k₂ * ϰ * (r₂_ₖ₋₁[k] - r₃_ₖ₋₁[k]) / ρ₂₃_ₖ₋₁
            a₁₂_ₖ₋₂ = k₃ * ϰ * (r₁_ₖ₋₂[k] - r₂_ₖ₋₂[k]) / ρ₁₂_ₖ₋₂
            a₁₃_ₖ₋₂ = k₃ * ϰ * (r₁_ₖ₋₂[k] - r₃_ₖ₋₂[k]) / ρ₁₃_ₖ₋₂
            a₂₃_ₖ₋₂ = k₃ * ϰ * (r₂_ₖ₋₂[k] - r₃_ₖ₋₂[k]) / ρ₂₃_ₖ₋₂

            a₂₁_ₖ = -a₁₂_ₖ
            a₃₁_ₖ = -a₁₃_ₖ
            a₃₂_ₖ = -a₂₃_ₖ
            a₂₁_ₖ₋₁ = -a₁₂_ₖ₋₁
            a₃₁_ₖ₋₁ = -a₁₃_ₖ₋₁
            a₃₂_ₖ₋₁ = -a₂₃_ₖ₋₁
            a₂₁_ₖ₋₂ = -a₁₂_ₖ₋₂
            a₃₁_ₖ₋₂ = -a₁₃_ₖ₋₂
            a₃₂_ₖ₋₂ = -a₂₃_ₖ₋₂

            r₁_ₖ₋₂[k] = r₁_ₖ₋₁[k]
            r₂_ₖ₋₂[k] = r₂_ₖ₋₁[k]
            r₃_ₖ₋₂[k] = r₃_ₖ₋₁[k]
            r₁_ₖ₋₁[k] = r₁[k]
            r₂_ₖ₋₁[k] = r₂[k]
            r₃_ₖ₋₁[k] = r₃[k]

            r₁[k] += k₁ * v₁[k] + k₂ * v₁_ₖ₋₁[k] + k₃ * v₁_ₖ₋₂[k]
            r₂[k] += k₁ * v₂[k] + k₂ * v₂_ₖ₋₁[k] + k₃ * v₂_ₖ₋₂[k]
            r₃[k] += k₁ * v₃[k] + k₂ * v₃_ₖ₋₁[k] + k₃ * v₃_ₖ₋₂[k]

            v₁_ₖ₋₂[k] = v₁_ₖ₋₁[k]
            v₂_ₖ₋₂[k] = v₂_ₖ₋₁[k]
            v₃_ₖ₋₂[k] = v₃_ₖ₋₁[k]
            v₁_ₖ₋₁[k] = v₁[k]
            v₂_ₖ₋₁[k] = v₂[k]
            v₃_ₖ₋₁[k] = v₃[k]

            v₁[k] += a₁₂_ₖ + a₁₃_ₖ + a₁₂_ₖ₋₁ + a₁₃_ₖ₋₁ + a₁₂_ₖ₋₂ + a₁₃_ₖ₋₂
            v₂[k] += a₂₁_ₖ + a₂₃_ₖ + a₂₁_ₖ₋₁ + a₂₃_ₖ₋₁ + a₂₁_ₖ₋₂ + a₂₃_ₖ₋₂
            v₃[k] += a₃₁_ₖ + a₃₂_ₖ + a₃₁_ₖ₋₁ + a₃₂_ₖ₋₁ + a₃₁_ₖ₋₂ + a₃₂_ₖ₋₂
        end
    end
    return r₁, r₂, r₃, v₁, v₂, v₃
end

println(" "^4, "> Integrating the 3-body problem...")

# Define the initial values of the position and velocity
r₀₁ = [0.97, -0.2431]
r₀₂ = [-0.97, 0.2431]
r₀₃ = [0.0, 0.0]
v₀₁ = [0.4662, 0.4324]
v₀₂ = [0.4662, 0.4324]
v₀₃ = [-0.9324, -0.8647]

# Calculate total energy
E₀ = 1 / 2 * (norm(v₀₁)^2 + norm(v₀₂)^2 + norm(v₀₃)^2) - ϰ^2 *
    (1 / norm(r₀₁ - r₀₂) + 1 / norm(r₀₃ - r₀₂) + 1 / norm(r₀₃ - r₀₁))

# Print the initial values
println(
    '\n',
    " "^6, "r₀₁: ", r₀₁[1], " ", r₀₁[2], '\n',
    " "^6, "r₀₂: ", r₀₂[1], " ", r₀₂[2], '\n',
    " "^6, "r₀₃: ", r₀₃[1], " ", r₀₃[2], '\n',
    " "^6, "v₀₁: ", v₀₁[1], " ", v₀₁[2], '\n',
    " "^6, "v₀₂: ", v₀₂[1], " ", v₀₂[2], '\n',
    " "^6, "v₀₃: ", v₀₃[1], " ", v₀₃[2], '\n',
    " "^6, "E₀: ", E₀,
)

# Define the output directories
const data_dir = "$(@__DIR__)/../data/3-body"
const files = joinpath.(
    data_dir,
    [
        "euler.dat",
        "ab2.dat",
        "ab3.dat",
    ]
)

# Prepare a list of methods
const methods = (euler, ab2, ab3)

# Integrate equations of motion using all available
# methods and write the values of position and velocity
# on the last step
function integrate(h::Vector{F}, n::Vector{I}) where {F <: AbstractFloat,I <: Unsigned}
    # Open the data files
    io = open.(files, "a")

    # Prepare results vectors

    r₁ = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    r₂ = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    r₃ = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    v₁ = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    v₂ = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    v₃ = [ Vector{Vector{F}}() for _ in 1:length(io) ]

    E = [ Vector{F}() for _ in 1:length(io) ]
    ΔE = [ Vector{F}() for _ in 1:length(io) ]

    # For each pair of parameters
    for i in eachindex(h)
        # Integrate and get the results of the last steps
        for j in eachindex(io)
            rᵢ₁, rᵢ₂, rᵢ₃, vᵢ₁, vᵢ₂, vᵢ₃ = methods[j](r₀₁, r₀₂, r₀₃, v₀₁, v₀₂, v₀₃, h[i], n[i])

            push!(r₁[j], rᵢ₁)
            push!(r₂[j], rᵢ₂)
            push!(r₃[j], rᵢ₃)
            push!(v₁[j], vᵢ₁)
            push!(v₂[j], vᵢ₂)
            push!(v₃[j], vᵢ₃)

            Eᵢ = 1 / 2 * (norm(vᵢ₁)^2 + norm(vᵢ₂)^2 + norm(vᵢ₃)^2) - ϰ^2 *
                (1 / norm(rᵢ₁ - rᵢ₂) + 1 / norm(rᵢ₃ - rᵢ₂) + 1 / norm(rᵢ₃ - rᵢ₁))
            ΔEᵢ = abs(E₀ - Eᵢ)

            push!(E[j], Eᵢ)
            push!(ΔE[j], ΔEᵢ)
        end
    end

    # Get the format string based on the index of the iteration
    function get_format(i)
        Printf.Format(
            "\$ 10^{-$(length("$(n[i])") - 1)} \$ " *
            "& \$ $(n[i] % 10 == 0 ? "10^$(length("$(n[i])") - 1)" : n[i]) \$ " *
            "& \$ %.14f \$ " *
            "& \$ %.14f \$ " *
            "\\\\\n"
        )
    end

    # Print the data in the files in the LaTeX format

    for j in eachindex(io)
        println(io[j], "# Positions #1")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), r₁[j][i][1], r₁[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Positions #2")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), r₂[j][i][1], r₂[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Positions #3")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), r₃[j][i][1], r₃[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Velocities #1")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i),  v₁[j][i][1], v₁[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Velocities #2")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i),  v₂[j][i][1], v₂[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Velocities #3")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i),  v₃[j][i][1], v₃[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Total energy")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), E[j][i], ΔE[j][i])
        end
        println(io[j])
        close(io[j])
    end
end

# Truncate the previous results
open.(files; truncate=true)

# Integrate n = 2^m iterations
integrate(
    [1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7],
    UInt.([ 10^i for i in 2:7 ]),
)

println()


text:
using HypothesisTests, UncertainData, Distributions

ts_vals = [UncertainValue(Normal, rand(Uniform(-10, 10)), rand()) for i = 1:10]
ts = UncertainDataset(xs)
#
# @show size(LjungBoxTest(ts))
# @show size(LjungBoxTest(ts, 10))
# @show size(BoxPierceTest(ts))
# @show size(BoxPiercexTest(ts, 10))
#
# @test pvalue(LjungBoxTest(ts)) isa Float64
# @test pvalue.(LjungBoxTest(ts, 10)) isa Vector{Float64}
#
# @test pvalue(BoxPierceTest(ts)) isa Float64
# @test pvalue.(BoxPierceTest(ts, 10)) isa Vector{Float64}


text:
using Test
using PowerModels
using NetDecOPF
using DualDecomposition
using Ipopt
const DD = DualDecomposition

sub_optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "warm_start_init_point" => "yes")
optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "warm_start_init_point" => "yes")

file = "../examples/case5.m"
data = parse_file(file)

# Partition network
partitions = metis_cluster(data, 2)

dn_model = decompose(data, partitions, ACRPowerModel, NetDecOPF.build_acopf_with_free_lines)
# dn_model = decompose(data, partitions, W_ACRModel, NetDecOPF.build_acopf_with_free_lines)

algo = init_DD_algo(dn_model)
set_subnet_optimizer!(dn_model, sub_optimizer)

# Change parameters (for example)
params = BM.Parameters()
BM.set_parameter(params, "maxiter", 3000)
BM.set_parameter(params, "ϵ_s", 1.e-5)

# Lagrange master method
LM = DD.BundleMaster(BM.ProximalMethod, optimizer, params)

DD.run!(algo, LM)

@test isapprox(DD.dual_objective_value(algo), 17551, rtol = 0.1)
# @show DD.dual_solution(algo)


text:
module poly
using LinearAlgebra

# Base.+ OK
# Base.- OK
# Base.* OK
# degree OK
# pow to a positive integer OK
# eval OK
# derivative OK
# integrate OK
# roots OK
# create polynomial from roots OK

struct Polynomial{N, T}
    coeff::NTuple{N, T}

    function Polynomial{N, T}(coeff::NTuple{N, T}) where {N, T}
        return new(coeff)
    end
end

function Polynomial(coeff::AbstractArray{T}) where {T}
    N = length(coeff)
    t = tuple(coeff...)
    return Polynomial{N, T}(t)
end

function Polynomial(coeff::Tuple)
    t = promote(coeff...)
    N = length(coeff)
    T = eltype(t)
    return Polynomial{N,T}(t)
    # return Polynomial(t)
end

function Polynomial(coeff::T) where {T <: Number}
    return Polynomial((coeff,))
end

function Base.eltype(p::Polynomial)
    return eltype(p.coeff)
end

function degree(p::Polynomial)
    return length(p.coeff)-1
end

function vector(p::Polynomial)
    N = degree(p)+1
    v = Vector{eltype(p)}(undef, N)
    for ii = 1:N
        v[ii] = p.coeff[ii]
    end
    return v
end

function Base.:+(p1::Polynomial, p2::Polynomial)
    N1 = degree(p1)+1
    N2 = degree(p2)+1
    N3 = max(N1, N2)
    T1 = eltype(p1)
    T2 = eltype(p2)
    T3 = promote_type(T1, T2)

    coeff = zeros(T3, N3)
    for ii = 1:N1
        coeff[ii] += p1.coeff[ii]
    end
    for ii = 1:N2
        coeff[ii] += p2.coeff[ii]
    end

    p3 = Polynomial(coeff)
    return p3
end

function Base.:-(p1::Polynomial, p2::Polynomial)
    N1 = degree(p1)+1
    N2 = degree(p2)+1
    N3 = max(N1, N2)
    T1 = eltype(p1)
    T2 = eltype(p2)
    T3 = promote_type(T1, T2)

    coeff = zeros(T3, N3)
    for ii = 1:N1
        coeff[ii] += p1.coeff[ii]
    end
    for ii = 1:N2
        coeff[ii] -= p2.coeff[ii]
    end

    p3 = Polynomial(coeff)
    return p3
end

function Base.:*(p1::Polynomial, p2::Polynomial)

    d1 = degree(p1)
    d2 = degree(p2)
    N1 = d1 + 1
    N2 = d2 + 1
    T1 = eltype(p1)
    T2 = eltype(p2)
    T = promote_type(T1, T2)
    N = d1 + d2 + 1
    v = Vector{T}(undef, N)

    if N1 <= N2
        v1 = p1.coeff
        v2 = p2.coeff
    else
        v1 = p2.coeff
        v2 = p1.coeff
        N1, N2 = N2, N1
    end

    for ii = 1:d1
        sum = 0.0
        for jj = 1:ii
            sum += v1[ii - jj + 1] * v2[jj]
        end
        v[ii] = sum
    end

    for ii = N1:-1:2
        sum = 0.0
        for jj = ii:N1
            sum += v1[jj] * v2[ii - jj + N2]
        end
        v[ii - N1 + N] = sum
    end

    for ii = N1:N2
        sum = 0.0
        for jj = 1:N1
            sum += v1[N1 - jj + 1] * v2[jj + ii - N1]
        end
        v[ii] = sum
    end

    q = Polynomial(v)
    return q
end

function Base.:^(p::Polynomial, n::Int)
    if n < 0
        error("This is not a polynomial any more")
    end
    if n == 0
        q = Polynomial(1.0)
        return q
    end

    q = p
    for ii = 1:(n-1)
        q = q * p
    end
    return q
end

function integrate(p::Polynomial, C = 0.0)
    # int x^m = x^(m+1) / (m+1)

    v = [C, vector(p)...]

    for ii = 2:length(v)
        v[ii] /= ii - 1
    end

    q = Polynomial(v)
    return q
end

function derivative(p::Polynomial)
    # der x^m = m * x^(m-1)

    v = vector(p)
    v = v[2:end]

    for ii = eachindex(v)
        v[ii] *= ii
    end

    q = Polynomial(v)
    return q
end

function highest_nonzero_degree(p::Polynomial; abs_tol = 0.0, rel_tol = eps(1.0))
    d = degree(p)
    max_degree = 0
    for ii = (d+1):-1:1
        if !isapprox(p.coeff[ii], zero(eltype(p)), atol = abs_tol, rtol = rel_tol)
            max_degree = ii - 1
            break
        end
    end

    return max_degree
end

function companion_matrix(p::Polynomial)
    hd = highest_nonzero_degree(p)
    m = zeros(hd, hd)
    for ii = 1:(hd-1)
        m[ii+1, ii] = 1.0
    end
    for ii = 1:hd
        m[1, ii] = - p.coeff[hd + 1 - ii] / p.coeff[hd+1]
    end

    return m
end

function roots(p::Polynomial)
    m = companion_matrix(p)
    F = eigen(m)
    z = F.values
    return z
end

function poly_from_roots(roots::AbstractArray; make_real::Bool = true)
    t = tuple(roots...)
    return poly_from_roots(t; make_real = make_real)
end

function poly_from_roots(roots::Tuple; make_real::Bool = true)
    L = length(roots)
    q = Polynomial(1.0)
    for ii = 1:L
        p = Polynomial( (-roots[ii], one(eltype(roots[1]))) )
        # q = q * p
        q = p * q
    end
    if make_real
        q_real_or_imag = Polynomial( map(real, q.coeff) )
    else
        q_real_or_imag = q
    end
    return q_real_or_imag
end

function poly_eval(p::Polynomial, x)
    return poly_eval(x, p.coeff)
end

function (p::Polynomial)(x)
    return poly_eval(p, x)
end

function poly_eval(x, a0, a1, a2, a3)
    return ((a3 .* x .+ a2) .* x + a1) .* x .+ a0
end

function poly_eval(x, a0, a1, a2)
    return (a2 .* x .+ a1) .* x + a0
end

function poly_eval(x, a0, a1)
    return a1 .* x .+ a0
end

function poly_eval(x::AbstractArray, p)
    # p array polynomial coeffiecents
    # p[1] + p[2] * x + p[3] * x^2 + ...
    y = zeros( length(x) )
    np = length(p)
    if np == 0
        return y
    end
    if np == 1
        # return y .* p[1]
        return p[1] .* ones(np)
    end

    y .= p[np]
    for ii = np:-1:2
        y = y .* x .+ p[ii-1]
    end

    return y
end

function poly_eval(x::T, p) where {T <: Number}
    y = zero(T)
    np = length(p)
    if np == 0
        return y
    end
    if np == 1
        return p[1]
    end

    y = p[np]
    for ii = np:-1:2
        y = y * x + p[ii-1]
    end

    return y
end

end

module ptest
using ..poly

# coeff = (1.0, 2.0, 3.0)
# # t_coeff = tuple(map(Float64, v_coeff)...)

# p_v = poly.Polynomial(v_coeff)
# p_t = poly.Polynomial(t_coeff)

coeff = (1.0, 2.0, 3.0, 5.0)
p = poly.Polynomial(coeff)
int_p = poly.integrate(p)
der_p = poly.derivative(p)

coeff1 = (3.0, 2.0, 4.0)
p1 = poly.Polynomial(coeff1)

coeff2 = (5.0, 6.0, 3.0, 7.0)
p2 = poly.Polynomial(coeff2)

# q = p1 * p2

# p1_sq = p1^2
# p2_sq = p2^2

# z = poly.roots(p1)
# q = poly.poly_from_roots(z)

z = [-3.0, 1.0, 5.0, 2.0 - im, 2.0 + im]
q = poly.poly_from_roots(z)
z_c = poly.roots(q)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Misc. Base and stdlib functions supports

import Base: cumsum, cumprod, sum, all, any
import Statistics: mean, std, var

const _tsmap  = Dict{Any,Dict{Number,Expr}}()  # timestamp map
const _colmap = Dict{Any,Dict{Number,Expr}}()  # colanmes map

"""
To handle :where clause
"""
_get_fname(sig::Expr) =
    (sig.head == :call) ?  sig.args[1] : _get_fname(sig.args[1])

macro _mapbase(sig::Expr, imp::Expr)
    fname = _get_fname(sig)

    # these default values are useful for reduction function
    ts = get(_tsmap, fname, Dict())
    dim1ts = get(ts, 1, :(timestamp(ta)[end]))
    dim2ts = get(ts, 2, :(timestamp(ta)))

    # these default values are useful for reduction function
    col = get(_colmap, fname, Dict())
    dim1col = get(col, 1, :(colnames(ta)))
    dim2col = get(col, 2, :([$(QuoteNode(fname))]))

    fbody = quote
        if dims == 1
            TimeArray($dim1ts, $imp, $dim1col, meta(ta))
        elseif dims == 2
            TimeArray($dim2ts, $imp, $dim2col, meta(ta))
        else
            throw(DimensionMismatch("dims should be 1 or 2"))
        end
    end

    doc = "   $sig"
    fdef = Expr(:function, sig, fbody)

    esc(quote
        $doc
        $fdef
    end)
end

# Cumulative functions
_tsmap[:cumsum] = Dict(1 => :(timestamp(ta)))
_colmap[:cumsum] = Dict(2 => :(colnames(ta)))
@_mapbase cumsum(ta::TimeArray; dims::Integer) cumsum(values(ta), dims = dims)
@_mapbase(cumsum(ta::TimeArray{T,1}, dims::Integer = 1) where{T},
          cumsum(values(ta), dims = dims))

_tsmap[:cumprod] = Dict(1 => :(timestamp(ta)))
_colmap[:cumprod] = Dict(2 => :(colnames(ta)))
@_mapbase cumprod(ta::TimeArray; dims::Integer) cumprod(values(ta), dims = dims)
@_mapbase(cumprod(ta::TimeArray{T,1}; dims::Integer = 1) where {T},
          cumprod(values(ta), dims = dims))

# Reduction functions
@_mapbase sum(ta::TimeArray; dims = 1) sum(values(ta), dims = dims)
@_mapbase all(ta::TimeArray; dims = 1) all(values(ta), dims = dims)
@_mapbase any(ta::TimeArray; dims = 1) any(values(ta), dims = dims)

@_mapbase mean(ta::TimeArray; dims = 1) mean(values(ta), dims = dims)
@_mapbase std(ta::TimeArray; dims = 1, kw...) std(values(ta); dims = dims, kw...)
@_mapbase var(ta::TimeArray; dims = 1, kw...) var(values(ta); dims = dims, kw...)


text:
endline = Sys.iswindows() ? "\r\n" : "\n"

@testset "formatted output" begin
    _, out = @grab_output @on_device @cuprintf("")
    @test out == ""

    _, out = @grab_output @on_device @cuprintf("Testing...\n")
    @test out == "Testing...$endline"

    # narrow integer
    _, out = @grab_output @on_device @cuprintf("Testing %d %d...\n", Int32(1), Int32(2))
    @test out == "Testing 1 2...$endline"

    # wide integer
    _, out = @grab_output if Sys.iswindows()
        @on_device @cuprintf("Testing %lld %lld...\n", Int64(1), Int64(2))
    else
        @on_device @cuprintf("Testing %ld %ld...\n", Int64(1), Int64(2))
    end
    @test out == "Testing 1 2...$endline"

    _, out = @grab_output @on_device begin
        @cuprintf("foo")
        @cuprintf("bar\n")
    end
    @test out == "foobar$endline"

    # c argument promotions
    function kernel(A)
        @cuprintf("%f %f\n", A[1], A[1])
        return
    end
    x = CuArray(ones(2, 2))
    _, out = @grab_output begin
        CUDA.@sync @cuda kernel(x)
    end
    @test out == "1.000000 1.000000$endline"
end

@testset "@cuprint" begin
    # basic @cuprint/@cuprintln

    _, out = @grab_output @on_device @cuprint("Hello, World\n")
    @test out == "Hello, World$endline"

    _, out = @grab_output @on_device @cuprintln("Hello, World")
    @test out == "Hello, World$endline"

    # argument interpolation (by the macro, so can use literals)

    _, out = @grab_output @on_device @cuprint("foobar")
    @test out == "foobar"

    _, out = @grab_output @on_device @cuprint(:foobar)
    @test out == "foobar"

    _, out = @grab_output @on_device @cuprint("foo", "bar")
    @test out == "foobar"

    _, out = @grab_output @on_device @cuprint("foobar ", 42)
    @test out == "foobar 42"

    _, out = @grab_output @on_device @cuprint("foobar $(42)")
    @test out == "foobar 42"

    _, out = @grab_output @on_device @cuprint("foobar $(4)", 2)
    @test out == "foobar 42"

    _, out = @grab_output @on_device @cuprint("foobar ", 4, "$(2)")
    @test out == "foobar 42"

    _, out = @grab_output @on_device @cuprint(42)
    @test out == "42"

    _, out = @grab_output @on_device @cuprint(4, 2)
    @test out == "42"

    # bug: @cuprintln failed to invokce @cuprint with endline in the case of interpolation
    _, out = @grab_output @on_device @cuprintln("foobar $(42)")
    @test out == "foobar 42$endline"

    # argument types

    # we're testing the generated functions now, so can't use literals
    function test_output(val, str)
        canary = rand(Int32) # if we mess up the main arg, this one will print wrong
        _, out = @grab_output @on_device @cuprint(val, " (", canary, ")")
        @test out == "$(str) ($(Int(canary)))"
    end

    for typ in (Int16, Int32, Int64, UInt16, UInt32, UInt64)
        test_output(typ(42), "42")
    end

    for typ in (Float32, Float64)
        test_output(typ(42), "42.000000")
    end

    test_output(Cchar('c'), "c")

    for typ in (Ptr{Cvoid}, Ptr{Int})
        ptr = convert(typ, Int(0x12345))
        test_output(ptr, Sys.iswindows() ? "0000000000012345" : "0x12345")
    end

    test_output(true, "1")
    test_output(false, "0")

    test_output((1,), "(1,)")
    test_output((1,2), "(1, 2)")
    test_output((1,2,3.), "(1, 2, 3.000000)")

    # escaping

    kernel1(val) = (@cuprint(val); nothing)
    _, out = @grab_output @on_device kernel1(42)
    @test out == "42"

    kernel2(val) = (@cuprintln(val); nothing)
    _, out = @grab_output @on_device kernel2(42)
    @test out == "42$endline"
end

@testset "@cushow" begin
    function kernel()
        seven_i32 = Int32(7)
        three_f64 = Float64(3)
        @cushow seven_i32
        @cushow three_f64
        @cushow 1f0 + 4f0
        return
    end

    _, out = @grab_output @on_device kernel()
    @test out == "seven_i32 = 7$(endline)three_f64 = 3.000000$(endline)1.0f0 + 4.0f0 = 5.000000$(endline)"
end

@testset "@cushow array pointers" begin
    function kernel()
        a = CuStaticSharedArray(Float32, 1)
        b = CuStaticSharedArray(Float32, 2)
        @cushow pointer(a) pointer(b)
        return
    end

    _, out = @grab_output @on_device kernel()
    @test occursin("pointer(a) = ", out)
    @test occursin("pointer(b) = ", out)
    @test occursin("= 0", out)  # 0x... on Linux, 0... on Windows
end