text:
@testset "LinearAlgebra" begin
    @testset "Constructor + basics" begin
        A = rand(ComplexF32, 12, 12)
        WS = HC2.MatrixWorkspace(A)
        @test WS isa HC2.MatrixWorkspace

        WS.factorized[] = true
        B = randn(ComplexF64, 12, 12)
        HC2.update!(WS, B)
        @test WS.A == B
        @test WS.A !== B
        @test WS.factorized[] == false

        # test indexing etc
        @test WS[2, 1] == B[2, 1]
        @test WS[2] == B[2]
        WS[2, 3] = 5.0
        @test WS[2, 3] == 5.0
        WS[5] = 2.32
        @test WS[5] == 2.32
        @test size(WS) == (12, 12)
        C = randn(ComplexF64, 12, 12)
        WS .= C
        @test WS.A == C

        #test struct array switch
        A = rand(ComplexF32, 30, 30)
        WS = HC2.MatrixWorkspace(A)
    end

    @testset "QR" begin
        A = randn(ComplexF64, 12, 7)
        MW = HC2.MatrixWorkspace(A)
        B = randn(ComplexF64, 12, 7)
        true_qr = LinearAlgebra.qrfactUnblocked!(copy(B))
        MW .= B

        HC2.updated!(MW)

        b = randn(ComplexF64, 12)
        b2 = copy(b)
        x = zeros(ComplexF64, 7)

        ldiv!(x, MW, b)

        @test norm(MW.qr.factors - true_qr.factors) / norm(MW.qr.factors) < 1e-14
        @test norm(b2 - b) == 0
        @test norm(x - (B \ b2)) / norm(x) < 1e-14
    end

    @testset "ldiv" begin
        for n in [3, 13, 31] # test that struct array also works
            A = randn(ComplexF64, n, n)
            b = randn(ComplexF64, n)
            x = zeros(ComplexF64, n)
            WS = HC2.MatrixWorkspace(A)
            HC2.update!(WS, A)
            ldiv!(x, WS, b)
            WS.factorized[] = false
            @test (@allocated ldiv!(x, WS, b)) == 0
            @test (lu(A) \ b) ≈ x rtol = cond(A) * 10
        end
    end

    @testset "Inf-norm estimator / cond" begin
        d_r = rand() * exp10.(range(-6, stop = 6, length = 6))
        D_R = diagm(d_r)
        d_l = rand() * exp10.(range(6, stop = -6, length = 6))
        D_L = diagm(d_l)

        A = randn(6, 6)

        B = A * inv(D_R)
        WB = HC2.MatrixWorkspace(B)
        @test 0.1 ≤ opnorm(inv(B), Inf) / HC2.inverse_inf_norm_est(WB) ≤ 10
        @test 0.1 ≤
              opnorm(inv(B * D_R), Inf) / HC2.inverse_inf_norm_est(WB, nothing, d_r) ≤
              10
        @test 0.1 ≤ opnorm(inv(D_R * B), Inf) / HC2.inverse_inf_norm_est(WB, d_r) ≤ 10

        @test 0.1 ≤ cond(B, Inf) / cond(WB) ≤ 10
        @test 0.1 ≤ cond(B * D_R, Inf) / cond(WB, nothing, d_r) ≤ 10
        @test 0.1 ≤ cond(D_R * B, Inf) / cond(WB, d_r, nothing) ≤ 10

        C = inv(D_L) * A
        WC = HC2.MatrixWorkspace(C)
        @test 0.1 ≤ opnorm(inv(C), Inf) / HC2.inverse_inf_norm_est(WC) ≤ 10
        @test 0.1 ≤ opnorm(inv(D_L * C), Inf) / HC2.inverse_inf_norm_est(WC, d_l) ≤ 10
        @test 0.1 ≤
              opnorm(inv(C * D_L), Inf) / HC2.inverse_inf_norm_est(WC, nothing, d_l) ≤
              10
        @test 0.1 ≤ cond(C, Inf) / cond(WC) ≤ 10
        @test 0.1 ≤ cond(C * D_L, Inf) / cond(WC, nothing, d_l) ≤ 10
        @test 0.1 ≤ cond(D_L * C, Inf) / cond(WC, d_l, nothing) ≤ 10

        D = inv(D_L) * A * inv(D_R)
        WD = HC2.MatrixWorkspace(D)
        @test 0.1 ≤ opnorm(inv(D), Inf) / HC2.inverse_inf_norm_est(WD) ≤ 10
        @test 0.1 ≤
              opnorm(inv(D_L * D * D_R), Inf) / HC2.inverse_inf_norm_est(WD, d_l, d_r) ≤
              10
        @test HC2.inverse_inf_norm_est(WD, d_r, d_l) >
              100 * HC2.inverse_inf_norm_est(WD, d_l, d_r)
        @test 0.1 ≤ cond(D, Inf) / cond(WD) ≤ 10
        @test 0.1 ≤ cond(D_L * D * D_R, Inf) / cond(WD, d_l, d_r) ≤ 10
    end

    @testset "Jacobian" begin
        A = randn(ComplexF64, 6, 6)
        b = randn(ComplexF64, 6)
        x̂ = similar(b)
        x = similar(b)
        JM = HC2.Jacobian(zeros(ComplexF64, 6, 6))
        @test JM isa HC2.Jacobian
        HC2.matrix(JM) .= A
        HC2.updated!(JM)
        ldiv!(x̂, JM, b)
        ldiv!(x, HC2.workspace(JM), b)
        @test x == x̂
    end
end


text:
export discount,
       GAEadvantages!,
       rollout!,
       mse,
       nstep_returns!

const AbsVecOfVec = AbstractVector{<:AbstractVector}

####
#### Discounted N-step returns
####

@inline discount(reward::Real, gamma::Real, t::Integer) = gamma^(t - 1) * reward

nstep_returns!(returns::AbsVec, rewards::AbsVec, gamma) =
    bootstrapped_nstep_returns!(returns, rewards, 0, gamma)

@propagate_inbounds function nstep_returns!(
    returns::AbsVecOfVec,
    rewards::AbsVecOfVec,
    gamma::Real,
)
    @boundscheck begin
        if size(returns) != size(rewards)
            throw(ArgumentError("returns and rewards must have same size"))
        end
    end

    for (rets, rews) in zip(returns, rewards)
        @inbounds nstep_returns!(rets, rews, gamma)
    end
    returns
end

@propagate_inbounds function nstep_returns!(returns::AbsMat, rewards::AbsMat, gamma::Real)
    @boundscheck begin
        if axes(returns) != axes(rewards)
            throw(ArgumentError("returns and rewards must have same axes"))
        end
    end

    @uviews returns rewards for k in axes(returns, 2)
        @inbounds nstep_returns!(view(returns, :, k), view(rewards, :, k), gamma)
    end
    returns
end

@propagate_inbounds function bootstrapped_nstep_returns!(
    returns::AbsVec{<:Real},
    rewards::AbsVec{<:Real},
    terminal_value::Real,
    gamma::Real,
)
    @boundscheck begin
        if axes(returns) != axes(rewards)
            throw(ArgumentError("returns and rewards must have same axes"))
        end
        require_one_based_indexing(returns, rewards)
    end

    DT = eltype(returns)
    gamma, ret_t = DT(gamma), DT(terminal_value)
    @inbounds for t in reverse(eachindex(returns))
        returns[t] = ret_t = rewards[t] + gamma * ret_t
    end
    returns
end

@propagate_inbounds function bootstrapped_nstep_returns!(
    returns::AbsVecOfVec,
    rewards::AbsVecOfVec,
    terminal_values::AbsVec,
    gamma::Real,
)
    @boundscheck begin
        if !(length(returns) == length(rewards) == length(terminal_values))
            throw(ArgumentError("returns, rewards, and terminal_values must have same length"))
        end
    end

    for (rets, rews, term) in zip(returns, rewards, terminal_values)
        @inbounds bootstrapped_nstep_returns!(rets, rews, term, gamma)
    end
    returns
end

@propagate_inbounds function bootstrapped_nstep_returns!(
    returns::AbsMat,
    rewards::AbsMat,
    terminal_values::AbsVec,
    gamma,
)
    @boundscheck begin
        if axes(returns) != axes(rewards) || axes(returns, 2) != axes(terminal_values, 1)
            throw(ArgumentError("returns and rewards must have same axes and share 2nd axis with terminal_values"))
        end
        require_one_based_indexing(returns, rewards, terminal_values)
    end

    @uviews returns rewards for k in axes(returns, 2)
        @inbounds bootstrapped_nstep_returns!(
            view(returns, :, k),
            view(rewards, :, k),
            terminal_values[k],
            gamma,
        )
    end
    returns
end

####
#### Generalized Advantage Estimation (https://arxiv.org/abs/1506.02438)
####

@propagate_inbounds function GAEadvantages!(
    advantages::AbsVec{<:Real},
    values::AbsVec{<:Real},
    rewards::AbsVec{<:Real},
    terminal_value::Real,
    gamma::Real,
    lambda::Real,
)
    @boundscheck begin
        if !(axes(advantages) == axes(values) == axes(rewards))
            throw(BoundsError("advantages, values, and rewards must have same axes"))
        end
        require_one_based_indexing(advantages, values, rewards)
    end
    DT = eltype(advantages)
    gamlam = DT(gamma * lambda)
    @inbounds advantages[end] = rewards[end] + gamma * terminal_value - values[end]
    @inbounds advsum = adv_t = advantages[end] # t=T
    @inbounds for t in reverse(eachindex(advantages)[1:end-1]) # t=(T-1):-1:1
        td_error = rewards[t] + gamma * values[t+1] - values[t]
        advantages[t] = adv_t = td_error + gamlam * adv_t
        advsum += adv_t
    end
    advantages
end

@propagate_inbounds function GAEadvantages!(
    advantages::AbsVecOfVec,
    values::AbsVecOfVec,
    rewards::AbsVecOfVec,
    terminal_values::AbsVec,
    gamma,
    lambda,
)
    @boundscheck if !(size(advantages) == size(values) == size(rewards) ==
                      size(terminal_values))
        throw(BoundsError("advantages, values, rewards, and terminal_values must have same size"))
    end
    for (a, v, r, term) in zip(advantages, values, rewards, terminal_values)
        @inbounds GAEadvantages!(a, v, r, term, gamma, lambda)
    end
    advantages
end

@propagate_inbounds function GAEadvantages!(
    advantages::AbsMat,
    values::AbsMat,
    rewards::AbsMat,
    terminal_values::AbsVec,
    gamma,
    lambda,
)
    @boundscheck begin
        if !(axes(advantages) == axes(values) == axes(rewards) && axes(
            advantages,
            2,
        ) == axes(terminal_values, 1))
            throw(ArgumentError("advantages, values, and rewards must have matching axes and share 2nd axis with terminal_values"))
        end
    end

    @uviews advantages values rewards for k in axes(advantages, 2)
        @inbounds GAEadvantages!(
            view(advantages, :, k),
            view(values, :, k),
            view(rewards, :, k),
            terminal_values[k],
            gamma,
            lambda,
        )
    end
    advantages
end

####
#### Miscellaneous
####

mse(ŷ::AbstractMatrix, y::AbstractMatrix) = sum((ŷ - y) .^ 2) * 1 // length(y)
mse(ŷ::AbstractVector, y::AbstractVector) = sum((ŷ - y) .^ 2) * 1 // length(y)

whiten(x::AbstractArray) = whiten(copy(x))
function whiten!(x::AbstractArray)
    μ, σ = mean(x), std(x)
    x .= (x .- μ) ./ (σ + eps(eltype(x)))
end

struct TimeFeatures{T}
    orders::Vector{T}
    coeffs::Vector{T}
    timestep::T
    function TimeFeatures{T}(orders::AbsVec, coeffs::AbsVec, timestep::Real) where {T<:AbstractFloat}
        if length(orders) != length(coeffs)
            throw(ArgumentError("length(orders) must equal length(coeffs)"))
        end
        timestep > 0 || throw(ArgumentError("timestep must be > 0"))

        orders = convert(Vector{T}, orders)
        coeffs = convert(Vector{T}, coeffs)
        timestep = convert(T, timestep)

        new{T}(orders, coeffs, timestep)
    end
end

function TimeFeatures(orders::AbsVec, coeffs::AbsVec, timestep::Real)
    TimeFeatures{Float64}(orders, coeffs, timestep)
end

function (op::TimeFeatures)(B::AbstractMatrix, A::AbstractMatrix)
    Base.require_one_based_indexing(B, A)
    if size(B, 1) != size(A, 1) + length(op.coeffs)
        throw(DimensionMismatch("size(B, 1) != size(A, 1) + length(op.coeffs)"))
    end
    size(B, 2) != size(A, 2) && throw(DimensionMismatch("size(B, 2) != size(A, 2)"))

    T = eltype(B)
    dimobs, N = size(A)
    dimfeats = length(op.coeffs)
    x = (0:(N - 1)) .* op.timestep

    r1 = 1:dimobs
    r2 = (dimobs + 1):(dimobs + dimfeats)

    for j = axes(B, 2)
        a = uview(A, :, j)
        b1 = uview(B, r1, j)
        b2 = uview(B, r2, j)
        copyto!(b1, a)
        for k = 1:dimfeats
            i = k + dimobs
            B[i, j] = op.coeffs[k] * x[j]^op.orders[k]
        end
    end
    B
end

@inline function (op::TimeFeatures{T})(A::AbstractMatrix) where {T}
    V = promote_type(T, eltype(A))
    B = Matrix{V}(undef, size(A, 1) + length(op.coeffs), size(A, 2))
    op(B, A)
end

# TODO
# 1) inner_size w/ BatchedVector
# 2) return BatchedVector{<:NestedView}
function (op::TimeFeatures{T})(A::AbsVec{<:AbsMat}) where {T}
    for a in A
        if size(a, 1) != size(first(A), 1)
            error("The first axis of each element in A must be equal")
        end
    end

    N = sum(a -> size(a, 2), A)
    dimobs = size(first(A), 1)

    B = Matrix{T}(undef, dimobs + length(op.coeffs), N)
    from = firstindex(B, 2)
    for a in A
        l = size(a, 2)
        to = from + l - 1
        b = view(B, :, from:to)
        op(b, a)
        from += l
    end
    B
end

@inline function (op::TimeFeatures{T})(A::AbsMat{<:Real}, ts::AbsVec{<:Integer}) where {T}
    Base.require_one_based_indexing(A, ts)
    if size(A, 2) != length(ts)
        throw(DimensionMismatch("size(A, 2) != length(ts)"))
    end

    dimobs, N = size(A)
    dimfeats = length(op.coeffs)
    B = Matrix{T}(undef, dimobs + dimfeats, N)

    r1 = 1:dimobs
    r2 = (dimobs + 1):(dimobs + dimfeats)

    for j = axes(B, 2)
        a = uview(A, :, j)
        b1 = uview(B, r1, j)
        b2 = uview(B, r2, j)
        t = (ts[j] - 1) * op.timestep
        copyto!(b1, a)
        for k = 1:dimfeats
            i = k + dimobs
            B[i, j] = op.coeffs[k] * t^op.orders[k]
        end
    end
    B
end

ext: jl


text:
export mseis!

function SeisData(U...)
  S = SeisData()
  for i = 1:length(U)
    Y = getindex(U,i)
    if typeof(Y) == SeisChannel
      push!(S, Y)
    elseif typeof(Y) <: GphysChannel
        push!(S, convert(SeisChannel, Y))
    elseif typeof(Y) == SeisData
      append!(S, Y)
    elseif typeof(Y) <: GphysData
      append!(S, convert(SeisData, Y))
    elseif typeof(Y) == SeisEvent
      append!(S, convert(SeisData, getfield(Y, :data)))
    else
      @warn(string("Tried to join incompatible type into SeisData at arg ", i, "; skipped."))
    end
  end
  return S
end

"""
    mseis!(S::SeisData, U...)

Merge multiple SeisData structures at once. The first argument (merge target)
must be a SeisData structure. Subsequent structures can be any type T <: Union{GphysData, GphysChannel, SeisEvent}.

    mseis!(C::GphysChannel, U...)

Merge all channels in U that match channel C into object C. To be merged, a
channel must match on fields `:id`, `:fs`, `:loc`, `:resp`, `:units`.

See also: `merge!`
"""
function mseis!(S1::SeisData, S...)
  U = Union{GphysData, GphysChannel, SeisEvent}
  L = Int64(length(S))
  (L == 0) && return
  (typeof(S1) == SeisData) || error("Target must be type SeisData!")
  for i in 1:L
    T = typeof(getindex(S, i))
    if (T <: U) == false
      @warn(string("Object of incompatible type passed to mseis! at ", i+1, "; skipped!"))
      continue
    end
    if T == SeisData
      append!(S1, getindex(S, i))
    elseif T == EventTraceData
      append!(S1, convert(SeisData, getindex(S, i)))
    elseif T == SeisChannel
      append!(S1, SeisData(getindex(S, i)))
    elseif T == EventChannel
      append!(S1, SeisData(convert(SeisChannel, getindex(S, i))))
    elseif T == SeisEvent
      append!(S1, convert(SeisData, getfield(getindex(S, i), :data)))
    end
  end
  merge!(S1)
  return S1
end

function mseis!(C::Y, S...) where Y<:GphysChannel
  U = Union{GphysData, GphysChannel, SeisEvent}
  for i in 1:length(S)
    X = S[i]
    T = typeof(X)

    # Only merge compatible types
    if (T <: U) == false
      @warn(string("Object of incompatible type passed to mchan at ", i+1, "; skipped!"))
      continue
    end

    # Only call merge_chan on channels with matching ID
    if T <: GphysChannel
      if channel_match(C, X, use_gain=false)
        merge!(C, X)
      end
    elseif T <: GphysData
      ID = getfield(X, :id)
      for i in 1:X.n
        D = getindex(X, i)
        if channel_match(C, D, use_gain=false)
          merge!(C, D)
        end
      end
    elseif T == SeisEvent
      X = X.data
      for i in 1:X.n
        D = getindex(X, i)
        if channel_match(C, D, use_gain=false)
          merge!(C, D)
        end
      end
    end
  end
  return nothing
end


text:
# functions to create parameters

dx_i/dt = r_i*X_i + sum(F_i(X_j)) - sum(F_j(X_i)) - m_i*X_i

dx_i/dt = fun1(X_i) + fun2(X_i) - fun3(X_i) - fun4(X_i)

function lorenz_added!(dx, x, p, t)
    dx[1] = sigma*(x[2]-x[3])
    dx[2] = x[2]*(rho - x[3]) - x[2]
    dx[3] = x[1]*x[2] - beta*x[3]
end

# dx is just a vector
# populate the vector via a loop through the functions where
# you look at each species
for i in species
    # get abundance from somewhere
    dx[i] = fun1(X_i) + fun2(X_i) - fun3(X_i) - fun4(X_i)
end

function sigma_fun(x::Array, sigma::Float)

    return (sigma*(x[2] - x[1]))

end

x0 = [1.0;0.0;0.0]

# dx1 = r1*x1 + sum(g12*x2, g13*x3) - sum(g21*x1, g31*x1) - m1*x1

x_metrics = [[2, 2], [1, 3], [1, 2.5]]
x_vals = [1, 2, 3]
x0 = [100, 200, 75]

"""
Function defines the density-independent growth of xi.
"""
function growth(r, i)
    return r*x[i]
end

"""
Determines for any focal x (i) and any other x (j) if the focal x will gain
biomass from j based on some arbitrary metric (metric defined in the x_metrics
vector), and if xi does gain from xj, then gives a g value (rate at which the
gain happens)
"""
function get_g(i, j)
    if x_metrics[i][2] > x_metrics[j][2]
        g = 0.5
    else
        g = 0
    end
end

"""
Figures out how much biomass is actually taken from any xj and given to a given
xi.
"""
function add(i, j)
    return get_g(i, j)*x[j]
end

"""
For a specific xi, gets all other x's that it can consume biomass from, gets the
biomass from each, and adds it to xi's biomass.
"""
function add_all(i)
    other_x = filter(x-> x ≠ i, x_vals)
    vals = []
    for j in other_x
        append!(vals, add(i, j))
    end

    return sum(vals)
end

"""
For any given xi, finds the other xj's that can remove biomass from xi, gets how
much biomass each xj will take from xi, and removes all of those values from xi.
"""
function remove_all(i)
    other_x = filter(x -> x ≠ i, x_vals)
    vals = []
    for j in other_x
        append!(vals, add(j, i))
    end

    return sum(vals)
end

function death(i)
    m = 1
    m*x[i]
end
death(i)
function eqns!(dx, x, p, t)
    dx = [[],[],[]]
    for i in x_vals
        dx[i] = growth(all_x[i][1], i) + add_all(i) - remove_all(i) - death(i)
    end
end

using DifferentialEquations
x0 = x
tspan = (0.0, 100.0)
prob = ODEProblem(eqns!, x0, tspan)
sol = solve(prob)

# dx1/dt = r1*x1 + x1(a12*x2) - a21*x1 - m1*x1
# dx2/dt = r2*x2 + x2(a21*x1) - a12*x2 - m2*x2

#dx/dt = indep_growth + dep_growth - dep_death - indep_death

function indep_growth(i, x)
    #define r
    if i == 1
        r = 1.5
        k = 5000
    else
        r = 1.9
        k = 2000
    end
    # return product of r and abundance of focal x
    return (r*x[i]*(1-(x[i]/k)))
end

function dep_growth(i, x)
    # determine if a12 or a21 should be used
    if i == 1 # use a12
        a = 0.010975
        k = 5000
        return x[i]*(a*x[2])*(1-(x[i]/k))
    else
        a = 0.01
        k = 2000
        return x[i]*(a*x[1])*(1-(x[i]/k))
    end
end

function dep_death(i, x)
    if i == 1
        a = 0.01
        return x[i]*a*x[2]
    else
        a = 0.010975
        return x[i]*a*x[1]
    end
    return death
end

function indep_death(i, x)
    if i == 1
        m = 1.4
    else
        m = 1.3
    end

    return m*x[i]
end

dx = zeros(2)
function eqns!(dx, x, p, t)
    for i in 1:2
        dx[i] = indep_growth(i, x) + dep_growth(i, x) - dep_death(i, x) - indep_death(i, x)
    end
end

r1 = 1.5
r2 = 1.9
a12 = 0.010975
a21 = 0.01
m1 = 1.4
m2 = 1.3
k1 = 5000
k2 = 2000
function eqns!(dx, x, p, t)
    dx[1] = r1*x[1]*(1-(x[1]/k1)) + x[1]*(a12*x[2])*(1-(x[1]/k1)) -
            x[1]*(a21*x[2]) - m1*x[1]
    dx[2] = r2*x[2]*(1-(x[2]/k2)) + x[2]*(a21*x[1])*(1-(x[2]/k2)) -
            x[2]*a12*x[1] - m2*x[2]
end

x0 = [100.0; 1000.0]
tspan = (0.0, 10.0)
prob = ODEProblem(eqns!, x0, tspan)
sol2 = solve(prob)
plot(sol2)

sol1[:] == sol2[:]

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
### Types

immutable SparseVector{Tv,Ti<:Integer} <: AbstractSparseVector{Tv,Ti}
    n::Int              # the number of elements
    nzind::Vector{Ti}   # the indices of nonzeros
    nzval::Vector{Tv}   # the values of nonzeros

    function SparseVector(n::Integer, nzind::Vector{Ti}, nzval::Vector{Tv})
        n >= 0 || throw(ArgumentError("The number of elements must be non-negative."))
        length(nzind) == length(nzval) ||
            throw(DimensionMismatch("The lengths of nzind and nzval are inconsistent."))
        new(convert(Int, n), nzind, nzval)
    end
end

SparseVector{Tv,Ti}(n::Integer, nzind::Vector{Ti}, nzval::Vector{Tv}) =
    SparseVector{Tv,Ti}(n, nzind, nzval)

### Basic properties

length(x::SparseVector) = x.n
size(x::SparseVector) = (x.n,)
nnz(x::SparseVector) = length(x.nzval)
countnz(x::SparseVector) = countnz(x.nzval)
nonzeros(x::SparseVector) = x.nzval
nonzeroinds(x::SparseVector) = x.nzind

### Construct empty sparse vector

sparsevector{T}(::Type{T}, len::Integer) = SparseVector(len, Int[], T[])

### Construction from lists of indices and values

function _sparsevector!{Ti<:Integer}(I::Vector{Ti}, V::Vector, len::Integer)
    # pre-condition: no duplicate indexes in I
    if !isempty(I)
        p = sortperm(I)
        permute!(I, p)
        permute!(V, p)
    end
    SparseVector(len, I, V)
end

function _sparsevector!{Tv,Ti<:Integer}(I::Vector{Ti}, V::Vector{Tv}, len::Integer, combine::BinaryOp)
    if !isempty(I)
        p = sortperm(I)
        permute!(I, p)
        permute!(V, p)
        m = length(I)

        # advances to the first non-zero element
        r = 1     # index of last processed entry
        while r <= m
            if V[r] == zero(Tv)
                r += 1
            else
                break
            end
        end
        r > m && SparseVector(len, Ti[], Tv[])

        # move r-th to l-th
        l = 1       # length of processed part
        i = I[r]    # row-index of current element
        if r > 1
            I[l] = i; V[l] = V[r]
        end

        # main loop
        while r < m
            r += 1
            i2 = I[r]
            if i2 == i  # accumulate r-th to the l-th entry
                V[l] = call(combine, V[l], V[r])
            else  # advance l, and move r-th to l-th
                pv = V[l]
                if pv != zero(Tv)
                    l += 1
                end
                i = i2
                if l < r
                    I[l] = i; V[l] = V[r]
                end
            end
        end
        if V[l] == zero(Tv)
            l -= 1
        end
        if l < m
            resize!(I, l)
            resize!(V, l)
        end
    end
    SparseVector(len, I, V)
end

function sparsevector{Tv,Ti<:Integer}(I::AbstractVector{Ti}, V::AbstractVector{Tv}, combine::BinaryOp)
    length(I) == length(V) ||
        throw(DimensionMismatch("The lengths of I and V are inconsistent."))
    len = 0
    for i in I
        i >= 1 || error("Index must be positive.")
        if i > len
            len = i
        end
    end
    _sparsevector!(_copy_convert(Ti, I), _copy_convert(Tv, V), len, combine)
end

function sparsevector{Tv,Ti<:Integer}(I::AbstractVector{Ti}, V::AbstractVector{Tv}, len::Integer, combine::BinaryOp)
    length(I) == length(V) ||
        throw(DimensionMismatch("The lengths of I and V are inconsistent."))
    maxi = convert(Ti, len)
    for i in I
        1 <= i <= maxi || error("An index is out of bound.")
    end
    _sparsevector!(_copy_convert(Ti, I), _copy_convert(Tv, V), len, combine)
end

sparsevector{Ti<:Integer}(I::AbstractVector{Ti}, V::AbstractVector) =
    sparsevector(I, V, AddFun())

sparsevector{Ti<:Integer}(I::AbstractVector{Ti}, V::AbstractVector, len::Integer) =
    sparsevector(I, V, len, AddFun())

sparsevector{Ti<:Integer}(I::AbstractVector{Ti}, v::Number, combine::BinaryOp) =
    sparsevector(I, fill(v, length(I)), combine)

sparsevector{Ti<:Integer}(I::AbstractVector{Ti}, v::Number, len::Integer, combine::BinaryOp) =
    sparsevector(I, fill(v, length(I)), len, combine)

sparsevector{Ti<:Integer}(I::AbstractVector{Ti}, v::Number) =
    sparsevector(I, v, AddFun())

sparsevector{Ti<:Integer}(I::AbstractVector{Ti}, v::Number, len::Integer) =
    sparsevector(I, v, len, AddFun())

### Construction from dictionary

function sparsevector{Tv,Ti<:Integer}(dict::Associative{Ti,Tv})
    m = length(dict)
    nzind = Array(Ti, m)
    nzval = Array(Tv, m)

    cnt = 0
    len = zero(Ti)
    for (k, v) in dict
        k >= 1 || error("Index must be positive.")
        if k > len
            len = k
        end
        if v != zero(v)
            cnt += 1
            @inbounds nzind[cnt] = k
            @inbounds nzval[cnt] = v
        end
    end
    resize!(nzind, cnt)
    resize!(nzval, cnt)
    _sparsevector!(nzind, nzval, len)
end

function sparsevector{Tv,Ti<:Integer}(dict::Associative{Ti,Tv}, len::Integer)
    m = length(dict)
    nzind = Array(Ti, m)
    nzval = Array(Tv, m)

    cnt = 0
    maxk = convert(Ti, len)
    for (k, v) in dict
        1 <= k <= maxk || error("An index (key) is out of bound.")
        if v != zero(v)
            cnt += 1
            @inbounds nzind[cnt] = k
            @inbounds nzval[cnt] = v
        end
    end
    resize!(nzind, cnt)
    resize!(nzval, cnt)
    _sparsevector!(nzind, nzval, len)
end

### Element access

function setindex!{Tv,Ti<:Integer}(x::SparseVector{Tv,Ti}, v::Tv, i::Ti)
    nzind = nonzeroinds(x)
    nzval = nonzeros(x)

    m = length(nzind)
    k = searchsortedfirst(nzind, i)
    if 1 <= k <= m && nzind[k] == i  # i found
        if v == zero(v)
            deleteat!(nzind, k)
            deleteat!(nzval, k)
        else
            nzval[k] = v
        end
    else  # i not found
        if v != zero(v)
            insert!(nzind, k, i)
            insert!(nzval, k, v)
        end
    end
    x
end

setindex!{Tv, Ti<:Integer}(x::SparseVector{Tv,Ti}, v, i::Integer) =
    setindex!(x, convert(Tv, v), convert(Ti, i))

### Conversion

# convert SparseMatrixCSC to SparseVector
function convert{Tv,Ti<:Integer}(::Type{SparseVector{Tv,Ti}}, s::SparseMatrixCSC{Tv,Ti})
    size(s, 2) == 1 || throw(ArgumentError("The input argument must have a single-column."))
    SparseVector(s.m, s.rowval, s.nzval)
end

convert{Tv,Ti}(::Type{SparseVector{Tv}}, s::SparseMatrixCSC{Tv,Ti}) =
    convert(SparseVector{Tv,Ti}, s)

convert{Tv,Ti}(::Type{SparseVector}, s::SparseMatrixCSC{Tv,Ti}) =
    convert(SparseVector{Tv,Ti}, s)

# convert Vector to SparseVector

function _dense2sparsevec{Tv}(s::Vector{Tv}, initcap::Int)
    # pre-condition: initcap > 0
    n = length(s)
    cap = initcap
    nzind = Array(Int, cap)
    nzval = Array(Tv, cap)
    c = 0
    @inbounds for i = 1:n
        v = s[i]
        if v != zero(v)
            if c >= cap
                cap *= 2
                resize!(nzind, cap)
                resize!(nzval, cap)
            end
            c += 1
            nzind[c] = i
            nzval[c] = v
        end
    end
    if c < cap
        resize!(nzind, c)
        resize!(nzval, c)
    end
    SparseVector(n, nzind, nzval)
end

convert{Tv}(::Type{SparseVector{Tv,Int}}, s::Vector{Tv}) =
    _dense2sparsevec(s, max(8, div(length(s), 8)))

convert{Tv}(::Type{SparseVector{Tv}}, s::Vector{Tv}) =
    convert(SparseVector{Tv,Int}, s)

convert{Tv}(::Type{SparseVector}, s::Vector{Tv}) =
    convert(SparseVector{Tv,Int}, s)

# convert between different types of SparseVector
convert{Tv,Ti,TvS,TiS}(::Type{SparseVector{Tv,Ti}}, s::SparseVector{TvS,TiS}) =
    SparseVector{Tv,Ti}(s.n, convert(Vector{Ti}, s.nzind), convert(Vector{Tv}, s.nzval))

convert{Tv,TvS,TiS}(::Type{SparseVector{Tv}}, s::SparseVector{TvS,TiS}) =
    SparseVector{Tv,TiS}(s.n, s.nzind, convert(Vector{Tv}, s.nzval))

### Rand Construction

@compat function sprand{T}(n::Integer, p::AbstractFloat, rfn::Function, ::Type{T})
    I = randsubseq(1:convert(Int, n), p)
    V = rfn(T, length(I))
    SparseVector(n, I, V)
end

@compat function sprand(n::Integer, p::AbstractFloat, rfn::Function)
    I = randsubseq(1:convert(Int, n), p)
    V = rfn(length(I))
    SparseVector(n, I, V)
end

@compat sprand{T}(n::Integer, p::AbstractFloat, ::Type{T}) = sprand(n, p, rand, T)
@compat sprand(n::Integer, p::AbstractFloat) = sprand(n, p, rand)
@compat sprandn(n::Integer, p::AbstractFloat) = sprand(n, p, randn)


text:
# ## Getting started

using MLJ, RDatasets, ScientificTypes, PrettyPrinting

@load LinearRegressor pkg=MLJLinearModels
@load RidgeRegressor pkg=MLJLinearModels
@load LassoRegressor pkg=MLJLinearModels

hitters = dataset("ISLR", "Hitters")
@show size(hitters)
names(hitters) |> pprint

# The target is `Salary`

y, X = unpack(hitters, ==(:Salary), col->true);

# It has missing values which we will just ignore:

no_miss = .!ismissing.(y)
y = collect(skipmissing(y))
X = X[no_miss, :]
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=424);

# Most features are currently encoded as integers but we will consider them as continuous

Xc = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))
scitype(Xc)

# There're a few features that are categorical which we'll one-hot-encode.

# ## Ridge pipeline
# ### Baseline
#
# Let's first fit a simple pipeline with a standardizer, a one-hot-encoder and a basic linear regression:

@pipeline RegPipe(std = Standardizer(),
                  hot = OneHotEncoder(),
                  reg = LinearRegressor())

model = RegPipe()
pipe  = machine(model, Xc, y)
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# ### Basic ridge
#
# Let's now swap the linear regressor for a ridge one without specifying the penalty (`1` by default):

pipe.model.reg = RidgeRegressor()
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# Ok that's a bit better but surely we can do better with an appropriate selection of the hyperparameter.

# ### Cross validating

# What penalty should you use? Let's do a simple CV to try  to find out:

r  = range(model, :(reg.lambda), lower=1e-2, upper=100_000, scale=:log10)
tm = TunedModel(model=model, ranges=r, tuning=Grid(resolution=50),
                resampling=CV(nfolds=3, rng=4141), measure=rms)
mtm = machine(tm, Xc, y)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.reg.lambda, sigdigits=4)

# right, and  with that we get:

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# ## Lasso pipeline
#
# Let's do the same as above but using a Lasso model and adjusting the range a bit:

mtm.model.model.reg = LassoRegressor()
mtm.model.ranges = range(model, :(reg.lambda), lower=500, upper=100_000, scale=:log10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.reg.lambda, sigdigits=4)

# Ok and let's see how that does:

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# Pretty good! and the parameters are reasonably sparse as expected:

coefs = mtm.fitresult.fitresult.machine.fitresult
round.(coefs, sigdigits=2)

# with around 50% sparsity:

sum(coefs .≈ 0) / length(coefs)

# ## Elastic net pipeline

@load ElasticNetRegressor pkg=MLJLinearModels

mtm.model.model.reg = ElasticNetRegressor()
mtm.model.ranges = [range(model, :(reg.lambda), lower=0.1, upper=100, scale=:log10),
                    range(model, :(reg.gamma),  lower=500, upper=10_000, scale=:log10)]
mtm.model.tuning = Grid(resolution=10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
@show round(best_mdl.reg.lambda, sigdigits=4)
@show round(best_mdl.reg.gamma, sigdigits=4)

# And it's not too bad in terms of accuracy either

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# But the simple ridge regression seems to work best here.

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
function setindex_no_overwrite!(
        dict::AbstractDict{<:K, <:V},
        key::K,
        new_value,
    ) where {K, V}
    if !haskey(dict, key)
        dict[key] = new_value
        return nothing
    end

    old_value = dict[key]
    msg = "Duplicate key found"
    @error msg key old_value new_value
    throw(ErrorException(msg))
end

function setindex_same_value!(
        dict::AbstractDict{<:K, <:V},
        key::K,
        new_value,
    ) where {K, V}
    if !haskey(dict, key)
        dict[key] = new_value
        return nothing
    end

    old_value = dict[key]
    if old_value !== new_value # we require egality
        msg = "Duplicate key found"
        @error msg key old_value new_value
        throw(ErrorException(msg))
    end

    return nothing
end


text:
###
# Format and filters

# Format (i.e. how files are packed together)

# Reader

support_format_7zip(archive::Reader) =
    @_la_call(archive_read_support_format_7zip, (Ptr{Cvoid},), archive)
support_format_all(archive::Reader) =
    @_la_call(archive_read_support_format_all, (Ptr{Cvoid},), archive)
support_format_ar(archive::Reader) =
    @_la_call(archive_read_support_format_ar, (Ptr{Cvoid},), archive)
support_format_by_code(archive::Reader, code) =
    @_la_call(archive_read_support_format_by_code,
              (Ptr{Cvoid}, Cint), archive, code)
support_format_cab(archive::Reader) =
    @_la_call(archive_read_support_format_cab, (Ptr{Cvoid},), archive)
support_format_cpio(archive::Reader) =
    @_la_call(archive_read_support_format_cpio, (Ptr{Cvoid},), archive)
support_format_empty(archive::Reader) =
    @_la_call(archive_read_support_format_empty, (Ptr{Cvoid},), archive)
support_format_gnutar(archive::Reader) =
    @_la_call(archive_read_support_format_gnutar, (Ptr{Cvoid},), archive)
support_format_iso9660(archive::Reader) =
    @_la_call(archive_read_support_format_iso9660, (Ptr{Cvoid},), archive)
support_format_lha(archive::Reader) =
    @_la_call(archive_read_support_format_lha, (Ptr{Cvoid},), archive)
support_format_mtree(archive::Reader) =
    @_la_call(archive_read_support_format_mtree, (Ptr{Cvoid},), archive)
support_format_rar(archive::Reader) =
    @_la_call(archive_read_support_format_rar, (Ptr{Cvoid},), archive)
support_format_rar5(archive::Reader) =
    @_la_call(archive_read_support_format_rar5, (Ptr{Cvoid},), archive)
support_format_raw(archive::Reader) =
    @_la_call(archive_read_support_format_raw, (Ptr{Cvoid},), archive)
support_format_tar(archive::Reader) =
    @_la_call(archive_read_support_format_tar, (Ptr{Cvoid},), archive)
support_format_xar(archive::Reader) =
    @_la_call(archive_read_support_format_xar, (Ptr{Cvoid},), archive)
support_format_warc(archive::Reader) =
    @_la_call(archive_read_support_format_warc, (Ptr{Cvoid},), archive)
support_format_zip(archive::Reader) =
    @_la_call(archive_read_support_format_zip, (Ptr{Cvoid},), archive)
support_format_zip_streamable(archive::Reader) =
    @_la_call(archive_read_support_format_zip_streamable, (Ptr{Cvoid},), archive)
support_format_zip_seekable(archive::Reader) =
    @_la_call(archive_read_support_format_zip_seekable, (Ptr{Cvoid},), archive)

const _format_ids = Dict{String,Cint}(
    "7zip"=>Format._7ZIP,
    "ar"=>Format.AR_BSD,
    "arbsd"=>Format.AR_BSD,
    "argnu"=>Format.AR_GNU,
    "arsvr4"=>Format.AR_GNU,
    "bsdtar"=>Format.TAR_PAX_RESTRICTED,
    "cd9660"=>Format.ISO9660,
    "cpio"=>Format.CPIO_POSIX,
    "gnutar"=>Format.TAR_GNUTAR,
    "iso"=>Format.ISO9660,
    "iso9660"=>Format.ISO9660,
    "mtree"=>Format.MTREE,
    "mtree-classic"=>Format.MTREE,
    "newc"=>Format.CPIO_SVR4_NOCRC,
    "odc"=>Format.CPIO_POSIX,
    "oldtar"=>Format.TAR,
    "pax"=>Format.TAR_PAX_INTERCHANGE,
    "paxr"=>Format.TAR_PAX_RESTRICTED,
    "posix"=>Format.TAR_PAX_INTERCHANGE,
    "raw"=>Format.RAW,
    "rpax"=>Format.TAR_PAX_RESTRICTED,
    "shar"=>Format.SHAR_BASE,
    "shardump"=>Format.SHAR_DUMP,
    "ustar"=>Format.TAR_USTAR,
    "v7tar"=>Format.TAR,
    "v7"=>Format.TAR,
    "warc"=>Format.WARC,
    "xar"=>Format.XAR,
    "zip"=>Format.ZIP,
)

# Functions to manually set the format and filters to be used. This is
# useful to bypass the bidding process when the format and filters to use
# is known in advance.
set_format(archive::Reader, fmt::Integer) =
    @_la_call(archive_read_set_format, (Ptr{Cvoid}, Cint), archive, fmt)
function set_format(archive::Reader, fmt::AbstractString)
    id = get(_format_ids, fmt, Cint(-1))
    if id < 0
        throw(ArchiveFatal("No such format '$fmt'"))
    end
    @_la_call(archive_read_set_format, (Ptr{Cvoid}, Cint), archive, id)
end
set_format_7zip(archive::Reader) = set_format(archive, Format._7ZIP)
set_format_ar_bsd(archive::Reader) = set_format(archive, Format.AR_BSD)
set_format_ar_gnu(archive::Reader) = set_format(archive, Format.AR_GNU)
const set_format_ar_svr4 = set_format_ar_gnu
set_format_cpio(archive::Reader) = set_format(archive, Format.CPIO_POSIX)
set_format_cpio_newc(archive::Reader) = set_format(archive, Format.CPIO_SVR4_NOCRC)
set_format_gnutar(archive::Reader) = set_format(archive, Format.TAR_GNUTAR)
set_format_iso9660(archive::Reader) = set_format(archive, Format.ISO9660)
set_format_mtree(archive::Reader) = set_format(archive, Format.MTREE)
set_format_mtree_classic(archive::Reader) = set_format(archive, Format.MTREE)
set_format_pax(archive::Reader) = set_format(archive, Format.TAR_PAX_INTERCHANGE)
set_format_pax_restricted(archive::Reader) = set_format(archive, Format.TAR_PAX_RESTRICTED)
set_format_raw(archive::Reader) = set_format(archive, Format.RAW)
set_format_shar(archive::Reader) = set_format(archive, Format.SHAR_BASE)
set_format_shar_dump(archive::Reader) = set_format(archive, Format.SHAR_DUMP)
set_format_ustar(archive::Reader) = set_format(archive, Format.TAR_USTAR)
set_format_v7tar(archive::Reader) = set_format(archive, Format.TAR)
set_format_warc(archive::Reader) = set_format(archive, Format.WARC)
set_format_xar(archive::Reader) = set_format(archive, Format.XAR)
set_format_zip(archive::Reader) = set_format(archive, Format.ZIP)

# Writer

"A convenience function to set the format based on the code or name."
set_format(archive::Writer, format_code::Integer) =
    @_la_call(archive_write_set_format, (Ptr{Cvoid}, Cint),
              archive, format_code)
set_format(archive::Writer, name::AbstractString) =
    @_la_call(archive_write_set_format_by_name, (Ptr{Cvoid}, Cstring),
              archive, name)

set_format_7zip(archive::Writer) =
    @_la_call(archive_write_set_format_7zip, (Ptr{Cvoid},), archive)
set_format_ar_bsd(archive::Writer) =
    @_la_call(archive_write_set_format_ar_bsd, (Ptr{Cvoid},), archive)
set_format_ar_gnu(archive::Writer) =
    @_la_call(archive_write_set_format_ar_svr4, (Ptr{Cvoid},), archive)
# const set_format_ar_svr4 = set_format_ar_gnu
set_format_cpio(archive::Writer) =
    @_la_call(archive_write_set_format_cpio, (Ptr{Cvoid},), archive)
set_format_cpio_newc(archive::Writer) =
    @_la_call(archive_write_set_format_cpio_newc, (Ptr{Cvoid},), archive)
set_format_gnutar(archive::Writer) =
    @_la_call(archive_write_set_format_gnutar, (Ptr{Cvoid},), archive)
set_format_iso9660(archive::Writer) =
    @_la_call(archive_write_set_format_iso9660, (Ptr{Cvoid},), archive)
set_format_mtree(archive::Writer) =
    @_la_call(archive_write_set_format_mtree, (Ptr{Cvoid},), archive)
set_format_mtree_classic(archive::Writer) =
    @_la_call(archive_write_set_format_mtree_classic, (Ptr{Cvoid},), archive)
set_format_pax(archive::Writer) =
    @_la_call(archive_write_set_format_pax, (Ptr{Cvoid},), archive)
set_format_pax_restricted(archive::Writer) =
    @_la_call(archive_write_set_format_pax_restricted, (Ptr{Cvoid},), archive)
set_format_raw(archive::Writer) =
    @_la_call(archive_write_set_format_raw, (Ptr{Cvoid},), archive)
set_format_shar(archive::Writer) =
    @_la_call(archive_write_set_format_shar, (Ptr{Cvoid},), archive)
set_format_shar_dump(archive::Writer) =
    @_la_call(archive_write_set_format_shar_dump, (Ptr{Cvoid},), archive)
set_format_ustar(archive::Writer) =
    @_la_call(archive_write_set_format_ustar, (Ptr{Cvoid},), archive)
set_format_v7tar(archive::Writer) =
    @_la_call(archive_write_set_format_v7tar, (Ptr{Cvoid},), archive)
set_format_warc(archive::Writer) =
    @_la_call(archive_write_set_format_warc, (Ptr{Cvoid},), archive)
set_format_xar(archive::Writer) =
    @_la_call(archive_write_set_format_xar, (Ptr{Cvoid},), archive)
set_format_zip(archive::Writer) =
    @_la_call(archive_write_set_format_zip, (Ptr{Cvoid},), archive)
zip_set_compression_deflate(archive::Writer) =
    @_la_call(archive_write_zip_set_compression_deflate, (Ptr{Cvoid},), archive)
zip_set_compression_store(archive::Writer) =
    @_la_call(archive_write_zip_set_compression_store, (Ptr{Cvoid},), archive)

# filter (i.e. compression and other post-processing on the packed files)

# Reader

support_filter_all(archive::Reader) =
    @_la_call(archive_read_support_filter_all, (Ptr{Cvoid},), archive)
support_filter_bzip2(archive::Reader) =
    @_la_call(archive_read_support_filter_bzip2, (Ptr{Cvoid},), archive)
support_filter_compress(archive::Reader) =
    @_la_call(archive_read_support_filter_compress, (Ptr{Cvoid},), archive)
support_filter_gzip(archive::Reader) =
    @_la_call(archive_read_support_filter_gzip, (Ptr{Cvoid},), archive)
support_filter_grzip(archive::Reader) =
    @_la_call(archive_read_support_filter_grzip, (Ptr{Cvoid},), archive)
support_filter_lrzip(archive::Reader) =
    @_la_call(archive_read_support_filter_lrzip, (Ptr{Cvoid},), archive)
support_filter_lz4(archive::Reader) =
    @_la_call(archive_read_support_filter_lz4, (Ptr{Cvoid},), archive)
support_filter_lzip(archive::Reader) =
    @_la_call(archive_read_support_filter_lzip, (Ptr{Cvoid},), archive)
support_filter_lzma(archive::Reader) =
    @_la_call(archive_read_support_filter_lzma, (Ptr{Cvoid},), archive)
support_filter_lzop(archive::Reader) =
    @_la_call(archive_read_support_filter_lzop, (Ptr{Cvoid},), archive)
support_filter_none(archive::Reader) =
    @_la_call(archive_read_support_filter_none, (Ptr{Cvoid},), archive)
"""
Data is fed through the specified external program before being dearchived.
Note that this disables automatic detection of the compression format, so it
makes no sense to specify this in conjunction with any other decompression
option.
"""
support_filter_program(archive::Reader, cmd) =
    @_la_call(archive_read_support_filter_program,
              (Ptr{Cvoid}, Cstring), archive, cmd)
support_filter_program_signature(archive::Reader, cmd, sig) =
    @_la_call(archive_read_support_filter_program_signature,
              (Ptr{Cvoid}, Cstring, Ptr{Cvoid}, Csize_t),
              archive, cmd, sig, sizeof(sig))
support_filter_rpm(archive::Reader) =
    @_la_call(archive_read_support_filter_rpm, (Ptr{Cvoid},), archive)
support_filter_uu(archive::Reader) =
    @_la_call(archive_read_support_filter_uu, (Ptr{Cvoid},), archive)
support_filter_xz(archive::Reader) =
    @_la_call(archive_read_support_filter_xz, (Ptr{Cvoid},), archive)
support_filter_zstd(archive::Reader) =
    @_la_call(archive_read_support_filter_zstd, (Ptr{Cvoid},), archive)

append_filter(archive::Reader, filter) =
    @_la_call(archive_read_append_filter, (Ptr{Cvoid}, Cint), archive, filter)
append_filter_program(archive::Reader, cmd) =
    @_la_call(archive_read_append_filter_program, (Ptr{Cvoid}, Cstring),
              archive, cmd)
append_filter_program_signature(archive::Reader, cmd, sig) =
    @_la_call(archive_read_append_filter_program_signature,
              (Ptr{Cvoid}, Cstring, Ptr{Cvoid}, Csize_t),
              archive, cmd, sig, sizeof(sig))

# Writer

"A convenience function to set the filter based on the code."
add_filter(archive::Writer, filter_code::Integer) =
    @_la_call(archive_write_add_filter, (Ptr{Cvoid}, Cint),
              archive, filter_code)
add_filter(archive::Writer, name::AbstractString) =
    @_la_call(archive_write_add_filter_by_name, (Ptr{Cvoid}, Cstring),
              archive, name)
add_filter_b64encode(archive::Writer) =
    @_la_call(archive_write_add_filter_b64encode, (Ptr{Cvoid},), archive)
add_filter_bzip2(archive::Writer) =
    @_la_call(archive_write_add_filter_bzip2, (Ptr{Cvoid},), archive)
add_filter_compress(archive::Writer) =
    @_la_call(archive_write_add_filter_compress, (Ptr{Cvoid},), archive)
add_filter_grzip(archive::Writer) =
    @_la_call(archive_write_add_filter_grzip, (Ptr{Cvoid},), archive)
add_filter_gzip(archive::Writer) =
    @_la_call(archive_write_add_filter_gzip, (Ptr{Cvoid},), archive)
add_filter_lrzip(archive::Writer) =
    @_la_call(archive_write_add_filter_lrzip, (Ptr{Cvoid},), archive)
add_filter_lz4(archive::Writer) =
    @_la_call(archive_write_add_filter_lz4, (Ptr{Cvoid},), archive)
add_filter_lzip(archive::Writer) =
    @_la_call(archive_write_add_filter_lzip, (Ptr{Cvoid},), archive)
add_filter_lzma(archive::Writer) =
    @_la_call(archive_write_add_filter_lzma, (Ptr{Cvoid},), archive)
add_filter_lzop(archive::Writer) =
    @_la_call(archive_write_add_filter_lzop, (Ptr{Cvoid},), archive)
add_filter_none(archive::Writer) =
    @_la_call(archive_write_add_filter_none, (Ptr{Cvoid},), archive)
add_filter_program(archive::Writer, cmd::AbstractString) =
    @_la_call(archive_write_add_filter_program,
              (Ptr{Cvoid}, Cstring), archive, cmd)
add_filter_uuencode(archive::Writer) =
    @_la_call(archive_write_add_filter_uuencode, (Ptr{Cvoid},), archive)
add_filter_xz(archive::Writer) =
    @_la_call(archive_write_add_filter_xz, (Ptr{Cvoid},), archive)
add_filter_zstd(archive::Writer) =
    @_la_call(archive_write_add_filter_zstd, (Ptr{Cvoid},), archive)

function set_format_filter_by_ext(archive::Writer, filename, def_ext=nothing)
    if def_ext === nothing
        @_la_call(archive_write_set_format_filter_by_ext,
                  (Ptr{Cvoid}, Cstring), archive, filename)
    else
        @_la_call(archive_write_set_format_filter_by_ext_def,
                  (Ptr{Cvoid}, Cstring, Cstring), archive, filename, def_ext)
    end
end


text:
module Plot

using Distributions
using KernelDensity

import PyPlot
const plt = PyPlot
export plt

# For Debugging
printKeys(x) = [print("$k  ") for k in keys(x)]

#=
plt.ioff() # turn off interactive graphics. turn on with plt.ion()
plt.ion() # turn on interactive graphics. turn on with plt.ion()
=#

function pyRange(n::Int)
  return collect(0:(n-1))
end

function plotZ(Z::Matrix; xticks=true, kw...)
  J, K = size(Z)
  cm = plt.cm_get_cmap(:Greys, 5)
  img = plt.imshow(Z, aspect="auto", vmin=0, vmax=1, cmap=cm; kw...);
  [ plt.axhline(y=i+.5, color="grey", linewidth=.5) for i in pyRange(J)];
  [ plt.axvline(x=i+.5, color="grey", linewidth=.5) for i in pyRange(K)];
  plt.yticks(pyRange(J), 1:J);
  if xticks
    plt.xticks(pyRange(K), 1:K, rotation=:vertical);
  end
  ax = plt.gca()
  ax[:tick_params](length=0)
  return ax
end

# Define BlueGreyRed ColorMap
let
  cm_name = :BlueGreyRed
  c = Matrix([[0, 0, 1] [.9, .9, .9] [1, 0, 0]]')
  cm = Plot.plt.ColorMap(c)
  Plot.plt.register_cmap(cm_name, cm)
end

function plotY(Y::Matrix; colorbar=true,
               cm=PyPlot.get_cmap(:BlueGreyRed, 9), vmin=-4, vmax=4, kw...)
  cm[:set_under](color=:blue)
  cm[:set_over](color=:red)
  cm[:set_bad](color=:black)
  img = plt.imshow(Y, aspect="auto", cmap=cm, vmin=vmin, vmax=vmax; kw...) # vmin, vmax
  if colorbar
    plt.plt[:colorbar]();
  end
  ax = plt.gca()
  ax[:tick_params](length=0)
  return img
end

"""
# Example

```julia
J, K = 20, 4
Y = randn(1000, J)
Y[:, 2] .= NaN
Z = Int.(randn(J, K) .> 0)
# Plot.plotY(Y)

Plot.plotYZ(Y, Z)
Plot.plt.savefig("yZ.pdf", bbox_inches="tight")
Plot.plt.close()
```
"""
function plotYZ(Y::Matrix, Z::Matrix, vmin=-4, vmax=4, y_ytick_fs=5, z_ytick_fs=5; kw...)
  J, K = size(Z)

  # Y
  plt.subplot2grid((10, 1), (0, 0), rowspan=7)
  plotY(Matrix(Y), colorbar=false, vmin=vmin, vmax=vmax);
  plt.yticks(fontsize=y_ytick_fs)
  plt.plt[:colorbar]()
  plt.xticks(pyRange(J), 1:J, rotation=:vertical);

  # Z
  plt.subplot2grid((10, 1), (8, 0), rowspan=3)
  ax = plotZ(Matrix(Z'))
  plt.xticks([])
  plt.yticks(fontsize=z_ytick_fs)
  # Add ticks to other axis
  ax = plt.gca()
  ax2 = ax[:twinx]()
  plt.plt[:colorbar](aspect=5)
  # Change default ticks! Comment this out to see the difference.
  ax2[:set_yticks](pyRange(K))
  plt.yticks((K-1) / K * pyRange(K) .+ .5, (1:K) / 10, fontsize=5)
  ax2[:tick_params](length=0)
end

# TODO: Plot in plot
# https://stackoverflow.com/questions/17458580/embedding-small-plots-inside-subplots-in-matplotlib
# https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html
function infig_position(ax, rect)
  # Get current figure
  fig = plt.gcf()

  box = ax[:get_position]()
  inax_position = ax[:transAxes][:transform](rect[1:2])
  transFigure = fig[:transFigure][:inverted]()
  infig_pos = transFigure[:transform](inax_position)
  width = box[:width] * rect[3]
  height = box[:height] * rect[4]
  xpos = infig_pos[1]
  ypos = infig_pos[2]

  return xpos, ypos, width, height
end

function plotPost(y, ax::PyPlot.PyObject; rect=[.7, .7, .3, .2], showTrace::Bool=true,
                  showAcc::Bool=true, digits_ci=2, c=:royalblue,
                  alpha_level=.05, useHistogram=false, normalizeHist=true, histBins=0,
                  accFontsize=6, annfs=6, hide_yaxis=true)
  # Get current figure
  fig = plt.gcf()

  if showTrace
    # Get box position / dimensions
    xpos, ypos, width, height = infig_position(ax, rect)

    subax = fig[:add_axes]([xpos, ypos, width, height])
    subax[:axis](:off)
    subax[:plot](1:length(y), y, c=:grey, linewidth=.5)
    accRate = length(unique(y)) / length(y)
    if showAcc
      subax[:set_title]("acc: $(accRate*100)%", fontsize=accFontsize)
    end
  end

  # plot histogram / density
  if useHistogram
    # Not really supported
    if histBins == 0
      ax[:hist](y, density=normalizeHist)
    else
      ax[:hist](y, density=normalizeHist, bins=histBins)
    end
  else
    if length(unique(y)) == 1
      ax[:annotate]("mean: $(round(mean(y), digits=digits_ci))", color=:red,
                    xy=(.8, .6), fontsize=annfs, xycoords="axes fraction")
      ax[:annotate]("SD: $(round(std(y), digits=digits_ci))",
                    xy=(.8, .4), fontsize=annfs, xycoords="axes fraction")
    else
      # Confidence Interval
      a_lower = alpha_level / 2
      a_upper = 1 - a_lower
      y_ci = quantile(y, [a_lower, a_upper])
      d = kde(y)
      # ax[:plot](d_ci[:x], d_ci[:dx], c=c)
      ax[:plot](d.x, d.density, c=c, alpha=.5)
      ax[:fill_between](d.x[y_ci[1] .< d.x .< y_ci[2]],
                        d.density[y_ci[1] .< d.x .< y_ci[2]], color=c)
      ax[:fill_between](d.x, d.density, color=c, alpha=.3)
      ax[:set_xticks](round.(y_ci, digits=digits_ci))
      y_mean = mean(y)
      ax[:vlines](y_mean, ymin=0, ymax=d.density[argmin(abs.(d.x .- y_mean))],
                  color=:red, linewidth=1)

      # Plot Annotations
      ax[:annotate]("mean: $(round(y_mean, digits=digits_ci))", color=:red,
                    xy=(.8, .6), fontsize=annfs, xycoords="axes fraction")
      ax[:annotate]("$(Int(round((1-alpha_level) * 100)))% CI", color=c,
                    xy=(.8, .5), fontsize=annfs, xycoords="axes fraction")
      ax[:annotate]("SD: $(round(std(y), digits=digits_ci))",
                    xy=(.8, .4), fontsize=annfs, xycoords="axes fraction")
    end
  end

  ax[:spines]["right"][:set_visible](false)
  ax[:spines]["top"][:set_visible](false)

  if hide_yaxis
    ax[:axes][:yaxis][:set_ticklabels]([])
    ax[:spines]["left"][:set_visible](false)
    ax[:tick_params](length=0, axis="y")
  end

  ax[:spines]["bottom"][:set_color]("grey")
end

function plotPost(y; kw...)
  fig, ax = plt.subplots()
  plotPost(y, ax; kw...)
  return fig, ax
end

"""
# Example

```julia
Y = cumsum(randn(1000, 4), dims=1)
Plot.plotPosts(Y)
```
"""
function plotPosts(Y; use_tight_layout=true,
                   fig_size_inches=missing,
                   fig_subplots_adjust=missing,
                   kw...)
  J = size(Y, 2)
  fig, ax = plt.subplots(J,J)
  if use_tight_layout
    fig[:tight_layout]()
  end

  if !ismissing(fig_size_inches)
    fig[:set_size_inches](fig_size_inches[1], fig_size_inches[2])
  end

  if !ismissing(fig_subplots_adjust)
    fig[:subplots_adjust](wspace=fig_subplots_adjust[1],
                          hspace=fig_subplots_adjust[2])
  end

  plotPosts(Y, fig, ax; kw...)
end

function plotPosts(Y, fig, ax; details=false, digits_cor=2, scatter_alpha=.7, scatter_s=.5,
                   warn=true, kw...)
  J = size(Y, 2)
  counter = 0
  for c in 1:J
    for r in 1:J
      counter += 1
      if r == c
        plotPost(Y[:, r], ax[counter]; kw...)
      elseif r > c
        ax[counter][:axis](:off)
        #xpos, ypos, w, h = infig_position(ax[counter], [.5, .5, .5, .5])
        cor_ycr = round(cor(Y[:, c], Y[:, r]), digits=digits_cor)
        if !isnan(cor_ycr)
          ax[counter][:text](0, .4, "r = $cor_ycr", fontsize=abs(cor_ycr)*10+10)
        else
          if warn
            @warn "Correlation was NaN."
          end
          ax[counter][:text](0, .4, "r = ???", fontsize=10)
        end
      else
        #ax[counter][:plot](Y[:, c], Y[:, r], c=:grey, linewidth=.5)
        ax[counter][:scatter](Y[:, c], Y[:, r], c=:grey, alpha=scatter_alpha, s=scatter_s)
        if !details
          ax[counter][:axis](:off)
        end
      end
    end
  end
  return fig, ax
end

end # Plot


text:
function Transformation()
    flip = node(:flip, (false, false, false))
    scale = node(:scale, Vec3f0(1))
    scale = lift(flip, scale) do f, s
        map((f, s)-> f ? -s : s, Vec(f), s)
    end
    translation, rotation, align = (
        node(:translation, Vec3f0(0)),
        node(:rotation, Quaternionf0(0, 0, 0, 1)),
        node(:align, Vec2f0(0))
    )
    trans = nothing
    model = map_once(scale, translation, rotation, align, flip) do s, o, q, a, flip
        parent = if trans !== nothing && isassigned(trans.parent)
            boundingbox(trans.parent[])
        else
            nothing
        end
        transformationmatrix(o, s, q, a, flip, parent)
    end
    trans = Transformation(
        translation,
        scale,
        rotation,
        model,
        flip,
        align,
        signal_convert(Node{Any}, identity)
    )
end

function Transformation(scene::SceneLike)
    flip = node(:flip, (false, false, false))
    scale = node(:scale, Vec3f0(1))
    translation, rotation, align = (
        node(:translation, Vec3f0(0)),
        node(:rotation, Quaternionf0(0, 0, 0, 1)),
        node(:align, Vec2f0(0))
    )
    pmodel = transformationmatrix(scene)
    trans = nothing
    model = map_once(scale, translation, rotation, align, pmodel, flip) do s, o, q, a, p, f
        bb = if trans !== nothing && isassigned(trans.parent)
            boundingbox(trans.parent[])
        else
            nothing
        end
        p * transformationmatrix(o, s, q, align, f, bb)
    end
    trans = Transformation(
        translation,
        scale,
        rotation,
        model,
        flip,
        align,
        signal_convert(Node{Any}, identity)
    )
    return trans
end

function translated(scene::Scene, translation...)
    tscene = Scene(scene, transformation = Transformation())
    transform!(tscene, translation...)
    tscene
end

function translated(scene::Scene; kw_args...)
    tscene = Scene(scene, transformation = Transformation())
    transform!(tscene; kw_args...)
    tscene
end

function transform!(
        scene::SceneLike;
        translation = Vec3f0(0),
        scale = Vec3f0(1),
        rotation = 0.0,
    )
    translate!(scene, translation)
    scale!(scene, scale)
    rotate!(scene, rotation)
end

transformation(t::Scene) = t.transformation
transformation(t::AbstractPlot) = t.transformation
transformation(t::Transformation) = t

scale(t::Transformable) = transformation(t).scale
scale!(t::Transformable, s) = (scale(t)[] = to_ndim(Vec3f0, Float32.(s), 1))
scale!(t::Transformable, xyz...) = scale!(t, xyz)

rotation(scene::Transformable) = transformation(scene).rotation
function rotate!(::Type{T}, scene::Transformable, q) where T
    rot = convert_attribute(q, key"rotation"())
    if T === Accum
        rot1 = rotation(scene)[]
        rotation(scene)[] = qmul(rot1, rot)
    elseif T == Absolute
        rotation(scene)[] = rot
    else
        error("Unknown transformation: $T")
    end
end

rotate!(::Type{T}, scene::Transformable, axis_rot...) where T = rotate!(T, scene, axis_rot)
rotate!(scene::Transformable, axis_rot...) = rotate!(Absolute, scene, axis_rot)
rotate!(scene::Transformable, axis_rot::Quaternion) = rotate!(Absolute, scene, axis_rot)
rotate!(scene::Transformable, axis_rot::AbstractFloat) = rotate!(Absolute, scene, axis_rot)

translation(scene::Transformable) = transformation(scene).translation

struct Accum end; struct Absolute end

function translate!(::Type{T}, scene::Transformable, t) where T
    offset = to_ndim(Vec3f0, Float32.(t), 0)
    if T === Accum
        translation(scene)[] = translation(scene)[] .+ offset
    elseif T === Absolute
        translation(scene)[] = offset
    else
        error("Unknown translation type: $T")
    end
end
translate!(scene::Transformable, xyz::VecTypes) = translate!(Absolute, scene, xyz)
translate!(scene::Transformable, xyz...) = translate!(Absolute, scene, xyz)
translate!(::Type{T}, scene::Transformable, xyz...) where T = translate!(T, scene, xyz)

function transform!(scene::Transformable, x::Tuple{Symbol, <: Number})
    plane, dimval = string(x[1]), Float32(x[2])
    if length(plane) != 2 || (!all(x-> x in ('x', 'y', 'z'), plane))
        error("plane needs to define a 2D plane in xyz. It should only contain 2 symbols out of (:x, :y, :z). Found: $plane")
    end
    if all(x-> x in ('x', 'y'), plane) # xy plane
        translate!(scene, 0, 0, dimval)
    elseif all(x-> x in ('x', 'z'), plane) # xz plane
        rotate!(scene, Vec3f0(1, 0, 0), 0.5pi)
        translate!(scene, 0, dimval, 0)
    else #yz plane
        q1 = qrotation(Vec3f0(1, 0, 0), -0.5pi)
        q2 = qrotation(Vec3f0(0, 0, 1), 0.5pi)
        rotate!(scene, q2 * q1)
        translate!(scene, dimval, 0, 0)
    end
    scene
end

transformationmatrix(x) = transformation(x).model


text:
struct LinearAnalysis
    options::Dict{Symbol, Any}
end

LinearAnalysis(; kwargs...) = LinearAnalysis(Dict{Symbol, Any}(kwargs))

function (l::LinearAnalysis)(le::Entry)
    return splitapply(le) do entry
        labels, mappings = map(getlabel, entry.mappings), map(getvalue, entry.mappings)
        x, y = mappings.positional
        x̂ = [x x]
        x̂[:, 2] .= 1
        a, b = x̂ \ y
        length = 100
        rg = range(extrema(x)...; length)
        m = arguments(rg, a .* rg .+ b)
        return Entry(
            AbstractPlotting.plottype(entry.plottype, Lines),
            map(Labeled, labels, m),
            entry.attributes
        )
    end
end

linear(; kwargs...) = Layer((LinearAnalysis(; kwargs...),))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

#=
test subjects:
- modal analysis, with mesh tie contact

Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ)cos(λᵢℓ)

1:  4.730040744862704
2:  7.853204624095838
3: 10.995607838001671

[1] De Silva, Clarence W. Vibration: fundamentals and practice. CRC press, 2006, p.355

Code Aster solution:
--------------------
numéro    fréquence (HZ)     norme d'erreur
    1       1.12946E+00        5.81018E-12
    2       1.13141E+00        6.33463E-12
    3       2.93779E+00        6.53408E-13
    4       2.94143E+00        5.43970E-13
    5       4.51684E+00        5.43252E-13
=#

comm_CA = """
DEBUT(PAR_LOT="NON")

MAIL = LIRE_MAILLAGE(FORMAT="MED", NOM_MED="CYLINDER_20_SPLITTED")

MO = AFFE_MODELE(
    MAILLAGE=MAIL,
    AFFE=_F(TOUT="OUI",
            PHENOMENE="MECANIQUE", MODELISATION="3D"))

MAT = DEFI_MATERIAU(
    ELAS=_F(E=50475.45, NU=0.3, RHO=1.0))

CHMAT = AFFE_MATERIAU(
    MAILLAGE=MAIL,
    AFFE=_F(TOUT="OUI", MATER=MAT))

BC1 = AFFE_CHAR_MECA(
    MODELE=MO,
    DDL_IMPO=(
        _F(GROUP_MA=("CYLINDER_20_1_FACE1"), DX=0, DY=0, DZ=0)))

BC2 = AFFE_CHAR_MECA(
    MODELE=MO,
    DDL_IMPO=(
        _F(GROUP_MA=("CYLINDER_20_2_FACE2"), DX=0, DY=0, DZ=0)))

# ESCL = SLAVE
# MAIT = MASTER
BC3 = AFFE_CHAR_MECA(
    MODELE=MO,
    LIAISON_MAIL=_F(
        GROUP_MA_ESCL="CYLINDER_20_1_FACE2",
        GROUP_MA_MAIT="CYLINDER_20_2"))

# assemble material stiffness matrix

RIGEL = CALC_MATR_ELEM(
    MODELE=MO,
    OPTION="RIGI_MECA",
    CHAM_MATER=CHMAT,
    CHARGE=(BC1, BC2, BC3))

NUMEDDL = NUME_DDL(
    MATR_RIGI=RIGEL)

RIGAS = ASSE_MATRICE(
    MATR_ELEM=RIGEL,
    NUME_DDL=NUMEDDL)

# assemble mass matrix

MASSEL = CALC_MATR_ELEM(
    MODELE=MO,
    OPTION="MASS_MECA",
    CHAM_MATER=CHMAT,
    CHARGE=(BC1, BC2, BC3))

MASSAS = ASSE_MATRICE(
    MATR_ELEM=MASSEL,
    NUME_DDL=NUMEDDL)

# modal analysis, without geometric stiffness

BRESU = CALC_MODES(
    MATR_RIGI=RIGAS,
    MATR_MASS=MASSAS,
    OPTION="BANDE",
    CALC_FREQ=_F(
        FREQ=(0.0, 5.0)))

# modal analysis, with geometric stiffness

BRESU = NORM_MODE(
    reuse=BRESU,
    MODE=BRESU,
    NORME="TRAN")

IMPR_RESU(
    MODELE=MO,
    FORMAT="RESULTAT",
    RESU=_F(RESULTAT=BRESU))

IMPR_RESU(
    FORMAT="MED",
    UNITE=80,
    RESU=_F(RESULTAT=BRESU))

FIN()
"""

@testset "splitted rod with tie contact" begin
    # CYLINDER_20_1_FACE1 -- CYLINDER_20_1_FACE2 -- CYLINDER_20_2_FACE_1 -- CYLINDER_20_2_FACE_2
    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, "CYLINDER_20_SPLITTED")
    body1 = Problem(mesh, Elasticity, "CYLINDER_20_1", 3)
    body2 = Problem(mesh, Elasticity, "CYLINDER_20_2", 3)
    for body in [body1, body2]
        update!(body.elements, "youngs modulus", 54475.45)
        update!(body.elements, "poissons ratio", 0.3)
        update!(body.elements, "density", 1.0)
    end
    bc1 = Problem(mesh, Dirichlet, "CYLINDER_20_1_FACE1", 3, "displacement")
    bc2 = Problem(mesh, Dirichlet, "CYLINDER_20_2_FACE2", 3, "displacement")
    for bc in [bc1, bc2]
        update!(bc.elements, "displacement 1", 0.0)
        update!(bc.elements, "displacement 2", 0.0)
        update!(bc.elements, "displacement 3", 0.0)
    end
    interface = Problem(Mortar, "interface between bodies", 3, "displacement")
    slave = create_elements(mesh, "CYLINDER_20_1_FACE2")
    master = create_elements(mesh, "CYLINDER_20_2_FACE1")
    update!(slave, "master elements", master)
    interface.elements = [slave; master]

    solver = Solver(Modal, body1, body2, bc1, bc2, interface)
    solver.properties.nev = 5
    solver.properties.which = :SM
    solver()
    freqs_jf = sqrt(solver.properties.eigvals)/(2*pi)

    freqs_ca = [1.12946E+00, 1.13141E+00, 2.93779E+00, 2.94143E+00, 4.51684E+00]
    freq_jf = freqs_jf[1]
    freq_ca = freqs_ca[1]
    rtol = norm(freq_jf - freq_ca)/max(freq_jf, freq_ca)
    info("rtol = $rtol")
    for (i, freq) in enumerate(freqs_jf)
        @printf "mode %i | freq JuliaFEM %8.3f | freq Code Aster %8.3f\n" i freqs_jf[i] freqs_ca[i]
    end
    if rtol > 1.0e-3
        outfile = tempname() * ".xmf"
        info("Something went wrong, results are saved to $outfile")
        result = XDMF()
        elems = [body1.elements; body2.elements]
        for (i, freq) in enumerate(freqs_jf)
            xdmf_new_result!(result, elems, freq)
            xdmf_save_field!(result, elems, freq, "displacement"; field_type="Vector")
        end
        xdmf_save!(result, outfile)
    end
    @test rtol < 0.05

end


text:
using LinearAlgebra
using Random
using SparseArrays
using Test

using Jutils.Elements
using Jutils.Functions
using Jutils.Integration
using Jutils.Mesh
using Jutils.Transforms
using Jutils.Topologies

const lineelt = Element(Simplex{1}(), 1)
const squareelt = Element(Tensor([Simplex{1}(), Simplex{1}()]), 1)

ev(func, pt, elt) = callable(func)(pt, elt)

@testset "Transforms" begin include("Transforms.jl") end

@testset "Functions" begin include("Functions.jl") end
@testset "Gradients" begin include("Gradients.jl") end

@testset "Optimization" begin include("Optimization.jl") end

@testset "Integration" begin include("Integration.jl") end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
###################################################
#
# Copyright 2016 SOASTA, Inc.
# Distributed under the terms of the MIT license
#
# File: mPulseAPI.jl
#
# Functions to communicate with the mPulse Query and Repository REST APIs
#
###################################################

__precompile__(true)

module mPulseAPI
# This should bind at compile time, so @__FILE__ is set to mPulseAPI.jl
const __module_dir = dirname(dirname(@__FILE__))

function readdoc(path::AbstractString...)
    docpath = joinpath(mPulseAPI.__module_dir, path...)
    if VERSION < v"0.5"
        return readall(docpath)
    elseif VERSION < v"0.7"
        return readstring(docpath)
    else
        return read(docpath, AbstractString)
    end
end

"""
$(replace(mPulseAPI.readdoc("README.md"), r"\n.*travis-ci\.org.*\n", ""))
"""
mPulseAPI

using Base.Dates
import Base: @__doc__

using Requests, LightXML, HttpCommon, Formatting, TimeZones

# Tells docgen.jl to document internal methods as well
const __document_internal = true

const default_SOASTAEndpoint = "https://mpulse.soasta.com/concerto"

function __init__()
    setEndpoints()

    global verbose = false
end

"""
Change the SOASTA API endpoint that we connect to.  The default is `$(mPulseAPI.default_SOASTAEndpoint)`

#### Example

```julia
mPulseAPI.setEndpoints("https://mpulse-alt.soasta.com/concerto")
```
"""
function setEndpoints(SOASTAEndpoint::AbstractString = default_SOASTAEndpoint)
    global mPulseEndpoint     = "$SOASTAEndpoint/mpulse/api/v2/"
    global RepositoryEndpoint = "$SOASTAEndpoint/services/rest"
    global RepositoryService  = "$RepositoryEndpoint/RepositoryService/v1"
    global TokenEndpoint      = "$RepositoryService/Tokens"
    global ObjectEndpoint     = "$RepositoryService/Objects"

    return (mPulseEndpoint, RepositoryService)
end

"""
Set verbosity of API calls.

If set to true, all URLs, headers and POST data will be printed to the console before making an API call.
"""
function setVerbose(vbs::Bool)
    global verbose = vbs
end

# Convenience method because the mPulse API is bad with dates
function iso8601ToDateTime(date::AbstractString)
    date = replace(date, r"(?:Z|\+[\d:]+)(\s|$)", "")
    return DateTime(date)
end

# Internal convenience function for retrieving object info
function getObjectInfo(token::AbstractString, objectType::AbstractString, objectID::Int64, name::AbstractString)

    if objectType == "alert"
        object = getRepositoryAlert(token, alertID=objectID, alertName = name)
    elseif objectType == "domain"
        object = getRepositoryDomain(token, domainID=objectID, appName = name)
    elseif objectType == "tenant"
        object = getRepositoryTenant(token, tenantID=objectID, name = name)
    elseif objectType == "statisticalmodel"
        object = getRepositoryStatModel(token, statModelID=objectID, statModelName = name)
    end

    return object
end

# Fix datatype of elements got from JSON because the mPulse API is bad in what it sends
function fixJSONDataType(value::Union{AbstractString, Void})
    if value == nothing || isa(value, Void)
        # nothing
    elseif value == "null"
        value = nothing
    elseif value == "true"
        value = true
    elseif value == "false"
        value = false
    elseif ismatch(r"^-?\d+$", value)
        value = parse(Int, value, 10)
    elseif ismatch(r"^-?\d+\.\d+$", value)
        value = parse(Float64, value)
    end

    return value
end

function readdocs(name::AbstractString, replacers=[]; indent=0)
    # read the file and strip out the newline at the end
    data = chomp( readdoc( "doc-snippets", name * ".md" ) )

    # If any placeholders have default values, pull them out and insert them into the replace
    # array if a value is not already set
    while( (defaults = match(r"\{(\d+)=(.*?)\}", data)) != nothing)
        id      = parse(Int, defaults.captures[1])
        default = defaults.captures[2]

        data = replace(data, Regex("\\{$id=(.*?)\\}"), "{$id}")

        if length(replacers) < id
            resize!(replacers, id)
            replacers[id] = default
        end
    end

    if length(replacers) > 0
        # If we have replacers, then replace all non-placeholders braces with {{}}
        data = replace(data, r"(\{[A-Za-z][\w,|]+\})", s"{\1}")

        # And run the whole thing through format
        try
            data = Formatting.format(data, replacers...)
        catch
            println(replacers)
            println(data)
            rethrow()
        end
    end

    if indent > 0
        data = replace(data, r"^"m, repeat(" ", indent))
    end

    return data
end

include(joinpath(dirname(@__FILE__), "exceptions.jl"))
include(joinpath(dirname(@__FILE__), "cache_utilities.jl"))

include(joinpath(dirname(@__FILE__), "RepositoryAPI.jl"))
include(joinpath(dirname(@__FILE__), "StatisticalModel.jl"))
include(joinpath(dirname(@__FILE__), "Alert.jl"))
include(joinpath(dirname(@__FILE__), "Domain.jl"))
include(joinpath(dirname(@__FILE__), "Tenant.jl"))
include(joinpath(dirname(@__FILE__), "Token.jl"))
include(joinpath(dirname(@__FILE__), "QueryAPI.jl"))

end


text:
function rand_matrix(n, symm=true, diag=4.0)
    A = randn(n,n)

    if symm
        A = 0.5*(A+A')
    end
    for i = 1:n
        A[i,i] += diag
    end

    A
end

simple_matrix(n) = [2.0 1.0; 1.0 2.0]

using Base.LinAlg.LAPACK.pttrf!
using Base.LinAlg.LAPACK.pttrs!

function exemplo_tri1(n)
    A = 2.0 * eye(n)
    for i = 1:(n-1)
        A[i,i+1] = 1.0
        A[i+1,i] = 1.0
    end

    D = 2.0 * ones(n)
    E = ones(n-1)

    x1 = randn(n)
    x2 = copy(x1)
    y1 = A\x1
    pttrf!(D, E)
    y2 = pttrs!(D, E, x2)

    maxabs(y1-y2)
end

function assemble!{T<:Number}(Ag::Array{T,2}, Ae, m)
    nloc = size(Ae,1)
    for i = 1:nloc
        ig = m[i]
        for k = 1:nloc
            kg = m[k]
            Ag[kg,ig] += Ae[k,i]
        end
    end
end

function exemplo_tri2(n=5)
    ne = n-1

    Ag = zeros(n,n)
    Atri = HPFEM.BBTriSym{Float64}(n, n)
    for e = 1:ne
        Ae = rand_matrix(2)
        m = [e,e+1]
        HPFEM.assemble!(Atri, Ae, m)
        assemble!(Ag, Ae, m)
    end

    x1 = randn(n)
    x2 = copy(x1)

    y1 = Ag\x1
    HPFEM.trf!(Atri)
    y2 = HPFEM.trs!(Atri, x2)

    maxabs(y2-y1)
end

function exemplo_tri3(n=5)

    ne = n-1

    Ag = zeros(n,n)
    Atri = HPFEM.BBTri{Float64}(n, n)
    for e = 1:ne
        Ae = rand_matrix(2, false)
        m = [e,e+1]
        HPFEM.assemble!(Atri, Ae, m)
        assemble!(Ag, Ae, m)
    end

    x1 = randn(n)
    x2 = copy(x1)

    y1 = Ag\x1
    HPFEM.trf!(Atri)
    y2 = HPFEM.trs!(Atri, x2)

    maxabs(x2 - y1)
    hcat(y1, y2, y2-y1)

end


text:
# Linear & Affine Transforms

### Transpose

immutable Transpose{T,A<:StridedMatrix{T}}
    mat::A
end

Transpose{T}(a::StridedMatrix{T}) = Transpose{T,typeof(a)}(a)

Base.size(a::Transpose) = ((m,n)=size(a.mat); (n,m))::(Int,Int)
Base.size(a::Transpose, i::Integer) = (i == 1 ? size(a.mat,2) :
                                       i == 2 ? size(a.mat,1) : 1)

Base.length(a::Transpose) = length(a.mat)
Base.ndims(a::Transpose) = 2

### Linear transformation

typealias LinearTransform{T<:Real} Union(T,
                                         StridedVector{T},
                                         StridedMatrix{T},
                                         Transpose{T})

typealias SimpleLinearTransform{T<:Real} Union(T, StridedVector{T})

_outdim(a::StridedVector) = length(a)
_outdim(a::StridedMatrix) = size(a,1)
_outdim(a::Transpose) = size(a.mat,2)

function transform!(y::StridedVector, a::Real, x::StridedVector)
    (length(y) == length(x)) || throw(DimensionMismatch("Mismatched dimensions."))
    for i = 1:length(x)
        @inbounds y[i] = a * x[i]
    end
    return y
end

transform!(a::Real, x::StridedVector) = scale!(x, a)

function transform!(y::StridedMatrix, a::Real, x::StridedMatrix)
    m, n = size(x)
    (size(y) == (m, n)) || throw(DimensionMismatch("Mismatched dimensions."))
    for j = 1:n
        xj = view(x,:,j)
        yj = view(y,:,j)
        for i = 1:m
            @inbounds yj[i] = a * xj[i]
        end
    end
    return y
end

transform!(a::Real, x::StridedMatrix) = scale!(x, a)

function transform!{T<:Real}(y::StridedVector, a::StridedVector{T}, x::StridedVector)
    (length(y) == length(a) == length(x)) || throw(DimensionMismatch("Mismatched dimensions."))
    for i = 1:length(x)
        @inbounds y[i] = a[i] * x[i]
    end
    return y
end

transform!{T<:Real}(a::StridedVector{T}, x::StridedVector) = multiply!(x, a)

function transform!{T<:Real}(y::StridedMatrix, a::StridedVector{T}, x::StridedMatrix)
    m, n = size(x)
    (size(y) == (m, n) && length(a) == m) || throw(DimensionMismatch("Mismatched dimensions."))
    for j = 1:n
        xj = view(x,:,j)
        yj = view(y,:,j)
        for i = 1:m
            @inbounds yj[i] = a[i] * xj[i]
        end
    end
    return y
end

transform!{T<:Real}(a::StridedVector{T}, x::StridedMatrix) = scale!(a, x)

transform!{T<:Real}(y::StridedVector, a::StridedMatrix{T}, x::StridedVector) = A_mul_B!(y, a, x)
transform!{T<:Real}(y::StridedMatrix, a::StridedMatrix{T}, x::StridedMatrix) = A_mul_B!(y, a, x)

transform!{T<:Real}(y::StridedVector, a::Transpose{T}, x::StridedVector) = At_mul_B!(y, a.mat, x)
transform!{T<:Real}(y::StridedMatrix, a::Transpose{T}, x::StridedMatrix) = At_mul_B!(y, a.mat, x)

transform(a::Real, x::StridedVector) = a * x
transform(a::Real, x::StridedMatrix) = a * x

transform{T<:Real}(a::LinearTransform{T}, x::StridedVector) =
    transform!(Array(T, _outdim(a)), a, x)

transform{T<:Real}(a::LinearTransform{T}, x::StridedMatrix) =
    transform!(Array(T, _outdim(a), size(x,2)), a, x)

### Affine Transformation: y <- a * x + b

immutable AffineTransform{T,A<:LinearTransform}
    a::A
    b::Vector{T}
end

AffineTransform{T<:Real}(a::T) = AffineTransform{T,T}(a,T[])
AffineTransform{T<:Real}(a::T, b::Vector{T}) = AffineTransform{T,T}(a, b)

AffineTransform{T<:Real}(a::LinearTransform{T}) = AffineTransform{T, typeof(a)}(a, T[])

function AffineTransform{T<:Real}(a::LinearTransform{T}, b::Vector{T})
    if !(isempty(b) || length(b) == _outdim(a))
        throw(DimensionMismatch("Mismatched dimensions."))
    end
    return AffineTransform{T, typeof(a)}(a, b)
end

hasshift(t::AffineTransform) = !isempty(t.b)

function transform!(y::StridedVector, t::AffineTransform, x::StridedVector)
    transform!(y, t.a, x)
    if hasshift(t)
        add!(y, t.b)
    end
    return y
end

function transform!(y::StridedMatrix, t::AffineTransform, x::StridedMatrix)
    transform!(y, t.a, x)
    if hasshift(t)
        badd!(y, t.b, 1)
    end
    return y
end

function transform!{T<:Real,A<:SimpleLinearTransform}(t::AffineTransform{T,A}, x::StridedVector)
    transform!(t.a, x)
    if hasshift(t)
        add!(x, t.b)
    end
    return x
end

function transform!{T<:Real,A<:SimpleLinearTransform}(t::AffineTransform{T,A}, x::StridedMatrix)
    transform!(t.a, x)
    if hasshift(t)
        badd!(x, t.b, 1)
    end
    return x
end

transform{T<:Real}(t::AffineTransform{T,T}, x::StridedVecOrMat) =
    transform!(Array(T, size(x)), t, x)

transform{T<:Real}(t::AffineTransform{T}, x::StridedVector) =
    transform!(Array(T, _outdim(t.a)), t, x)

transform{T<:Real}(t::AffineTransform{T}, x::StridedMatrix) =
    transform!(Array(T, _outdim(t.a), size(x,2)), t, x)


text:
module Colors_2021

# Write your package code here.
using Markdown
using Pkg

md"""
## Load these Packages First - Needed for Images Manipulation
Let's use the Images.jl package to load an image and see what we get
"""

md"""
	Add These Packages ["Images", "ImageIO", "ImageMagick"]
"""
begin
	Pkg.add(["Images", "ImageIO", "ImageMagick", "TestImages"])
	using Images
end

md"""
## Open image, from URL
"""
# defines a variable called `url_variable`
# whose value is a string (written inside `"`):

Url_Variable = "https://emahesha.files.wordpress.com/2008/08/9indian_flag.jpg"

md"""
    download to a local file
"""
Indi_Flag_file = download(Url_Variable, "Flag.jpg")  # download to a local file

md"""
    Load the Image from Local Stored Image
"""
Tricolor = load(Indi_Flag_file)

end


text:
using MusicModels.MusicTransformer
using NoteSequences: SeqNote

@testset "MidiPerformanceEncoder" begin
    @testset "Encoding/Decoding" begin
        encoder = MidiPerformanceEncoder(100, 32, 21, 108, true)

        ns = NoteSequence()
        append!(ns.notes, [
            SeqNote(21, 127, 0, 220, 1, 0),
            SeqNote(50, 1, 110, 220, 1, 0),
            SeqNote(108, 100, 220, 440, 1, 0),
            SeqNote(80, 100, 330, 440, 1, 0),
        ])

        expected = [
            310,
            3,
            203,
            279,
            32,
            203,
            91,
            120,
            303,
            90,
            203,
            62,
            203,
            178,
            150,
            2
        ]

        @test expected == encode_notesequence(ns, encoder)
    end

    @testset "Misc" begin
        encoder = MidiPerformanceEncoder(100, 32, 21, 108, true)

        @test encode_notesequence(NoteSequence(), encoder) == [2]
        @test encoder.vocab_size == 310
        @test encoder.num_reserved_tokens == 2
    end
end

@testset "TextMelodyEncoder" begin
    @testset "Encoding/Decoding" begin
        encoder = TextMelodyEncoder(21, 108, 10)

        ns = NoteSequence()
        append!(ns.notes, [
            SeqNote(21, 127, 0, 220, 1, 0),
            SeqNote(107, 100, 220, 440, 1, 0),
            SeqNote(80, 100, 440, 660, 1, 0),
        ])

        expected = [
            5,
            3,
            3,
            3,
            4,
            91,
            3,
            3,
            3,
            4,
            64,
            3,
            3,
            3,
            3
        ]

        @test expected == encode_notesequence(ns, encoder)
    end

    @testset "Misc" begin
        encoder = TextMelodyEncoder(21, 109, 100)

        @test encode_notesequence(NoteSequence(), encoder) == Int[]
        @test encoder.vocab_size == 92
        @test encoder.num_reserved_tokens == 2
    end
end


text:
######### Stan batch program example  ###########

using Stan, Mamba

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

const dyes ="
data {
  int BATCHES;
  int SAMPLES;
  real y[BATCHES, SAMPLES];
  // vector[SAMPLES] y[BATCHES];
}

parameters {
  real<lower=0> tau_between;
  real<lower=0> tau_within;
  real theta;
  real mu[BATCHES];
}

transformed parameters {
  real sigma_between;
  real sigma_within;
  sigma_between <- 1/sqrt(tau_between);
  sigma_within <- 1/sqrt(tau_within);
}

model {
  theta ~ normal(0.0, 1E5);
  tau_between ~ gamma(.001, .001);
  tau_within ~ gamma(.001, .001);

  mu ~ normal(theta, sigma_between);
  for (n in 1:BATCHES)
    y[n] ~ normal(mu[n], sigma_within);
}

generated quantities {
  real sigmasq_between;
  real sigmasq_within;

  sigmasq_between <- 1 / tau_between;
  sigmasq_within <- 1 / tau_within;
}
"

const dyesdata = [
  Dict("BATCHES" => 6,
    "SAMPLES" => 5,
    "y" => reshape([
      [1545, 1540, 1595, 1445, 1595];
      [1520, 1440, 1555, 1550, 1440];
      [1630, 1455, 1440, 1490, 1605];
      [1595, 1515, 1450, 1520, 1560];
      [1510, 1465, 1635, 1480, 1580];
      [1495, 1560, 1545, 1625, 1445]
    ], 6, 5)
  )
]

stanmodel = Stanmodel(name="dyes", model=dyes);
@time sim1 = stan(stanmodel, dyesdata, ProjDir, CmdStanDir=CMDSTAN_HOME)

nodesubset = ["theta", "mu.1", "mu.2", "mu.3", "mu.4", "mu.5", "mu.6", "sigma_between", "sigma_within"]
## Subset Sampler Output
sim = sim1[1:size(sim1, 1), nodesubset, 1:size(sim1, 3)]

describe(sim)
println()

## Brooks, Gelman and Rubin Convergence Diagnostic
try
  gelmandiag(sim1, mpsrf=true, transform=true) |> display
catch e
  #println(e)
  gelmandiag(sim, mpsrf=false, transform=true) |> display
end

## Geweke Convergence Diagnostic
gewekediag(sim) |> display

## Highest Posterior Density Intervals
hpd(sim) |> display

## Cross-Correlations
cor(sim) |> display

## Lag-Autocorrelations
autocor(sim) |> display

## Deviance Information Criterion
#dic(sim) |> display

## Plotting

p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)

# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
        for i in 1:4
          isfile("$(stanmodel.name)-summaryplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-summaryplot-$(i).svg"`)
        end
      end : println()

end # cd


text:
module DayNNN
export dayNNN

function parseLine(line)
    line
end

function readInput(inputfile="inputs/inputNNN.txt")
    lines = readlines(inputfile)
    return lines
    # return [parseLine(line) for line in eachline(inputfile)]
end

function dayNNN(inputdata = readInput("inputs/inputNNN.txt"))
    return (dayNNNPart1(inputdata), dayNNNPart2(inputdata))
end

function dayNNNPart1(values)

end

function dayNNNPart2(values)

end

end #module


text:
include("../src/geom_module.jl")

module meshtest
using PyPlot
PyPlot.pygui(true)
using ..geom

rectangle = geom.Rectangle(-1.0, -1.0, 1.0, 1.0)
circle = geom.Circle(0.0, 0.0, 0.4)

distance_fcn = p -> geom.difference_signed_distance_function(
    geom.signed_distance_function(rectangle, p),
    geom.signed_distance_function(circle, p) )

src_point     = [0.0 0]
h_0           = 0.15
h_inf         = 0.5
spread_factor = 1.0
node_spacing_fcn(p) = geom.point_source_node_spacing(p, src_point, h_0, h_inf, spread_factor)
# node_spacing_fcn(p) = geom.uniform_node_spacing(p)

initial_spacing 	= 0.1
bounding_box 		= geom.Rectangle(-1.0, -1.0, 1.0, 1.0)
fixed_points 		= [
    -1.0 -1.0;
    +1.0 -1.0;
    +1.0 +1.0;
    -1.0 +1.0] # zeros(0,2)
max_iteration 		= 20

options = geom.MeshGenerationOptions(
    distance_fcn,
	node_spacing_fcn,
	initial_spacing,
	bounding_box,
	fixed_points;
	max_iteration = max_iteration)

(point_array, triangle_index, delmesh) = geom.create_mesh(options);

n_points = size(point_array, 1)

px = point_array[:,1]
py = point_array[:,2]
Lt = size(triangle_index,1)

PyPlot.figure()
PyPlot.grid()

PyPlot.plot( px, py,
    marker = :., linestyle = :none, markersize = 10)

for tt = 1:size(triangle_index,1)
    t = triangle_index[tt,:]
    tx = px[ [t[1], t[2], t[3], t[1]]]
    ty = py[ [t[1], t[2], t[3], t[1]]]
    PyPlot.plot(tx, ty)
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#does any data processign that is needed
#generally called from ComputeNCCS since it is cacheable
function prepsensitivityregressions!(data::NCCSData;
    toosmalltolog::Float64 = TOO_SMALL_TO_LOG,
    period::Symbol = :fisyr, group=:ein,)

  #logs if possible or returns missing
  plog(f::MFloat64)::MFloat64 = ismissing(f) || f ≤ toosmalltolog ? missing : log(f)

  #cleans up a common-size ratio
  cleanratio(f::MFloat64) = (!ismissing(f)) && isfinite(f) && f ≤ 1.0  && f > 0.0 ? f : missing

  #log our focal variables
  #data.df[:lcontributions] = (plog).(data.df[:contributions])
  #data.df[:lprogramexpense] = (plog).(data.df[:programexpense])
  #data.df[:lprogramexpense_avg] = (plog).(data.df[:programexpense_avg])

  #taking the average for control purposes only with valid contribution fields
  #data.df[:lregnetassets] = ((c::MFloat64, a::MFloat64)->
  #  ismissing(c) ? missing : a).(data.df[:contributions], data.df[:lnetassets])
  sort!(data.df, [group, period])
  lagwithin!(data, :netassets, sorted=true)
  #  lagwithin!(data, :totexp, sorted=true)

  #scale everythign to lagged net assets
  data.df[:pcontributions] = data.df[:contributions] ./ data.df[:Lnetassets]
  data.df[:pcontributions] .= (cleanratio).(data.df[:pcontributions])
  #data.df[:lpcontributions] = (plog).(data.df[:pcontributions])

  data.df[:nonzeroprogramexpense] = (f->
    (!ismissing(f)) && (f>0.0) ? f : missing ).(data.df[:programexpense])

  data.df[:pprogramexpense] = data.df[:programexpense] ./ data.df[:totexp]
    #data.df[:mgmtexpense] ./ data.df[:Lnetassets]
    #(data.df[:totexp] .- data.df[:nonzeroprogramexpense]) ./ data.df[:Lnetassets]
    #data.df[:programexpense] ./ data.df[:totexp] #./ data.df[:Lnetassets]
  data.df[:pprogramexpense] .= (cleanratio).(data.df[:pprogramexpense])
  #data.df[:lpprogramexpense] = (plog).(data.df[:pprogramexpense])

  #get the average expense
  fieldstoaggregate::Vector{Symbol} = [:pprogramexpense]
  aggdf::DataFrame = aggregatefielddf(data.df[[:ein;
    fieldstoaggregate]], rejoin=false)
  unique!(aggdf)

  #the setdiff below is to avoid duplicates
  fieldstonotmerge::Vector{Symbol} = (s->Symbol(s,"_avg")).(fieldstoaggregate)
  data.df = join(data.df[setdiff(names(data.df), fieldstonotmerge)],
    aggdf, on=[:ein], kind=:left)

  #resort after the join
  sort!(data.df, [group, period])
  #lag our focal variables
  #lagwithin!(data, :lcontributions)
  lagwithin!(data, :pcontributions, sorted=true)
  #lagwithin!(data, :lpcontributions, sorted=true)
  lagwithin!(data, :publicsupport, sorted=true)
  #lagwithin!(data, :lprogramexpense, sorted=true)
  lagwithin!(data, :pprogramexpense, sorted=true)
  #lagwithin!(data, :lpprogramexpense, sorted=true)
  lagwithin!(data, :lreturn3yr, sorted=true)
  lagwithin!(data, :lreturnnetbmsp5003yr, sorted=true)
  lagwithin!(data, :lreturnnetbmff33yr, sorted=true)

  createquantiles!(data.df, :pcontributions)
  createquantiles!(data.df, :pprogramexpense)
  createquantiles!(data.df, :lreturn3yr)
  createquantiles!(data.df, :lreturnnetbmsp5003yr)
  createquantiles!(data.df, :lreturnnetbmff33yr)

  differencewithin!(data, :pcontributions, sorted=true)
  createquantiles!(data.df, :Dpcontributions)

  #println(describe(data.df[[:ein, :lpprogramexpense, :lpprogramexpense_avg, :Llpprogramexpense]]))

  #generate the interactions
  data.df[:LpublicsupportXLpcont] = data.df[:Lpublicsupport] .* data.df[:Lpcontributions]
  data.df[:publicsupport_avgXLpcont] = data.df[:publicsupport_avg] .* data.df[:Lpcontributions]

  data.df[:LpprogramexpenseXLpcont] = data.df[:Lpprogramexpense] .* data.df[:Lpcontributions]
  data.df[:pprogramexpense_avgXLpcont] = data.df[:pprogramexpense_avg] .* data.df[:Lpcontributions]

  data.df[:publicsupportXpcont] = data.df[:publicsupport] .* data.df[:pcontributions]
  data.df[:publicsupport_avgXpcont] = data.df[:publicsupport_avg] .* data.df[:pcontributions]

  data.df[:pprogramexpenseXpcont] = data.df[:pprogramexpense] .* data.df[:pcontributions]
  data.df[:pprogramexpense_avgXpcont] = data.df[:pprogramexpense_avg] .* data.df[:pcontributions]

  #I think regressions work better with symbols instead of strings
  rename!(data.df, :ein=>:ein_old)
  data.df[:ein] = (Symbol).(data.df[:ein_old])
  deletecols!(data.df, :ein_old)

  return nothing
end

#makes a scatterplot of betas
function plotβ!(df::DataFrame;
  xfield::Symbol = :pprogramexpense_avg,
  yfield::Symbol = :b1lRet3yrNetSP500Con,#:b1lRet3yrCon,
  suffix::String = "",
  outputpath::String = GRAPH_PATH,
  shortnameindex::Dict = SHORT_NAME_INDEX,
  computequantiles::Bool = true
  )

  local plotnames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local title::String

  if computequantiles
    df = deepcopy(df)
    createquantiles!(df, yfield, groupfields = nothing)
    #    createquantiles!(df, xfield, groupfields = nothing)
    #    xfield = Symbol(:q, xfield)
    yfield = Symbol(:q, yfield)
    push!(plotnames, "betaplot-$(shortnameindex[yfield])V$(shortnameindex[xfield])$suffix")
    title = "Return-Contrib. Beta ($(shortnameindex[yfield]))
      vs. Prog. Expense Percentile ($(shortnameindex[xfield])$suffix)"

  else
    push!(plotnames, "betaplot-$(shortnameindex[yfield])V$(shortnameindex[xfield])$suffix")
    title = "Return-Contribution Beta ($(shortnameindex[yfield]))
      vs. Program Expense ($(shortnameindex[xfield])$suffix)"
  end

  push!(plots, plot(df[completecases(df[[xfield,yfield]]),:],
          x=xfield, y = yfield,
          Guide.title("Return-Contribution Beta ($(shortnameindex[yfield]))
            vs. Program Expense ($(shortnameindex[xfield])$suffix)"),
          Guide.xlabel("Governance ($(shortnameindex[xfield]))"),
          Guide.ylabel("density$(shortnameindex[yfield])"),
          Geom.hexbin(xbincount=25, ybincount=25),
          #Geom.point,
          #Coord.cartesian(xmin=0.0, xmax=1.0)
          ))

  for i ∈ 1:length(plotnames)
    draw(PDF("$outputpath\\$(plotnames[i]).pdf", 9inch, 7inch), plots[i])
  end

  return nothing
end

#primarily generates the regeression table
function olssensitivitypanel(df::AbstractDataFrame;
  decimals::Int=DECIMALS,
  suffix::String = "",
  displayedvars::Union{Vector{Symbol},Nothing} = nothing,
  contentRowNames::Union{Vector{String},Nothing}=nothing,
  displayedvarsscaling::Union{Vector{Float64},Nothing} = nothing,
  regressionfunction::Function = olssensitivitypanelreg,
  titleCaption::String = "TBD",
  outputpath::String = TABLE_PATH,
  errorfunction::T where T<:Union{Function, Vector{Function}} = getClustered!,
  scalingvalue::Float64 = 100.0)

  #first run the regressions
  local specs::FMSpecs = regressionfunction(df)

  #apply the defaults based on the sepcs
  isnothing(displayedvars) &&  (displayedvars = setdiff(unique([specs.xnames...;]), [:intercept]))
  isnothing(contentRowNames) &&  (contentRowNames = (string).(displayedvars))
  isnothing(displayedvarsscaling) && (displayedvarsscaling = fill(scalingvalue, length(displayedvars)))

  local descRowNames::Vector{String} = [
    "Nonprofit(NP) FE",
    "Year FE",
    "NP Clustered",
    "\$R^{2}\$ (within)",
    "N"]

  #allocate space for the descriptive rows
  local descContent::Vector{Vector{String}} =
    ((i::Int)->Vector{String}(undef, specs.N[])).(1:length(descRowNames))

    #this builds the descriptive rows. There is Probably a better way to do this,
  #but its fairly specific to the project.
  for i ∈ 1:specs.N[]
    r::Int = 0
    descContent[r+=1][i] = "$(specs.withinspecs[i]==:ein ||
        occursin("ein", string(specs.xspecs[i])) ? "X" : "")"
    descContent[r+=1][i] = "$(specs.withinspecs[i]==:fisyr ||
        occursin("fisyr", string(specs.xspecs[i])) ? "X" : "")"
    descContent[r+=1][i] = "$(specs.clusterspecs[i]==:ein ? "X" : "")"
    descContent[r+=1][i] =
      "$(num2Str(getR²(specs.results[i], adjusted=false),  decimals))"
    descContent[r+=1][i] = "$(specs.results[i].N)"
  end

  #NOTE delete this if the console gets too cluttered
  #display(descContent)

  local colNames::Vector{Vector{String}} = [
    ["return", "SP500 Adj", "FF3 Adj"],
    specs.specnames]

  local widthColNames::Vector{Vector{Int}} = [[2,2,2], ones(Int, specs.N[])]

  local alignmentColNames::Vector{Vector{String}} =
    [fill("c",length(widthColNames[1])), fill("r", specs.N[])]

  #checks based on the dimensions of the table parameters
  @assert (minimum((sum).(widthColNames)) == specs.N[] &&
    maximum((sum).(widthColNames)) == specs.N[] &&
    length(contentRowNames) == length(displayedvars)
    )

  local tabletext::String = texTable(specs.results,
    errorfunction#=getNeweyWestFunc(5) (lm::CTLM)->getNeweyWestSlow(lm, 5)=#,
    displayedvars,
    titleCaption = titleCaption,
    colNames = colNames,
    widthColNames = widthColNames,
    contentRowNames = contentRowNames,
    descRowNames = descRowNames,
    descContent = descContent,
    decimalDigits = 2,
    scaling= displayedvarsscaling,
    stars=true,
    starStrings = OVERRIDE_STAR_STRINGS,
    #clearMem = USE_AGGRESSIVE_GC,
    caption = "to be written",
    nakedTable=true)

  writeNakedTable(tabletext, path=outputpath, outName="olspanel$suffix.tex")
  println("Regressions for olspanel$suffix.tex complete and recorded.")

  return tabletext
end

function prefilterforolssensitivity(df::DataFrame;
  minpointspluskforβ::Int = MIN_POINTS_PLUS_K_FOR_BETA,
  essentialfields::Vector{Symbol}=[#=:Llcontributions,=# :Lpcontributions])::DataFrame

  local subdf::SubDataFrame
  local ssubdf::SubDataFrame

  #make sure we ahve enough unique points for each organization
  df[:useforregressions] = trues(size(df,1))

  checkfields::Vector{Symbol} = [:lreturn; essentialfields; :useforregressions]

  #loop over all unique eins
  for subdf ∈ groupby(df, :ein)

    ssubdf = view(subdf, completecases(subdf[checkfields]), checkfields)
    subN = minimum((f->size(unique(ssubdf[f]),1)).(essentialfields))

    #make sure we have enough points for the regression
    if subN < minpointspluskforβ + 2 #add two for the existing dof
      ssubdf[:useforregressions] .= false
    end
  end

  return df[df[:useforregressions],:]
end

function makeolspaneltables(data::NCCSData; rerunprep::Bool = true)

  #NOTE: This is also called in the compute function
  #if nothing is changed it does not need to be run
  rerunprep && prepsensitivityregressions!(data)

  #can disable the below

  local displayedvars::Vector{Symbol}

  #this gets repeated a lot so we can save time doing it now
  regdf::DataFrame = prefilterforolssensitivity(data.df)

  #expense regressions
  displayedvars = [:LagContribAvgExp, :LagContribExp, :LagContrib]
  olssensitivitypanel(regdf, suffix = "programexpense",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olssensitivitypanelregPE,
    titleCaption = "Program Expense vs Contribution Return Sensitivity")

  displayedvars = [:ContribAvgExp, :ContribExp, :Contrib]
  olssensitivitypanel(regdf, suffix = "programexpensehist",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olssensitivitypanelregPEhist,
    titleCaption = "Governance vs Contribution Lagged 3YReturn Sensitivity")
  #NOTE: Failed proxy: Public Support
  #=displayedvars = [:LagContribAvgSup, :LagContribSup, :LagContrib]
  olssensitivitypanel(regdf, suffix = "publicsupport",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olssensitivitypanelregPS,
    titleCaption = "PublicSupport vs Contribution Return Sensitivity")=#

  #Univariate Regressions
  displayedvars = [:intercept, :LagContrib, :LagExp]
  olssensitivitypanel(regdf, suffix = "univariate",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olscontributionreg,
    titleCaption = "Program Expense vs Contribution Return Sensitivity")

  #pick the standard error functions
  nwerror::Function = getNeweyWestFunc(3)
  mwerror::Function = getWhiteΣ!
  errorfunctions::Vector{Function} = [mwerror, nwerror, mwerror, nwerror, mwerror, nwerror]
  olssensitivitypanel(regdf, suffix = "2pass",
    #displayedvars = nothing,
    #contentRowNames = nothing,
    regressionfunction = olssensitivity2pass,
    errorfunction = errorfunctions,
    titleCaption = "Effect of program expense on sensitivity to expected returns")
  olssensitivitypanel(regdf,
    suffix = "2passhist",
    #displayedvars = nothing,
    #contentRowNames = nothing,
    regressionfunction = olssensitivity2passhist,
    errorfunction = errorfunctions,
    titleCaption = "Effect of program expense on sensitivity to historical returns")

  tab::CrossTabTable = CrossTabTable(regdf, tabfield=:qDpcontributions)
  writecrosstabtable(tab)

end

function olssensitivitytables(data::NCCSData)
  makeolspaneltables(data)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text: using ExprSearch, POMDPs

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
@testset "mixed complementarity problem" begin
    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    myfunc(x) = M*x + q

    n = 4
    lb = zeros(n)
    ub = 100*ones(n)

    status, z, f = solveMCP(myfunc, lb, ub)
    @show status
    @show z
    @show f

    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    println("-------------------------------------------------------")

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    myfunc(x) = M*x + q

    n = 4
    lb = zeros(n)
    ub = 100*ones(n)

    options(convergence_tolerance=1e-2, output=:no, time_limit=3600)

    status, z, f = solveMCP(myfunc, lb, ub)

    @show status
    @show z
    @show f

    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    println("-------------------------------------------------------")

    function elemfunc(x)
        val = similar(x)
        val[1] = -x[3]-x[4] + q[1]
        val[2] = x[3] -2x[4] + q[2]
        val[3] = x[1]-x[2]+2x[3]-2x[4] + q[3]
        val[4] = x[1]+2x[2]-2x[3]+4x[4] + q[4]
        return val
    end

    n = 4
    lb = zeros(n)
    ub = 100*ones(n)

    var_name = ["first var", "second var", "third var", "fourth var"]
    con_name = ["func 1", "func 2", "func 3", "func 4"]

    options(convergence_tolerance=1e-2, output=:yes, time_limit=3600)

    jacfunc(x) = M

    status, z, f = solveMCP(elemfunc, lb, ub)
    status, z, f = solveMCP(elemfunc, lb, ub, var_name)
    status, z, f = solveMCP(elemfunc, lb, ub, var_name, con_name)
    status, z, f = solveMCP(elemfunc, jacfunc, lb, ub)
    status, z, f = solveMCP(elemfunc, jacfunc, lb, ub, var_name)
    status, z, f = solveMCP(elemfunc, jacfunc, lb, ub, var_name, con_name)

    @show status
    @show z
    @show f

    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    function test_in_local_scope()
        # Verify that we can solve MCPs in local scope. Surprisingly, this is
        # is relevant because it affects the way closures are generated. To be
        # specific, you can do the following in global scope:
        #
        # julia> y = [1]
        # julia> cfunction(x -> x + y[1], Int, (Int,))
        #
        # but running the same code inside a function will fail with:
        #   ERROR: closures are not yet c-callable

        M = [0  0 -1 -1 ;
             0  0  1 -2 ;
             1 -1  2 -2 ;
             1  2 -2  4 ]

        q = [2; 2; -2; -6]

        myfunc(x) = M*x + q

        n = 4
        lb = zeros(n)
        ub = 100*ones(n)

        options(convergence_tolerance=1e-2, output=:yes, time_limit=3600, lemke_start=:first, nms_searchtype=:line)

        status, z, f = solveMCP(myfunc, lb, ub)

        @show status
        @show z
        @show f

        @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
        @test status == :Solved
    end

    test_in_local_scope()
end