text:
"""
`@generated function nfcalls(::Val{N}, f)`

Effectively generates a loop of functions calls, but fully unrolled and
therefore type-stable:
```{julia}
f(Val(1))
f(Val(2))
f(Val(3))
# ...
f(Val(N))
```
"""
@generated function nfcalls(::Val{N}, f) where {N}
   code = Expr[]
   for n = 1:N
      push!(code, :(f(Val($n))))
   end
   quote
      $(Expr(:block, code...))
      return nothing
   end
end

"""
`@generated function valnmapreduce(::Val{N}, v, f)`

Generates a map-reduce like code, with fully unrolled loop which makes this
type-stable,
```{julia}
begin
   v += f(Val(1))
   v += f(Val(2))
   # ...
   v += f(Val(N))
   return v
end
```
"""
@generated function valnmapreduce(::Val{N}, v, f) where {N}
   code = Expr[]
   for n = 1:N
      push!(code, :(v += f(Val($n))))
   end
   quote
      $(Expr(:block, code...))
      return v
   end
end

# ----------------------------------------------------------------------
# sparse matrix multiplication with weaker type restrictions

function _my_mul!(C::AbstractVector, A::SparseMatrixCSC, B::AbstractVector)
   A.n == length(B) || throw(DimensionMismatch())
   A.m == length(C) || throw(DimensionMismatch())
   nzv = A.nzval
   rv = A.rowval
   fill!(C, zero(eltype(C)))
   @inbounds for col = 1:A.n
      b = B[col]
      for j = A.colptr[col]:(A.colptr[col + 1] - 1)
         C[rv[j]] += nzv[j] * b
      end
   end
   return C
end

# # gradients
#
# site_alloc_dB(basis::PIBasis, Rs::AbstractVector, args...) =
#    site_alloc_dB(basis, length(Rs))
#
# site_alloc_dB(basis::PIBasis, nmax::Integer) =
#       zeros( JVec{fltype(basis)}, maximum(length.(basis.inner)), nmax )
#
# alloc_temp_d(basis::PIBasis, Rs::AbstractVector, args...) =
#    alloc_temp_d(basis, length(Rs))
#
# alloc_temp_d(basis::PIBasis, nmax::Integer) =
#       (
#         A = alloc_B(basis.basis1p, nmax),
#         dA = alloc_dB(basis.basis1p, nmax),
#         tmp_basis1p = alloc_temp(basis.basis1p, nmax),
#         tmpd_basis1p = alloc_temp_d(basis.basis1p, nmax),
#         tmpd_inner = alloc_temp_d(basis.evaluator, basis, nmax)
#       )
#
#
# function evaluate_d!(AA, dAA, tmpd, basis::PIBasis, Rs, Zs, z0)
#    iz0 = z2i(basis, z0)
#    AAview = @view AA[basis.inner[iz0].AAindices]
#    dAAview = @view dAA[basis.inner[iz0].AAindices, 1:length(Rs)]
#    site_evaluate_d!(AAview, dAAview, tmpd, basis, Rs, Zs, z0)
#    return dAA
# end
#
# function site_evaluate_d!(AA, dAA, tmpd, basis::PIBasis, Rs, Zs, z0)
#    iz0 = z2i(basis, z0)
#    # precompute the 1-p basis and its derivatives
#    evaluate_d!(tmpd.A, tmpd.dA, tmpd.tmpd_basis1p, basis.basis1p, Rs, Zs, z0)
#    site_evaluate_d!(AA, dAA, tmpd, basis.inner[iz0], basis.evaluator, tmpd.A, tmpd.dA)
#    return nothing
# end
#
#
# #--------------------------------------------------------
# # ------- specialised code for the standard evaluator
#
# alloc_temp(::StandardEvaluator, basis::PIBasis) = nothing
#
# alloc_temp_d(::StandardEvaluator, basis::PIBasis, args...) = nothing
#
# function site_evaluate_d!(AA, dAA, tmpd, inner::InnerPIBasis,
#                           _ev::StandardEvaluator, A, dA)
#    # evaluate the AA basis
#    evaluate!(AA, nothing, inner, _ev, tmpd.A)
#    # loop over all neighbours
#    for j = 1:size(dA, 2)
#       # write the gradients into the correct slice of the dAA matrix
#       for iAA = 1:length(inner)
#          @inbounds dAA[iAA, j] = grad_AAi_Rj(iAA, j, inner, A, dA)
#       end
#    end
#    return nothing
# end
#
#
#
# """
# Compute ∂∏A_a / ∂A_b = ∏_{a ≂̸ b} A_a
# """
# function grad_AAi_Ab(iAA, b, inner, A)
#    g = one(eltype(A))
#    @inbounds for a = 1:inner.orders[iAA]
#       if a != b
#          g *= A[inner.iAA2iA[iAA, a]]
#       end
#    end
#    return g
# end
#
# @doc raw"""
# Compuate ∂∏A_a / ∂Rⱼ:
# ```math
#    \frac{\partial \prod_a A_a}{\partial \bm r_j}
#    =
#    \sum_b \frac{\prod_{a} A_a}{\partial A_b} \cdot \frac{\partial A_b}{\partial \bm r_j}
#    =
#    \sum_b \prod_{a \neq b} A_a \cdot \frac{\partial \phi_b}{\partial \bm r_j}
# ```
# """
# function grad_AAi_Rj(iAA, j, inner, A, dA)
#    g = zero(eltype(dA))
#    @inbounds for b = 1:inner.orders[iAA] # interaction order
#       # A_{k_b} = A[iA]
#       iA = inner.iAA2iA[iAA, b]
#       # dAAi_dAb = ∂(∏A_{n_a}) / ∂A_{n_b}
#       dAAi_dAb = grad_AAi_Ab(iAA, b, inner, A)
#       g += dAAi_dAb * dA[iA, j]
#    end
#    return g
# end
#
#
# # --------------------------------------------------------
# #   Evaluation code for DAG evaluator
#
#
# alloc_temp(::DAGEvaluator, basis::PIBasis) =
#    zeros(fltype(basis), maximum(inner.dag.numstore for inner in basis.inner))
#
#
# function alloc_temp_d(::DAGEvaluator, basis::PIBasis,  nmax::Integer)
#    maxstore = maximum(inner.dag.numstore for inner in basis.inner)
#    return (
#       AA = zeros(fltype(basis), maxstore),
#       dAA = zeros(JVec{fltype(basis)}, maxstore, nmax)
#    )
# end
#
# function evaluate!(AA, tmp, basis::InnerPIBasis, ::DAGEvaluator, A)
#    AAdag = tmp.tmp_inner
#    fill!(AAdag, 1)
#    traverse_fwd!(AAdag, basis.dag, A,
#                  (idx, AAval) -> if idx > 0; AA[idx] = AAval; end)
#    return AA
# end
#
#
#
# function site_evaluate_d!(AA, dAA, tmpd, inner::InnerPIBasis,
#                           ::DAGEvaluator, A, dA)
#    dag = inner.dag
#    nodes = dag.nodes
#    idxs = dag.vals
#    AAtmp = tmpd.tmpd_inner.AA
#    dAAtmp = tmpd.tmpd_inner.dAA
#    numneigs = size(dA, 2)
#    # --- manual bounds checking
#    @assert size(dAA, 2) >= numneigs
#    @assert size(dAAtmp, 2) >= numneigs
#    @assert length(AAtmp) >= inner.dag.numstore
#    @assert size(dAAtmp, 1) >= inner.dag.numstore
#    # ------------------
#
#    @inline function _copyrow!(A_, iA, B_, iB)
#       @simd for j = 1:numneigs
#          @inbounds A_[iA,j] = B_[iB,j]
#       end
#    end
#
#    @inline function _fwdgrad!(dAA_, i, n1, AA1, n2, AA2)
#       @simd for j = 1:numneigs
#          @inbounds dAA_[i, j] = AA1 * dAAtmp[n2,j] + AA2 * dAAtmp[n1,j]
#       end
#    end
#
#    @inbounds for i = 1:dag.num1
#       idx = idxs[i]
#       AAtmp[i] = A[i]
#       _copyrow!(dAAtmp, i, dA, i)
#       if idx > 0
#          AA[idx] = AAtmp[i]
#          _copyrow!(dAA, idx, dAAtmp, i)
#          # @. dAA[idx,:] = (@view dAAtmp[i,:])
#       end
#    end
#
#    @inbounds for i = (dag.num1+1):dag.numstore
#       n1, n2 = nodes[i]
#       idx = idxs[i]
#       AA1, AA2 = AAtmp[n1], AAtmp[n2]
#       AAtmp[i] = AA1 * AA2
#       _fwdgrad!(dAAtmp, i, n1, AA1, n2, AA2)
#       # @. dAAtmp[i, :] = AA1 * (@view dAAtmp[n2,:]) + AA2 * (@view dAAtmp[n1,:])
#       if idx > 0
#          AA[idx] = AAtmp[i]
#          _copyrow!(dAA, idx, dAAtmp, i)
#          # @. dAA[idx, :] .= (@view dAAtmp[i, :])
#       end
#    end
#
#    @inbounds for i = (dag.numstore+1):length(dag)
#       idx = idxs[i]
#       n1, n2 = nodes[i]
#       # @assert idx > 0
#       AA1, AA2 = AAtmp[n1], AAtmp[n2]
#       AA[idx] = AA1 * AA2
#       _fwdgrad!(dAA, idx, n1, AA1, n2, AA2)
#       # @simd for j = 1:numneigs
#       #    @inbounds dAA[idx, j] = AA1 * dAAtmp[n2,j] + AA2 * dAAtmp[n1,j]
#       # end
#    end
#
#    return nothing
# end

   # ----- compute the adjoints dAA / dA
   for iAA = 1:length(basis)
      dAAdA[1] = 1
      AAfwd = A[iAA2iA[iAA, 1]]
      ord = orders[iAA]
      for a = 2:ord
         dAAdA[a] = AAfwd
         AAfwd *= A[iAA2iA[iAA, a]]
      end
      AA[iAA] = AAfwd
      AAbwd = A[iAA2iA[iAA, ord]]
      for a = ord-1:-1:1
         dAAdA[a] *= AAbwd
         AAbwd *= A[iAA2iA[iAA, a]]
      end
      for a = 1:ord
         W[iAA2iA[iAA, a]] += dAAdA[a]
      end
   end


text:
"""
    baggr_oob(object::Baggr, X, Y; score = rmsep)
Compute the out-of-bag (OOB) error after bagging a regression model.
* `object` : Output of a bagging.
* `X` : X-data used in the bagging.
* `Y` : Y-data used in the bagging.
* `score` : Function computing the prediction error (default: msep).

## References

Breiman, L., 1996. Bagging predictors. Mach Learn 24, 123–140.
https://doi.org/10.1007/BF00058655

Breiman, L., 2001. Random Forests. Machine Learning 45, 5–32.
https://doi.org/10.1023/A:1010933404324

Genuer, R., 2010. Forêts aléatoires : aspects théoriques,
sélection de variables et applications. PhD Thesis. Université Paris Sud - Paris XI.

Gey, S., 2002. Bornes de risque, détection de ruptures, boosting :
trois thèmes statistiques autour de CART en régression (These de doctorat).
Paris 11. http://www.theses.fr/2002PA112245
"""
function baggr_oob(object::Baggr, X, Y; score = rmsep)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    n = size(X, 1)
    q = size(Y, 2)
    rep = length(object.fm)
    pred_oob = similar(X, n, q)
    k = zeros(n)
    @inbounds for i = 1:n
        zpred = fill(0., 1, q)
        zk = 0
        @inbounds for j = 1:rep
            soob = object.s_oob[j]
            if in(i, soob)
                zpred .+= predict(object.fm[j], X[i:i, vcol(object.s_col, j)]).pred
                zk += 1
            end
        end
        pred_oob[i, :] = zpred / zk
        k[i] = zk
    end
    scor = score(pred_oob, Y)
    (scor = scor, pred_oob, k)
end

"""
    baggr_vi(object::Baggr, X, Y; score = rmsep)
Compute variance importances (with permutation method) after bagging a regression model.
* `object` : Output of a bagging.
* `X` : X-data that was used in the model bagging.
* `Y` : Y-data that was used in the model bagging.

Variances importances are computed
by permuting sucessively each column of the out-of-bag (X_OOB),
and by looking at the effect on the MSEP.

## References

Breiman, L., 2001. Random Forests. Machine Learning 45, 5–32.
https://doi.org/10.1023/A:1010933404324

Genuer, R., 2010. Forêts aléatoires : aspects théoriques,
sélection de variables et applications. PhD Thesis. Université Paris Sud - Paris XI.

Gey, S., 2002. Bornes de risque, détection de ruptures, boosting :
trois thèmes statistiques autour de CART en régression. PhD Thesis.
Univ. Paris 11. http://www.theses.fr/2002PA112245
"""
function baggr_vi(object::Baggr, X, Y; score = rmsep)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    p = size(X, 2)
    q = size(Y, 2)
    rep = length(object.s_oob)
    ncol = size(object.s_col, 1)
    scol = similar(object.s_col, ncol)
    res = similar(X, p, rep, q)
    @inbounds for i = 1:rep
        soob = object.s_oob[i]  # soob has a variable length
        m = length(soob)
        scol .= vcol(object.s_col, i)
        zpred = predict(object.fm[i], X[soob, scol]).pred
        zY = @view(Y[soob, :]) ;
        zscore = score(zpred, zY)
        zX = similar(X[soob, :])
        @inbounds for j = 1:p
            zX .= X[soob, :]
            s = sample(1:m, m, replace = false)
            zX[:, j] .= zX[s, j]
            zpred .= predict(object.fm[i], zX[:, scol]).pred
            zscore_perm = score(zpred, zY)
            res[j, i, :] = zscore_perm .- zscore
        end
    end
    imp = reshape(mean(res, dims = 2), p, q)
    (imp = imp, res = res)
end

ext: jl


text:
# op1 is the first operator to be applied to a vector.
# compose_and_simplify(op1, op2, op3, ops...) = Base.afoldl(compose_and_simplify, compose_and_simplify(compose_and_simplify(op1,op2),op3), ops...)
compose_and_simplify(op1::AbstractOperator, op2::AbstractOperator, op3::AbstractOperator, ops::AbstractOperator...) =
    _compose_and_simplify(op1::AbstractOperator, op2::AbstractOperator, op3::AbstractOperator, ops...)

function _compose_and_simplify(ops::AbstractOperator...)
    ops = collect(ops)
    unchanged = false
    while !unchanged
        unchanged = true
        i = 1
        if length(ops) == 1
            break;
        end
        while i <= length(ops)-1
            changed, result = _compose_and_simplify(ops[i], ops[i+1])
            unchanged = !changed
            if changed
                if result isa AbstractOperator
                    ops[i] = result
                    deleteat!(ops, i+1)
                else
                    ops[i] = result[1]
                    ops[i+1] = result[2]
                end
            end
            i += 1
        end
    end
    tuple(ops...)
end

cas_src(op1, op2) = src(op1)
cas_dest(op1, op2) = dest(op2)

assert_compatible_compose_src_dest(op1, op2) = iscompatible(dest(op1), src(op2))

iterate(x::AbstractOperator) = (x, nothing)
iterate(x::AbstractOperator, ::Any) = nothing

function compose_and_simplify(op1::DictionaryOperator, op2::DictionaryOperator)
    assert_compatible_compose_src_dest(op1, op2)
    result = unsafe_compose_and_simplify(op1, op2)
    if result isa Tuple && result[1] isa Bool
        result[2]
    else
        tuple(result...)
    end
end

function _compose_and_simplify(op1::DictionaryOperator, op2::DictionaryOperator)
    assert_compatible_compose_src_dest(op1, op2)
    result = unsafe_compose_and_simplify(op1, op2)
    if result isa Tuple && result[1] isa Bool
        false, result[2]
    else
        true, result
    end
end

unsafe_compose_and_simplify(op1, op2) =
    unsafe_compose_and_simplify1(op1, op2)

unsafe_compose_and_simplify1(op1::DictionaryOperator, op2) =
        unsafe_compose_and_simplify2(op1, op2)

unsafe_compose_and_simplify2(op1, op2::DictionaryOperator) =
        default_unsafe_compose_and_simplify(op1, op2)

default_unsafe_compose_and_simplify(op1, op2) = (false, (op1, op2))

unsafe_compose_and_simplify(op1::IdentityOperator, op2::IdentityOperator) = op1
unsafe_compose_and_simplify1(op1::IdentityOperator, op2) = op2
unsafe_compose_and_simplify2(op1, op2::IdentityOperator) = op1

# Move scalars to last place
function unsafe_compose_and_simplify1(op1::ScalingOperator{T}, op2) where {T}
    if scalar(op1) == 1
        op2
    else
        op2, ScalingOperator{T}(dest(op2), scalar(op1))
    end
end

# Combine scalars to one / not type stable
function unsafe_compose_and_simplify(op1::ScalingOperator, op2::ScalingOperator)
    T = promote_type(eltype(op1),eltype(op2))
    z = scalar(op1)*scalar(op2)
    if z ≈ 1 && iscompatible(cas_src(op1,op2), cas_dest(op1,op2))
        IdentityOperator{T}(cas_src(op1,op2), cas_dest(op1,op2))
    else
        ScalingOperator{T}(cas_src(op1,op2), cas_dest(op1,op2), z)
    end
end

for (OP) in (:DiagonalOperator, :CirculantOperator)
    @eval begin
        unsafe_compose_and_simplify(op1::$OP, op2::$OP) = ArrayOperator(_checked_mul(unsafe_matrix(op2), unsafe_matrix(op1)), cas_src(op1,op2), cas_dest(op1,op2))
    end
end

function unsafe_compose_and_simplify(op1::IndexExtension, op2::IndexRestriction)
    if subindices(op1) == subindices(op2)
        IdentityOperator(cas_src(op1,op2), cas_dest(op1,op2))
    else
        default_unsafe_compose_and_simplify(op1, op2)
    end
end

unsafe_compose_and_simplify(op1::DiagonalOperator, op2::ScalingOperator) =
    DiagonalOperator(cas_src(op1,op2),cas_dest(op1,op2),Diagonal(diag(op1)*scalar(op2)))

unsafe_compose_and_simplify(op1::ScalingOperator, op2::DiagonalOperator) =
    DiagonalOperator(cas_src(op1,op2),cas_dest(op1,op2),Diagonal(diag(op2)*scalar(op1)))


text:
abstract type NetworkParams{T<:Real} <: AbstractParams end
abstract type TwoPortParams{T<:Real} <: NetworkParams{T} end

check_row_column(data::Array{T, 2}) where {T<:Number} =
    begin nr, nc = size(data)
        (nr == nc)? nr: error("Number of rows and columns doesn't match")
    end

"""
    Sparams(nPort, data) <: NetworkParams
Scattering parameters for microwave network.
"""
mutable struct Sparams{T<:Real} <: NetworkParams{T}
    nPort::Int
    data::Matrix{Complex{T}}
    Sparams(data::Matrix{Complex{T}}) where {T<:Real} =
        new{T}(check_row_column(data), data)
end
Sparams(data::Matrix{T}) where {T<:Real} = Sparams(complex(data))
Sparams(P::Vector{Matrix{T}}) where {T<:Number} = [Sparams(p) for p in P]

"""
    Yparams(nPort, data) <: NetworkParams
Admittance parameters for microwave network
"""
mutable struct Yparams{T<:Real} <: NetworkParams{T}
    nPort::Int
    data::Matrix{Complex{T}}
    Yparams(data::Matrix{Complex{T}}) where {T<:Real} =
        new{T}(check_row_column(data), data)
end
Yparams(data::Matrix{T}) where {T<:Real} = Yparams(complex(data))
Yparams(P::Vector{Matrix{T}}) where {T<:Number} = [Yparams(p) for p in P]

"""
    Zparams(nPort, data) <: NetworkParams
Impedance-parameters for microwave network
"""
mutable struct Zparams{T<:Real} <: NetworkParams{T}
    nPort::Int
    data::Matrix{Complex{T}}
    Zparams(data::Matrix{Complex{T}}) where {T<:Real} =
        new{T}(check_row_column(data), data)
end
Zparams(data::Matrix{T}) where {T<:Real} = Zparams(complex(data))
Zparams(P::Vector{Matrix{T}}) where {T<:Number} = [Zparams(p) for p in P]

"""
    ABCDparams(nPort, data) <: TwoPortParams
Transfer-parameters for microwave network
"""
mutable struct ABCDparams{T<:Real} <: TwoPortParams{T}
    nPort::Int
    data::Matrix{Complex{T}}
    ABCDparams(data::Matrix{Complex{T}}) where {T<:Real} = begin
        n = check_row_column(data)
        (n == 2)? new{T}(n, data): error("ABCDparams are defined only for 2×2 matrices")
    end
end
ABCDparams(data::Matrix{T}) where {T<:Real} = ABCDparams(complex(data))
ABCDparams(P::Vector{Matrix{T}}) where {T<:Number} = [ABCDparams(p) for p in P]

"""
    show(io::IO, params::NetworkParams)
Pretty-printing of `NetworkParams`
"""
function show(io::IO, params::NetworkParams)
    write(io, "$(params.nPort)-port $(typeof(params))\n")
    for nr in 1:(params.nPort)
        for nc in 1:(params.nPort)
            write(io, "$((params.data)[nr, nc])\t")
        end
        write(io, "\n")
    end
end


text:
# belief
type RevealedGraph <: Belief
    revealed_edges::Set{Edge}
    revealed_adj::Vector{Vector{Int}}
    revealed_nodes::Set{Int}
    known_labels::Dict{Int,Int} # constant
    probed::Set{Int}
    # highest_degree_not_probed::Int # -1 if unknown
    budget::Int
    p_intra::Float64
    p_inter::Float64
    n_comms::Int
    RevealedGraph() = new()
end
n_nodes(r::RevealedGraph) = length(r.revealed_adj)

create_belief(p::CommunityAffiliationPOMDP) = RevealedGraph()

function initial_belief(p::CommunityAffiliationPOMDP, r::RevealedGraph=RevealedGraph())
    r.revealed_edges = copy(p.known_edges)
    r.revealed_nodes = copy(p.known_nodes)
    r.known_labels = copy(p.known_labels)
    r.probed = Set{Int}()
    r.budget = p.budget
    r.p_intra = p.p_intra
    r.p_inter = p.p_inter
    r.n_comms = p.true_network.n_comms
    r.revealed_adj = [Array(Int,0) for i in 1:length(p)]
    for ed in r.revealed_edges
        push!(r.revealed_adj[ed.first], ed.second)
        push!(r.revealed_adj[ed.second], ed.first)
    end
    node_vector = collect(r.revealed_nodes)
    # r.highest_degree_not_probed = node_vector[indmax([degree(r,n) for n in node_vector])]
    return r
end

function clone!(new::RevealedGraph, b::RevealedGraph)
    new.revealed_edges = copy(b.revealed_edges)
    new.revealed_adj = copy(b.revealed_adj)
    new.revealed_nodes = copy(b.revealed_nodes)
    new.known_labels = b.known_labels
    new.probed = copy(b.probed)
    # new.highest_degree_not_probed = b.highest_degree_not_probed
    new.budget = b.budget
    new.p_intra = b.p_intra
    new.p_inter = b.p_inter
    new.n_comms = b.n_comms
end

function add_edge!(g::RevealedGraph, i::Int, j::Int)
    ed = Edge(i,j)
    if !(ed in g.revealed_edges || reverse(ed) in g.revealed_edges)
        push!(g.revealed_edges, ed)
        push!(g.revealed_adj[ed.first], ed.second)
        push!(g.revealed_adj[ed.second], ed.first)
    end
end

function degree(g::RevealedGraph, i::Int)
    return length(g.revealed_adj[i])
end

function neighbors(g::RevealedGraph, i::Int)
    return g.revealed_adj[i]
end

# right now this generates a whole graph. should it generate something partial?
function Base.rand!(rng::AbstractRNG, s::CommunityNetwork, d::RevealedGraph)
    # fill in edges for probed nodes
    s.graph = Graph(n_nodes(d))
    s.labels = Array(Int, n_nodes(d))
    s.n_comms = d.n_comms
    s.budget = d.budget

    # assign node affiliation
    for n in 1:n_nodes(d)
        if length(d.revealed_adj[n]) > 0
            adj_labeled = filter(x->haskey(d.known_labels,x), d.revealed_adj[n])
            n_labeled = length(adj_labeled)
            if n_labeled > 0
                adj_labels = [d.known_labels[n] for n in adj_labeled]
                # use heuristic bayes rule approximation to guess label
                rel_p = Array(Float64, d.n_comms)
                for c in 1:d.n_comms
                    n_intra = sum(adj_labels.==c)
                    rel_p[c] = d.p_inter^n_intra*d.p_inter^(n_labeled-n_intra)
                end
                # select from the probability vector
                x = rand(rng)*sum(rel_p); i = 1; p_sum = rel_p[1];
                while p_sum < x
                    i += 1; p_sum += rel_p[i];
                end
                s.labels[n] = i
                continue
            end
        end
        # if we have not continued by now, assign randomly
        s.labels[n] = ceil(rand(rng)*d.n_comms)
    end

    # connect based on probabilities
    # XXX speed up by moving if statements outside?
    for i in 1:n_nodes(d)
        if !(i in d.probed)
            for j in i+1:n_nodes(d)
                if !(j in d.probed)
                    if s.labels[i] == s.labels[j]
                        if rand(rng) < d.p_intra
                            add_edge!(s.graph, i, j)
                        end
                    else
                        if rand(rng) < d.p_inter
                            add_edge!(s.graph, i, j)
                        end
                    end
                end
            end
        end
    end

    # make sure to add edges we know about
    for ed in d.revealed_edges
        if !has_edge(s.graph, ed)
            add_edge!(s.graph, ed)
        end
    end

    return s
end

function belief(p::CommunityAffiliationPOMDP, b::RevealedGraph, a::ProbeNode, o::Neighborhood, d=RevealedGraph())

    clone!(d, b)
    push!(d.probed, a.node)
    d.budget = d.budget-1
    for n in o.neighbors
        add_edge!(d, a.node, n)
        push!(d.revealed_nodes, n)
    end

    return d
end

function belief(p::CommunityAffiliationPOMDP, b::RevealedGraph, a::GuessAffiliation, o::Neighborhood, d=RevealedGraph())
    clone!(d,b)
    d.budget -= 1
    return d
end


text:
using Documenter
using PowerSimulations

makedocs(
    sitename = "PowerSimulations.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [PowerSimulations],
    strict = true,
    authors = "Jose Daniel Lara, Clayton Barrows and Dheepak Krishnamurthy",
    pages = Any[ "Introduction" => "index.md",
                 #"Quick Start Guide" => "qs_guide.md",
                 "Operation Model" => "man/op_problem.md",
                 "API" => Any["PowerSimulations" => "api/PowerSimulations.md"]
    ],

)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

deploydocs(
    repo = "github.com/NREL/PowerSimulations.jl.git",
    branch = "gh-pages",
    target = "build",
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    make = nothing,
)


text:
struct PSMSDSTProblem
    numjobs::Integer
    nummachines::Integer
    numfamilies::Integer
    durations::Vector{Integer}
    families::Vector{Integer}
    setuptimes::Matrix{Integer}
    lowerbound::Integer
end

"""
    parsepsmsdst(path)

Parse a PSMSDST problem from a text file into a PSMSDSTProblem struct.

Will error if the file does not exists or if the job indices are presented
in an unexpected format.

# Examples
```julia-repl
julia> parsepsmsdst("PSMSDSTdata/PSMSDST1.txt")
PSMSDSTProblem(20, 2, 3, Integer[18, 7, 12, 11, 7, 18, 19, 19, ...
```
"""
function parsepsmsdst(path::String)::PSMSDSTProblem

    # Read the raw file as text
    text = open(path,"r") do f
        read(f,String)
    end

    # Split each line
    lines = split(text,r"\n|\r\n")

    # Remove empty lines
    lines = filter(x->x!="",lines)

    # Split each line on a whitespace separator
    segmented = map(line->split(strip(line),r"(\s+)"),lines)

    # Convert to numbers
    numerical = map(line->map(segment->parse(Int,segment),line),segmented)

    # Extract data from the first line
    numjobs, nummachines, numfamilies = numerical[1]

    # Extract the data for each job
    jobsinfo = vcat(transpose.(numerical[2:1+numjobs])...)

    if jobsinfo[:,1] != collect(1:numjobs)
        error("Job indices have an unexpected format")
    end

    durations = jobsinfo[:,2]
    families = jobsinfo[:,3]

    # Extract the setup times
    setuptimes = vcat(transpose.(numerical[end-numfamilies+1:end])...)

    inevitablesetups = 0
    if numfamilies > nummachines
        inevitablesetups = numfamilies - nummachines
    end
    lowestsetup = typemax(Int64)
    for i in 1:numfamilies
        for j in 1:numfamilies
            if i != j
                if setuptimes[i,j] < lowestsetup
                    lowestsetup = setuptimes[i,j]
                end
            end
        end
    end

    lowerbound = round(Int64,(sum(durations)/nummachines)
        + inevitablesetups*lowestsetup)

    # Construct the result
    PSMSDSTProblem(
        numjobs,
        nummachines,
        numfamilies,
        durations,
        families,
        setuptimes,
        lowerbound
    )
end

"""
    isfeasible(problem,solution)

Determine if a given solution is feasible for the PSMSDST problem.

The i-th elementh of the solution vector corresponds to the list of jobs
(in order) that are to be processed by machine i

Returns a tuple where the first element determines if the solution is feasible
and, when not feasible, the second element contains the reason.

# Examples
```julia-repl
julia> testproblem = PSMSDSTProblem(5,3,2,[6,9,4,2,7],[2,1,2,2,1,2],[0 4; 3 0])
PSMSDSTProblem(5, 3, 2, Integer[6, 9, 4, 2, 7], Integer[2, 1, 2, 2, 1, 2]...
julia> testsolution = Array{Int}.([[2,5],[3,4,1]])
2-element Vector{Vector{Int64}}:...
julia> feasible = isfeasible(testproblem, testsolution)
(false,"The problem has 3 machine but the solution only lists 2 machines")
```
"""
function isfeasible(problem::PSMSDSTProblem,solution::Vector{Vector{Int}}
    )::Tuple{Bool,Union{Nothing,String}}

    # Verify solution length
    if length(solution) != problem.nummachines
        return false, "The problem has $(problem.nummachines) machines but" *
        " the solution only lists $(length(solution)) machines"
    end

    # Create sets for counting
    programmedjobs = Set{Int}()
    availablejobs = Set{Int}(1:problem.numjobs)

    for machine in solution
        for job in machine
            if job in programmedjobs
                return false, "Job $job is programmed twice."
            end
            if !(job in availablejobs)
                return false, "Job $job is not listed in the problem."
            end
            push!(programmedjobs,job)
            delete!(availablejobs,job)
        end
    end

    # Determine if there are jobs that were not programmed
    if !isempty(availablejobs)
        return false, "jobs $availablejobs were not programmed."
    end

    # Looks like it is feasible!
    return true, nothing
end

"""
    evaluatesolution(problem,solution[,verifyfeasibility=true])

Go throught a proposed solution and calculate the start time of each job,
the finish time of each machine, and the finish time of the last job.

If ``verifyfeasibility`` is set to true, this function will also verify
that the solution is feasible and error if it isn't.

# Examples
```julia-repl
julia> starttimes, finishtimes, lastfinish = evaluatesolution(problem,solution)
Vector{Vector{Int64}}, Vector{Int64}, Int64...
```
"""
function evaluatesolution(problem::PSMSDSTProblem,solution::Vector{Vector{Int}}
    ;verifyfeasibility::Bool=false)::Tuple{Vector{Vector{Int}},Vector{Int},Int,Float64}

    # Verify feasibility when asked
    if verifyfeasibility
        feasible, reason = isfeasible(problem,solution)
        if not feasible
            error(reason)
        end
    end

    # Initialize start times and finish times to 0
    starttimes = [ [ 0 for (_,_) in enumerate(jobs) ]
                        for (_,jobs) in enumerate(solution) ]

    finishtimes = [ 0 for (_,_) in enumerate(solution) ]

    # Compute metrics for each machine
    for (machineid,jobs) in enumerate(solution)
        currtime = 0
        for (scheduledorder,jobid) in enumerate(jobs)
            starttime = currtime
            if scheduledorder != 1
                # Find the setup time considering the previous processed job
                thisfamily = problem.families[jobid]
                prevfamily = problem.families[jobs[scheduledorder-1]]
                setuptime = problem.setuptimes[prevfamily,thisfamily]
                starttime += setuptime
            end
            starttimes[machineid][scheduledorder] = starttime
            finishtime = starttime + problem.durations[jobid]

            # Increment the current time
            currtime = finishtime
        end

        # The finish time is the last current time
        finishtimes[machineid] = currtime
    end

    # The maximum job finish time is the maximum finish time over all machines.
    maxtime = maximum(finishtimes)

    relative = ((maxtime - problem.lowerbound) / problem.lowerbound)*100

    return starttimes, finishtimes, maxtime, relative
end

"""
    presolve(solver)

Solves a small problem using the provided solver. This ensures that the solver
function is pre-compiled and no performance is lost during evaluation.

Note that this method is not perfect, as Julia will only compile the code that
is reached by solving this problem, and not necessarily the entire thing.
"""
function presolve(solver::Function)
    testproblem = PSMSDSTProblem(5,3,2,[6,9,4,2,7],[1,2,2,1,2],[0 4; 3 0],100)
    solver(testproblem)
    return
end

function __test()
    parsed = parsepsmsdst("PSMSDSTdata/PSMSDST1.txt")
    @info "Parse result" parsed

    testproblem = PSMSDSTProblem(5,3,2,[6,9,4,2,7],[1,2,2,1,2],[0 4; 3 0],100)
    @info "Test problem" testproblem

    testsolution = Array{Int}.([[2,5],[3,1,4],[]])
    feasible = isfeasible(testproblem, testsolution)
    @info "Test solution" testsolution feasible

    starttimes, finishtimes, loss = evaluatesolution(testproblem,testsolution)
    @info "Evaluation" starttimes finishtimes loss
end

if abspath(PROGRAM_FILE) == @__FILE__
    using Logging
    __test()
end

max_stars_repo_name: mkualquiera/huepl

max_issues_repo_name: mkualquiera/huepl

max_forks_repo_name: mkualquiera/huepl


text:
module Tableaus

    using ..Config
    using ..CommonFunctions
    using ..BasisFunctions
    using ..Quadratures
    using ..Integrators
    using ..Utils

    export getCoefficientsGLRK,
           getCoefficientsGLRK1, getCoefficientsGLRK2, getCoefficientsGLRK3,
           getCoefficientsGLRK4, getCoefficientsGLRK5, getCoefficientsGLRK6

    include("tableaus/coefficients_glrk.jl")

    export getCoefficientsLobIII2,  getCoefficientsLobIII3,  getCoefficientsLobIII4,
           getCoefficientsLobIIIA2, getCoefficientsLobIIIA3, getCoefficientsLobIIIA4,
           getCoefficientsLobIIIB2, getCoefficientsLobIIIB3, getCoefficientsLobIIIB4,
           getCoefficientsLobIIIC2, getCoefficientsLobIIIC3, getCoefficientsLobIIIC4,
           getCoefficientsLobIIID2, getCoefficientsLobIIID3, getCoefficientsLobIIID4,
           getCoefficientsLobIIIE2, getCoefficientsLobIIIE3, getCoefficientsLobIIIE4,
           getCoefficientsLobIIIF2, getCoefficientsLobIIIF3, getCoefficientsLobIIIF4,
           getCoefficientsLobIIIG2, getCoefficientsLobIIIG3, getCoefficientsLobIIIG4

    export getCoefficientsLobIII,  getCoefficientsLobIIIA, getCoefficientsLobIIIB,
           getCoefficientsLobIIIC, getCoefficientsLobIIID, getCoefficientsLobIIIE

    include("tableaus/coefficients_lob.jl")

    export getCoefficientsRadIIA2, getCoefficientsRadIIA3

    include("tableaus/coefficients_rad.jl")

    export getCoefficientsPGLRK, getTableauPGLRK

    include("tableaus/coefficients_pglrk.jl")

    export getCoefficientsSRK3

    include("tableaus/coefficients_srk.jl")

    include("tableaus/tableaus_erk.jl")

    export getTableauExplicitEuler, getTableauExplicitMidpoint, getTableauHeun,
           getTableauKutta, getTableauERK4, getTableauERK438, getTableauVerner

    include("tableaus/tableaus_dirk.jl")

    export getTableauCrouzeix

    include("tableaus/tableaus_firk.jl")

    export getTableauImplicitEuler, getTableauImplicitMidpoint,
           getTableauGLRK,
           getTableauLobIIIA2, getTableauLobIIIA3, getTableauLobIIIA4,
           getTableauLobIIIB2, getTableauLobIIIB3, getTableauLobIIIB4,
           getTableauLobIIIC2, getTableauLobIIIC3, getTableauLobIIIC4,
           getTableauLobIIID2, getTableauLobIIID3, getTableauLobIIID4,
           getTableauLobIIIE2, getTableauLobIIIE3, getTableauLobIIIE4,
           getTableauLobIIIF2, getTableauLobIIIF3, getTableauLobIIIF4,
           getTableauLobIIIG2, getTableauLobIIIG3, getTableauLobIIIG4,
           getTableauRadIIA2,  getTableauRadIIA3,
           getTableauSRK3

    include("tableaus/tableaus_sirk.jl")
    include("tableaus/tableaus_siprk.jl")
    include("tableaus/tableaus_sisprk.jl")

    export  getTableauStochasticGLRK, getTableauStochasticDIRK
    export  getTableauStochasticStormerVerlet, getTableauStochasticSymplecticEuler
    export  getTableauStochasticLobIIIABD2, getTableauModifiedStochasticStormerVerlet

    include("tableaus/tableaus_serk.jl")

    export getTableauPlaten, getTableauBurrageR2, getTableauBurrageCL
    export getTableauBurrageE1, getTableauBurrageG5, getTableauStochasticHeun
    export getTableauStochasticEuler

    include("tableaus/tableaus_werk.jl")

    export getTableauRosslerRS1, getTableauRosslerRS2

    include("tableaus/tableaus_wirk.jl")

    export getTableauSRKw1, getTableauSRKw2

    include("tableaus/tableaus_eprk.jl")

    export getTableauSymplecticEulerA, getTableauSymplecticEulerB,
           getTableauLobattoIIIAIIIB2, getTableauLobattoIIIBIIIA2

    export getTableauSPARKGLRK

    include("tableaus/tableaus_spark.jl")

    export getTableauVPGLRK,
           getTableauVPLobIIIA2, getTableauVPLobIIIA3, getTableauVPLobIIIA4,
           getTableauVPLobIIIB2, getTableauVPLobIIIB3, getTableauVPLobIIIB4,
           getTableauVPLobIIIC2, getTableauVPLobIIIC3, getTableauVPLobIIIC4,
           getTableauVPLobIIID2, getTableauVPLobIIID3, getTableauVPLobIIID4,
           getTableauVPLobIIIE2, getTableauVPLobIIIE3, getTableauVPLobIIIE4,
           getTableauVPLobIIIF2, getTableauVPLobIIIF3, getTableauVPLobIIIF4,
           getTableauVPLobIIIG2, getTableauVPLobIIIG3, getTableauVPLobIIIG4,
           getTableauVPSRK3,
           getTableauVPLobIIIAIIIA2, getTableauVPLobIIIAIIIA3, getTableauVPLobIIIAIIIA4,
           getTableauVPRadIIAIIA2, getTableauVPRadIIAIIA3

    include("tableaus/tableaus_vprk.jl")

    export getTableauSymplecticProjection,
           getTableauLobIIIAIIIB2pSymplectic,
           getTableauLobIIIAIIIB3pSymplectic,
           getTableauLobIIIAIIIB4pSymplectic,
           getTableauGLRKpSymplectic

    include("tableaus/tableaus_vpark.jl")

    export getTableauVSPARKMidpointProjection,
           getTableauVSPARKSymmetricProjection,
           getTableauVSPARKSymplecticProjection,
           getTableauVSPARKGLRKpMidpoint,
           getTableauVSPARKGLRKpSymmetric,
           getTableauVSPARKGLRKpSymplectic,
           getTableauVSPARKLobIIIAIIIB2pSymmetric,
           getTableauVSPARKLobIIIAIIIB3pSymmetric,
           getTableauVSPARKLobIIIAIIIB4pSymmetric

    include("tableaus/tableaus_vspark.jl")

    export getTableauVSPARKLobIIIAB,
           getTableauVSPARKLobIIIC,
           getTableauVSPARKLobIIID,
           getTableauVSPARKLobIIIE,
           getTableauVSPARKGLRKLobIIIAB,
           getTableauVSPARKGLRKLobIIIC,
           getTableauVSPARKGLRKLobIIID,
           getTableauVSPARKGLRKLobIIIE

    include("tableaus/tableaus_vspark_psecondary.jl")

    export getTableauHPARK, getTableauHPARKGLRK,
           getTableauHPARKLobIIIAIIIB2,
           getTableauHPARKLobIIIAIIIB3,
           getTableauHPARKLobIIIAIIIB4

    include("tableaus/tableaus_hpark.jl")

    export getTableauHSPARKSymmetricProjection,
           getTableauHSPARKGLRKpSymmetric,
           getTableauHSPARKLobIIIAIIIB2pSymmetric,
           getTableauHSPARKLobIIIAIIIB3pSymmetric,
           getTableauHSPARKLobIIIAIIIB4pSymmetric

    include("tableaus/tableaus_hspark.jl")

    export getTableauSLRKLobIIIAB,
           getTableauSLRKLobIIIC,
           getTableauSLRKLobIIID,
           getTableauSLRKLobIIIE

    include("tableaus/tableaus_slrk.jl")

    export getTableauLieA, getTableauLieB, getTableauStrang,
           getTableauMcLachlan2, getTableauMcLachlan4,
           getTableauTripleJump, getTableauSuzukiFractal

    include("tableaus/tableaus_splitting.jl")

    function __init__()
        add_config(:tab_compensated_summation, true)
    end

end


text:
# Functions related to creating homogenised equivalent medium tensors

"""
    pitl(d1, vp1, vs1, rho1, d2, vp2, vs2, rho2) -> C, rho

Return the effective 6x6 Voigt matrix `C` and density `rho` for a medium defined
by two periodic layers where each layer `i` is defined by:

`di`:   The proportion of the total medium

`vpi`:  P-wave velocity (m/s)

`vsi`:  S-wave velocity (m/s)

`rhoi`: Density (kg/m^3)

`C` is density-normalised.
"""
function pitl(d1, vp1, vs1, rho1, d2, vp2, vs2, rho2)
    C = zero(EC)
    # Lamé parameters from velocities
    m1 = rho1*vs1^2
    m2 = rho2*vs2^2
    l1 = rho1*vp1^2 - 2m1
    l2 = rho2*vp2^2 - 2m2
    # Time-saving terms
    l1p2m1 = l1 + 2*m1
    l2p2m2 = l2 + 2*m2
    # D term, p. 785
    D = (d1 + d2)*(d1*l2p2m2 + d2*l1p2m1)
    # Eq. (7)
    @inbounds begin
        C[1,1] = ((d1+d2)^2*l1p2m1*l2p2m2 + 4*d1*d2*(m1 - m2)*((l1 + m1) - (l2 + m2)))/D
        C[2,2] = C[1,1]
        C[1,2] = C[2,1] = ((d1 + d2)^2*l1*l2 + 2*(l1*d1 + l2*d2)*(m2*d1 + m1*d2))/D
        C[1,3] = ((d1 + d2)*(l1*d1*l2p2m2 + l2*d2*l1p2m1))/D
        C[3,1] = C[2,3] = C[1,3]
        C[3,2] = C[2,3]
        C[3,3] = ((d1 + d2)^2*l1p2m1*l2p2m2)/D
        C[4,4] = C[5,5] = (d1 + d2)*m1*m2/(d1*m2 + d2*m1)
        C[6,6] = (m1*d1 + m2*d2)/(d1 + d2)
    end
    # Normalise back to average density
    rho = (d1*rho1 + d2*rho2)/(d1 + d2)
    C ./= rho
    return C, rho
end

"""
    tandon_and_weng(vp, vs, rho, del, c, vpi, vsi, rhoi) -> C, rho

Return the effective elastic constants `C` and density `rho` of a medium with matrix
velocities `vp` and `vs` (m/s) and density `rho` (kg/m^3), and inclusion
velocities `vpi` and `vsi`, density `rhoi`.  The symmetry axis is parallel to x1.

`del` is the aspect ration of spheroidal inclusions: <1 oblate, >1 prolate

`c` is the volume fraction of inclusions (0 <= c <= 1).
"""
function tandon_and_weng(vp, vs, rho, del, c, vpi, vsi, rhoi)
    #This implementation is based on the Fortran code by Mike Kendall
    del > 0 || error("CIJ.tandon_and_weng: `del` must be > 0.")
    del == 1 && error("CIJ.tandon_and_weng: Theory not valid for `del == 1`")
    0 <= c <= 1 || error("CIJ.tandon_and_weng: `c` must be in range 0 - 1")
    vp == vpi && vs == vsi && rho == rhoi && return CIJ.iso(vp=vp, vs=vs), rho

    C = zero(EC)
    #  weighted average density
    rho_out = (1 - c)*rho + c*rhoi

    amu  = vs^2*rho
    amui = vsi^2*rhoi
    alam = vp^2*rho - 2*amu
    alami = vpi^2*rhoi - 2*amui
    bmi = alami + amui*2/3
    bmps = alam + amu
    #  Young's modulus for matrix
    E0 = amu*(3*alam + 2*amu)/(alam + amu)
    #  Poisson's ratio of the matrix.
    anu = alam/(2*(alam + amu))

    #  Some time saving terms
    t1 = del^2 - 1
    t2 = 1 - anu
    t3 = 1 - 2*anu
    t4 = 3 * del^2
    t5 = 1 - del^2

    # D1, D2 and D3 from Tandon and Weng (1984) (just before equation (18)).
    D1 = 1 + 2*(amui - amu)/(alami - alam)
    D2 = (alam + 2*amu)/(alami - alam)
    D3 = alam/(alami - alam)

    # g and g' terms (appendix of Tandon and Weng 1984). g is for prolate spheroidal
    # inclusions (del>1), whilst g' is for disc-like (oblate) inclusions (del<1).
    if (del >= 1)
        acshdel = log(del + sqrt(t1))
        g = (del*sqrt(t1) - acshdel)*del/sqrt(t1^3)
    else
        # g' below
        g = (acos(del) - del*sqrt(t5))*del/sqrt(t5^3) ;
    end

    # Eshelby's Sijkl tensor (appendix of Tandon and Weng 1984).
    s11 = (t3 + (t4-1)/t1 - (t3 + t4/t1)*g)/(2*t2)
    s22 = (t4/(t1*2) + (t3 - 9/(4*t1))*g)/(4*t2)
    s33 = s22
    s23 = (del^2/(2*t1) - (t3 + 3/(4*t1))*g)/(4*t2)
    s32 = s23
    s21 = (-2*del*del/t1 + (t4/t1 - t3)*g)/(4*t2)
    s31 = s21
    s12 = (-(t3 + 1/t1) + (t3 + 3/(2*t1))*g)/(2*t2)
    s13 = s12
    s44 = (del*del/(2*t1) + (t3 - 3/(4*t1))*g)/(4*t2)
    s66 = (t3 - (t1+2)/t1 - (t3 - 3*(t1+2)/t1)*g/2)/(4*t2)
    s55 = s66

    # Tandon and Weng's B terms (after equation 17).
    B1 = c*D1 + D2 + (1-c)*(D1*s11 + 2*s21)
    B2 = c + D3 + (1-c)*(D1*s12 + s22 + s23)
    B3 = c + D3 + (1-c)*(s11 + (1+D1)*s21)
    B4 = c*D1 + D2 + (1-c)*(s12 + D1*s22 + s23)
    B5 = c + D3 + (1-c)*(s12 + s22 + D1*s23)

    # Tandon and Weng's A terms (after equation 20).
    A1 = D1*(B4 + B5) - 2*B2
    A2 = (1 + D1)*B2 - (B4 + B5)
    A3 = B1 - D1*B3
    A4 = (1 + D1)*B1 - 2*B3
    A5 = (1 - D1)/(B4 - B5)
    A = 2*B2*B3 - B1*(B4+B5)

    # Tandon and Weng (1984) equations (25) (28) (31) (32)
    E11 = E0/(1+c*(A1+2*anu*A2)/A)
    E22 = E0/(1+c*(-2*anu*A3 + (1-anu)*A4 + (1+anu)*A5*A)/(2*A))
    amu12 = amu*(1 + c/(amu/(amui-amu) + 2*(1-c)*s66))
    amu23 = amu*(1 + c/(amu/(amui-amu) + 2*(1-c)*s44))

    # Sayers equation (36)
    anu31 = anu - c*(anu*(A1+2*anu*A2)+(A3-anu*A4)) / (A + c*(A1+2*anu*A2))

    # T&W equation (36)
    #     aK12 term; bmps=plane strain bulk modulus
    anum = (1 + anu)*(1 - 2*anu)
    denom = 1 - anu*(1+2*anu31) + c*(2*(anu31-anu)*A3 + (1-anu*(1+2*anu31))*A4)/A
    aK23 = bmps*anum/denom
    anu12tst = E11/E22 - (1/amu23 + 1/aK23)*E11/4

    # Cij - Sayers' (1992) equations (24)-(29).
    # Conversion
    C[2,2] = amu23 + aK23
    C[3,3] = C[2,2]
    C[1,1] = E11 + 4*anu12tst*aK23
    C[2,3] = -amu23 + aK23
    C[1,2] = 2*anu31*aK23
    C[1,3] = C[1,2]
    C[5,5] = amu12
    C[6,6] = C[5,5]
    C[4,4] = (C[2,2] - C[2,3])/2

    # Fill out matrix by symmetry
    CIJ.symm!(C)

    # apply density normalisation
    C ./= rho_out

    C, rho_out
end

"""
    hudson(vp, vs, rho, del, ϵ, vpi, vsi, rhoi) -> C, rho_bulk

Return the effective elastic constants `C` and density `rho_bulk` of a medium with
matrix isotropic velocities `vp` and `vs`, and density `rho`, and inclusion
velocities `vpi` and `vsi`, and density `rhoi`, where the crack density is `ϵ`
and the aspect ratio of ellipsoidal inclusions is `del`.

The theory is valid when `ϵ` ≪ 1, where `ϵ` is given by ``N a^2/\nu`` and ``N`` is
the number of cracks of radius ``a`` in a volume ``\nu``.

The formulation is according to Hudson (1982), as given in Crampin (1984).
"""
function hudson(vp, vs, rho, del, ϵ, vpi, vsi, rhoi)
    # error("`hudson` has not been tested yet")
    ϵ > 0.1 && warn("Theory of Hudson (1982) only valid for `ϵ` < 0.1, but using ϵ = $ϵ")
    λ, μ, κ = lame(vp, vs, rho)
    λ′, μ′, κ′ = lame(vpi, vsi, rhoi)
    K = (κ′ + 4/3*μ′)/(π*del*μ) * (λ + 2μ)/(λ + μ)
    M = 4μ′/(π*del*μ) * (λ + 2μ)/(3λ + 4μ)
    U₁ = 4/3*(λ + 2μ)/(λ + μ)/(1 + K)
    U₃ = 16/3*(λ + 2μ)/(3λ + 4μ)/(1 + M)
    q = 15*(λ/μ)^2 + 28*(λ/μ) + 28
    X = 2μ*(3λ + 8μ)/(λ + 2μ)
    c⁰ = iso(lam=λ, mu=μ)
    # First-order correction from Hudson (1981)
    c¹ = -ϵ./μ .* @SMatrix [ (λ + 2μ)^2*U₁  λ*(λ + 2μ)*U₁  λ*(λ + 2μ)*U₁ 0    0      0
                             λ*(λ + 2μ)*U₁     λ^2*U₁         λ^2*U₁     0    0      0
                             λ*(λ + 2μ)*U₁     λ^2*U₁         λ^2*U₁     0    0      0
                                   0             0              0        0    0      0
                                   0             0              0        0  μ^2*U₃   0
                                   0             0              0        0    0    μ^2*U₃]
    # Second-order correction from Hudson (1982)
    c² = ϵ.^2 ./ 15 .* @SMatrix [ (λ + 2μ)*q*U₁^2       λ*q*U₁^2             λ*q*U₁^2       0   0      0
                                     λ*q*U₁^2     λ^2*q/(λ + 2μ)*U₁^2  λ^2*q/(λ + 2μ)*U₁^2  0   0      0
                                     λ*q*U₁^2     λ^2*q/(λ + 2μ)*U₁^2  λ^2*q/(λ + 2μ)*U₁^2  0   0      0
                                        0                  0                    0           0   0      0
                                        0                  0                    0           0 X*U₃^2   0
                                        0                  0                    0           0   0    X*U₃^2]
    rho_bulk = (1 - ϵ)*rho + ϵ*rhoi
    c = EC{DEFAULT_FLOAT}((c⁰ .+ c¹ .+ c²)./rho_bulk)
    c, rho_bulk
end

"""
    grechka_cracks!(C, ξ, ϕ=0) -> C

Add dry cracks using the theory of Grechka (2007) [1] to a tensor `C` which has VTI symmetry
about the 3-axis.  Cracks have fracture density `ξ` and strike `ϕ`°, measured from the
1-axis towards the negative 2-axis.  `C` is updated in-place.

    grechka_cracks(C, ξ, ϕ=0) -> C′

Copying version of `grechka_cracks!`.

Both forms take `ξ` and `ϕ` as both scalars and `AbstractArray`s.  In the latter case,
multiple fracture sets are added, and are assumed to be non-interacting.  This is valid for
'small' ξ only.

N.B.: No check is made on the input constants `C`.

#### References

1. Vladimir Grechka (2007). Multiple cracks in VTI rocks: Effective properties and fracture
   characterization.  Geophysics, 72, D81–D91.  doi:10.1190/1.2751500
"""
function grechka_cracks!(C, ξ, ϕ=zero(eltype(C)))
    C .= CIJ.rot3(C, 0, 0, -ϕ)
    # Excess normal and tangential compliances of crack
    Bn = 4/3*ξ*C[2,2]/C[6,6]/(C[2,2] - C[6,6])
    Bth = 16/3*ξ*C[2,2]/C[6,6]/(3C[2,2] - 2C[6,6])
    Btv = 16/3*ξ*C[3,3]/C[5,5]/(3C[3,3] - 2C[5,5])
    # Update compliance after Eschelby (1957)
    S = C2S!(C)
    S[2,2] += Bn
    S[4,4] += Btv
    S[6,6] += Bth
    C .= CIJ.rot3(S2C!(S), 0, 0, ϕ)
end
grechka_cracks(C, ξ, ϕ=zero(eltype(C))) = grechka_cracks!(deepcopy(C), ξ, ϕ)

function grechka_cracks!(C, ξ::AbstractArray, ϕ::AbstractArray=zeros(ξ, eltype(C)))
    length(ξ) == length(ϕ) || throw(ArgumentError("Lengths of ξ and ϕ must be the same"))
    for i in 1:length(ξ)
        grechka_cracks!(C, ξ[i], ϕ[i])
    end
    C
end
grechka_cracks(C, ξ::AbstractVector, ϕ::AbstractVector=zeros(ξ, eltype(C))) =
    grechka_cracks!(deepcopy(C), ξ, ϕ)

@doc (@doc grechka_cracks!) grechka_cracks


text:
# util.jl

"""
Return a vector with reversed indices, if `p = reverse_indices(a)` then
`p[i] = a[-i]`. The result is an OffsetArray.
"""
function reverse_indices(a::AbstractVector)
    firstindex = first(linearindices(a))
    lastindex = last(linearindices(a))
    OffsetArray([a[i] for i in lastindex:-1:firstindex], -lastindex:-firstindex)
end


text:
module GuillotineCallback

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export guillotine_callback, get_involved_boxes_flbx, find_all_guillotine_cuts_exceeding_limit

#Finds a guillotine cut that is larger than the allowed amount
function guillotine_callback(cb_data, m)

    max_cut = guillotine_ds_var * callback_value(cb_data, m[:maxHeight])

    fixed_x = []
    fixed_y = []

    #Check if there exists a cut higher than the allowed cuts
    for i in 1:length(boxes) - 1

        if length(fixed_x) > 0 || length(fixed_y) > 0
            print("Fixed X ", fixed_x, "\n")
            print("Fixed Y ", fixed_y, "\n")
        end

        print("Round ", i, "\n")

        #Find out where the cut is and pass information onwards
        if (callback_value(cb_data, m[:cut_flbx_left_i][i]) + 0.001 >= max_cut && check_if_fixed_x(cb_data, m, 1, i, fixed_x))
            print("Nr 1 \n")
            get_involved_boxes(cb_data, m, 1, i)
            x_value = callback_value(cb_data, m[:flbx_front_left_bot_x][i])
            push!(fixed_x, x_value)
        end
        if (callback_value(cb_data, m[:cut_flbx_j_plus_left_i][i]) + 0.001 >= max_cut && check_if_fixed_x(cb_data, m, 2, i, fixed_x))
            print("Nr 2 \n")
            get_involved_boxes(cb_data, m, 2, i)
            x_value = callback_value(cb_data, m[:flbx_front_left_bot_x][i])
            push!(fixed_x, x_value)
        end
        if (callback_value(cb_data, m[:cut_flbx_i_plus_right][i]) + 0.001 >= max_cut && check_if_fixed_x(cb_data, m, 3, i, fixed_x))
            print("Nr 3 \n")
            get_involved_boxes(cb_data, m, 3, i)
            x_value = callback_value(cb_data, m[:flbx_front_left_bot_x][i]) + callback_value(cb_data, m[:box_x_axis_cover][i])
            push!(fixed_x, x_value)
        end
        if (callback_value(cb_data, m[:cut_flbx_ij_plus_right][i]) + 0.001 >= max_cut && check_if_fixed_x(cb_data, m, 4, i, fixed_x))
            print("Nr 4 \n")
            get_involved_boxes(cb_data, m, 4, i)
            x_value = callback_value(cb_data, m[:flbx_front_left_bot_x][i]) + callback_value(cb_data, m[:box_x_axis_cover][i])
            push!(fixed_x, x_value)
        end
        if (callback_value(cb_data, m[:cut_flby_right][i]) + 0.001 >= max_cut && check_if_fixed_y(cb_data, m, 5, i, fixed_y))
            print("Nr 5 \n")
            get_involved_boxes(cb_data, m, 5, i)
            y_value = callback_value(cb_data, m[:flby_front_left_bot_y][i])
            push!(fixed_y, y_value)
        end
        if (callback_value(cb_data, m[:cut_flby_j_plus_right][i]) + 0.001 >= max_cut && check_if_fixed_y(cb_data, m, 6, i, fixed_y))
            print("Nr 6 \n")
            get_involved_boxes(cb_data, m, 6, i)
            y_value = callback_value(cb_data, m[:flby_front_left_bot_y][i])
            push!(fixed_y, y_value)
        end
        if (callback_value(cb_data, m[:cut_flby_i_plus_left][i]) + 0.001 >= max_cut && check_if_fixed_y(cb_data, m, 7, i, fixed_y))
            print("Nr 7 \n")
            get_involved_boxes(cb_data, m, 7, i)
            y_value = callback_value(cb_data, m[:flby_front_left_bot_y][i]) + callback_value(cb_data, m[:box_y_axis_cover][i])
            push!(fixed_y, y_value)
        end
        if (callback_value(cb_data, m[:cut_flby_ij_plus_left][i]) + 0.001 >= max_cut && check_if_fixed_y(cb_data, m, 8, i, fixed_y))
            print("Nr 8 \n")
            get_involved_boxes(cb_data, m, 8, i)
            y_value = callback_value(cb_data, m[:flby_front_left_bot_y][i]) + callback_value(cb_data, m[:box_y_axis_cover][i])
            push!(fixed_y, y_value)
        end
    end

    return
end

function check_if_fixed_x(cb_data, m, pos, i, fixed_x)
    fixed = false
    if pos == 1 || pos == 2
        for fixed_pos in 1:length(fixed_x)
            if (callback_value(cb_data, m[:flbx_front_left_bot_x][i]) - 0.0001 <=
                fixed_x[fixed_pos] <=
                callback_value(cb_data, m[:flbx_front_left_bot_x][i]) + 0.0001)
                fixed = true
            end
        end
    elseif pos == 3 || pos == 4
        for fixed_pos in 1:length(fixed_x)
            if (callback_value(cb_data, m[:flbx_front_left_bot_x][i]) + callback_value(cb_data, m[:box_x_axis_cover][i]) - 0.0001 <=
                fixed_x[fixed_pos] <=
                callback_value(cb_data, m[:flbx_front_left_bot_x][i]) + callback_value(cb_data, m[:box_x_axis_cover][i]) + 0.0001)
                fixed = true
            end
        end
    end
    return ~fixed
end

function check_if_fixed_y(cb_data, m, pos, i, fixed_y)
    fixed = false
    if pos == 5 || pos == 6
        for fixed_pos in 1:length(fixed_y)
            if (callback_value(cb_data, m[:flby_front_left_bot_y][i]) - 0.0001 <=
                fixed_y[fixed_pos] <=
                callback_value(cb_data, m[:flby_front_left_bot_y][i]) + 0.0001)
                fixed = true
            end
        end
    elseif pos == 7 || pos == 8
        for fixed_pos in 1:length(fixed_y)
            if (callback_value(cb_data, m[:flby_front_left_bot_y][i]) + callback_value(cb_data, m[:box_y_axis_cover][i]) - 0.0001 <=
                fixed_y[fixed_pos] <=
                callback_value(cb_data, m[:flby_front_left_bot_y][i]) + callback_value(cb_data, m[:box_y_axis_cover][i]) + 0.0001)
                fixed = true
            end
        end
    end
    return ~fixed
end

#Takes a box number to check what boxes are involved in guillotine
function get_involved_boxes(cb_data, m, pos, i)
    involved = []
    push!(involved, i)
    # Check the position of the cut
    for j in 1:length(boxes) - 1
        if i!=j
            if (pos == 1)
                if (callback_value(cb_data, m[:guillotine_flbx][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 2)
                if (callback_value(cb_data, m[:guillotine_flbx_j_plus][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 3)
                if (callback_value(cb_data, m[:guillotine_flbx_i_plus][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 4)
                if (callback_value(cb_data, m[:guillotine_flbx_ij_plus][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 5)
                if (callback_value(cb_data, m[:guillotine_flby][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 6)
                if (callback_value(cb_data, m[:guillotine_flby_j_plus][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 7)
                if (callback_value(cb_data, m[:guillotine_flby_i_plus][i, j]) == 1)
                    push!(involved, j)
                end
            elseif (pos == 8)
                if (callback_value(cb_data, m[:guillotine_flby_ij_plus][i, j]) == 1)
                    push!(involved, j)
                end
            end
            break
        end
    end

    if (length(involved) >= 2)
        create_constraint(cb_data, m, involved, pos)
    end
end

function create_constraint(cb_data, m, involved, pos)
    print(involved)

    max_height = 0
    max_height_box = 0
    snd_max_height = 0
    snd_max_height_box = 0

    for i in 1:length(involved)
        height_of_i = callback_value(cb_data, m[:flbz_front_left_bot_z][involved[i]]) + callback_value(cb_data, m[:hob_height_of_box][involved[i]])
        if height_of_i >= max_height
            #update second highest
            snd_max_height = max_height
            snd_max_height_box = max_height_box

            #update highest
            max_height = height_of_i
            max_height_box = involved[i]
        elseif height_of_i >= snd_max_height
            #update second highest
            snd_max_height = height_of_i
            snd_max_height_box = involved[i]
        end
    end

    for i in 1:3
        print('\n')
        print(i)
        print('\n')
        @show callback_value(cb_data, m[:flbx_front_left_bot_x][i])
        @show callback_value(cb_data, m[:flby_front_left_bot_y][i])
        @show callback_value(cb_data, m[:flbz_front_left_bot_z][i])
    end

    print('\n')
    print("First ", max_height, " ", max_height_box, "\n")
    print("Second ", snd_max_height, " ", snd_max_height_box, "\n")

    #If same size as sentinel, then the new constraints can't apply.
    if ~((callback_value(cb_data, m[:lob_length_of_box][max_height_box]) >=
        (1 - guillotine_fix_var) * callback_value(cb_data, m[:lob_length_of_box][length(boxes)])) &&
        (callback_value(cb_data, m[:wob_width_of_box][max_height_box]) >=
        (1 - guillotine_fix_var) * callback_value(cb_data, m[:wob_width_of_box][length(boxes)])))

        #check if lowest dimension still same as width of sentinel
        if ~((callback_value(cb_data, m[:lob_length_of_box][max_height_box]) >=
            (1 - guillotine_fix_var) * callback_value(cb_data, m[:lob_length_of_box][length(boxes)])) &&
            (callback_value(cb_data, m[:wob_width_of_box][max_height_box]) >=
            (1 - guillotine_fix_var) * callback_value(cb_data, m[:lob_length_of_box][length(boxes)])))
            if pos == 1
                con = @build_constraint(
                    m[:flbx_front_left_bot_x][max_height_box] -
                    m[:flbx_front_left_bot_x][snd_max_height_box] >= guillotine_fix_var
                    )
                print("added constraint 1 \n")
            elseif pos == 2
                con = @build_constraint(
                    (m[:flbx_front_left_bot_x][max_height_box] + m[:box_x_axis_cover][max_height_box]) -
                    m[:flbx_front_left_bot_x][snd_max_height_box] >= guillotine_fix_var
                    )
                print("added constraint 2 \n")
            elseif pos == 3
                con = @build_constraint(
                    (m[:flbx_front_left_bot_x][snd_max_height_box] + m[:box_x_axis_cover][snd_max_height_box]) -
                    m[:flbx_front_left_bot_x][max_height_box] >= guillotine_fix_var
                    )
                print("added constraint 3 \n")
            elseif pos == 4
                con = @build_constraint(
                    (m[:flbx_front_left_bot_x][snd_max_height_box] + m[:box_x_axis_cover][snd_max_height_box]) -
                    (m[:flbx_front_left_bot_x][max_height_box] + m[:box_x_axis_cover][max_height_box]) >= guillotine_fix_var
                    )
                print("added constraint 4 \n")
            end
            if pos == 1 || pos == 2 || pos == 3 || pos == 4
                MOI.submit(m, MOI.LazyConstraint(cb_data), con)
            end
        end
        if ~((callback_value(cb_data, m[:lob_length_of_box][max_height_box]) >=
            (1 - guillotine_fix_var) * callback_value(cb_data, m[:wob_width_of_box][length(boxes)])) &&
            (callback_value(cb_data, m[:wob_width_of_box][max_height_box]) >=
            (1 - guillotine_fix_var) * callback_value(cb_data, m[:wob_width_of_box][length(boxes)])))
            if pos == 5
                con = @build_constraint(
                    m[:flby_front_left_bot_y][max_height_box] -
                    m[:flby_front_left_bot_y][snd_max_height_box] >= guillotine_fix_var
                    )
                print("added constraint 5 \n")
            elseif pos == 6
                con = @build_constraint(
                    (m[:flby_front_left_bot_y][max_height_box] + m[:box_y_axis_cover][max_height_box]) -
                    m[:flby_front_left_bot_y][snd_max_height_box] >= guillotine_fix_var
                    )
                print("added constraint 6 \n")
            elseif pos == 7
                con = @build_constraint(
                    (m[:flby_front_left_bot_y][snd_max_height_box] + m[:box_y_axis_cover][snd_max_height_box]) -
                    m[:flby_front_left_bot_y][max_height_box] >= guillotine_fix_var
                    )
                print("added constraint 7 \n")
            elseif pos == 8
                con = @build_constraint(
                    (m[:flby_front_left_bot_y][snd_max_height_box] + m[:box_y_axis_cover][snd_max_height_box]) -
                    (m[:flby_front_left_bot_y][max_height_box] + m[:box_y_axis_cover][max_height_box]) >= guillotine_fix_var
                    )
                print("added constraint 8 \n")
            end
            if pos == 5 || pos == 6 || pos == 7 || pos == 8
                MOI.submit(m, MOI.LazyConstraint(cb_data), con)
            end
        end
    end

end

end


text:
# This file is part of Kpax3. License is MIT.

function traceR(fileroot::String;
                maxlag::Int=50)
  fpR = open(string(fileroot, "_row_partition.bin"), "r")

  tmp = zeros(Int, 1)

  read!(fpR, tmp)
  n = tmp[1]

  read!(fpR, tmp)
  read!(fpR, tmp)
  N = tmp[1]

  if N - maxlag < 20
    close(fpR)
    throw(KDomainError(string("Not enough samples to estimate ",
                              "autocorrelation. Increase chain length or ",
                              "reduce value of 'maxlag'.")))
  end

  # at any given point we store maxlag + 1 partitions
  kset = zeros(Int, maxlag + 1)
  pset = zeros(Int, n, maxlag + 1)

  # read partitions into these variables
  k = zeros(Int, 1)
  R = zeros(Int, n)
  p = zeros(Int, n)
  v = zeros(Int, n)

  # output
  nl = zeros(Float64, maxlag)
  avgd = zeros(Float64, maxlag)
  entropy = fill(log(n), N)

  T = 1
  t = 1
  s = 1
  l = 1
  z = 0.0

  # the first maxlag + 1 behaves differently
  while !eof(fpR) && T <= maxlag
    read!(fpR, k)
    read!(fpR, R)

    copyto!(p, normalizepartition(R, n))

    fill!(v, 0.0)
    for i in 1:n
      v[p[i]] += 1
    end

    # compute (Shannon) partition entropy: -sum((v / n) * log(v / n))
    z = 0.0
    for g in 1:k[1]
      z += v[g] * log(v[g])
    end
    entropy[T] -= z / n

    kset[t] = k[1]
    copyto!(pset, 1 + n * (t - 1), p, 1, n)

    # compute all the lags starting from this point going backward
    l = 1
    s = t - l
    while s > 0
      nl[l] += 1
      avgd[l] += (jaccard(pset[:, t], kset[t],
                          pset[:, s], kset[s], n) - avgd[l]) / nl[l]
      l += 1
      s -= 1
    end

    t = (t <= maxlag) ? t + 1 : 1
    T += 1
  end

  while !eof(fpR)
    read!(fpR, k)
    read!(fpR, R)

    copyto!(p, normalizepartition(R, n))

    fill!(v, 0.0)
    for i in 1:n
      v[p[i]] += 1
    end

    # compute (Shannon) partition entropy: -sum((v / n) * log(v / n))
    z = 0.0
    for g in 1:k[1]
      z += v[g] * log(v[g])
    end
    entropy[T] -= z / n

    kset[t] = k[1]
    copyto!(pset, 1 + n * (t - 1), p, 1, n)

    # compute all the lags starting from this point going backward
    l = 1
    s = t - l
    while s > 0
      nl[l] += 1
      avgd[l] += (jaccard(pset[:, t], kset[t],
                          pset[:, s], kset[s], n) - avgd[l]) / nl[l]
      l += 1
      s -= 1
    end

    s = maxlag + 1
    while s > t
      nl[l] += 1
      avgd[l] += (jaccard(pset[:, t], kset[t],
                          pset[:, s], kset[s], n) - avgd[l]) / nl[l]
      l += 1
      s -= 1
    end

    t = (t <= maxlag) ? t + 1 : 1
    T += 1
  end

  close(fpR)

  if T != N + 1
    warn(string("Expecting ", N, " simulations but instead found ", T-1, "."))
  end

  (entropy, avgd)
end

function traceC(fileroot::String;
                maxlag::Int=50)
  fpC = open(string(fileroot, "_col_partition.bin"), "r")

  tmp = zeros(Int, 1)

  read!(fpC, tmp)
  n = tmp[1]

  read!(fpC, tmp)
  m = tmp[1]

  read!(fpC, tmp)
  N = tmp[1]

  if N - maxlag < 20
    close(fpC)
    throw(KDomainError(string("Not enough samples to estimate ",
                              "autocorrelation. Increase chain length or ",
                              "reduce value of 'maxlag'.")))
  end

  # at any given point we store maxlag + 1 variables
  cset = zeros(UInt8, m, maxlag + 1)

  # read column classification into these variable
  C = zeros(UInt8, m)

  # output
  avgd = zeros(Float64, maxlag)
  entropy = fill(log(m), N)

  # temporary variables
  ct = zeros(Float64, 3)
  nl = zeros(Float64, maxlag)

  T = 1
  t = 1
  s = 1
  l = 1

  # the first maxlag + 1 behaves differently
  while !eof(fpC) && T <= maxlag
    readbytes!(fpC, C, m)

    # compute (Shannon) entropy of C: -sum((ct / m) * log(ct / m))
    fill!(ct, 0.0)
    for b in 1:m
      ct[C[b]] += 1
    end

    if ct[1] > 0.0
      ct[1] *= log(ct[1])
    end

    if ct[2] > 0.0
      ct[2] *= log(ct[2])
    end

    if ct[3] > 0.0
      ct[3] *= log(ct[3])
    end

    entropy[T] -= (ct[1] + ct[2] + ct[3]) / m

    copyto!(cset, 1 + m * (t - 1), C, 1, m)

    # compute all the lags starting from this point going backward
    l = 1
    s = t - l
    while s > 0
      nl[l] += 1
      avgd[l] += (Distances.hamming(cset[:, t], cset[:, s]) / m -
                  avgd[l]) / nl[l]
      l += 1
      s -= 1
    end

    t = (t <= maxlag) ? t + 1 : 1
    T += 1
  end

  while !eof(fpC)
    readbytes!(fpC, C, m)

    # compute entropy of C: -sum((ct / m) * log(ct / m))
    fill!(ct, 0.0)
    for b in 1:m
      ct[C[b]] += 1
    end

    if ct[1] > 0.0
      ct[1] *= log(ct[1])
    end

    if ct[2] > 0.0
      ct[2] *= log(ct[2])
    end

    if ct[3] > 0.0
      ct[3] *= log(ct[3])
    end

    entropy[T] -= (ct[1] + ct[2] + ct[3]) / m

    copyto!(cset, 1 + m * (t - 1), C, 1, m)

    # compute all the lags starting from this point going backward
    l = 1
    s = t - l
    while s > 0
      nl[l] += 1
      avgd[l] += (Distances.hamming(cset[:, t], cset[:, s]) / m -
                  avgd[l]) / nl[l]
      l += 1
      s -= 1
    end

    s = maxlag + 1
    while s > t
      nl[l] += 1
      avgd[l] += (Distances.hamming(cset[:, t], cset[:, s]) / m -
                  avgd[l]) / nl[l]
      l += 1
      s -= 1
    end

    t = (t <= maxlag) ? t + 1 : 1
    T += 1
  end

  close(fpC)

  if T != N + 1
    warn(string("Expecting ", N, " simulations but instead found ", T-1, "."))
  end

  (entropy, avgd)
end

# Initial monotone sequence estimator (IMSE) of Monte Carlo variance
# Geyer C. J. (1992) Practical Markov Chain Monte Carlo.
# Statistical Science, 7 (4), pp 473-483
function imsevar(ac::Vector{Float64},
                 nsim::Int)
  if length(ac) < 2
    throw(KInputError("'maxlag' value is less than 1."))
  end

  u = div(length(ac) - 2, 2)

  s = ac[1] + ac[2]
  old = s
  cur = 0.0
  counter = 0
  for m in 1:u
    counter += 1
    cur = ac[2 * m + 1] + ac[2 * m + 2]

    if cur <= 0.0
      break
    end

    s += (cur <= old) ? cur : old
    old = cur
  end

  if counter == u
    @warn string("'maxlag' value of ", length(ac) - 1, " might not be enough ",
                 "for estimating MCMC variance. Try increasing its value.")
  end

  (-ac[1] + 2 * s) / nsim
end

function ess(variid::Float64,
             varmcmc::Float64,
             nsim::Int)
  floor(Int, nsim * variid / varmcmc)
end


text:
import LinearAlgebra: inv, lu, lutype, copy_oftype, \
using Base: @propagate_inbounds

struct SpectralArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::OffsetArray{T,N,AA}
end

@forward SpectralArray.parent (Base.getindex, Base.setindex, Base.setindex!,
    Base.size, Base.iterate, Base.axes, Base.show, Base.in, Base.strides, Base.elsize,
    Base.dataids)

const SpectralVector{T,AA<:AbstractArray} = SpectralArray{T,1,AA}

"""Alias for `SpectralArray{T,1}`"""
SpectralVector(A::AbstractVector) = SpectralArray(OffsetArray(A, 0:(length(A)-1)))

# UnitRange slicing yields the parent array
@propagate_inbounds Base.getindex(x::SpectralArray, r::Vararg{UnitRange}) =
    (x.parent[r...])

# other slicing yields a SpectralArray
@propagate_inbounds Base.getindex(x::SpectralArray, r::Vararg{AbstractRange}) =
    SpectralArray(x.parent[r...])

"""
    SpectralArray(A::AbstractArray, [ranges])

A renamed OffsetArray. By default, it produces a 0-indexed array.
"""
function SpectralArray(A::AbstractArray{T,N}) where {T,N}
    SpectralArray{T,N,typeof(A)}(OffsetArray(A, map(x->-1, size(A))))
end

function SpectralArray(A::AA, offsets) where {T,N,AA<:AbstractArray{T,N}}
    SpectralArray{T,N,AA}(OffsetArray(A, offsets))
end

function SpectralVector(A::AA, offsets) where {T,AA<:AbstractArray{T,1}}
    SpectralArray{T,1,AA}(OffsetArray(A, offsets))
end

@inline SpectralArray(A::AbstractArray, inds::Vararg) = SpectralArray(OffsetArray(A, inds))

Base.parent(A::SpectralArray) = parent(A.parent)

Base.IndexStyle(::Type{OA}) where {OA<:SpectralArray} = IndexStyle(parenttype(OA))
parenttype(::Type{SpectralArray{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::SpectralArray) = parenttype(typeof(A))

function Base.similar(A::SpectralArray)
    return SpectralArray(similar(parent(A)), A.parent.offsets)
end
# function similar_SA(A::AbstractArray, ::Type{T}, inds::Tuple) where T
#     B = similar(A, T, map(OffsetArray._indexlength, inds))
#     return SpectralArray(OffsetArray(B, map(OffsetArray._offset, axes(B), inds)))
# end
Base.BroadcastStyle(::Type{<:SpectralArray}) = Broadcast.ArrayStyle{SpectralArray}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{SpectralArray}}, ::Type{ElType}) where ElType
    SpectralArray(similar(Array{ElType}, axes(bc)))
end
Broadcast.broadcast_unalias(dest::SpectralArray, src::SpectralArray) =
    parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src)

"""
    spectralzeros(size1, size2, ...)
    spectralzeros(range1, range2, ...)

Utility function for generating a SpectralArray by passing arguments of
ranges or sizes, just like zeros.
"""
spectralzeros(r::Vararg{Int}) = SpectralArray(zeros(r), map(x->-1, r))
spectralzeros(r::Vararg{AbstractRange}) = SpectralArray(zeros(r))

"""
    spectralones(size1::Int, size2::Int, ...)
    spectralones(range1::AbstractRange, range2::AbstractRange, ...)

Utility function for generating a SpectralArray by passing arguments of
ranges or sizes, just like ones.
"""
spectralones(r::Vararg{Int}) = SpectralArray(ones(r), map(x->-1, r))
spectralones(r::Vararg{AbstractRange}) = SpectralArray(ones(r))

function LinearAlgebra.inv(A::SpectralArray{T,2}) where T
    return SpectralArray(inv(parent(A)), A.parent.offsets)
end

# function LinearAlgebra.lu(A::SpectralArray{T,2}, pivot::Union{Val{false}, Val{true}}=Val(true);
#             check::Bool = true) where T
#     S = lutype(T)
#     pA = parent(A)
#     F = lu!(copy_oftype(pA, S), pivot; check = check)
#     return SpectralFactorization(F, A.parent.offsets)
# end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

for (funcr,funcc,func2,typ) in ((:test_csric0_real,:test_csric0_complex,:test_csric02,:CudaSparseMatrixCSR),
                                (:test_cscic0_real,:test_cscic0_complex,:test_cscic02,:CudaSparseMatrixCSC))
    @eval begin
        function $funcr(elty)
            A = rand(elty,m,m)
            A = A + transpose(A)
            A += m * eye(elty,m)
            d_A = Symmetric($typ(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N','S','U',d_A,'O')
            d_B = CUSPARSE.ic0('N','S',d_A,info,'O')
            h_A = to_host(d_B)
            Ac = sparse(full(cholfact(A)))
            h_A = transpose(h_A) * h_A
            @test_approx_eq(h_A.rowval,Ac.rowval)
            @test reduce(&,isfinite(h_A.nzval))
        end

        function $funcc(elty)
            A = rand(elty,m,m)
            A = A + ctranspose(A)
            A += m * eye(elty,m)
            d_A = Hermitian($typ(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N','H','U',d_A,'O')
            d_B = CUSPARSE.ic0('N','H',d_A,info,'O')
            h_A = to_host(d_B)
            Ac = sparse(full(cholfact(A)))
            h_A = ctranspose(h_A) * h_A
            @test_approx_eq(h_A.rowval,Ac.rowval)
            @test reduce(&,isfinite(h_A.nzval))
        end

        function $func2(elty)
            A = rand(elty,m,m)
            A += transpose(A)
            A += m * eye(elty,m)
            d_A = $typ(sparse(tril(A)))
            d_B = CUSPARSE.ic02(d_A,'O')
            h_A = to_host(d_B)
            Ac = sparse(full(cholfact(A)))
            h_A = transpose(h_A) * h_A
            @test_approx_eq(h_A.rowval,Ac.rowval)
            @test reduce(&,isfinite(h_A.nzval))
            A = rand(elty,m,n)
            d_A = $typ(sparse(tril(A)))
            @test_throws(DimensionMismatch,CUSPARSE.ic02(d_A,'O'))
        end
    end
end

types = [Float32,Float64]
for elty in types
    tic()
    test_csric0_real(elty)
    test_cscic0_real(elty)
    println("ic0 took ", toq(), " for ", elty)
    tic()
    test_csric02(elty)
    test_cscic02(elty)
    println("ic02 took ", toq(), " for ", elty)
end

types = [Complex64,Complex128]
for elty in types
    tic()
    test_csric0_complex(elty)
    test_cscic0_complex(elty)
    println("ic0 took ", toq(), " for ", elty)
    tic()
    test_csric02(elty)
    test_cscic02(elty)
    println("ic02 took ", toq(), " for ", elty)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DdosSettings2 <: SwaggerModel
    ddosCustomPolicy::Any # spec type: Union{ Nothing, SubResource } # spec name: ddosCustomPolicy
    protectionCoverage::Any # spec type: Union{ Nothing, String } # spec name: protectionCoverage
    protectedIP::Any # spec type: Union{ Nothing, Bool } # spec name: protectedIP

    function DdosSettings2(;ddosCustomPolicy=nothing, protectionCoverage=nothing, protectedIP=nothing)
        o = new()
        validate_property(DdosSettings2, Symbol("ddosCustomPolicy"), ddosCustomPolicy)
        setfield!(o, Symbol("ddosCustomPolicy"), ddosCustomPolicy)
        validate_property(DdosSettings2, Symbol("protectionCoverage"), protectionCoverage)
        setfield!(o, Symbol("protectionCoverage"), protectionCoverage)
        validate_property(DdosSettings2, Symbol("protectedIP"), protectedIP)
        setfield!(o, Symbol("protectedIP"), protectedIP)
        o
    end
end # type DdosSettings2

const _property_map_DdosSettings2 = Dict{Symbol,Symbol}(Symbol("ddosCustomPolicy")=>Symbol("ddosCustomPolicy"), Symbol("protectionCoverage")=>Symbol("protectionCoverage"), Symbol("protectedIP")=>Symbol("protectedIP"))
const _property_types_DdosSettings2 = Dict{Symbol,String}(Symbol("ddosCustomPolicy")=>"SubResource", Symbol("protectionCoverage")=>"String", Symbol("protectedIP")=>"Bool")
Base.propertynames(::Type{ DdosSettings2 }) = collect(keys(_property_map_DdosSettings2))
Swagger.property_type(::Type{ DdosSettings2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DdosSettings2[name]))}
Swagger.field_name(::Type{ DdosSettings2 }, property_name::Symbol) =  _property_map_DdosSettings2[property_name]

const _allowed_DdosSettings2_protectionCoverage = ["Basic", "Standard"]

function check_required(o::DdosSettings2)
    true
end

function validate_property(::Type{ DdosSettings2 }, name::Symbol, val)
    if name === Symbol("protectionCoverage")
        Swagger.validate_param(name, "DdosSettings2", :enum, val, _allowed_DdosSettings2_protectionCoverage)
    end
end


text:
using TypeInspections
const TI = TypeInspections
using Test

struct Foo{T}
    a::T
end

struct Bar{T}
    b::T
end

struct NestedFoo{TA,TB}
    a::TA
    b::TB
end

@testset "Apply" begin
    @test isempty(TI.apply(x->!TI.isconcretetype_params(x), Foo((1,))))
    @test TI.apply(x->!TI.isconcretetype_params(x), Foo(Any[1]))[1][1] == Vector{Any}

    @test isempty(TI.apply(x->!TI.isconcretetype_params(x), NestedFoo(Foo((1,)), Foo((1,)))))
    @test TI.apply(x->!TI.isconcretetype_params(x), NestedFoo(Foo(Any[1]), Foo((1,))))[1][1] == Vector{Any}

    @test isempty(TI.apply(x->!TI.isconcretetype_params(x), NestedFoo(Bar(Foo((1,))), Bar(Foo((1,))))))
    @test TI.apply(x->!TI.isconcretetype_params(x), NestedFoo(Bar(Foo(Any[1])), Bar(Foo((1,)))))[1][1] == Vector{Any}
end

@testset "warntype_param" begin
    TI.warntype_param(Foo((1,)))
    TI.warntype_param(Foo(Any[1]))

    TI.warntype_param(NestedFoo(Foo((1,)), Foo((1,))))
    TI.warntype_param(NestedFoo(Foo(Any[1]), Foo((1,))))

    TI.warntype_param(NestedFoo(Bar(Foo((1,))), Bar(Foo((1,)))))
    TI.warntype_param(NestedFoo(Bar(Foo(Any[1])), Bar(Foo((1,)))))
end


text:
# Compare with cached output
cachedout = joinpath(testdir, "data")
differentfiles = AbstractString[]
const creator_producer = r"(Creator|Producer)"
for output in readdir(cachedout)
    cached = open(readlines, joinpath(cachedout, output))
    genned = open(readlines, joinpath(testdir, output))
    same = (n=length(cached)) == length(genned)
    if same
        lsame = Bool[cached[i] == genned[i] for i = 1:n]
        if !all(lsame)
            for idx in find(lsame.==false)
                # Don't worry about lines that are due to
                # Creator/Producer (e.g., Cairo versions)
                if !isempty(search(cached[idx], creator_producer))
                    lsame[idx] = true
                end
            end
        end
        same = same & all(lsame)
    end
    if !same
        push!(differentfiles, output)
    else #Delete generated file
        rm(joinpath(testdir, output))
    end
end

#Print out which files differ and their diffs
if length(differentfiles)>0
    #Capture diffs
    diffs = map(
        output -> output * ":\n" *
            readstring(ignorestatus(`diff $(joinpath(cachedout, output)) $(joinpath(testdir, output))`)) *
            "\n\n",
        differentfiles)
    error(string("Generated output differs from cached test output:\n",
        join(differentfiles, "\n"), "\n\n", join(diffs, "\n")))
end


text:
export
    LineSegment,
    parallel

struct LineSegment
    A::VecE2
    B::VecE2
end
Base.:+(seg::LineSegment, V::VecE2) = LineSegment(seg.A + V, seg.B + V)
Base.:-(seg::LineSegment, V::VecE2) = LineSegment(seg.A - V, seg.B - V)
Base.convert(::Type{Line}, seg::LineSegment) = Line(seg.A, seg.B)
Base.convert(::Type{LineSegment}, line::Line) = LineSegment(line.A, line.B)

get_polar_angle(seg::LineSegment) = mod2pi(atan(seg.B.y - seg.A.y, seg.B.x - seg.A.x))

"""
The distance between the line segment and the point P
"""
function get_distance(seg::LineSegment, P::VecE2)

    ab = seg.B - seg.A
    pb = P - seg.A

    denom = normsquared(ab)
    if denom == 0.0
        return 0.0
    end

    r = (ab⋅pb)/denom

    if r ≤ 0.0
        LinearAlgebra.norm(P - seg.A)
    elseif r ≥ 1.0
        LinearAlgebra.norm(P - seg.B)
    else
        LinearAlgebra.norm(P - (seg.A + r*ab))
    end
end

"""
What side of the line you are on, based on A → B
"""
get_side(seg::LineSegment, p::VecE2) = sign((seg.B.x-seg.A.x) * (p.y-seg.A.y) - (seg.B.y-seg.A.y) * (p.x-seg.A.x))

"""
The angular distance between the two line segments
"""
function angledist(segA::LineSegment, segB::LineSegment)
    u = segA.B - segA.A
    v = segB.B - segB.A
    sqdenom = (u⋅u)*(v⋅v)
    if isapprox(sqdenom, 0.0, atol=1e-10)
        return NaN
    end
    return acos((u⋅v) / sqrt(sqdenom))
end

"""
True if the two segments are parallel
"""
function parallel(segA::LineSegment, segB::LineSegment, ε::Float64=1e-10)
    θ = angledist(segA, segB)
    return isapprox(θ, 0.0, atol=ε)
end

"""
Given P colinear with seg, the function checks if
point P lies on the line segment.
"""
function on_segment(P::VecE2, seg::LineSegment)
    return P.x ≤ max(seg.A.x, seg.B.x) && P.x ≥ min(seg.A.x, seg.B.x) &&
           P.y ≤ max(seg.A.y, seg.B.y) && P.y ≥ min(seg.A.y, seg.B.y)
end
on_segment(A::VecE2, P::VecE2, B::VecE2) = on_segment(P, LineSegment(A,B))

"""
returns true if line segments segP and segQ intersect
"""
function intersects(segP::LineSegment, segQ::LineSegment)

    p1 = segP.A
    q1 = segP.B
    p2 = segQ.A
    q2 = segQ.B

    # Find the four orientations needed for general and
    # special cases
    o1 = orientation(p1, q1, p2)
    o2 = orientation(p1, q1, q2)
    o3 = orientation(p2, q2, p1)
    o4 = orientation(p2, q2, q1)

    # Special Cases
    if (o1 != o2 && o3 != o4) || # General case
       (o1 == 0 && on_segment(p1, p2, q1)) || # p1, q1 and p2 are colinear and p2 lies on segment p1q1
       (o2 == 0 && on_segment(p1, q2, q1)) || # p1, q1 and p2 are colinear and q2 lies on segment p1q1
       (o3 == 0 && on_segment(p2, p1, q2)) || # p2, q2 and p1 are colinear and p1 lies on segment p2q2
       (o4 == 0 && on_segment(p2, q1, q2))    # p2, q2 and q1 are colinear and q1 lies on segment p2q2

        return true
    end
    return false # Doesn't fall in any of the above cases
end


text:
# approximate <=
@test _leq(2e-15, 1e-15) && _leq(1e-15, 2e-15)
@test _leq(1//100, 1//99) && _leq(1//100, 0.099)

# approximate <= with ztol
@test !_leq(2e-15, 1e-15, ztol=1e-15) && _leq(1e-15, 2e-15, ztol=1e-15)

# approximate >=
@test _geq(2e-15, 1e-15) && _geq(1e-15, 2e-15)

# approximate >= with ztol
@test _geq(2e-15, 1e-15, ztol=1e-15) && !_geq(1e-15, 2e-15, ztol=1e-15)

# default absolute zero tolerance for rational
_ztol(eltype(1/100)) == zero(Rational{Int})

# default absolute zero tolerance for FP
_ztol(eltype(0.01)) == sqrt(eps(Float64))

# approximately zero tests
@test isapproxzero(0//1)
@test isapproxzero(1e-8) && !isapproxzero(1e-8, ztol=1e-10)

# approximate equality
@test !_isapprox(1//1, 1//1 + 1//1000000000000)
@test _isapprox(2e-15, 1e-15)
@test !_isapprox(2e-6, 1e-6)

# approximate equality for arrays
for N in [Float32, Float64]
    @test _isapprox(N[1, 2, 3], N[1, 2, 3])
    @test _isapprox(N[1 0; 0 1], N[1 0; 0 1])
end


text:
using VarianceComponentsHDFE
using Test

#run_full_test = get(ENV, "VCHDFE_FULL_TEST", "1") == "1" ? true : false
#include("prepare_test_data.jl")

#@testset "VarianceComponentsHDFE.jl" begin
#    include("test_matrices.jl")
#    run_full_test && include("fulltest.jl")
#end


text:
#ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер.
#Робот - в произвольной клетке поля.
#РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.

function moves!(r,side,num)
    for i = 1:num
        if !ismarker(r)
            move!(r,side)
        end
    end
end
function spiral(r)
    global num
    moves!(r,Nord,num)
    moves!(r,Ost,num)
    num += 1
    moves!(r,Sud,num)
    moves!(r,West,num)
    num += 1
end

function mfind(r)
    global num = 1
    while !ismarker(r)
        spiral(r)
    end
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
using DiffEqBiological
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

sir_model = @reaction_network sir_rn begin
  0.5/1000, s + i --> 2i
  0.25, i --> r
end

tmax = 40.0
tspan = (0.0,tmax);

δt = 0.1
t = 0:δt:tmax;

u0 = [990,10,0]; # S,I,R

Random.seed!(1234);

prob = DiscreteProblem(u0,tspan)

prob_jump = JumpProblem(prob,Direct(),sir_model)

sol_jump = solve(prob_jump,SSAStepper());

out_jump = sol_jump(t);

df_jump = DataFrame(out_jump')
df_jump[!,:t] = out_jump.t;

@df df_jump plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_jump,FunctionMap())

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


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

using FEMBase, FEMBeam, Test, LinearAlgebra, SparseArrays

@testset "FEMBeam.jl" begin

    @testset "test beam 2d" begin
        include("test_beam2d.jl")
    end

    @testset "test beam 3D stiffness" begin
        include("test_beam3d_ex1.jl")
    end

    @testset "test beam 3D mass matrix" begin
        include("test_beam3d_mass_matrix.jl")
    end

    @testset "test supports" begin
        include("test_supports.jl")
    end

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

end


text:
type Test
  a::Int
  b::Int
  showTest::Function
  function Test()
    self = new()
    self.a = 1
    self.b = 2
    self.showTest = function()
      self.a, self.b
    end

    return self
  end
end

t1 = Test()
println(t1.showTest())


text:
using VectorizedArgMinMax
using Test

@testset "VectorizedArgMinMax.jl" begin
    @testset "vargmax.jl" begin
        n = 10
        A = rand(n, n, n);
        @test argmax(A) == vargmax(A)
        @test argmax(A, dims=1) == vargmax(A, dims=1)
        @test argmax(A, dims=2) == vargmax(A, dims=2)
        @test argmax(A, dims=3) == vargmax(A, dims=3)
    end
    @testset "vfindmax.jl" begin
        n = 10
        A = rand(n, n, n);
        @test findmax(A) == vfindmax(A)
        @test findmax(A, dims=1) == vfindmax(A, dims=1)
        @test findmax(A, dims=2) == vfindmax(A, dims=2)
        @test findmax(A, dims=3) == vfindmax(A, dims=3)
    end
end


text:
using NetworkEpidemics
using LightGraphs
using Plots
using Random
using Statistics
using LinearAlgebra
using GraphRecipes

cd("metaplex-networks/old")

include("networks/empirical_networks.jl")

Random.seed!(2021);

N = 1000
M = 20

h = contiguous_usa()
h = london_transport()
h = erdos_renyi(M, 0.3)

M = nv(h)
N = M*1000
graphplot(h)

begin
    p = sortperm(degree(h))
    eig = eigen(Array(normalized_laplacian(h)))
    V = Array{Float64,2}(undef, nv(h), nv(h))
    for i in 1:nv(h)
        V[:,i] .= eig.vectors[p,i]
    end
    # for each node, find the eigenvector whose entry for that node is highest in magnitude
    eigenmode = Vector{Int}(undef, nv(h))
    for i in 1:nv(h)
        eigenmode[i] = findmax(abs.(eig.vectors[i,:]))[2]
    end
end

histogram(degree(h), label="")
scatter!(degree_histogram(h), label="")

scatter(eig.values, label="")

heatmap(V', color=cgrad(:balance))

bar(eigenmode, label="")

β = 0.39
γ = 0.4
D = 0.1

critical_k = (N / M) * γ / β

ks = fill(50, M)
ks_ = copy(ks)
ks_[1] = 400

gs = [random_regular_graph(N, k) for k in ks]
mpx = HeterogeneousMetaplex(gs, h, D, SIS(β,γ))

mp = HeterogeneousMetapopulation(h, N, ks, D, SIS(β,γ))
mp_ = HeterogeneousMetapopulation(h, N, ks_, D, SIS(β,γ))

x0_i = fill(1, N)
x0_i[1:10] .= 2

x0_μ = [div(i-1,div(N,M))+1 for i in 1:N]
x0_μ = [mod(i,M)+1 for i in 1:N]

x0 = Array{Int,2}(undef, M, 2)
x0[2:M,:] .= [1000 0]
x0[1,:] .= [980, 20]

tmax = 100.0
nmax = 1000000
nsims = 100
nbins = 200

#ts_mpx, u_mpx = gillespie(mpx, [x0_i,x0_μ], nmax=nmax, tmax=tmax)
ts_mp, u_mp = gillespie(mp, x0, nmax=nmax, tmax=tmax)

u_mp

#plot(ts_mpx, u_mpx[2], label="")
plot(ts_mp, u_mp[2], label="")

#ts_av_mpx, u_av_mpx = average(mpx, [x0_i, x0_μ], nmax=nmax, tmax=tmax, nbins=nbins, nsims=nsims)
#ts_mf_mpx, u_mf_mpx = meanfield(mpx, [x0_i,x0_μ], tmax=tmax)

ts_av_mp, u_av_mp = average(mp, x0, nmax=nmax, tmax=tmax, nbins=nbins, nsims=nsims, progressbar=true)
ts_mf_mp, u_mf_mp = meanfield(mp, x0, tmax=200.0)
ts_mf_mp_, u_mf_mp_ = meanfield(mp_, x0, tmax=200.0)
#plot(ts_av_mpx, sum(u_av_mpx[2], dims=2)/N, label="Average", ylims=(0.0,1.0))
#plot!(ts_mf_mpx, sum(u_mf_mpx[2], dims=2)/N,  label="Meanfield")

plot(ts_av_mp, sum(u_av_mp[2], dims=2)/N, label="Average")
plot!(ts_mf_mp, sum(u_mf_mp[2], dims=2)/N, label="Meanfield")

plot(ts_mf_mp, sum(u_mf_mp[2], dims=2)/N, label="Meanfield", legend=:topleft)
plot!(ts_mf_mp_, sum(u_mf_mp_[2], dims=2)/N, label="Meanfield", legend=:topleft)

plot(ts_av_mp, u_av_mp[2], label="")
plot(ts_mf_mp, u_mf_mp[2], label="")

#plot(ts_av_mpx, u_av_mpx[2], label="")
#plot(ts_mf_mpx, u_mf_mpx[2], label="")

v = inv(eig.vectors)
ψ_mf_mp = copy(u_mf_mp[2])
for i in 1:size(u_mf_mp)[1]
   ψ_mf_mp[i,:] .= v * u_mf_mp[2][i,:]
end
plot(ts_mf_mp, ψ_mf_mp, label="")

plot(ts_mf_mp, (eig.vectors * u_mf_mp[2]')', label="")

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
#!/usr/bin/env julia

using Luxor, Colors

using Test

using Random
Random.seed!(42)

function poly_areas()
    # test area calculations
    @test polyarea(ngon(O, 100, 4, 0, vertices=true)) ≈ 20_000
    @test polyarea(ngonside(O, 100, 4, 0, vertices=true)) < polyarea(ngonside(O, 100, 5, 0, vertices=true))
    @test polyarea(ngonside(O, 100, 3, 0, vertices=true)) ≈ 2500sqrt(3)
end

function simple_polys()
    gsave()
    sethue(0, 0, 0)
    vstep = 50
    hstep = 50
    x = 0
    y = 0
    translate(-currentwidth/2 + 50, -currentheight/2 + 50)
    for sides in 3:12
        ngon(x, y, 20, sides, 0, :stroke)
        x += hstep
        if x > currentwidth
            x = 0
            y += vstep
        end
    end
    translate(0, currentheight - 150)
    x = 0
    for sides in 3:12
        ngon(x, y, 20, sides, 0, :fill)
        x += hstep
        if x > currentwidth
            x = 0
            y += vstep
        end
    end
    grestore()

    @layer begin
        translate(BoundingBox()[1] + (0, 130))
        pgon = ngonside(O, 20, 5, vertices=true)
        for i in 1:10
            polymove!(pgon, O, Point(50, 0))
            poly(pgon, :stroke, close=true)
            polyrotate!(pgon, rand() * π, center=polycentroid(pgon))
            polyscale!(pgon, isodd(i) ? 1.5 : 1/1.5, center=polycentroid(pgon))
        end
    end
end

function hex_mixtures()
    gsave()
    vstep = 50
    hstep = 50
    for y in -currentheight:vstep:currentheight
        for x in -currentwidth:hstep:currentwidth
            setopacity(0.2)
            setline(rand()*10)
            sethue(backgroundcolors[rand(1:end)])
            gsave()
                 translate(x, y)
                 rotate(rand() * pi)
                 # point version
                 poly(hexagon(rand() * 5, rand() * 5, 15 + rand() * 5), :fill, close=true)
                 setopacity(0.1)
                 poly(hexagon(rand() * 10, rand() * 5, 15 + rand() * 5))
                 setopacity(0.3)
                 randomhue()
                 poly([
                        Point(rand() * -100, rand() * 100),
                        Point(rand() * -100, rand() * 100),
                        Point(rand() * -100, rand() * 100),
                        Point(rand() * -100, rand() * 100)
                     ], :fill)
            grestore()
            sethue(backgroundcolors[rand(1:end)])
            gsave()
                 translate(x, y)
                 ngon(x, y, 50, rand(3:13), 0, :stroke)
                 rotate(rand() * pi)
                 poly(hexagon(rand() * 5, rand() * 5, 15 + rand() * 5), :stroke)
                 ngon(x, y, 50, rand(3:13), 0, :fill)
                 rect(0, 0, hstep-5, vstep-5, :stroke)
                 ngon(x, y, 50, rand(3:13), 0, :none)
            grestore()
        end
    end
    grestore()
end

hexagon(x, y, size) = [Point(x + size * cos(2pi/6 * i), y + size * sin(2pi/6 * i)) for i in 1:6]

function polygon_test(fname)
    global currentwidth = 595 # pts
    global currentheight = 842 # pts
    Drawing(currentwidth, currentheight, fname)

    origin()
    background("antiquewhite")

    # test polyareas
    poly_areas()

    global foregroundcolors = diverging_palette(230, 280, 200, s = 0.99, b=0.8)
    global backgroundcolors = diverging_palette(200, 260, 280, s = 0.8, b=0.5)
    setline(2)
    # fill, then clip to heptagon
    setcolor("lightcyan")
    setline(3)
    ngon(0, 0, 270, 7, 0, :fillpreserve) # fill it
    sethue("orange")
    strokepreserve()                     # stroke it
    clip()                               # then use to clip
    hex_mixtures()
    clipreset()
    simple_polys()
    @test finish() == true
    println("...finished test: output in $(fname)")
end

polygon_test("polygon-test.pdf")


text:
using Flux
using Flux.Tracker
using Flux.Optimise: ADAM
using Gym: reset!, make, step!, render!
using DataStructures: CircularBuffer

#-----------------Environment Setup--------------------------------------------#
env = make("CartPole-v0", :human_pane)
reset!(env)
ACTION_SPACE = 2
OBSERVATION_SPACE = 4
REWARD_THRESHOLD = 200
MAX_EPISODES = 5000

#---------------------Model Architechture--------------------------------------#
model = Chain(
    Dense(OBSERVATION_SPACE, 128, relu),
    Dense(128, ACTION_SPACE),
    softmax
)
opt = ADAM(3f-4)
const γ = 9f-1

#W1, b1, W2, b2 = params(model).order

#-------------------- Functions -----------------------------------------------#

include("forwardpass.jl")
include("backwardpass.jl")
include("utils.jl")

function episode()
    current_state = reset!(env)
    rewards = Float32[]
    log_probs = []

    terminal = false

    while !terminal
        action, log_prob = execpolicy(model, current_state)
        next_state, reward, terminal, _ = step!(env, action)

        push!(rewards, reward)
        push!(log_probs, log_prob)

        current_state = next_state
    end
    return rewards, log_probs
end

function learnpolicy()
    eps_rewards = Float32[]
    mean_rewards = Float32[]

    for eps=1:MAX_EPISODES
        rewards, log_probs = episode()
        policyupdate!(rewards, log_probs)
        eps_reward = sum(rewards)
        push!(eps_rewards, eps_reward)

        mean_reward = eps < 10 ? mean(eps_reward) : mean(eps_rewards[end-9:end])
        push!(mean_rewards, mean_reward)
        println("Episode $eps = $eps_reward | Mean of the last 10 episodes = $mean_reward")
        mean_reward ≥ REWARD_THRESHOLD && break
    end

    plotresults(mean_rewards, eps_rewards)
end

function policydemo()
    current_state = reset!(env)
    terminal = false
    for frame_num ∈ 1:500
        action, _ = execpolicy(model, current_state)
        next_state, reward, terminal, _ = step!(env, action)

        render!(env)

        terminal && (next_state = reset!(env))

        current_state = next_state
        sleep(0.001)
    end
end

# Driver code
function REINFORCE()
    learnpolicy()
    policydemo()
end

REINFORCE()


text:
using Environment
using Test

@testset "Environment" begin

include("space.jl")
include("gymenv.jl")

end


text:
# ------------------------------------------------------------------
# Code to simulate logit data and estimate nonparametrically
# Plots estimates of own and cross-price elasticities
# Written by James Brand
# ------------------------------------------------------------------
using Statistics, NPDemand
using RCall, DataFrames
@rlibrary ggplot2

J=2; # # of products
T =2000; # # of markets
beta = -0.4; # price coefficient
sdxi = 0.15; # standard deviation of xi

S = 200;
G = 20; # size of grid on which to evaluate price elasticities
elast_own = zeros(S,G);
elast_cross = zeros(S,G);
elast_own_dist = zeros(S,T);
elastTrue = zeros(S,G);

s, p, z, xi = simulate_logit(J,T, beta, sdxi);
p_points = range(quantile(p[:,1],.25),stop = quantile(p[:,1],.75),length = G);
p_points = convert.(Float64, p_points)

# ------------------------------------------------------
# Set options for estimation and elasticity calculation
# ------------------------------------------------------
included = ones(J,J);         # Matrix of ones and zeros, specifying which products are substitutes.
                                # Each row corresponds to an inverse demand function.
                                # Set element (j,k) to 1 if k is a substitute for j
                                # Mostly helpful if you're incorporating model selection prior
                                # to estimation. Otherwise set to ones(J).

bernO = 2*ones(J,1);        # Order of Bernstein Polynomial
iv=0;                       # Order of IV Polynomial = (bernO + iv)
constrained = false;            # Monotonicity Constraint (experience says this can help but much slower in current build)
x = zeros(size(p));       # No exogenous product characteristics
trueS = true;                    # Evaluate at true market shares or not
own = [1,1];                # [derivative of j, with respect to price of k]
cross = [1,2];

# ------------------------------------------------------
# Simulation
# ------------------------------------------------------
for si = 1:1:S
    # Returns market shares, prices, instruments, and the market demand shock, respectively
    s, p, z  = simulate_logit(J,T, beta, sdxi);
    df = toDataFrame(s,p,z);
    # Estimate demand nonparametrically
        # If you want to include an additional covariate in all demand
        # functions, add an additional argument "marketvars" after included. If it is an
        # additional product characteristic, marketvars should be T x J
    inv_sigma, designs = inverse_demand(df; included = included);

    # Calculate price elasticities
    prices = median(p).*ones(G,J);
    prices[:,1] = p_points;

    elast, Jacobians, share_vec = price_elasticity(inv_sigma, df, prices; included = included,
        whichProducts = own, trueS = false)

    trueelast = beta.*p_points.*(1 .-  share_vec[:,1])

    elast_own[si,:] = elast;
    elastTrue[si,:] = trueelast;
end

# Now generate plots using RCall and ggplot
elast025 = zeros(G,1)
elast50 = zeros(G,1)
elast975 = zeros(G,1)
for i = 1:G
    elast025[i] = quantile(elast_own[:,i], 0.25)
    elast50[i] = quantile(elast_own[:,i], 0.5)
    elast975[i] = quantile(elast_own[:,i], 0.975)
end
elast025 = dropdims(elast025,dims = 2)
elast50 = dropdims(elast50,dims = 2)
elast975 = dropdims(elast975,dims = 2)

df = DataFrame(p = p_points, e025 = elast025, e50 = elast50, e975 = elast975)

# Plot 95% interval of price elasticities along interquartile range
ggplot(df, aes(x=:p, y=:e50)) + geom_line() + geom_line(aes(y=:e975), color = "gray", linetype = "dashed") +
    geom_line(aes(y=:e025), color = "gray", linetype = "dashed") +
    xlab("Price") + ylab("Own-Elasticity") + theme_light()

# ------------------------------------------------------
# ------------------------------------------------------
# You can also calculate price elasticities at the realized prices and market shares
        # Jacobians2 is particularly helpful here. This is an array, each element of which is the Jacobian
        # of demand with respect to prices. This can be used to calculate markups under many models
        # of competition.
s, p, z  = simulate_logit(J,T, beta, sdxi);
df = toDataFrame(s,p,z);
inv_sigma, designs = inverse_demand(df; included = included);
elast2, Jacobians2 = price_elasticity(inv_sigma, df, p, included = included; deltas = -1 .* p);

df2 = DataFrame(Estimate = elast2, True = beta.*p[:,1].*(1 .- s[:,1]))
ggplot(df2, aes(x=:True))+
 geom_density(aes(x=:Estimate), color = "gray", linetype = "dashed") +
 geom_density(aes(x=:True), color = "black") + xlab("Elasticity") +
 ylab("Density")