text:
using Visualize, ModernGL
using Visualize.GLRasterization
using FileIO, Images, GeometryTypes

include(Visualize.dir("src", "julia","mesh.jl"))
catmesh = normalmesh(load(Pkg.dir("GLVisualize", "assets", "cat.obj")))
catmesh = Base.view(reinterpret(MeshRender.Vertex, catmesh.parent), catmesh.indexes[1])

proj = perspectiveprojection(42f0, 1f0, 0.1f0, 10f0)
view = lookat(Vec3f0(1.5), Vec3f0(0.5), Vec3f0(0, -1, 0))
solid_color = Vec3f0(0.0, 0.7, 1.0)

light = MeshRender.Light(
    Vec3f0(10),
    Vec3f0(0.0),
    Vec3f0(0.9),
    1.0f0,
    Vec3f0(0.2),
    1.0f0
)
shading = MeshRender.Shading(
    Vec3f0(0.0),
    Vec3f0(0.3),
    8.0f0
)
program_args = (light, shading, (solid_color, proj, view))
events = WindowEvents()
window = events[Window]

uniforms = map(program_args) do x
    UniformBuffer(x)
end;
vbo = VertexArray(catmesh)

draw_cat = GLRasterizer(
    vbo, uniforms,
    MeshRender.vertexshader, MeshRender.fragmentshader
)

glEnable(GL_DEPTH_TEST)
glEnable(GL_CULL_FACE)
glClearColor(1, 1, 1, 0)
GLAbstraction.enabletransparency()

while isopen(window)
    GLWindow.poll_glfw()
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    draw_cat(vbo, uniforms)
    GLWindow.swapbuffers(window)
end
GLFW.DestroyWindow(window)


text:
###
### LongSequence specific specializations of src/biosequence/transformations.jl
###

"""
    resize!(seq, size, [force::Bool])

Resize a biological sequence `seq`, to a given `size`. Does not resize the underlying data
array unless the new size does not fit. If `force`, always resize underlying data array.
"""
function Base.resize!(seq::LongSequence{A}, size::Integer, force::Bool=false) where {A}
    if size < 0
        throw(ArgumentError("size must be non-negative"))
    elseif seq.shared
        # May seem wasteful to orphan here, but if you resize! a shared seq,
        # you're probably going to write to it anyway.
        return _orphan!(seq, size)
    else
        if force | (seq_data_len(A, size) > seq_data_len(A, length(seq)))
            resize!(seq.data, seq_data_len(A, size))
        end
        seq.part = 1:size
        return seq
    end
end

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

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

    return seq
end

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

"""
    reverse!(seq::LongSequence)

Reverse a biological sequence `seq` in place.
"""
Base.reverse!(seq::LongSequence{<:Alphabet}) = _reverse!(orphan!(seq), BitsPerSymbol(seq))

"""
    reverse(seq::LongSequence)

Create reversed copy of a biological sequence.
"""
Base.reverse(seq::LongSequence{<:Alphabet}) = _reverse(orphan!(seq), BitsPerSymbol(seq))

# Fast path for non-inplace reversion
@inline function _reverse(seq::LongSequence{A}, B::BT) where {A <: Alphabet,
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    cp = LongSequence{A}(unsigned(length(seq)))
    reverse_data_copy!(identity, cp.data, seq.data, B)
    return zero_offset!(cp)
end

_reverse(seq::LongSequence{<:Alphabet}, ::BitsPerSymbol) = reverse!(copy(seq))

# Generic fallback
function _reverse!(seq::LongSequence{<:Alphabet}, ::BitsPerSymbol)
    i, j = 1, lastindex(seq)
    @inbounds while i < j
        seq[i], seq[j] = seq[j], seq[i]
        i += 1
        j -= 1
    end
    return seq
end

@inline function _reverse!(seq::LongSequence{<:Alphabet}, B::BT) where {
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    reverse_data!(identity, seq.data, B)
    return zero_offset!(seq)
end

# Reversion of chunk bits may have left-shifted data in chunks, this function right shifts
# all chunks by up to 63 bits. Only works on orphan sequences.
# This is written so it SIMD parallelizes - careful with changes
@inline function zero_offset!(seq::LongSequence{A}) where A <: Alphabet
    lshift = offset(bitindex(seq, last(seq.part)) + bits_per_symbol(A()))
    rshift = 64 - lshift
    len = length(seq.data)
    @inbounds if !iszero(lshift)
        this = seq.data[1]
        for i in 1:len-1
            next = seq.data[i+1]
            seq.data[i] = (this >>> (unsigned(rshift) & 63)) | (next << (unsigned(lshift) & 63))
            this = next
        end
        seq.data[len] >>>= (unsigned(rshift) & 63)
    end
    return seq
end

# Reverse chunks in data vector and each symbol within a chunk. Chunks may have nonzero
# offset after use, so use zero_offset!
@inline function reverse_data!(pred, data::Vector{UInt64}, B::BT) where {
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    len = length(data)
    @inbounds @simd ivdep for i in 1:len >>> 1
        data[i], data[len-i+1] = pred(reversebits(data[len-i+1], B)), pred(reversebits(data[i], B))
    end
    @inbounds if isodd(len)
        data[len >>> 1 + 1] = pred(reversebits(data[len >>> 1 + 1], B))
    end
end

@inline function reverse_data_copy!(pred, dst::Vector{UInt64}, src::Vector{UInt64}, B::BT) where {
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    len = length(dst)
    @inbounds @simd for i in eachindex(dst)
        dst[i] = pred(reversebits(src[len - i + 1], B))
    end
end

"""
    complement!(seq)

Make a complement sequence of `seq` in place.
"""
function complement!(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet}
    orphan!(seq)
    next = index(firstbitindex(seq))
    stop = index(lastbitindex(seq))
    seqdata = seq.data
    @inbounds for i in index(firstbitindex(seq)):index(lastbitindex(seq))
        seqdata[i] = complement_bitpar(seqdata[i], Alphabet(seq))
    end
    return seq
end

function reverse_complement!(seq::LongSequence{<:NucleicAcidAlphabet})
    orphan!(seq)
    pred = x -> complement_bitpar(x, Alphabet(seq))
    reverse_data!(pred, seq.data, BitsPerSymbol(seq))
    return zero_offset!(seq)
end

function reverse_complement(seq::LongSequence{<:NucleicAcidAlphabet})
    cp = typeof(seq)(unsigned(length(seq)))
    pred = x -> complement_bitpar(x, Alphabet(seq))
    reverse_data_copy!(pred, cp.data, seq.data, BitsPerSymbol(seq))
    return zero_offset!(cp)
end

###
### Shuffle
###

function Random.shuffle!(seq::LongSequence)
    orphan!(seq) # TODO: Is this call to orphan nessecery, given setindex should call `orphan!` for us?
    # Fisher-Yates shuffle
    for i in 1:lastindex(seq) - 1
        j = rand(i:lastindex(seq))
        seq[i], seq[j] = seq[j], seq[i]
    end
    return seq
end

ext: jl


text:
"""
    compare_summary(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1}; kwargs...)

# Arguments
- df_lst::Array{DataFrame,1}: List of DataFrames to compare.
    These must all share the same column names.
- `inds::Array{Symbol,1}`: List of indicators that describe each DataFrame and track which
    values/keys are present in each DataFrame. There must be an equal number of input
    DataFrames and indicators.

# Keyword Argument
- `tol::Float64 = 1E-6`: Tolerance used when determining whether values are equal.
    Default values is `1E-6`.
"""
function compare_summary(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1}; tol = 1E-6)
    df_lst = copy.(df_lst)
    N = length(df_lst)

    # Rename columns to indicate which values go with which data set.
    val_0 = intersect(find_oftype.(df_lst, AbstractFloat)...)
    vals = [Symbol.(val_0, :_, ind) for ind in inds]
    cols = setdiff(intersect(propertynames.(df_lst)...), val_0)

    # Print warning if attempting to compare multiple values at once. We're not there yet.
    if length(val_0) > 1
        val_0 = ensurearray(val_0[1])
        @warn("compare_summary currently supports comparing one column of values/DataFrame.")
    end

    # Make all keys lowercase to focus on value comparisons.
    df_lst = [edit_with(df, [Rename.(val_0, val); Drop.(val,0.0,"==")])
        for (df, val) in zip(df_lst, vals)]
    vals = [vals...;]

    # Join all dataframes.
    df = df_lst[1]
    [df = outerjoin(df, df_lst[ii], on = cols) for ii in 2:N]
    [df[!,ind] .= .!ismissing.(df[:,val]) for (ind, val) in zip(inds, vals)]

    # Are all keys equal/present in the DataFrame?
    df[!,:equal_keys] .= prod.(eachrow(df[:,inds]))

    # Are there discrepancies between PRESENT values (within the specified tolerance)?
    # All values in a row x will be considered "equal" if (max(x) - x_i) / mean(x) < tol
    df_comp = (maximum.(skipmissing.(eachrow(df[:,vals]))) .- df[:,vals]) ./
        Statistics.mean.(skipmissing.(eachrow(df[:,vals])))

    df[!,:equal_values] .= all.(skipmissing.(eachrow(df_comp .< tol)))
    df[all.(eachrow(.|(ismissing.(df[:,vals]), df[:,vals].==0))), :equal_values] .= true
    df = df[:,[cols; sort(vals); sort(inds); [:equal_keys, :equal_values]]]

    ii = df[:,:equal_keys] .* df[:,:equal_values]
    df = df[.!ii,:]

    return df
end

"""
    compare_values(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1})

# Arguments
- df_lst::Array{DataFrame,1}: List of DataFrames to compare.
    These must all share the same column names.
- `inds::Array{Symbol,1}`: List of indicators that describe each DataFrame and track which
    values/keys are present in each DataFrame. There must be an equal number of input
    DataFrames and indicators.

# Keyword Argument
- `tol::Float64 = 1E-6`: Tolerance used when determining whether values are equal.
    Default values is `1E-6`.
"""
function compare_values(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1}; tol = 1E-6)
    df_lst = copy.(df_lst)
    df = compare_summary(copy.(df_lst), inds; tol = tol)
    df = df[.!df[:,:equal_values],:]

    size(df,1) > 0 && @warn("Inconsistent values:", df)
    return df
end

"""
    compare_keys(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1})

# Arguments
- df_lst::Array{DataFrame,1}: List of DataFrames to compare.
    These must all share the same column names.
- `inds::Array{Symbol,1}`: List of indicators that describe each DataFrame and track which
    values/keys are present in each DataFrame. There must be an equal number of input
    DataFrames and indicators.
"""
function compare_keys(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1})
    df_lst = copy.(df_lst)
    N = length(inds)
    cols = intersect(find_oftype.(df_lst, Not(AbstractFloat))...)

    # Determine whether to consider case when comparing keys. Only consider case if there
    # instances of the same keys with differing cases in the same DataFrame.
    d_unique = Dict(col => Dict(inds[ii] => sort(unique(df_lst[ii][:,col]))
        for ii in 1:N) for col in cols)
    d_lower = Dict(col => Dict(inds[ii] => lowercase.(d_unique[col][inds[ii]])
        for ii in 1:N) for col in cols);
    CHECKCASE = Dict(col => any(length.(unique.(values(d_lower[col]))) .!==
        length.(values(d_unique[col]))) for col in cols)

    d_all = Dict(col => CHECKCASE[col] ? sort(unique([values(d_unique[col])...;])) :
        sort(unique([values(d_lower[col])...;])) for col in cols)

    df = DataFrame()

    for col in cols
        df_temp = DataFrame(key = fill(col, size(d_all[col])))
        d_check = CHECKCASE[col] ? d_unique[col] : d_lower[col]

        [df_temp[!,ind] = [v in d_check[ind] ? d_unique[col][ind][v .== d_check[ind]][1] :
            missing for v in d_all[col]] for ind in inds]
        df_temp = unique(df_temp[length.(unique.(eachrow(df_temp[:,inds]))) .> 1, :])

        df = [df; df_temp]
    end

    size(df,1) > 0 && @warn("Inconsistent keys:", df)
    return df
end

"""
    benchmark!(d_summ::Dict, k::Symbol, d_bench::Dict, d_calc::Dict;
"""
function benchmark!(d_summ::Dict, k::Symbol, d_bench::Dict, d_calc::Dict;
    tol = 1E-3, small = 1E-7)

    !(k in collect(keys(d_bench))) && return

    df_calc = copy(d_calc[k])
    df_bench = ensurenames(copy(d_bench[k]), propertynames(df_calc))

    # Remove very small numbers. These might be zero or missing in the other DataFrame,
    # and we're not splitting hairs here.
    df_calc = df_calc[abs.(df_calc[:,:value] .> small), :]
    df_bench = df_bench[abs.(df_bench[:,:value] .> small), :]

    println("  Comparing keys and values for ", k)
    df_comp = compare_summary([df_calc, df_bench], [:calc,:bench]; tol = tol)

    k == :utd && (df_comp = edit_with(df_comp, Drop(:yr,2002,"<")))
    k == :utd_new && (df_comp = edit_with(df_comp, Drop(:yr,2002,"<")))

    # If the dataframes are in agreement, store this value as "true".
    # Otherwise, store the comparison dataframe rows that are not in agreement.
    d_summ[k] = size(df_comp,1) == 0 ? true : df_comp
    return d_summ
end

"""
    verify_over(df::DataFrame, col::Any; tol = 1E-6)
"""
function verify_over(df::DataFrame, col::Any; tol = 1E-6)
    df = combine_over(df, col)
    df = df[(df[:,:value] .- 1.0) .> tol, :]
    return size(df,1) == 0 ? true : df
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using GraphicalModelLearning
using Ipopt
using JuMP

using Test
using Random

import DelimitedFiles: readdlm
import LinearAlgebra: diag

const SOLVER = with_optimizer(Ipopt.Optimizer, print_level=0)

include("common.jl")

@testset "GraphicalModelLearning" begin

@testset "factor graphs" begin
    for (name, gm) in gms
        matrix = convert(Array{Float64,2}, gm)
        gm2 = FactorGraph(matrix)
        for key in keys(gm)
            @test isapprox(gm[key], gm2[key])
            if length(key) == 1
                @test isapprox(gm[key], matrix[key..., key...])
            else
                @test isapprox(gm[key], matrix[key...])
            end
        end
    end
end

@testset "gibbs sampler" begin
    for (name, gm) in gms
        gm_tmp = deepcopy(gm)
        gm_tmp.order = 3
        Random.seed!(0) # fix random number generator
        samples = sample(gm_tmp, gibbs_test_samples)
        #base_samples = readdlm("data/$(name)_samples.csv", ',', Int64)
        base_samples = readdlm("data/$(name)_samples.csv", ',')
        @test isapprox(samples, base_samples)
    end
end

@testset "gibbs sampler, 2nd order" begin
    for (name, gm) in gms
        Random.seed!(0) # fix random number generator
        samples = sample(gm, gibbs_test_samples)
        base_samples = readdlm("data/$(name)_samples.csv", ',')
        #println(name)
        #println(base_samples)
        #println(samples)
        #println(abs.(base_samples-samples))
        @test isapprox(samples, base_samples)
    end
end

@testset "gibbs sampler with replicates" begin
    for (name, gm) in gms
        Random.seed!(0) # fix random number generator
        samples_set = sample(gm, gibbs_test_samples, 3)
        for samples in samples_set
            @test sum(samples[:,1]) == gibbs_test_samples
        end
    end
end

@testset "inverse ising formulations" begin

    for (form_name, formulation) in formulations
        @testset "  $(form_name)" begin
            for (name, gm) in gms
                samples = readdlm("data/$(name)_samples.csv", ',')
                Random.seed!(0) # fix random number generator
                #learned_gm = inverse_ising(samples, method=formulation)
                learned_gm = learn(samples, formulation)
                base_learned_gm = readdlm("data/$(name)_$(form_name)_learned.csv", ',')
                #println(maximum(abs.(learned_gm - base_learned_gm)))
                @test isapprox(learned_gm, base_learned_gm)
            end
        end
    end

    samples = readdlm("data/mvt_samples.csv", ',')

    Random.seed!(0) # fix random number generator
    learned_gm_rise = learn(samples, RISE(0.2, false), NLP(SOLVER))
    base_learned_gm = readdlm("data/mvt_RISE_learned.csv", ',')
    #println(abs.(learned_gm_rise - base_learned_gm))
    @test isapprox(learned_gm_rise, base_learned_gm)

    Random.seed!(0) # fix random number generator
    learned_gm_logrise = learn(samples, logRISE(0.2, false), NLP(SOLVER))
    base_learned_gm = readdlm("data/mvt_logRISE_learned.csv", ',')
    #println(abs.(learned_gm_logrise - base_learned_gm))
    @test isapprox(learned_gm_logrise, base_learned_gm)

    Random.seed!(0) # fix random number generator
    learned_gm_rple = learn(samples, RPLE(0.2, false), NLP(SOLVER))
    base_learned_gm = readdlm("data/mvt_RPLE_learned.csv", ',')
    #println(abs.(learned_gm_rple - base_learned_gm))
    @test isapprox(learned_gm_rple, base_learned_gm)
end

accuracy_tests = [
    AccuracyTest(RISE(),     1000, 0.15)
    AccuracyTest(logRISE(),  1000, 0.15)
    AccuracyTest(RPLE(),     1000, 0.15)
    AccuracyTest(RISE(),    10000, 0.05)
    AccuracyTest(logRISE(), 10000, 0.05)
    AccuracyTest(RPLE(),    10000, 0.05)
]

Random.seed!(0) # fix random number generator
@testset "learned model accuracy" begin
    for act in accuracy_tests
        @testset "  $(act.formulation) $(act.samples) $(act.threshold)" begin
            for (gm_name, gm) in gms
                sample_histo = sample(gm, act.samples)
                #learned_gm = inverse_ising(sample_histo, method=act.formulation)
                learned_gm = learn(sample_histo, act.formulation)
                max_error = maximum(abs.(convert(Array{Float64,2}, gm) - learned_gm))
                @test max_error <= act.threshold
            end
        end
    end
end

@testset "inverse multi-body formulations" begin

    for (name, gm) in gms
        samples = readdlm("data/$(name)_samples.csv", ',')
        Random.seed!(0) # fix random number generator
        learned_ising = learn(samples, RISE(0.2, false))
        learned_two_body = learn(samples, multiRISE(0.2, false, 2))

        learned_ising_dict = convert(Dict, learned_ising)
        #println(learned_ising_dict)
        #println(learned_two_body)

        @test length(learned_ising_dict) == length(learned_two_body)
        for (key, value) in learned_ising_dict
            @test isapprox(learned_two_body[key], value, atol=1e-7)
        end
    end

    samples = readdlm("data/mvt_samples.csv", ',')

    rand(0) # fix random number generator
    learned_ising = learn(samples, RISE(0.2, false), NLP(SOLVER))
    learned_two_body = learn(samples, multiRISE(0.2, false, 2), NLP(SOLVER))

    learned_ising_dict = convert(Dict, learned_ising)
    @test length(learned_ising_dict) == length(learned_two_body)
    for (key, value) in learned_ising_dict
        @test isapprox(learned_two_body[key], value)
    end

    for (name, gm) in gms
        gm_tmp = deepcopy(gm)
        gm_tmp.order = 4
        Random.seed!(0) # fix random number generator
        samples = sample(gm_tmp, 10000)

        learned_gm = learn(samples, multiRISE(0.0, false, 4))

        for (key, value) in gm_tmp
            #println(learned_gm[key], value)
            @test isapprox(learned_gm[key], value, atol = 0.15)
        end

        samples = sample(learned_gm, 10000)
        learned_gm2 = learn(samples, multiRISE(0.0, false, 4))

        for (key, value) in gm_tmp
            #println(learned_gm2[key], " ", value)
            # this is a bug, the learned_gm2 value should not be twice as large
            @test isapprox(learned_gm2[key]/2.0, value, atol = 0.15)
        end

    end
end

Random.seed!(0) # fix random number generator
@testset "docs example" begin
    model = FactorGraph([0.0 0.1 0.2; 0.1 0.0 0.3; 0.2 0.3 0.0])
    samples = sample(model, 100000)
    learned = learn(samples)

    err = abs.(convert(Array{Float64,2}, model) - learned)
    @test maximum(err) <= 0.01
end

end


text:
"""
Partial out variables

### Arguments
* `f` : Formula,
* `df` : AbstractDataFrame
* `add_mean` : should intial mean added to the returned variable
* `weight` : Symbol for weight variables.
* `subset` : AbstractVector{Bool} for subsample
* `maxiter` : Maximum number of iterations
* `tol` : tolerance
* `method` : A symbol for the method. Default is :lsmr (akin to conjugate gradient descent). Other choices are :qrfact and :cholfact (factorization methods)

### Returns
* `::DataFrame` : a dataframe with as many columns as there are dependent variables and as many rows as the original dataframe.

### Details
`partial_out` returns the residuals of a set of variables after regressing them on a set of regressors. The syntax is similar to `reg` - but it accepts multiple dependent variables. It returns a dataframe with as many columns as there are dependent variables and as many rows as the original dataframe.
The regression model is estimated on only the rows where *none* of the dependent variables is missing. With the option `add_mean = true`, the mean of the initial variable is added to the residuals.

### Examples
```julia
using  RDatasets, DataFrames, FixedEffectModels, Gadfly
df = dataset("datasets", "iris")
result = partial_out(SepalWidth + SepalLength ~ 1|> Species, df, add_mean = true)
plot(
   layer(result, x="SepalWidth", y="SepalLength", Stat.binmean(n=10), Geom.point),
   layer(result, x="SepalWidth", y="SepalLength", Geom.smooth(method=:lm))
)
```
"""

function partial_out(f::Formula, df::AbstractDataFrame;
                     add_mean = false, weight::Union{Symbol, Void} = nothing,
                     maxiter::Integer = 10000, tol::Real = 1e-8,
                     method::Symbol = :lsmr)

    rf = deepcopy(f)
    (has_absorb, absorb_formula, absorb_terms,
        has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose!(rf)
    if has_iv
        error("partial_out does not support instrumental variables")
    end
    rt = Terms(rf)
    has_weight = weight != nothing
    xf = Formula(nothing, rf.rhs)
    xt = Terms(xf)

    # create a dataframe without missing values & negative weights
    vars = allvars(rf)
    absorb_vars = allvars(absorb_formula)
    all_vars = vcat(vars, absorb_vars)
    all_vars = unique(convert(Vector{Symbol}, all_vars))
    esample = complete_cases(df[all_vars])
    if weight != nothing
        esample &= isnaorneg(df[weight])
    end
    subdf = df[esample, all_vars]
    all_except_absorb_vars = unique(convert(Vector{Symbol}, vars))
    for v in all_except_absorb_vars
        dropUnusedLevels!(subdf[v])
    end

    # Compute weight vector
    sqrtw = get_weight(df, esample, weight)

    # initialize iterations & converged
    iterations = Int[]
    converged = Bool[]

    # Build fixedeffects, an array of AbtractFixedEffects
    if has_absorb
        fixedeffects = FixedEffect(subdf, absorb_terms, sqrtw)
        # in case there is any intercept fe, remove the intercept
        if any([typeof(f.interaction) <: Ones for f in fixedeffects])
            xt.intercept = false
        end
        pfe = FixedEffectProblem(fixedeffects, Val{method})
    else
        pfe = nothing
    end

    # Compute residualized Y
    yf = Formula(nothing, rf.lhs)
    yt = Terms(yf)
    yt.intercept = false
    mfY = ModelFrame2(yt, subdf, esample)
    Y = ModelMatrix(mfY).m
    broadcast!(*, Y, Y, sqrtw)
    if add_mean
        m = mean(Y, 1)
    end
    residualize!(Y, pfe, iterations, converged, maxiter = maxiter, tol = tol)

    # Compute residualized X
    xvars = allvars(xf)
    if length(xvars) > 0 || xt.intercept
        if length(xvars) > 0
            mf = ModelFrame2(xt, subdf, esample)
            X = ModelMatrix(mf).m
        else
            X = fill(one(Float64), (size(subdf, 1), 1))
        end
        broadcast!(*, X, X, sqrtw)
        residualize!(X, pfe, iterations, converged, maxiter = maxiter, tol = tol)
    end

    # Compute residuals
    if length(xvars) > 0 || xt.intercept
        residuals = Y - X * (X \ Y)
    else
        residuals = Y
    end

    # rescale residuals
    broadcast!(/, residuals,  residuals, sqrtw)
    if add_mean
        broadcast!(+, residuals, m, residuals)
    end

    # Return a dataframe
    yvars = convert(Vector{Symbol}, map(string, yt.eterms))
    out = DataFrame()
    j = 0
    for y in yvars
        j += 1
        out[y] = DataArray(Float64, size(df, 1))
        out[esample, y] = residuals[:, j]
    end
    return(out)
end


text:
module test_julia_io_read

using Test

buf = IOBuffer()

println(buf, "abc")
seekstart(buf)
@test "abc\n" == String(read(buf))
@test "" == String(read(buf))

println(buf, "xyz")
seekstart(buf)
@test "abc\nxyz\n" == String(read(buf))

end # module test_julia_io_read

module test_julia_io_skipchars

using Test

buf = IOBuffer()

println(buf, "\r\nabc--")
seekstart(buf)
@test position(buf) == 0

skipchars(==(Char(UInt8('\r'))), buf)
skipchars(==('\n'), buf)

@test "abc" == String(readuntil(buf, "--"))

seekend(buf)
@test eof(buf)

end # module test_julia_io_skipchars


text:
using DashBootstrapComponents, DashHtmlComponents

button_groups = html_div([
    dbc_buttongroup(
        [dbc_button("Left"), dbc_button("Middle"), dbc_button("Right")],
        size="lg",
        className="mr-1",
    ),
    dbc_buttongroup(
        [dbc_button("Left"), dbc_button("Middle"), dbc_button("Right")],
        size="md",
        className="mr-1",
    ),
    dbc_buttongroup(
        [dbc_button("Left"), dbc_button("Middle"), dbc_button("Right")],
        size="sm",
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""
Note: you can pass `result` and `block` as same pointers
"""
function AESEncryptBlock!(result::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key) where {B<:AESBlock}
	cache = gen_cache(key)
	AESEncryptBlock!(result, block, key, cache)
end

function AESEncryptBlock(block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key, cache) where {B<:AESBlock}
	result = zeros(UInt8, 16)
	AESEncryptBlock!(result, block, key, cache)
	result
end

function AESEncryptBlock(block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key) where {B<:AESBlock}
	result = zeros(UInt8, 16)
	cache = gen_cache(key)
	AESEncryptBlock!(result, block, key, cache)
	result
end

function AESEncryptBlock!(result::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key::AES128Key, cache::CipherCache) where {B<:AESBlock}
	current = result
	copyto!(current, block)
	K = cache.K
	tmp = cache.tmp

	# Copy Key
	for i in 1:4
		for j in 1:4
			K[i][j] = key[4i+j-4]
		end
	end
	# Start Rounds
	for i in 1:11
		if i > 1
			#GenRoundKey
			copyto!(tmp, K[4])
			for j in 1:4
				tmp[j] = SBOX[tmp[j]+1]
			end
			tmp1 = tmp[1]
			tmp[1] = tmp[2]
			tmp[2] = tmp[3]
			tmp[3] = tmp[4]
			tmp[4] = tmp1
			K[1] .= K[1] .⊻ tmp
			K[1][1] = K[1][1] ⊻ RCON[i - 1]
			K[2] .= K[2] .⊻ K[1]
			K[3] .= K[3] .⊻ K[2]
			K[4] .= K[4] .⊻ K[3]
			# SubBytes
			for j in 1:16
				current[j] = SBOX[current[j]+1]
			end
			# ShiftRows
			tmp1 = current[2]
			current[2] = current[6]
			current[6] = current[10]
			current[10] = current[14]
			current[14] = tmp1

			tmp1 = current[3]
			current[3] = current[11]
			current[11] = tmp1
			tmp1 = current[7]
			current[7] = current[15]
			current[15] = tmp1

			tmp1 = current[16]
			current[16] = current[12]
			current[12] = current[8]
			current[8] = current[4]
			current[4] = tmp1
			# MixColumns
			if i != 11
				for j in 1:4
					tmp[1] = aes_mul(current[4j-3], 0x2) ⊻ aes_mul(current[4j-2], 0x3) ⊻ aes_mul(current[4j-1], 0x1) ⊻ aes_mul(current[4j], 0x1)
					tmp[2] = aes_mul(current[4j-3], 0x1) ⊻ aes_mul(current[4j-2], 0x2) ⊻ aes_mul(current[4j-1], 0x3) ⊻ aes_mul(current[4j], 0x1)
					tmp[3] = aes_mul(current[4j-3], 0x1) ⊻ aes_mul(current[4j-2], 0x1) ⊻ aes_mul(current[4j-1], 0x2) ⊻ aes_mul(current[4j], 0x3)
					tmp[4] = aes_mul(current[4j-3], 0x3) ⊻ aes_mul(current[4j-2], 0x1) ⊻ aes_mul(current[4j-1], 0x1) ⊻ aes_mul(current[4j], 0x2)
					current[4j-3] = tmp[1]
					current[4j-2] = tmp[2]
					current[4j-1] = tmp[3]
					current[4j  ] = tmp[4]
				end
			end
		end
		# AddRoundKey
		for j in 1:4
			for k in 1:4
				current[4j+k-4] = current[4j+k-4] ⊻ K[j][k]
			end
		end
	end
	current
end

function AESEncryptBlock!(result::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, block::Union{Array{UInt8, 1}, SubArray{UInt8}, B}, key::AES192Key, cache::CipherCache) where {B<:AESBlock}
	current = result
	copyto!(current, block)
	K = cache.K
	tmp = cache.tmp

	# Copy Key
	for i in 1:6 #DIFF
		for j in 1:4
			K[i][j] = key[4i+j-4]
		end
	end
	# Start Rounds
	rid = 1
	for i in 1:13 #DIFF
		if i > 1
			#GenRoundKey
			if (i % 3 == 2) || (i % 3 == 1)
				copyto!(tmp, K[6]) #DIFF
				for j in 1:4
					tmp[j] = SBOX[tmp[j]+1]
				end
				tmp1 = tmp[1]
				tmp[1] = tmp[2]
				tmp[2] = tmp[3]
				tmp[3] = tmp[4]
				tmp[4] = tmp1
				K[1] .= K[1] .⊻ tmp
				K[1][1] = K[1][1] ⊻ RCON[rid]
				rid += 1
				K[2] .= K[2] .⊻ K[1]
			end
			if (i % 3 == 0) || (i % 3 == 1)
				K[3] .= K[3] .⊻ K[2]
				K[4] .= K[4] .⊻ K[3]
				K[5] .= K[5] .⊻ K[4] #DIFF
				K[6] .= K[6] .⊻ K[5] #DIFF
			end
			# SubBytes
			for j in 1:16
				current[j] = SBOX[current[j]+1]
			end
			# ShiftRows
			tmp1 = current[2]
			current[2] = current[6]
			current[6] = current[10]
			current[10] = current[14]
			current[14] = tmp1

			tmp1 = current[3]
			current[3] = current[11]
			current[11] = tmp1
			tmp1 = current[7]
			current[7] = current[15]
			current[15] = tmp1

			tmp1 = current[16]
			current[16] = current[12]
			current[12] = current[8]
			current[8] = current[4]
			current[4] = tmp1
			# MixColumns
			if i != 13 #DIFF
				for j in 1:4
					tmp[1] = aes_mul(current[4j-3], 0x2) ⊻ aes_mul(current[4j-2], 0x3) ⊻ aes_mul(current[4j-1], 0x1) ⊻ aes_mul(current[4j], 0x1)
					tmp[2] = aes_mul(current[4j-3], 0x1) ⊻ aes_mul(current[4j-2], 0x2) ⊻ aes_mul(current[4j-1], 0x3) ⊻ aes_mul(current[4j], 0x1)
					tmp[3] = aes_mul(current[4j-3], 0x1) ⊻ aes_mul(current[4j-2], 0x1) ⊻ aes_mul(current[4j-1], 0x2) ⊻ aes_mul(current[4j], 0x3)
					tmp[4] = aes_mul(current[4j-3], 0x3) ⊻ aes_mul(current[4j-2], 0x1) ⊻ aes_mul(current[4j-1], 0x1) ⊻ aes_mul(current[4j], 0x2)
					current[4j-3] = tmp[1]
					current[4j-2] = tmp[2]
					current[4j-1] = tmp[3]
					current[4j  ] = tmp[4]
				end
			end
		end
		# AddRoundKey
		if i % 3 == 1
			for j in 1:4
				for k in 1:4
					current[4j+k-4] = current[4j+k-4] ⊻ K[j][k]
				end
			end
		elseif i % 3 == 2
			for j in 1:4
				for k in 1:4
					current[4j+k-4] = current[4j+k-4] ⊻ K[(j + 3) % 6 + 1][k]
				end
			end
		else
			for j in 1:4
				for k in 1:4
					current[4j+k-4] = current[4j+k-4] ⊻ K[j + 2][k]
				end
			end
		end
	end
	current
end


text:
# Partial Pose3 constraints

type PartialPriorRollPitchZ <: IncrementalInference.FunctorSingleton
  rp::Distributions.MvNormal
  z::Distributions.Normal
  partial::Tuple
  PartialPriorRollPitchZ() = new()
  PartialPriorRollPitchZ(rp::MvNormal,z::Normal) = new(rp, z, (3,4,5))
  # PartialPriorRollPitchZ(x1,x2,x3,x4) = new(MvNormal(x1,x2), Normal(x3,x4), (3,4,5))
  # PartialPriorRollPitchZ(rpz::PackedPartialPriorRollPitchZ) = new(MvNormal(d.rpmu, reshapeVec2Mat(d.rpsig,2)),
                                                    # Normal(rpz.zmu, rpz.zsig) )
end
function getSample(pprz::PartialPriorRollPitchZ, N::Int=1)
  return ([rand(pprz.z,N)';rand(pprz.rp,N)], )
end

type PackedPartialPriorRollPitchZ <: IncrementalInference.PackedInferenceType
  rpmu::Vector{Float64}
  rpsig::Vector{Float64}
  zmu::Float64
  zsig::Float64
  PackedPartialPriorRollPitchZ() = new()
  PackedPartialPriorRollPitchZ(x1,x2,x3,x4) = new(x1,x2,x3,x4)
  PackedPartialPriorRollPitchZ(rp::MvNormal,z::Normal) = new(rp.μ, rp.Σ.mat[:], z.μ, z.σ)
  # PackedPartialPriorRollPitchZ(rpz::PartialPriorRollPitchZ) = new(rpz.rp.μ, rpz.rp.Σ.mat[:], rpz.z.μ, rpz.z.σ)
end
function convert(::Type{PartialPriorRollPitchZ}, d::PackedPartialPriorRollPitchZ)
  PartialPriorRollPitchZ( MvNormal(d.rpmu, reshapeVec2Mat(d.rpsig,2)),Normal(d.zmu, d.zsig)  )
end
function convert(::Type{PackedPartialPriorRollPitchZ}, d::PartialPriorRollPitchZ)
  PackedPartialPriorRollPitchZ( d.rp, d.z )
end

function compare(a::PartialPriorRollPitchZ, b::PartialPriorRollPitchZ; tol::Float64=1e-10)
  TP = true
  TP = TP && norm(a.rp.μ-b.rp.μ) < tol
  TP = TP && norm(a.rp.Σ.mat[:]-b.rp.Σ.mat[:]) < tol
  TP = TP && norm(a.z.μ-b.z.μ) < tol
  TP = TP && norm(a.z.σ-b.z.σ) < tol
  TP = TP && norm(collect(a.partial)-collect(b.partial)) < tol
  return TP
end

# Partial pairwise constraint between poses X,Y,Yaw
# ------------------------------------------------------------------------------

type PartialPose3XYYaw <: RoME.BetweenPoses
  xyy::Distributions.MvNormal
  partial::Tuple
  PartialPose3XYYaw() = new()
  PartialPose3XYYaw(xyy::MvNormal) = new(xyy, (1,2,6))
end
function getSample(pxyy::PartialPose3XYYaw, N::Int=1)
  return (rand(pxyy.xyy,N), )
end
function (pxyy::PartialPose3XYYaw)(res::Array{Float64},
                            idx::Int,
                            meas::Tuple{Array{Float64,2}},
                            wXi::Array{Float64,2},
                            wXj::Array{Float64,2}  )
  #
  wXjhat = SE2(wXi[[1;2;6],idx])*SE2(meas[1][:,idx]) #*SE2(pp2.Zij[:,1])*SE2(meas[1][:,idx])
  jXjhat = SE2(wXj[[1;2;6],idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

type PackedPartialPose3XYYaw <: IncrementalInference.PackedInferenceType
  vecZij::Array{Float64,1} # 3translations, 3rotation
  vecCov::Array{Float64,1}
  PackedPartialPose3XYYaw() = new()
  PackedPartialPose3XYYaw(x1::Vector{Float64}, x2::Array{Float64,2}) = new(x1, x2[:])
end
function convert(::Type{PartialPose3XYYaw}, d::PackedPartialPose3XYYaw)
  return PartialPose3XYYaw( Distributions.MvNormal(d.vecZij,
               reshapeVec2Mat(d.vecCov, 3))  )
end
function convert(::Type{PackedPartialPose3XYYaw}, d::PartialPose3XYYaw)
  return PackedPartialPose3XYYaw(d.xyy.μ, d.xyy.Σ.mat )
end

function compare(a::PartialPose3XYYaw, b::PartialPose3XYYaw; tol::Float64=1e-10)
  TP = true
  TP = TP && norm(a.xyy.μ-b.xyy.μ) < tol
  TP = TP && norm(a.xyy.Σ.mat[:]-b.xyy.Σ.mat[:]) < tol
  TP = TP && norm(collect(a.partial)-collect(b.partial)) < tol
  return TP
end

type PartialPose3XYYawNH <: IncrementalInference.FunctorPairwiseNH
  xyy::Distributions.MvNormal
  partial::Tuple
  nullhypothesis::Distributions.Categorical
  PartialPose3XYYawNH() = new()
  PartialPose3XYYawNH(xyy::MvNormal, vh::Vector{Float64}) = new(xyy, (1,2,6),  Distributions.Categorical(vh))
end
function getSample(pxyy::PartialPose3XYYawNH, N::Int=1)
  return (rand(pxyy.xyy,N), )
end
function (pxyy::PartialPose3XYYawNH)(res::Array{Float64},
                            idx::Int,
                            meas::Tuple{Array{Float64,2}},
                            wXi::Array{Float64,2},
                            wXj::Array{Float64,2}  )
  #
  wXjhat = SE2(wXi[[1;2;6],idx])*SE2(meas[1][:,idx]) #*SE2(pp2.Zij[:,1])*SE2(meas[1][:,idx])
  jXjhat = SE2(wXj[[1;2;6],idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

type PackedPartialPose3XYYawNH <: IncrementalInference.PackedInferenceType
  vecZij::Array{Float64,1} # 3translations, 3rotation
  vecCov::Array{Float64,1}
  nullhypothesis::Vector{Float64}
  PackedPartialPose3XYYawNH() = new()
  PackedPartialPose3XYYawNH(x1::Vector{Float64}, x2::Vector{Float64},x3::Vector{Float64}) = new(x1, x2, x3)
  PackedPartialPose3XYYawNH(x1::Vector{Float64}, x2::Array{Float64,2},x3::Vector{Float64}) = new(x1, x2[:], x3)
end
function convert(::Type{PartialPose3XYYawNH}, d::PackedPartialPose3XYYawNH)
  return PartialPose3XYYawNH( Distributions.MvNormal(d.vecZij,
               reshapeVec2Mat(d.vecCov, 3)), d.nullhypothesis  )
end
function convert(::Type{PackedPartialPose3XYYawNH}, d::PartialPose3XYYawNH)
  return PackedPartialPose3XYYawNH(d.xyy.μ, d.xyy.Σ.mat, d.nullhypothesis.p )
end

function compare(a::PartialPose3XYYawNH, b::PartialPose3XYYawNH; tol::Float64=1e-10)
  TP = true
  TP = TP && norm(a.xyy.μ-b.xyy.μ) < tol
  TP = TP && norm(a.xyy.Σ.mat[:]-b.xyy.Σ.mat[:]) < tol
  TP = TP && norm(collect(a.partial)-collect(b.partial)) < tol
  TP = TP && norm(a.nullhypothesis.p-b.nullhypothesis.p) < tol
  return TP
end

#


text:
using StatisticalRethinking, Parameters

struct globe_toss{TO <: AbstractVector}
  # Observations of water
  w::TO
  # Number of tosses
  n::Int
end

function (problem::globe_toss)(theta)
  @unpack w, n = problem
  ll = 0.0
  ll += sum(log.(pdf.(Beta(1, 1), theta[1])))
  ll += sum(log.(pdf.(Binomial(n[1], theta[1]), w)))
  -ll
end

n = 9        # Total number of tosses
w = [6]     # Water in n tosses

llf = globe_toss(w, n);
theta = [0.5]
llf(theta)

# Compute the MAP (maximum_a_posteriori) estimate

x0 = [0.5]
lower = [0.0]
upper = [1.0]

(qmap, opt) = quap(llf, x0, lower, upper)


text:
module Burgers

using DataInterpolations

export FunArr
export FunArr2D
export evolve
export get_timeslice
export strip_boundary_conditions
export integrate_difference1D
export absolute_difference
export array_smoother
export SlopeCalc
export minmod
export linear
export constant

@enum SlopeCalc minmod linear constant

mutable struct FunArr2D{T, U}
	xmin::T
	xmax::T
	dx::T
	Nx::Int64

	tmin::T
	tmax::T
	dt::T
	Nt::Int64

	xdat::Vector{T}
	tdat::Vector{T}
	u2dat::Array{U, 2}

	function FunArr2D{T, U}(_xdat, t_duration, _Nt) where {T, U}
		# test for appropriate sizes
		if length(_xdat) < 2
			error("The first dimension of the array must be at least 2.")
		end
		if _Nt < 2
			error("The second dimension of the array must be at least 2.")
		end

		# the duration of elapsed time must be positive
		if t_duration <= 0.0
			error("The duration of time that passes must be positive.")
		end

		_xmin = _xdat[1]
		_xmax = _xdat[end]
		_dx   = _xdat[2] - _xdat[1]
		_Nx   = length(_xdat)
		_tmin = 0.0
		_tmax = t_duration
		_dt = (_tmax - _tmin)/(_Nt - 1)
		_tdat = range(0.0, stop = t_duration, length = _Nt)
		_u2dat = Array{U, 2}(undef, _Nx, _Nt)

		new(_xmin, _xmax, _dx, _Nx, _tmin, _tmax, _dt, _Nt, _xdat, _tdat, _u2dat)
	end
end

mutable struct FunArr{T, U}
	xmin::T            # x-domain minimum
	xmax::T            # x-domain maximum
	xmin_g::T          # x-domain minimum with ghost
	xmax_g::T          # x-domain maximum with ghost
	Ncells::Int64      # size of the array held, excluding ghost cells
	Nghost::Int64      # number of ghost elements at EACH END of the array
	Nsize::Int64       # the total size of the array, with ghost cells included
	dx::T              # step size of z-data
	xdat::Vector{T}    # domain data
	udat::Vector{U}    # range data

	function FunArr{T, U}(func::Function, _Ncells, _Nghost) where {T, U}
		# func : function at which the xdat array is evaluated
		#      : must have the signature func(::Vector{T})::Vector{T}
		_xmin   = 0.0
		_xmax   = 1.0
		_dx     = (_xmax - _xmin)/(_Ncells)
		_xmin_g = _xmin + T(_dx*(0.5 - _Nghost))
		_xmax_g = _xmax + T(_dx*(_Nghost - 0.5))
		_Nsize  = _Ncells + 2*_Nghost
		_xdat   = range(_xmin_g, stop = _xmax_g, length = _Nsize)
		_udat   = func(_xdat)
		new(_xmin, _xmax, _xmin_g, _xmax_g, _Ncells, _Nghost, _Nsize, _dx, _xdat, _udat)
	end

	function FunArr{T, U}(_udat, _Ncells, _Nghost) where {T, U}
		_xmin   = 0.0
		_xmax   = 1.0
		_dx     = (_xmax - _xmin)/(_Ncells)
		_xmin_g = _xmin + _dx*(0.5 - T(_Nghost))
		_xmax_g = _xmax + _dx*(T(_Nghost) - 0.5)
		_Nsize  = _Ncells + 2*_Nghost
		_xdat   = range(_xmin_g, stop = _xmax_g, length = _Nsize)
		new(_xmin, _xmax, _xmin_g, _xmax_g, _Ncells, _Nghost, _Nsize, _dx, _xdat, _udat)
	end
end

function Base.:+(a::FunArr{T, U}, b::FunArr{T, U})::FunArr{T, U} where {T, U}
	# addition for two FunArr{T, U} instances
	if (a.Ncells != b.Ncells) || (a.Nghost != b.Nghost)
		error("Two FunArr instances must hold vectors of the same size for addition.")
	end

	_udat = a.udat + b.udat
	FunArr{T, U}(_udat, b.Ncells, b.Nghost)
end

function Base.:+(a::Number, b::FunArr{T, U})::FunArr{T, U} where {T, U}
	# addition for two FunArr{T, U} instances
	_udat = a .+ b.udat
	FunArr{T, U}(_udat, b.Ncells, b.Nghost)
end

function Base.:-(a::FunArr{T, U}, b::FunArr{T, U})::FunArr{T, U} where {T, U}
	# subtraction for two FunArr{T, U} instances
	if (a.Ncells != b.Ncells) || (a.Nghost != b.Nghost)
		error("Two FunArr instances must hold vectors of the same size for subtraction.")
	end

	_udat = a.udat - b.udat
	FunArr{T, U}(_udat, b.Ncells, b.Nghost)
end

function Base.:*(a::FunArr{T, U}, b::Number)::FunArr{T, U} where {T, U}
	# multiplication of a FunArr{T, U} instance by a Number
	_udat = a.udat * b
	FunArr{T, U}(_udat, a.Ncells, a.Nghost)
end

function Base.:*(a::Number, b::FunArr{T, U})::FunArr{T, U} where {T, U}
	# multiplication of a FunArr{T, U} instance by a Number
	_udat = a * b.udat
	FunArr{T, U}(_udat, b.Ncells, b.Nghost)
end

function Base.:*(a::FunArr{T, U}, b::FunArr{T, U})::FunArr{T, U} where {T, U}
	# multiplication of a FunArr{T, U} instance by a Number
	_udat = a.udat .* b.udat
	FunArr{T, U}(_udat, b.Ncells, b.Nghost)
end

function Base.:/(a::FunArr{T, U}, b::FunArr{T, U})::FunArr{T, U} where {T, U}
	# multiplication of a FunArr{T, U} instance by a Number
	_udat = a.udat ./ b.udat
	FunArr{T, U}(_udat, b.Ncells, b.Nghost)
end

function Base.:/(a::FunArr{T, U}, b::Number)::FunArr{T, U} where {T, U}
	# multiplication of a FunArr{T, U} instance by a Number
	_udat = a.udat / b
	FunArr{T, U}(_udat, a.Ncells, a.Nghost)
end

###############################################################################
###############################################################################
###############################################################################

function throw_hasNaN(uarr::Vector{T}, msg) where {T}
	if sum(isnan.(uarr)) > 0
		println(uarr)
		error(msg)
	end
end

function physFlux(ufun::Vector{T})::Vector{T} where {T}
	# returns the value of F in the PDE:
	# (dt)u + (dx)F = 0
	return 0.5 * ufun .* ufun
end

function flux_naive(uarrP::Vector{T}, uarrM::Vector{T})::Vector{T} where {T}
	# returns the uncorrected average flux function at a boundary between cells
	return 0.5*(physFlux(uarrP) + physFlux(uarrM))
end

function flux_roe(uarrP::Vector{T}, uarrM::Vector{T})::Vector{T} where {T}
	# returns the average flux function at a boundary between cells
	# with the Roe's method correction, using an arithmetic average
	Fluxfun_naive = 0.5*(physFlux(uarrP) + physFlux(uarrM))
	uavgfun       = 0.5*abs.(uarrP + uarrM)
	ustepfun      = uarrM - uarrP
	return Fluxfun_naive - 0.5 * uavgfun .* ustepfun
end

function get_minmod_plusminus(uarr::Vector{T}) where {T}
	# uses the min-mod function to calculate the interface
	N = length(uarr)

	throw_hasNaN(uarr, "uarr")

	# the forward difference slopes
	sF      = Vector{T}(undef, N)
	sF[end] = uarr[1] - uarr[end]
	for i in 1:(N-1)
		sF[i] = uarr[i+1] - uarr[i]
	end

	# the backward difference slopes
	sB      = Vector{T}(undef, N)
	sB[1]   = uarr[1] - uarr[end]
	for i in 2:N
		sB[i] = uarr[i] - uarr[i-1]
	end

	# minmod-applied slope
	s = Vector{T}(undef, N)
	for i in 1:N
		if sF[i]*sB[i] < 0.0
			s[i] = 0.0
		elseif abs(sF[i]) < abs(sB[i])
			s[i] = sF[i]
		elseif abs(sB[i]) <= abs(sF[i])
			s[i] = sB[i]
		else
			error("This condition shouldn't happen, but just in case.")
		end
	end

	uarr_plus  = uarr + 0.5*s
	uarr_minus = uarr - 0.5*s

	return (uarr_plus, uarr_minus)
end

function get_linear_plusminus(uarr::Vector{T}) where {T}
	# "linear interpolation to interface"
	uarr_plus  = copy(uarr)
	uarr_minus = copy(uarr)
	uarr_plus[1:end-1] = 0.5*(uarr[2:end] + uarr[1:end-1])
	uarr_minus[2:end]  = 0.5*(uarr[2:end] + uarr[1:end-1])
	return (uarr_plus, uarr_minus)
end

function get_constant_plusminus(uarr::Vector{T}) where {T}
	# "linear interpolation to interface"
	uarr_plus  = copy(uarr)
	uarr_minus = copy(uarr)
	return (uarr_plus, uarr_minus)
end

function get_plusminus(uarr::Vector{T}, sc::SlopeCalc) where {T}
	# choose which method of calculating uarr_plus and uarr_minus to use
	if     sc == minmod
		return get_minmod_plusminus(uarr)
	elseif sc == linear
		return get_linear_plusminus(uarr)
	elseif sc == constant
		return get_constant_plusminus(uarr)
	else
		error("Invalid SlopeCalc instance used.")
	end
end

function udot(ufun::FunArr{T, U}, dx::T, sc::SlopeCalc)::FunArr{T, U} where {T, U}
	# used in the calculation of each step of RK2
	# finding ufun_plus and ufun_minus
	uarr_plus, uarr_minus = get_plusminus(ufun.udat, sc)

	Farr = flux_roe(uarr_plus[1:end-1], uarr_minus[2:end])
	Farr = vcat([0], Farr, [0])
	udot_arr = (Farr[1:end-1] - Farr[2:end])/dx
	udot_fun = 1.0*ufun
	udot_fun.udat = udot_arr

	throw_hasNaN(udot_arr, "udot_arr")

	return udot_fun
end

function step_RK2(ufun::FunArr{T, U}, dt::T, dx::T, sc::SlopeCalc)::FunArr{T, U} where {T, U}
	# performs a Runge-Kutta 2 step in time

	# first half-step of RK2
	ufun_half = ufun + 0.5*dt*udot(ufun, dx, sc)
	ufun_half = apply_boundary_conditions(ufun_half)
	throw_hasNaN(ufun_half.udat, "ufun_half.udat")

	# second half-step of RK2
	ufun_full = ufun + dt*udot(ufun_half, dx, sc)
	ufun_full = apply_boundary_conditions(ufun_full)
	throw_hasNaN(ufun_full.udat, "ufun_full.udat")

	return ufun_full
end

function strip_boundary_conditions(ufun::FunArr{T, U})::Vector{T} where {T, U}
	# removes the ghost cells from the ends of a FunArr instance
	Nghost = ufun.Nghost
	uarr   = ufun.udat[1+Nghost:end-Nghost]
	return uarr
end

function apply_boundary_conditions(ufun::FunArr{T, U})::FunArr{T, U} where {T, U}
	# copies ghost cells from one end to the other
	# assuming periodic boundary conditions
	uarr_bc = copy(ufun.udat)
	Nghost  = ufun.Nghost
	uarr_bc[1:Nghost]         = ufun.udat[end-2*Nghost+1:end-Nghost]
	uarr_bc[end-Nghost+1:end] = ufun.udat[Nghost+1:2*Nghost]

	ufun_bc = 1.0*ufun  # copy instead of point? IDK
	ufun_bc.udat[1:Nghost]         = ufun.udat[end-2*Nghost+1:end-Nghost]
	ufun_bc.udat[end-Nghost+1:end] = ufun.udat[Nghost+1:2*Nghost]
	return ufun_bc
end

function evolve(ufun::FunArr{T, U}, Nt::I, duration::T, sc::SlopeCalc, itstop::I) where {T, U, I<:Int}
	# ufun  : initial positions at t = 0
	# Nt    : total number of time steps
	# duration : total amount of time elapsed
	# returns : FunArr2D instance for the u(x, t) grid, and an array of energies

	dx = ufun.dx          # position step
	dt = duration/(Nt-1)  # time step

	Nghost = ufun.Nghost
	xarr  =  ufun.xdat[1+Nghost:end-Nghost]
	u2fun = FunArr2D{T, U}(xarr, duration, Nt)

	u2fun.u2dat[:,1] = strip_boundary_conditions(ufun)

	# propagate forward in time (need to make)
	ufun_next = ufun
	for it in 2:Nt
		ufun_next = step_RK2(ufun_next, dt, dx, sc)
		u2fun.u2dat[:,it] = strip_boundary_conditions(ufun_next)
		if it == itstop
			break
		end
	end

	return u2fun
end

function nearest_lower_index(N::I, xmin::T, xmax::T, x::T) where {T, I<:Int}
	# takes care of the rand() = 0.0 corner case
	if x == xmin
		x += eps(1.0)
	elseif x == xmax
		x -= eps(1.0)
	end
	if !(xmax > x > xmin)
		error("x = $x")
	end
	ratio = (x - xmin)/(xmax - xmin)
	lower_index = ceil(typeof(N), (N-1)*ratio)
	return lower_index
end

function get_timeslice(u2fun::FunArr2D{T, U}, t::T, Ncells::I, Nghost::I)::FunArr{T, U} where {T, U, I<:Int}
	# accepts a u2fun containing a time-evolved grid, and a time-step
	# returns the function u(x) evaluated at the timeslice it

	if !(u2fun.tmin <= t <= u2fun.tmax)
		error("Time out of bounds for this FunArr2D instance.")
	end

	NtL    = nearest_lower_index(u2fun.Nt, u2fun.tmin, u2fun.tmax, t)
	uarr1  = copy(u2fun.u2dat[:,NtL])
	uarr2  = copy(u2fun.u2dat[:,NtL+1])
	tratio = (t - u2fun.tdat[NtL])/u2fun.dt
	udat   = uarr1 .+ (uarr2 .- uarr1)*tratio
	ufun   = FunArr{T, U}(udat, Ncells, 0)

	return ufun
end

###############################################################################
# FUNCTIONS FOR CALCULATING L2 DIFFERENCE BETWEEN FUNCTIONS
###############################################################################

#function absolute_difference(ufun1::FunArr{T, U}, ufun2::FunArr{T, U}, Ndiv::I)::FunArr{T, U} where {T, U, I<:Int}
function absolute_difference(ufun1::FunArr{T, U}, ufun2::FunArr{T, U}, xdata::Vector{T})::FunArr{T, U} where {T, U, I<:Int}
	# returns the absolute difference between the linearly-extrapolated function ufun1 and ufun2
	@assert ufun1.xmin == ufun2.xmin
	@assert ufun1.xmax == ufun2.xmax

	xmin = ufun1.xmin
	xmax = ufun1.xmax
	Nghost = ufun1.Nghost

	interp1 = CubicSpline(ufun1.udat, ufun1.xdat)
	interp2 = CubicSpline(ufun2.udat, ufun2.xdat)

	uarr_diff = Vector{T}(undef, length(xdata))
	for (i, x) in enumerate(xdata)
		u1 = interp1(x)
		u2 = interp2(x)
		if (u1 == nothing) || (u2 == nothing)
			uarr_diff[i] = 0.0
		else
			uarr_diff[i] = abs(u2 - u1)
		end
	end

	return FunArr{T, U}(uarr_diff, length(xdata), 0)
end

function array_smoother(ufun::FunArr{T, U}, Nsmooth::I)::FunArr{T, U} where {T, U, I<:Int}
	# returns an array ufun_smooth where ufun_smooth[i] = avg(ufun, (i-Nsmooth):(i+Nsmooth))

	Nsize = ufun.Nsize
	uarr_smooth = Vector{T}(undef, Nsize)

	for i in 1:Nsize
		if ((i - Nsmooth) < 1) || ((i + Nsmooth) > Nsize)
			uarr_smooth[i] = ufun.udat[i]
			continue
		end

		avg = sum(ufun.udat[i-Nsmooth:i+Nsmooth])/(1 + 2*Nsmooth)
		uarr_smooth[i] = avg
	end

	ufun_smooth = FunArr{T, U}(uarr_smooth, Nsize, 0)
end

end # module

#function interp_regular1D(ufun::FunArr{T, U}, x::T) where {T, U}
#	# performs linear interpolation on a function discretized on a 1D regular grid
#	NxL = nearest_lower_index(ufun.Nsize, ufun.xmin, ufun.xmax, x)
#	x1 = ufun.xmin + (NxL-1)*ufun.dx
#	x2 = ufun.xmin + NxL*ufun.dx
#	u1 = ufun.udat[NxL]
#	u2 = ufun.udat[NxL+1]
#	slope = (u2 - u1)/ufun.dx
#
#	interp_value = u1 + slope*(x-x1)
#	return interp_value
#end
#
#function integrate_difference1D(ufun1::FunArr{T, U}, ufun2::FunArr{T, U}, Ndiv::I) where {T, U, I<:Int}
#	@assert ufun1.xmin == ufun2.xmin
#	@assert ufun1.xmax == ufun2.xmax
#
#	xmin = ufun1.xmin
#	xmax = ufun1.xmax
#	dx   = (xmax-xmin)/T(Ndiv-1)
#
#	total = 0.0
#	for i in 1:(Ndiv-1)
#		x = (i-0.5)*dx
#		term_fun1 = interp_regular1D(ufun1, x)
#		term_fun2 = interp_regular1D(ufun2, x)
#		total    += dx * (term_fun1 - term_fun2)^2
#	end
#
#	return total
#end


text:
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("kplot object has correct attributes") do

    p = kplot(percentchange(cl))

    context("plot attributes are correct") do
        @fact typeof(p.x1)   --> Winston.HalfAxisX
        @fact p.attr[:title] --> "AAPL Close"
    end
end


text:
include("rw_1d.jl")

nsteps = length(ARGS)>0 ? parse(ARGS[1]) : 1000
prob = length(ARGS)>1 ? parse(ARGS[1]) : 0.5

result = randomwalk(nsteps,prob)

for (t,x) in enumerate(result)
    println("$(t-1) $x")
end

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
export solve_JuMP, JuMP_from_SDP_Problem

# for complex to real conversion, see klorel/mathprogcomplex a78593860b662fc1c09dbd3910656b84425729f7

"""
  (primalobj, dualobj) = solve_JuMP(problem, solver, primal, dual; debug, logname, printlog, msk_maxtime, sol_info, optsense)

Calls any JuMP interfaced SDP solver on `problem::SDP_Problem`. Returns the primal and dual objectives if possible.

TODO: update arguments
*Arguments* :
- `problem::SDP_Problem`
- `primal::SortedDict{Tuple{String,String,String}, Float64}`: primal solution `x`,
- `dual::SortedDict{Tuple{String, String, String}, Float64}` : dual solution `s`,
- `debug`: default is false, dump Mosek loaded problem
- `printlog` : if true (default), Mosek will write its log to the console,
- `logname` : if specified, Mosek log will be written to the given file name,
- `msk_maxtime` : Mosek max computation time, in seconds. Default -1 means no limit,
- `sol_info` : Information on problem and solution status upon termination of solve,
- `optsense` : default is `:Max`.

**Note**:
- Mosek expects lower triangular terms of the coefficient matrices. Hence diagonal or non-diagonal terms will not be scaled.
"""
function solve_JuMP(problem::SDP_Problem, solver::T,
                                            primal::SortedDict{Tuple{String,String,String}, Float64},
                                            dual::SortedDict{Tuple{String, String, String}, Float64};
                                            debug = false,
                                            logname = "",
                                            printlog = true,
                                            msk_maxtime = -1,            # Default -1 means no time limit
                                            sol_info = OrderedDict(),
                                            optsense = :Max) where T<:MathProgBase.AbstractMathProgSolver

    empty!(primal)
    empty!(dual)

    m = JuMP_from_SDP_Problem(problem, solver)

    JuMP.setobjectivesense(m, optsense)

    JuMP.solve(m)

    debug && MathProgBase.writeproblem(m.internalModel, "myfile.jtask")
    debug && mv("myfile.jtask", "solve_JuMP_task.json", remove_destination = true)

    primobj = JuMP.getobjectivevalue(m)

    ## TODO: - build up primal, dual solutions...
    ##       - have functions signature converge

    return primobj, primobj
end

function solve_JuMP(problem::SDP_Problem, solver::Symbol,
                                            primal::SortedDict{Tuple{String,String,String}, Float64},
                                            dual::SortedDict{Tuple{String, String, String}, Float64};
                                            debug = false,
                                            logname = "",
                                            printlog = true,
                                            msk_maxtime = -1,            # Default -1 means no time limit
                                            sol_info = OrderedDict(),
                                            optsense = :Max)

    empty!(primal)
    empty!(dual)

    @assert solver in OrderedSet([:MosekSolver, :SCSSolver, :CSDPSolver])
    if solver == :MosekSolver
        options = Any[]

        !printlog && push!(options, (:MSK_IPAR_LOG, 0))

        mysolver = Mosek.MosekSolver(options)

    elseif solver == :SCSSolver
        options = Any[]

        !printlog && push!(options, (:verbose, 0))

        mysolver = SCS.SCSSolver(options)

    elseif solver == :CSDPSolver
        options = Dict{Symbol, Any}()

        !printlog && (options[:printlevel] = 0)

        mysolver = CSDP.CSDPSolver(options)
    end

    m = JuMP_from_SDP_Problem(problem, mysolver)

    JuMP.setobjectivesense(m, optsense)

    JuMP.solve(m)

    debug && MathProgBase.writeproblem(m.internalModel, "myfile.jtask")
    debug && mv("myfile.jtask", "solve_JuMP_task.json", remove_destination = true)

    primobj = JuMP.getobjectivevalue(m)

    ## TODO: - build up primal, dual solutions...
    ##       - have functions signature converge

    return primobj, primobj
end

"""
    m = JuMP_from_SDP_Problem(sdp_pb::SDP_Problem, mysolver::T) where T<:MathProgBase.AbstractMathProgSolver

Given a JuMP interfaced SDP solver, convert the input `SDP_Problem` into a JuMP model `m`.
"""
function JuMP_from_SDP_Problem(sdp_pb::SDP_Problem, mysolver::T) where T<:MathProgBase.AbstractMathProgSolver
    m = JuMP.Model(solver = mysolver)

    ## Variables
    Zi = Dict{String, Array{JuMP.Variable, 2}}()
    for (blockname, block) in sdp_pb.name_to_sdpblock
        n = length(block.var_to_id)
        var = JuMP.@variable(m, [1:n,1:n], SDP, basename=blockname)
        Zi[blockname] = var
    end

    varlin = Dict()
    for varname in keys(sdp_pb.scalvar_to_id)
        varlin[varname] = JuMP.@variable(m, basename=varname)
    end

    ## Constraints

    # Gathering constraint and objective keys
    ctr_names = SortedSet([k[1] for k in keys(sdp_pb.matrices)])
    union!(ctr_names, [k[1] for k in keys(sdp_pb.linear)])
    union!(ctr_names, [k for k in keys(sdp_pb.cst_ctr)])
    ctr_names = SortedSet(collect(ctr_names))

    objective_keys = SortedSet()
    for ctr_name in ctr_names
        if ctr_name[1] == "1" && ctr_name[2] == "1"
            delete!(ctr_names, ctr_name)
            push!(objective_keys, ctr_name)
        end
    end

    obj = 0
    bodys = Dict{Tuple{String, String, String}, JuMP.GenericAffExpr{Float64,JuMP.Variable}}([ctr_name=>0 for ctr_name in ctr_names])

    # Building constraints and objective bodys
    for ((ctr_name, blockname, var1, var2), f_αβ) in sdp_pb.matrices
        block = sdp_pb.name_to_sdpblock[blockname]
        n = length(block.var_to_id)
        i = block.var_to_id[var1]
        j = block.var_to_id[var2]

        Ai = sparse([i], [j], [1], n, n)
        if (i != j) #&& (ctr_name ∉ objective_keys)
            Ai = sparse([i; j], [j; i], [1.; 1.], n, n)
        end

        if ctr_name in ctr_names
            bodys[ctr_name] += vecdot( Zi[blockname], Ai ) * f_αβ
        else
            @assert ctr_name in objective_keys
            obj += vecdot( Zi[blockname], Ai ) * f_αβ
        end
    end

    for ((ctr_name, varname), f_αβ) in sdp_pb.linear
        !haskey(bodys, ctr_name) && (bodys[ctr_name] = 0)

        if ctr_name in ctr_names
            bodys[ctr_name] += varlin[varname] * f_αβ
        else
            @assert ctr_name in objective_keys
            obj += varlin[varname] * f_αβ
        end
    end

    for (ctr_name, f_αβ) in sdp_pb.cst_ctr
        if ctr_name in ctr_names
            bodys[ctr_name] += f_αβ
        else
            @assert ctr_name in objective_keys
            obj += f_αβ
        end
    end

    # Adding constraints and objective to JuMP model
    JuMP.@objective(m, Max, obj)
    JuMP.@constraintref ctrs[1:length(ctr_names)]

    for (i, ctr_name) in enumerate(ctr_names)
        ctrs[i] = JuMP.@constraint(m, bodys[ctr_name] == 0)
    end

    return m

end


text:
using Base: front, tail, setindex

@inbounds begin
    _tfill(t::Tuple{}, x::Unsigned) = iszero(x) ? () : throw(InexactError)
    _tfill(t::NTuple{N,U}, x::Unsigned) where {N, U<:BitUnsigned} =
        (convert(U, x & typemax(U)), _tfill(tail(t), x >> fullshift(U))...)

    _tzero(t::Tuple{}) = ()
    _tzero(t::NTuple{N,U}) where {N, U<:BitUnsigned} = (zero(U), _tzero(tail(t))...)

    _tone(t::Tuple{}) = ()
    _tone(t::NTuple{N,U}) where {N, U<:BitUnsigned} = (one(U), _tzero(tail(t))...)

    _tcmp(t1::Tuple{}, t2::Tuple{}) = 0
    _tcmp(t1::NTuple{N,U}, t2::NTuple{N,U}) where {N, U<:BitUnsigned} =
        t1[N] == t2[N] ? _tcmp(front(t1), front(t2)) : cmp(t1[N], t2[N])

    _tfullleftshift(::Tuple{}) = ()
    _tfullleftshift(t::NTuple{N,U}) where {N, U<:BitUnsigned} = (zero(U), front(t)...)
    _tleftshift(::Tuple{}, n::Unsigned) = ()
    _tleftshift(::Tuple{}, n::Unsigned, carry::BitUnsigned) = ()
    function _tleftshift(t::NTuple{N,U}, n::Unsigned) where {N, U<:BitUnsigned}
        while n >= fullshift(U)
            t = _tfullleftshift(t)
            n -= fullshift(U)
        end
        n == 0 ? t : (t[1] << n, _tleftshift(tail(t), n, t[1] >> (fullshift(U)-n))...)
    end
    _tleftshift(t::NTuple{N,U}, n::Unsigned, carry::U) where {N, U<:BitUnsigned} =
        ((t[1] << n) | carry, _tleftshift(tail(t), n, t[1] >> (fullshift(U)-n))...)

    _tfullrightshift(::Tuple{}) = ()
    _tfullrightshift(t::NTuple{N,U}) where {N, U<:BitUnsigned} = (tail(t)..., zero(U))
    _trightshift(::Tuple{}, n::Unsigned) = ()
    _trightshift(::Tuple{}, n::Unsigned, carry::BitUnsigned) = ()
    function _trightshift(t::NTuple{N,U}, n::Unsigned) where {N, U<:BitUnsigned}
        while n >= fullshift(U)
            t = _tfullrightshift(t)
            n -= fullshift(U)
        end
        n == 0 ? t : (_trightshift(front(t), n, t[N] << (fullshift(U)-n))..., t[N] >> n)
    end
    _trightshift(t::NTuple{N,U}, n::Unsigned, carry::U) where {N, U<:BitUnsigned} =
        (_trightshift(front(t), n, t[N] << (fullshift(U)-n))..., carry | (t[N] >> n))

    _tand(::Tuple{}, ::Tuple{}) = ()
    _tand(t1::NTuple{N,U}, t2::NTuple{N,U}) where {N, U<:BitUnsigned} =
        (t1[1] & t2[1], _tand(tail(t1), tail(t2))...)

    _tor(::Tuple{}, ::Tuple{}) = ()
    _tor(t1::NTuple{N,U}, t2::NTuple{N,U}) where {N, U<:BitUnsigned} =
        (t1[1] | t2[1], _tor(tail(t1), tail(t2))...)

    _tnot(::Tuple{}) = ()
    _tnot(t::Tuple{N,U}) where {N, U<:BitUnsigned} = (~t[1], _tnot(tail(t))...)

    _tadd(::Tuple{}, ::Tuple{}) = ()
    _tadd(::Tuple{}, ::Tuple{}, carry::BitUnsigned) = ()
    function _tadd(t1::NTuple{N,U}, t2::NTuple{N,U}, carry::U = zero(U)) where {N, U<:BitUnsigned}
        s, ncarry = _add(t1[1], t2[1])
        if iszero(carry)
            carry = ncarry
        elseif iszero(ncarry)
            s, carry = _add(s, carry)
        else
            s += carry # cannot overflow
            carry = ncarry
        end
        (s, _tadd(tail(t1), tail(t2), carry)...)
    end

    _tsub(::Tuple{}, ::Tuple{}) = ()
    _tsub(::Tuple{}, ::Tuple{}, borrow::BitUnsigned) = ()
    function _tsub(t1::NTuple{N,U}, t2::NTuple{N,U}, borrow::U = zero(U)) where {N, U<:BitUnsigned}
        s, nborrow = _sub(t1[1], t2[1])
        if iszero(borrow)
            borrow = nborrow
        elseif iszero(nborrow)
            s, borrow = _sub(s, borrow)
        else
            s -= borrow # cannot overflow
            borrow = nborrow
        end
        (s, _tsub(tail(t1), tail(t2), borrow)...)
    end

    _tmul(::Tuple{}, ::Tuple{}) = ()
    _tmul(::Tuple{}, a::U, carry::U = zero(U)) where {U<:BitUnsigned} = ()
    function _tmul(t::NTuple{N,U}, a::U, carry::U = zero(U)) where {N, U<:BitUnsigned}0
        t1, carry1 = _mul(t[1], a) # carry1 <= typemax - 2
        t1, carry2 = _add(t1, carry) # carry2 <= 1
        carry = carry1 + carry2 # cannot overflow
        (t1, _tmul(tail(t), a, carry)...)
    end
    function _tmul(t1::NTuple{N,U}, t2::NTuple{N,U}) where {N, U<:BitUnsigned}
        (all(iszero, t1) || all(iszero, t2)) && return _tzero(t1)
        _tadd(_tmul(t1, t2[1]), (zero(U), _tmul(front(t1), tail(t2))...))
    end

    function _tdivrem(num::NTuple{N,U}, den::U) where {N, U<:BitUnsigned}
        iszero(den) && throw(DivideError())
        if iszero(num[N])
            qa, ra = _tdivrem(front(num), den)
            return (qa..., zero(U)), (ra..., zero(U))
        elseif num[N] >= den
            qN, rN = divrem(num[N], den)
            qb, rb = _tdivrem2((front(num)..., rN), den)
            return (qb..., qN), (rb..., zero(U))
        else
            qc, rc = _tdivrem2(num, den)
            return (qc..., zero(U)), (rc..., zero(U))
        end
    end
    function _tdivrem2(num::NTuple{N,U}, den::U) where {N, U<:BitUnsigned}
        q, r = _divrem(num[N-1], num[N], den)
        tq, tr = _tdivrem2((front(front(num))..., r), den)
        return (tq..., q), (tr..., zero(U))
    end
    _tdivrem2(num::Tuple{U,U}, den::U) where {U<:BitUnsigned} = map(tuple, _divrem(num[1], num[2], den))

    _tdivrem(num::Tuple{}, den::Tuple{})  = (), ()
    function _tdivrem(num::NTuple{N,U}, den::NTuple{N,U}) where {N, U<:BitUnsigned}
        all(iszero, tail(den)) && return _tdivrem(num, den[1])
        c = _tcmp(num, den)
        c < 0 && return (_tzero(num), num)
        c == 0 && return (_tone(num), _tzero(num))

        n = N
        while iszero(den[n])
            n-=1
        end
        s = unsigned(leading_zeros(den[n]))
        v = _tleftshift(den, s)
        u = _tleftshift((num..., zero(U)), s)
        m = N+1
        while iszero(u[m])
            m-=1
        end
        if u[m] > v[n]
            m+=1
        end
        m = m - n - 1

        q = _tzero(v)
        for j = m:-1:0
            q̂ = _estimateq(u[j+n-1],u[j+n],u[j+n+1], v[n-1], v[n])
            vs = _tleftshift((v..., zero(U)), unsigned(j*fullshift(U)))
            p = _tmul(vs, q̂)
            if _tcmp(p, u) > 0
                q̂ -= one(q̂)
                p = _tsub(p, vs)
            end
            u = _tsub(u, p)
            q = Base.setindex(q, q̂, j+1)
        end
        r = _trightshift(front(u), s)
        return q, r
    end
end


text:
module TestGrouping
    using Base.Test
    using DataFrames

    df = DataFrame(a=rep(1:4, 2), b=rep(2:-1:1, 4), c=randn(8))
    #df[6, :a] = NA
    #df[7, :b] = NA

    cols = [:a, :b]

    f(df) = DataFrame(cmax = maximum(df[:c]))

    sdf = sort(df, cols=cols)
    bdf = by(df, cols, f)

    @test bdf[cols] == unique(sdf[cols])

    byf = by(df, :a, df -> DataFrame(bsum = sum(df[:b])))

    @test all(T -> T <: AbstractVector, map(typeof, colwise([sum], df)))
    @test all(T -> T <: AbstractVector, map(typeof, colwise(sum, df)))

    gd = groupby(df, cols)
    ga = map(f, gd)

    @test bdf == combine(ga)

    g(df) = DataFrame(cmax1 = df[:cmax] + 1)
    h(df) = g(f(df))

    @test combine(map(h, gd)) == combine(map(g, ga))
end


text:
module GainPatternsTex

# package code goes here
using GainPatterns
export plot
import PGFPlots: PolarAxis, Plots, save, Axis
export PolarAxis, save, Axis

###########################################################################
# PLOTTING
###########################################################################
# Returns the PolarAxis object containing a plot of the gain pattern
function plot(gp::GainPattern; title=nothing, ymin::Real=0.0, ymax=nothing, showsamples::Bool=false, lastleg::Bool=true, style=nothing, degrees::Bool=false)
	plot([gp], title=title, ymin=ymin, ymax=ymax, showsamples=showsamples, lastleg=lastleg, styles=[style],degrees=degrees)
end

function plot(gp_array::Vector{GainPattern}; title=nothing,ymin::Real=0.0, ymax=nothing, showsamples::Bool=false, lastleg::Bool=true, legendentries=nothing, colors=nothing, styles=nothing, degrees::Bool=false)

	# Create an array with length of angles for each gain pattern
	# Create an array with minimum mean gain for each pattern
	# Initialize the minimum_gain to the lowest of these minimum gains
	num_gp = length(gp_array)
	mingain_array = zeros(num_gp)
	nangles_array = zeros(Int, num_gp)
	for i = 1:num_gp
		nangles_array[i] = length(gp_array[i].angles)
		mingain_array[i] = minimum(gp_array[i].meangains)
	end
	mingain = minimum(mingain_array)

	# Currently, only plot samples for one gain pattern
	# TODO: Don't allow legend entries for showing samples
	#  Or do, but it has to be done well...
	if showsamples && (num_gp == 1)
		gp = gp_array[1]
		for i = 1:nangles_array[1]
			tempmin = minimum(gp.samples[i])
			mingain = (tempmin < mingain ? tempmin : mingain)
		end

		ymin = min(ymin, mingain, 0.0)
		if typeof(ymax) != Void
			ymax -= ymin
			emsg = "ymax is smaller than smallest gain, 0, and specified ymin"
			ymax < ymin ? error(emsg) : nothing
		end

		le = legendentries == nothing ? nothing : legendentries[1]
		st = styles == nothing ? nothing : styles[1]
		gain_plot = plotgains(gp.angles, gp.meangains, ymin, lastleg, le, st)
		plot_array = plotsamples(gp.angles, gp.samples, ymin)
		push!(plot_array, gain_plot)
		pa = PolarAxis(plot_array, ymax=ymax, yticklabel="{\\pgfmathparse{$ymin+\\tick} \\pgfmathprintnumber{\\pgfmathresult}}")
	else
		ymin = min(ymin, mingain, 0.0)
		if typeof(ymax) != Void
			ymax -= ymin
			emsg = "ymax is smaller than smallest gain, 0, and specified ymin"
			ymax < ymin ? error(emsg) : nothing
		end

		# legendentries, styles must be indexable
		if legendentries == nothing
			legendentries = Array(Void, num_gp)
		end
		if styles == nothing
			styles = Array(Void, num_gp)
		end

		# Do we include the degrees
		if degrees
			xtl = "{\$\\pgfmathprintnumber{\\tick}^{\\circ}\$}"
		else
			xtl = nothing
		end

		plot_array = Array(Plots.Linear, num_gp)
		for i = 1:num_gp
			gp = gp_array[i]
			plot_array[i] = plotgains(gp.angles, gp.meangains, ymin, lastleg, legendentries[i], styles[i])
		end
		#xticklabel=$\pgfmathprintnumber{\tick}^\circ$
		pa = PolarAxis(plot_array, title=title, ymax=ymax, yticklabel="{\\pgfmathparse{$ymin+\\tick} \\pgfmathprintnumber{\\pgfmathresult}}", xticklabel=xtl)

	end

	# Return the polar axis object
	return pa

end

# Creates Plots.Linear object (from PGFPlots) with gains vs angles
# Can handle negative values, which radial plots normally cannot.
#
# ymin = minimum y (gain or radial) value. If it is positive, it is ignored.
#  If it is less than the minimum value of gains, also ignored.
function plotgains{T1<:Real,T2<:Real}(angles::Vector{T1}, gains::Vector{T2}, ymin::Real, lastleg::Bool, legendentry, style)

	# Make copies before we make some changes
	# Remember to shift gains by the minimum value to be plotted
	plot_gains = copy(gains) - ymin
	plot_angles = copy(angles)

	# Last point must be same as first point to complete the plot
	# If not, it will be missing a section between last point and first
	if (angles[1] != angles[end]) && lastleg
		push!(plot_angles, plot_angles[1])
		push!(plot_gains, plot_gains[1])
	end

	# Create a linear plot and return it
	#return Plots.Linear(plot_angles, plot_gains, mark="none", style="red,thick")
	return Plots.Linear(plot_angles, plot_gains, mark="none", legendentry=legendentry, style=style)
end

# Returns an array of plots
function plotsamples{T1<:Real,T2<:Real}(angles::Vector{T1}, samples::Vector{Vector{T2}}, ymin::Real)

	# Create an array of plots
	num_angles = length(angles)
	plot_array = Array(Plots.Linear, num_angles)

	# Create a linear plot of the samples for each angle
	for i = 1:num_angles
		# check that
		num_samples = length(samples[i])
		temp_ones = Float64[]
		temp_samples = Float64[]
		for j = 1:num_samples
			if validgain(samples[i][j])
				push!(temp_ones, 1)
				push!(temp_samples, samples[i][j]-ymin)
			end
		end
		plot_array[i] = Plots.Linear(angles[i]*temp_ones, temp_samples, mark="x", style="blue,smooth")
	end

	# Return the array of plots
	return plot_array
end
end # module


text:
module parallel
using Distributed, SharedArrays
using LinearAlgebra, Interpolations, Wavelets, FFTW
using ContinuousWavelets
using SpecialFunctions, LinearAlgebra
using HDF5, JLD
using Plots
# todo remove nScales controls
# the overall container defines a few types we need here
import ScatteringTransform: stParallel, Scattered, ScatteredFull, ScatteredOut, depth, eltypes

include("subsampling.jl")
export sizes, bsplineType, bilinearType
include("modifiedTransforms.jl")
include("Utils.jl")
export calcuateSizes, calculateThinStSizes, createFFTPlans, remoteMultiply,
    createRemoteFFTPlan, computeAllWavelets, plotAllWavelets, getQ
include("nonlinearities.jl")
export spInverse, aTanh, Tanh, ReLU, piecewiseLinear, plInverse
include("transform.jl")
export st, transformMidLayer!, transformFinalLayer!
include("pathMethods.jl")
export pathToThinIndex, MatrixAggrigator, plotCoordinate,
    reshapeFlattened, numberSkipped, logabs, maxPooling, numScales,
    incrementKeeper, numInLayer
include("postProcessing.jl")
export logabs, ReLU, MatrixAggrigator, reshapeFlattened,
    loadSyntheticMatFile, transformFolder, flatten, roll
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

include("./examples/3d/twocubegrids.jl")

# generation of the chain 3-complex of two intersecting cell 2-complexes in 3-space
W, copEV, copFE, copCF = twocubegrids(2,2,2);
# generation of arrays of simplicial complexes: by body, by face, and by edges
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)
# generation and explosion of meshes by body
meshes = GL.GLExplode(V,CVs[2:end],3,3,3,99,.5);
GL.VIEW(meshes);
meshes = GL.GLExplode(V,FVs[2:end],3,3,3,99,.5);
# visualization
GL.VIEW(meshes);


text:
function RealVector()
    return RealVector(())
end

function add(obj::RealVector, arg0::RealVector)
    return jcall(obj, "add", RealVector, (RealVector,), arg0)
end

function append(obj::RealVector, arg0::RealVector)
    return jcall(obj, "append", RealVector, (RealVector,), arg0)
end

function append(obj::RealVector, arg0::jdouble)
    return jcall(obj, "append", RealVector, (jdouble,), arg0)
end

function copy(obj::RealVector)
    return jcall(obj, "copy", RealVector, ())
end

function ebe_divide(obj::RealVector, arg0::RealVector)
    return jcall(obj, "ebeDivide", RealVector, (RealVector,), arg0)
end

function ebe_multiply(obj::RealVector, arg0::RealVector)
    return jcall(obj, "ebeMultiply", RealVector, (RealVector,), arg0)
end

function equals(obj::RealVector, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_dimension(obj::RealVector)
    return jcall(obj, "getDimension", jint, ())
end

function get_distance(obj::RealVector, arg0::RealVector)
    return jcall(obj, "getDistance", jdouble, (RealVector,), arg0)
end

function get_entry(obj::RealVector, arg0::jint)
    return jcall(obj, "getEntry", jdouble, (jint,), arg0)
end

function get_l1_distance(obj::RealVector, arg0::RealVector)
    return jcall(obj, "getL1Distance", jdouble, (RealVector,), arg0)
end

function get_l_inf_distance(obj::RealVector, arg0::RealVector)
    return jcall(obj, "getLInfDistance", jdouble, (RealVector,), arg0)
end

function get_sub_vector(obj::RealVector, arg0::jint, arg1::jint)
    return jcall(obj, "getSubVector", RealVector, (jint, jint), arg0, arg1)
end

function hash_code(obj::RealVector)
    return jcall(obj, "hashCode", jint, ())
end

function is_infinite(obj::RealVector)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::RealVector)
    return jcall(obj, "isNaN", jboolean, ())
end

function map_add_to_self(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapAddToSelf", RealVector, (jdouble,), arg0)
end

function set(obj::RealVector, arg0::jdouble)
    return jcall(obj, "set", void, (jdouble,), arg0)
end

function set_entry(obj::RealVector, arg0::jint, arg1::jdouble)
    return jcall(obj, "setEntry", void, (jint, jdouble), arg0, arg1)
end

function set_sub_vector(obj::RealVector, arg0::jint, arg1::RealVector)
    return jcall(obj, "setSubVector", void, (jint, RealVector), arg0, arg1)
end

function subtract(obj::RealVector, arg0::RealVector)
    return jcall(obj, "subtract", RealVector, (RealVector,), arg0)
end

function to_array(obj::RealVector)
    return jcall(obj, "toArray", Vector{jdouble}, ())
end


text:
module CommonCrawl

using Compat
using AWS
using AWS.S3
using GZip
using HTTPClient
using URIParser
using PublicSuffix

import Base.open

export CrawlCorpus, clear_cache, segments, archives, open, read_entry, read_entries
export HTTPHeaderStats, analyze_header, CorpusStats, analyze_corpus, merge_corpus_stats, print_corpus_stats

if isless(Base.VERSION, v"0.4.0-")
import Base.split
split{T<:String}(str::T, splitter; limit::Integer=0, keep::Bool=true) = split(str, splitter, limit, keep)
end

include("parse.jl")
include("analyze.jl")

end


text:
module MemoryBuffers

using ArrayInterface
using ArrayInterface: AbstractDevice, CPUPointer, CPUTuple, static_length, StaticInt
using Base.Broadcast: Broadcasted
using Base: tail, @propagate_inbounds, AbstractBroadcasted
using LinearAlgebra

include("layouts.jl")
include("buffers.jl")

""" allocate """
function allocate(::Union{CPUPointer,CPUTuple}, x::DynamicMemory{T}) where {T}
    return Vector{T}(undef, length(x))
end
function allocate(::Union{CPUPointer,CPUTuple}, x::StaticMemory{T,L}) where {T,L}
    return StaticBuffer{T,L}()
end
function allocate(::Union{CPUPointer,CPUTuple}, x::FixedMemory{T}) where {T}
    return MutableBuffer{T}(length(x))
end

# FIXME - what should structured_buffer produce
function allocate(d::AbstractDevice, x::ArrayLayout)
    return structured_buffer(allocate(memory_layout(x)), axes(x))
end

""" dereference """
dereference(::CPUTuple, x::StaticBuffer) = x
function dereference(::CPUTuple, x::MutableBuffer)
    return ImmutableBuffer(materialize(CPUTuple(), x.data), x.length)
end
dereference(::CPUPointer, x::StaticBuffer) = x
dereference(::CPUPointer, x::MutableBuffer) = x

# splits (buffer, pointer) up and assigns them symbols say they can be passed to gc_preserve
@generated function preserve(f, args::Tuple{Vararg{Tuple{Any,Any},N}}) where {N}
    blk = Expr(:block)
    buffers = Vecotr{Symbol}(undef, N)
    ptrs = Vecotr{Symbol}(undef, N)
    for i in 1:N
        bufsym = buffers[i] = gensym(:b, i)
        ptrsym = pointers[i] = gensym(:p, i)
        ptrcall = :(@inbounds(args[$i]))
        push!(blk.args, Expr(:(=), Expr(:tuple, bufsym, ptrsym), ptrcall))
    end
    call = append!(Expr(:call, :f), ptrs)
    push!(blk.args, esc(Expr(:gc_preserve, call, buffers...)))
    return blk
end

# FIXME this needs to be figured out with structured_buffer
buffer_pointer(x) = x, pointer(x)

###
### broadcast interface
###
@inline function instantiate(bc::Broadcasted)
    return materialize!(allocate(device(bc), combine_layouts(bc)), bc)
end

# this method needs to change to facilitate iteration or `f` needs to be initialized to be
# in charge of iterating for this to broadcast
function materialize!(dst, bc::Broadcasted)
    preserve(bc.f, buffer_pointer(dst, bc.args...))
    return dereference(dst)
end

###
### indexing interface
###
function unsafe_get_collection(A::AbstractArray{T}, inds) where {T}
    axs = to_axes(A, inds)
    dest = allocate(MemoryLayout{T}(axs))
    if map(Base.unsafe_length, axes(dest)) == map(Base.unsafe_length, axs)
        preserve(buffer_pointer(dest, A)) do d, a
            _unsafe_getindex!(d, a, inds...)
        end
    else
        Base.throw_checksize_error(dest, axs)
    end
    return dereference(dest)
end

end