text:
#=======================================,2)
 Fully Independent Training Conditional
=========================================#

"""
    Positive Definite Matrix for Fully Independent Training Conditional approximation.
"""
mutable struct FullyIndepPDMat{T,M<:AbstractMatrix,PD<:AbstractPDMat{T},M2<:AbstractMatrix{T},DIAG<:Diagonal} <: SparsePDMat{T}
    inducing::M
    ΣQR_PD::PD
    Kuu::PD
    Kuf::M2
    Λ::DIAG
end
size(a::FullyIndepPDMat) = (size(a.Kuf,2), size(a.Kuf,2))
size(a::FullyIndepPDMat, d::Int) = size(a.Kuf,2)
"""
    We have
        Σ ≈ Kuf' Kuu⁻¹ Kuf + Λ
    where Λ is a diagonal matrix (here given by σ²I + diag(Kff - Qff))
    The Woodbury matrix identity gives
        (A+UCV)⁻¹ = A⁻¹ - A⁻¹ U(C⁻¹ + V A⁻¹ U)⁻¹ V A⁻¹
    which we use here with
        A ← Λ
        U = Kuf'
        V = Kuf
        C = Kuu⁻¹
    which gives
        Σ⁻¹ = Λ⁻¹ - Λ⁻² Kuf'(Kuu + Kuf Λ⁻¹ Kuf')⁻¹ Kuf
            = Λ⁻¹ - Λ⁻² Kuf'(        ΣQR       )⁻¹ Kuf
"""
function \(a::FullyIndepPDMat, x)
    Lk = whiten(a.ΣQR_PD, a.Kuf)
    return a.Λ \ (x .- Lk' * (Lk * (a.Λ \ x)))
end

"""
    trinvAB(A::AbstractPDMat, B)

    Computes tr(A⁻¹ B).
"""
function trinvAB(A::AbstractPDMat, B)
    tr(A \ B)
end
"""
    trinvAB(A::FullyIndepPDMat, B::Diagonal)

    Computes tr(A⁻¹ B) efficiently under the FITC approximation:
        Σ ≈ Kuf' Kuu⁻¹ Kuf + Λ

    Derivation:
        tr(Σ⁻¹ B) = tr[ (Λ⁻¹ - Λ⁻² Kuf'(Kuu + Kuf Λ⁻¹ Kuf')⁻¹ Kuf) B ]
                  = tr(Λ⁻¹ B) - tr( Λ⁻¹ Kuf' ΣQR⁻¹ Kuf Λ⁻¹ B )
                  = tr(Λ⁻¹ B) - tr( Λ⁻¹ Kuf' ΣQR^(-1/2) ΣQR^(-1/2) Kuf Λ⁻¹ B )
                                                        ╰────────────────╯
                                                            ≡ L
                  = tr(Λ⁻¹ B) - tr(L'*L*B)
                  = tr(Λ⁻¹ B) - dot(L, L*B)

    See also: [`\`](@ref).
"""
function trinvAB(A::FullyIndepPDMat, B::Diagonal)
    Λ = A.Λ
    L = whiten(A.ΣQR_PD, A.Kuf * inv(Λ))
    return tr(Λ \ B) - dot(L, L * B)
end

"""
    The matrix determinant lemma states that
        logdet(A+UWV') = logdet(W⁻¹ + V'A⁻¹U) + logdet(W) + logdet(A)
    So for
        Σ ≈ Kuf' Kuu⁻¹ Kuf + Λ
        logdet(Σ) = logdet(Kuu + Kuf Λ⁻¹ Kuf')       + logdet(Kuu⁻¹) + logdet(Λ)
                  = logdet(        ΣQR             ) - logdet(Kuu)   + logdet(Λ)
"""
logdet(a::FullyIndepPDMat) = logdet(a.ΣQR_PD) - logdet(a.Kuu) + logdet(a.Λ)#sum(log.(a.Λ))
function Base.Matrix(a::FullyIndepPDMat)
    Lk = whiten(a.Kuu, a.Kuf)
    Σ = Lk'Lk + a.Λ
    nobs = size(Σ,1)
    return Σ
end

function wrap_cK(cK::FullyIndepPDMat, inducing, ΣQR_PD, Kuu, Kuf, Λ::Diagonal)
    FullyIndepPDMat(inducing, ΣQR_PD, Kuu, Kuf, Λ)
end
function wrap_cK(cK::FullyIndepPDMat, inducing, ΣQR_PD, Kuu, Kuf, Λ::Vector)
    FullyIndepPDMat(inducing, ΣQR_PD, Kuu, Kuf, Diagonal(Λ))
end
"""
    tr(a::FullyIndepPDMat)

    Trace of the FITC approximation to the covariance matrix:

    tr(Σ) = tr(Kuf' Kuu⁻¹ Kuf + Λ)
          = tr(Kuf' Kuu⁻¹ Kuf) + tr(Λ)
          = tr(Kuf' Kuu^{-1/2} Kuu^{-1/2} Kuf) + tr(Λ)
                              ╰──────────────╯
                                 ≡  Lk
          = dot(Lk, Lk) + sum(diag(Λ))
"""
function LinearAlgebra.tr(a::FullyIndepPDMat)
    Lk = whiten(a.Kuu, a.Kuf)
    return tr(a.Λ) + dot(Lk, Lk)
end

"""
    Fully Independent Training Conditional (FITC) covariance strategy.
"""
struct FullyIndepStrat{M<:AbstractMatrix} <: SparseStrategy
    inducing::M
end
SubsetOfRegsStrategy(fitc::FullyIndepStrat) = SubsetOfRegsStrategy(fitc.inducing)
DeterminTrainCondStrat(fitc::FullyIndepStrat) = DeterminTrainCondStrat(fitc.inducing)

function alloc_cK(covstrat::FullyIndepStrat, nobs)
    # The objects that need to be allocated are very similar
    # to SoR, so we'll use that as a starting point:
    SoR = SubsetOfRegsStrategy(covstrat)
    cK_SoR = alloc_cK(SoR, nobs)
    # Additionally, we need to store the diagonal corrections.
    Λ = Diagonal(Vector{Float64}(undef, nobs))

    cK_FITC = FullyIndepPDMat(
        covstrat.inducing,
        cK_SoR.ΣQR_PD,
        cK_SoR.Kuu,
        cK_SoR.Kuf,
        Λ)
    return cK_FITC
end
function update_cK!(cK::FullyIndepPDMat, X::AbstractMatrix, kernel::Kernel,
                    logNoise::Real, kerneldata::KernelData, covstrat::FullyIndepStrat)
    inducing = covstrat.inducing
    Kuu = cK.Kuu
    Kuubuffer = mat(Kuu)
    cov!(Kuubuffer, kernel, inducing)
    Kuubuffer, chol = make_posdef!(Kuubuffer, cholfactors(cK.Kuu))
    Kuu_PD = wrap_cK(cK.Kuu, Kuubuffer, chol)
    Kuf = cov!(cK.Kuf, kernel, inducing, X)
    Kfu = Kuf'

    dim, nobs = size(X)
    Kdiag = [cov_ij(kernel, X, X, EmptyData(), i, i, dim) for i in 1:nobs]
    Qdiag = [invquad(Kuu_PD, Kuf[:,i]) for i in 1:nobs]
    Λ = Diagonal(exp(2*logNoise) .+ Kdiag .- Qdiag)

    ΣQR = Kuf * (Λ \ Kfu) + Kuu
    LinearAlgebra.copytri!(ΣQR, 'U')

    ΣQR, chol = make_posdef!(ΣQR, cholfactors(cK.ΣQR_PD))
    ΣQR_PD = wrap_cK(cK.ΣQR_PD, ΣQR, chol)
    return wrap_cK(cK, inducing, ΣQR_PD, Kuu_PD, Kuf, Λ)
end

#==========================================
  Log-likelihood gradients
===========================================#
function init_precompute(covstrat::FullyIndepStrat, X, y, kernel)
    # here we can re-use the Subset of Regressors pre-computations
    SoR = SubsetOfRegsStrategy(covstrat)
    return init_precompute(SoR, X, y, kernel)
end

"""
dmll_kern!(dmll::AbstractVector, kernel::Kernel, X::AbstractMatrix, cK::AbstractPDMat, kerneldata::KernelData,
                    alpha::AbstractVector, Kuu, Kuf, Kuu⁻¹Kuf, Kuu⁻¹KufΣ⁻¹y, Σ⁻¹Kfu, ∂Kuu, ∂Kuf,
                    covstrat::FullyIndepStrat)
Derivative of the log likelihood under the Fully Independent Training Conditional (FITC) approximation.

Helpful reference: Vanhatalo, Jarno, and Aki Vehtari.
                   "Sparse log Gaussian processes via MCMC for spatial epidemiology."
                   In Gaussian processes in practice, pp. 73-89. 2007.

Generally, for a multivariate normal with zero mean
    ∂logp(Y|θ) = 1/2 y' Σ⁻¹ ∂Σ Σ⁻¹ y - 1/2 tr(Σ⁻¹ ∂Σ)
                    ╰───────────────╯     ╰──────────╯
                           `V`                 `T`

where Σ = Kff + σ²I.

Notation: `f` is the observations, `u` is the inducing points.
          ∂X stands for ∂X/∂θ, where θ is the kernel hyperparameters.

In the case of the FITC approximation, we have
    Σ = Λ + Qff
    The second component gives
    ∂(Qff) = ∂(Kfu Kuu⁻¹ Kuf)
    which is used by the gradient function for the subset of regressors approximation,
    and so I don't repeat it here.

    the ith element of diag(Kff-Qff) is
    Λi = Kii - Qii + σ²
       = Kii - Kui' Kuu⁻¹ Kui + σ²
    ∂Λi = ∂Kii - Kui' ∂(Kuu⁻¹) Kui - 2 ∂Kui' Kuu⁻¹ Kui
        = ∂Kii + Kui' Kuu⁻¹ ∂(Kuu) Kuu⁻¹ Kui - 2 ∂Kui' Kuu⁻¹ Kui
"""
function dmll_kern!(dmll::AbstractVector, kernel::Kernel, X::AbstractMatrix,
                    cK::AbstractPDMat, kerneldata::SparseKernelData,
                    alpha::AbstractVector, Kuu, Kuf, Kuu⁻¹Kuf, Kuu⁻¹KufΣ⁻¹y, Σ⁻¹Kfu, ∂Kuu, ∂Kuf,
                    covstrat::FullyIndepStrat)
    # first compute the SoR component
    SoR = SubsetOfRegsStrategy(covstrat)
    dmll_kern!(dmll, kernel, X, cK, kerneldata, alpha,
               Kuu, Kuf, Kuu⁻¹Kuf, Kuu⁻¹KufΣ⁻¹y, Σ⁻¹Kfu, ∂Kuu, ∂Kuf,
               SoR)
    nparams = num_params(kernel)
    dim, nobs = size(X)
    inducing = covstrat.inducing
    for iparam in 1:nparams
        # TODO: the grad_slice! calls here are redundant with the ones in the dmll_kern! call above
        grad_slice!(∂Kuu, kernel, inducing, inducing, kerneldata.Kuu, iparam)
        grad_slice!(∂Kuf, kernel, inducing, X       , kerneldata.Kux1, iparam)

        ∂Λ = Diagonal([
               (dKij_dθp(kernel, X, X, EmptyData(), i, i, iparam, dim) # ∂Kii
                + dot(Kuu⁻¹Kuf[:,i], ∂Kuu * Kuu⁻¹Kuf[:,i])  # Kui' Kuu⁻¹ ∂(Kuu) Kuu⁻¹ Kui
                - 2 * dot(∂Kuf[:,i], Kuu⁻¹Kuf[:,i]))        # -2 ∂Kui' Kuu⁻¹ Kui
               for i in 1:nobs
               ])
        V = dot(alpha, ∂Λ * alpha)
        T = trinvAB(cK, ∂Λ)
        # # FOR DEBUG ONLY
        # Talt = tr(cK \ ∂Λ) # inefficient
        # @show T, Talt
        # @assert isapprox(T, Talt, atol=1e-5)
        # # END DEBUG

        dmll[iparam] += (V-T)/2
    end
    return dmll
end
function dmll_kern!(dmll::AbstractVector, gp::GPBase, precomp::SoRPrecompute, covstrat::FullyIndepStrat)
    return dmll_kern!(dmll, gp.kernel, gp.x, gp.cK, gp.data, gp.alpha,
                      gp.cK.Kuu, gp.cK.Kuf,
                      precomp.Kuu⁻¹Kuf, precomp.Kuu⁻¹KufΣ⁻¹y, precomp.Σ⁻¹Kfu,
                      precomp.∂Kuu, precomp.∂Kuf,
                      covstrat)
end

"""
    dmll_noise(gp::GPE, precomp::SoRPrecompute, covstrat::FullyIndepStrat)

∂logp(Y|θ) = 1/2 y' Σ⁻¹ ∂Σ Σ⁻¹ y - 1/2 tr(Σ⁻¹ ∂Σ)

∂Σ = I for derivative wrt σ², so
∂logp(Y|θ) = 1/2 y' Σ⁻¹ Σ⁻¹ y - 1/2 tr(Σ⁻¹)
            = 1/2[ dot(α,α) - tr(Σ⁻¹) ]

We have:
    Σ⁻¹ = Λ⁻¹ - Λ⁻² Kuf'(        ΣQR       )⁻¹ Kuf
Use the identity tr(A'A) = dot(A,A) to get:
    tr(Σ⁻¹) = tr(Λ⁻¹) - dot(Lk, Lk) .
where
    Lk ≡ ΣQR^(-1/2) Kuf Λ⁻¹
"""
function dmll_noise(gp::GPE, precomp::SoRPrecompute, covstrat::FullyIndepStrat)
    nobs = gp.nobs
    cK = gp.cK
    Λ = cK.Λ
    Lk = whiten(cK.ΣQR_PD, cK.Kuf) * inv(Λ)
    # # DEBUG
    # noiseT = sum(1 ./ Λ) - dot(Lk, Lk)
    # Talt = tr(cK \ Matrix(1.0*I, nobs, nobs))
    # @show noiseT, Talt # should be same
    # # END DEBUG
    return exp(2*gp.logNoise) * ( # Jacobian
        dot(gp.alpha, gp.alpha)
        - tr(inv(Λ)) # sum(1 ./ Λ)
        + dot(Lk, Lk)
        )
end

function get_alpha_u(Ktrain::FullyIndepPDMat, xtrain::AbstractMatrix, ytrain::AbstractVector, meanf::Mean)
    ΣQR_PD = Ktrain.ΣQR_PD
    Kuf = Ktrain.Kuf
    meantrain = mean(meanf, xtrain)
    Λ = Ktrain.Λ
    alpha_u = ΣQR_PD \ (Kuf * (Λ \ (ytrain-meantrain)) )
    return alpha_u
end

"""
predictMVN(xpred::AbstractMatrix, xtrain::AbstractMatrix, ytrain::AbstractVector,
                    kernel::Kernel, meanf::Mean, logNoise::Real,
                    alpha::AbstractVector,
                    covstrat::FullyIndepStrat, Ktrain::FullyIndepPDMat)
    See Quiñonero-Candela and Rasmussen 2005, equations 24b.
    Some derivations can be found below that are not spelled out in the paper.

    Notation: Qab = Kau Kuu⁻¹ Kub
              ΣQR = Kuu + σ⁻² Kuf Kuf'

              x: prediction (test) locations
              f: training (observed) locations
              u: inducing point locations

    We have
        Σ ≈ Kuf' Kuu⁻¹ Kuf + Λ
    By Woodbury
        Σ⁻¹ = Λ⁻¹ - Λ⁻² Kuf'(Kuu + Kuf Λ⁻¹ Kuf')⁻¹ Kuf
            = Λ⁻¹ - Λ⁻² Kuf'(        ΣQR       )⁻¹ Kuf

    The predictive mean can be derived (assuming zero mean function for simplicity)
    μ = Qxf (Qff + Λ)⁻¹ y
      = Kxu Kuu⁻¹ Kuf [Λ⁻¹ - Λ⁻² Kuf' ΣQR⁻¹ Kuf] y   # see Woodbury formula above.
      = Kxu Kuu⁻¹ [ΣQR - Kuf Λ⁻¹ Kfu] ΣQR⁻¹ Kuf Λ⁻¹ y # factoring out common terms
      = Kxu Kuu⁻¹ [Kuu] ΣQR⁻¹ Kuf Λ⁻¹ y               # using definition of ΣQR
      = Kxu ΣQR⁻¹ Kuf Λ⁻¹ y                           # matches equation 16b

    Similarly for the posterior predictive covariance:
    Σ = Σxx - Qxf (Qff + Λ)⁻¹ Qxf'
      = Σxx - Kxu ΣQR⁻¹ Kuf Λ⁻¹ Qxf'                # substituting result from μ
      = Σxx - Kxu ΣQR⁻¹  Kuf Λ⁻¹ Kfu Kuu⁻¹ Kux      # definition of Qxf
      = Σxx - Kxu ΣQR⁻¹ (ΣQR - Kuu) Kuu⁻¹ Kux       # using definition of ΣQR
      = Σxx - Kxu Kuu⁻¹ Kux + Kxu ΣQR⁻¹ Kux         # expanding
      = Σxx - Qxx           + Kxu ΣQR⁻¹ Kux         # definition of Qxx
"""
function predictMVN(xpred::AbstractMatrix,
                    xtrain::AbstractMatrix, ytrain::AbstractVector,
                    kernel::Kernel, meanf::Mean,
                    alpha::AbstractVector,
                    covstrat::FullyIndepStrat, Ktrain::FullyIndepPDMat)
    DTC = DeterminTrainCondStrat(covstrat)
    μ_DTC, Σ_DTC = predictMVN(xpred, xtrain, ytrain, kernel, meanf, alpha, DTC, Ktrain)
    return μ_DTC, Σ_DTC
end

function FITC(x::AbstractMatrix, inducing::AbstractMatrix, y::AbstractVector, mean::Mean, kernel::Kernel, logNoise::Real)
    covstrat = FullyIndepStrat(inducing)
    GPE(x, y, mean, kernel, logNoise, covstrat)
end


text:
using CUDA
using ExaPF
using KernelAbstractions
using Test

import ExaPF: PowerSystem

@testset "Powerflow solver" begin
    datafile = joinpath(dirname(@__FILE__), "..", "data", "case14.raw")
    pf = PowerSystem.PowerNetwork(datafile)
    # Parameters
    npartitions = 8
    tolerance = 1e-6
    if has_cuda_gpu()
        DEVICES = [CPU(), CUDADevice()]
    else
        DEVICES = [CPU()]
    end

    @testset "Deport computation on device $device" for device in DEVICES
        polar = PolarForm(pf, device)
        jac = ExaPF._state_jacobian(polar)
        precond = ExaPF.LinearSolvers.BlockJacobiPreconditioner(jac, npartitions, device)
        # Retrieve initial state of network
        x0 = ExaPF.initial(polar, State())
        uk = ExaPF.initial(polar, Control())
        p = ExaPF.initial(polar, Parameters())

        @testset "Powerflow solver $(LinSolver)" for LinSolver in ExaPF.list_solvers(device)
            algo = LinSolver(precond)
            xk = copy(x0)
            nlp = ExaPF.ReducedSpaceEvaluator(polar, xk, uk, p;
                                              ε_tol=tolerance, linear_solver=algo)
            convergence = ExaPF.update!(nlp, uk; verbose_level=ExaPF.VERBOSE_LEVEL_NONE)
            @test convergence.has_converged
            @test convergence.norm_residuals < tolerance
            @test convergence.n_iterations == 2
        end
    end
end

ext: jl


text:
# TODO: make sure this is safe
# multiple sharing?  other writes before reads?
function finddif(r::Net, n, nsparse)
    dif0 = nothing
    if (!isa(r.op[n], Par) &&
        !r.toincr[n])
        for i=n+1:length(r.op)
            if (isassigned(r.dif0, i) &&
                in(n, r.inputs[i]) &&
                overwrites(r.op[i]) &&
                !r.toincr[i])
                dif0 = r.dif0[i]
                break
            end
        end
    end
    if dif0 == nothing
        dif0 = newarray(gpu(), nsparse, eltype(r.out0[n]), size(r.out0[n]))
    end
    return dif0
end

function findtmp(r::Net, n, nsparse)
    tmp = nothing
    for i=n+1:length(r.op)
        if (isassigned(r.tmp, i) &&
            size(r.tmp[i]) == size(r.dif0[i]) &&
            issparse(r.tmp[i]) == nsparse)
            tmp = r.tmp[i]
            break
        end
    end
    if tmp == nothing
        tmp = newarray(gpu(), nsparse, eltype(r.dif0[n]), size(r.dif0[n]))
    end
    return tmp
end

function initback(r::Net, dy, dx...; seq=false, a...)
    @assert dx == () || length(dx) == ninputs(r)
    N = length(r.op)
    for n=1:N; r.toincr[n] = (r.multi[n] || (seq && isa(r.op[n], Par))); end
    isassigned(r.dif0, 1) || initback0(r, dy, dx...; a...)
    fill!(r.dif, nothing)
    for n=1:length(r.op)                                # TODO-OPTIMIZATION
        isassigned(r.dif0, i) && fill!(r.dif0[i], 0)
    end
end

function initback0(r::Net, dy, dx...; a...)
    N = length(r.op)
    dxback = falses(N)
    if !isempty(dx)                                     # TODO: what if this changes after net init?
        lastinput = 0
        for n=1:N
            isa(r.op[n], Input) || continue
            dx[lastinput += 1] == nothing || (dxback[n] = true)
        end
    end
    for n=N:-1:1
        nsparse = difsparse(r, dy, n)
        r.toback[n] || dxback[n] || continue
        # r.dif0[n] = finddif(r, n, nsparse)  # TODO-OPTIMIZATION
        r.dif0[n] = newarray(gpu(), nsparse, eltype(r.out0[n]), size(r.out0[n]))
        if r.toincr[n]    # TODO: what if this changes after init?
            # r.tmp[n] = findtmp(r, n, nsparse) # TODO-OPTIMIZATION
            r.tmp[n] = newarray(gpu(), nsparse, eltype(r.out0[n]), size(r.out0[n]))
        end
    end
end

function difsparse(r::Net, dy, n)                       # TODO: test this, compare with old initback
    r.toincr[n] && return false
    n == length(r.op) && return issparse(dy)
    for i=1:length(r.op)
        if (isa(r.op[i], Dot) &&
            n == r.inputs[i][1] &&
            issparse(r.out0[r.inputs[i][2]]))
            return true
        end
    end
    return false
end

### DEAD CODE

# # instead of changing the ops, record the eltype in net and use it when generating stuff.
#     # for o in r.op
#     #     t = eltype(o)
#     #     if t == nothing
#     #         continue
#     #     elseif it == nothing
#     #         it = t
#     #     elseif t != it
#     #         seteltype!(o, it)
#     #     end
#     # end

# #####

# function netinit(r::Net)
#     initout(r, inputs...)
#     initdif(r, inputs...)
#     inittmp(r, inputs...)
# end

# """

# initout(r::Net, inputs...) initializes r.out[n] which points to the
# array holding the output of op[n].  r.out0[n] holds the actual array,
# allowing r.out[n] point to stack elements when necessary during back
# calculation.  We optimize memory use by sharing arrays when we can.

# """
# function initout(r::Net)
#     N = length(r.op)
#     index = zeros(Int,N)          			# index==0 represents need for new register
#     for n=1:N
#         overwrites(r.op[n]) || continue                 # see if we can reuse our input for output
#         r.tosave[n] && continue                         # a saved register should only be written by op[n]
#         i = r.inputs[n][1]                              # see if we can overwrite the first input
#         r.tosave[i] && continue                         # do not overwrite if you are going to save for back
#         ow = true; k = n                                # TODO: checkout other inputs if this one fails
#         while true                                      # see if anybody else uses i before its next update
#             k = mod1(k+1, N)                            # TODO: This is suboptimal, gives 13 regs for LSTM
#             in(i, r.inputs[k]) && (ow = false; break)   # Should look for existing regs no longer used
#             k == i && break                             # this is nontrivial because the sizes are unknown
#             k == N && i > N && break                    # but can possibly be inferred
#         end                                             # fencepost check: N==1:OK, i==n:OK, i>N:OK
#         ow && (index[n]=i)
#     end
#     r.out0 = cell(N)
#     for n=1:N                                           # index==0 ops get their own registers
#         index[n]==0 || continue
#         r.out0[n] = Any[]
#     end
#     for n=1:N                                             # other ops will overwrite existing registers
#         index[n]==0 && continue
#         k = index[n]
#         while index[k]!=0; k=index[k]; end
#         r.out0[n] = r.out0[k]
#     end
#     r.out = fill!(cell(N), nothing)                       # use nothing to represent identity element: 0 for add, 1 for mul
# end

# """
# initdif(r::Net) initializes r.dif[n], which is the loss gradient of the last output of op[n].
# r.dif0[n] holds the actual array, r.dif[n] can point to this or be 'nothing' representing the zero matrix.
# """
# function initdif(r::Net)
#     # TODO: this is not optimal for add: no need for two copies of dy when neither input is overwriting.
#     # TODO: tmp register to replace dif1 and serve other needs

#     index = zeros(Int,nops(r))                                  # index==0 represents need for new register
#     for n=1:nops(r)                                             # find out if back(op[n]) can overwrite dif[n]
#         overwrites(r.op[n]) || continue                         # overwrites means potentially both forw x<-y and back dy<-dx
#         r.toincr[n] && continue                                  # don't share dif[n] with multi-output
#         n == nops(r) && continue                                # don't share last output, dy may come in sparse
#         for i in r.inputs[n]
#             r.toincr[i] && continue                              # don't share dif[n] with multi-output
#             index[i]==0 || error("should not happen")
#             index[i]=n; break                                   # op[n] will overwrite dif[n] to get dif[i]
#         end
#     end
#     r.dif0 = cell(nops(r))
#     for n=1:nops(r)
#         index[n]==0 || continue
#         r.dif0[n] = Any[]
#     end
#     for n=1:nops(r)
#         index[n]==0 && continue
#         k = index[n]
#         while index[k]!=0; k=index[k]; end
#         r.dif0[n] = r.dif0[k]
#     end
#     r.dif = fill!(cell(nops(r)), nothing)
# end

# """
# inittmp(r::Net) initializes r.tmp[n] which is an extra array for:
# (1) incremental updates if r.toincr[n].
# (2) for forw/back ops that need extra space (TODO)
# """
# function inittmp(r::Net)
#     # TODO: this is inefficient, could use a single tmp for each size.
#     r.tmp = cell(nops(r))
#     for n=1:length(r.tmp)
#         r.tmp[n] = r.toincr[n] ? Any[] : nothing
#     end
# end

# ### DEAD CODE:

#     # initop(r, a...)
#     # initinputs(r, a...)
#     # @assert length(r.op)==length(r.inputs)
#     # initninputs(r)
#     # initparams(r)
# # function initstack(r::Net)
# #     r.stack = Any[]
# #     r.sp = 0
# # end

# #    initstack(r)

# # r.op[n] is the n'th operation in the net
# # The user specifies the operations in the MLP constructor arguments
# # If an argument is another MLP, its operations are spliced in

# # function initop(r::Net, a...)
# #     r.op = Op[]
# #     for ai in a
# #         isa(ai,Tuple) && (ai=ai[1])
# #         isa(ai,Op) ?  push!(r.op, ai) :
# #         isa(ai,Net)   ?  append!(r.op, ai.op) :
# #         error("Bad op: $ai")
# #     end
# # end

# # r.inputs[n] is an array of k indices for the inputs of the n'th op.
# # k is typically 1 but can be 2 for Add2 and Mul2 operations
# # index i<=nops(r) indicates output of r.op[i], i>nops(r) indicates network inputs
# # By default an op takes the results of the previous k op outputs
# # (or network inputs for initial ops)
# # The user can override this by providing a tuple argument for an op
# # The first element of the tuple is the op, the rest are user indices for inputs
# # userindex j>0 indicates output of userarg[j], j<=0 indicates network input 1-j.

# # function initinputs(r::Net, a...)
# #     newindex = Array(Int, length(a))
# #     lastindex = 0
# #     for i=1:length(a)
# #         ai = isa(a[i],Tuple) ? a[i][1] : a[i]
# #         lastindex += (isa(ai,Op) ? 1 : length(ai.op))
# #         newindex[i] = lastindex
# #     end
# #     r.inputs = Any[]
# #     for i=1:length(a)
# #         ai = isa(a[i],Tuple) ? a[i][1] : a[i]
# #         bi = isa(a[i],Tuple) ? a[i][2:end] : ((i-ninputs(ai)):(i-1))
# #         length(bi) == ninputs(ai) || error("Wrong number of inputs for $i:$(typeof(ai))")
# #         if isa(ai, Op)
# #             push!(r.inputs, map(j->(j>0 ? newindex[j] : nops(r)+1-j), Int[bi...]))
# #         else
# #             j0 = length(r.inputs)
# #             for aii in ai.inputs
# #                 push!(r.inputs, map(j->(j<=nops(ai) ? j+j0 : (j=bi[j-nops(ai)]; j>0 ? newindex[j] : nops(r)+1-j)), aii))
# #             end
# #         end
# #     end
# # end

# # r.ninputs is the number of inputs the whole MLP expects
# # indices i>nops(r) in r.inputs refer to network inputs

# # function initninputs(r::Net)
# #     n = 0
# #     for ai in r.inputs
# #         for aij in ai
# #             aij - nops(r) > n && (n = aij - nops(r))
# #         end
# #     end
# #     r.ninputs = n
# # end

# # r.params points to all op parameters

# # function initparams(r::Net)
# #     r.params = Any[]
# #     for o in r.op
# #         if isa(o,par) && isupdated(o)
# #             push!(r.params, o)
# #         end
# #     end
# # end

# # isupdated(o)=true # TODO: should be false for rand and const

#     #     @assert !haskey(dict, target)			# each target has to be unique
#     #     for x in args
#     #         @assert haskey(dict, x)                     #
#     #     end
#     #     if func == :input && inputs != nothing          # subroutines take their input registers from caller
#     #         @assert length(inputs) >= nextinput
#     #         dict[target] = inputs[nextinput]
#     #         nextinput += 1
#     #     end
#     #     op = eval(Expr(:call, func, pars...))
#     # end
#     #

# # initforw(r::Net,x...) is called at the beginning of a sequence or
# # before processing a stand-alone item.  It is not called between
# # elements of a sequence.  It allocates and/or resizes r.out0.

# # function initforw2(r::Net, inputs...; keepstate=false, a...)
# #     r.dbg && display((:initforw0,keepstate,vecnorm0(r.out),vecnorm0(r.stack[1:r.sp])))
# #     r.sp == 0 || error("Stack corruption")
# #     length(inputs) == ninputs(r) || error("Wrong number of inputs")
# #     out = fill!(cell(length(r.out0)), nothing)          # TODO: (minor) can we get rid of alloc
# #     lastinput = 0
# #     while findfirst(out,nothing) > 0                            # allow multiple passes for size inference
# #         nnothing = count(x->(x==nothing), out)
# #         for n = 1:nops(r)
# #             if out[n] != nothing
# #                 continue
# #             elseif isa(r.op[n], Input)
# #                 out[n] = initarray(r.out0, n, inputs[lastinput += 1])
# #             elseif isempty(r.inputs[n])                         # par, rnd, con
# #                 s = psize(r, n)                                 # TODO: write psize inference
# #                 s == nothing && continue                        # may happen first pass
# #                 out[n] = initarray(r.out0, n, inputs[1], s; dense=true)
# #             else
# #                 s = ysize(r.op[n], out[r.inputs[n]]...)
# #                 s == nothing && continue # may happen with recurrent connections
# #                 out[n] = initarray(r.out0, n, out[r.inputs[n]][1], s; dense=true)
# #             end
# #         end
# #         nnothing == count(x->(x==nothing), out) && error("Cannot determine size of array $(findfirst(out,nothing))")
# #     end
# #     # We recover or reset r.out:
# #     keepstate ? copy!(r.out, r.out0) : fill!(r.out, nothing) # TODO: r.out[i] = arr(r.out0[i]) elsewhere..
# #     r.dbg && display((:initforw1,keepstate,vecnorm0(r.out),vecnorm0(r.stack[1:r.sp])))
# # end

# # function initarray(a, i, x, dims=size(x); dense=false)
# #     (display((:initarray0,summary(a[i]),i,summary(x),dims,dense,summary(a)));println())
# #     if isempty(a[i])
# #         oldai = a[i]
# #         if !dense && issparse(x)
# #             a[i] = spzeros(eltype(x), dims...)
# #             gpu() && (a[i] = CudaSparseMatrixCSC(a[i]))
# #         else
# #             a[i] = fill!(KUdense(gpu()?CudaArray:Array, eltype(x), dims), 0)
# #         end
# #         for j=1:length(a); a[j]===oldai && (a[j]=a[i]); end # preserve array sharing
# #     elseif eltype(a[i]) != eltype(x)
# #         error("Element type mismatch")
# #     elseif size(a[i]) != dims
# #         warn("Resizing $(size(a[i]))->$dims")
# #         fill!(resize!(a[i], dims), 0) # TODO: this does not work for sparse
# #     end
# #     (display((:initarray1,summary(a[i]),i,summary(x),dims,dense,summary(a)));println())
# #     return a[i]
# # end

#             # elseif in(nothing, dims[r.inputs[n]])
#             #     continue
#             # else
#             # end

#             # if isempty(r.inputs[n])                     # happens with par, rnd, con, input
#             #     elseif prod(size(r.op[n])) > 0
#             #         dims[n] = size(r.op[n])
#             #     end
#             # end
#             # if isa(op, Input)
#             #     dims[n] = size(inputs[lastinput += 1])
#             # elseif isa(op, Par) && prod(op.dims) > 0
#             #     dims[n] = op.dims
#             # elseif isa(op, Dot)
#             #     (i1,i2) = r.inputs[n]
#             #     (d1,d2) = dims[i1,i2]
#             #     if d1 == nothing || d2 == nothing
#             #         continue
#             #     end
#             #     m1 = prod(d1[1:end-1])                   # treat [d...,dn] as [prod(d...), dn] matrix
#             #     n1 = d1[end]
#             #     m2 = prod(d2[1:end-1])
#             #     n2 = d2[end]
#             #     if m1 == 0 || n2 == 0
#             #         continue
#             #     end
#             #     if n1 == m2 == 0
#             #         continue
#             #     elseif n1 == 0
#             #         setdims!(r.op[i1], (d1[1:end-1]..., m2))
#             #     elseif m2 == 0
#             #         length(d2) == 2 || error("Cannot infer size of tensor")
#             #         setdims!(r.op[i2], (n1, d2[2]))
#             #     elseif n1 == m2
#             #         dims[n] = (m1,n2)
#             #     else
#             #         error("Dot size mismatch")
#             #     end
#             # end
#             # # input, par, dot, add, mul, actf, loss, conv, pool

# function initback0(r::Net, dy; a...)
#     N = length(r.op)
#     firstwrite = zeros(Int, N)                          # dif[n] is read at t=n, find out when it is first written.
#     # TODO: dif[N] is special!
#     for n=N:-1:1
#         r.toback[n] || continue
#         k = n
#         while true
#             in(n, r.inputs[k]) && (firstwrite[n] = k; break)
#             k = mod1(k-1, N)
#         end
#     end
#     return firstwrite
#     for n=N:-1:1
#         nsparse = false # TODO
#         r.toback[n] || continue
#         r.dif[n] = finddif(r, n, nsparse)
#         r.toincr[n] && (r.tmp[n] = findtmp(r, n, nsparse))
#     end
# end

# function finddif(r::Net, n, nsparse)
#     isa(r.op[n], Par) && return                         # par difs should persist across iterations
#     N = length(r.op)
#     k = n; nw1 = 0                                      # dif[n] is read at t=n, find out when it is first written.
#     while true
#         in(n, r.inputs[k]) && (nw1 = k; break)          # it could potentially be written at t=n right after reading
#         k = mod1(k-1, N)
#     end
#     @assert nw1 > 0                                     # dif[n] will be busy from nw1 to n going back cyclical
#     free = nothing                                      # we need dif[i] that is free during this period
#     for i = n+1:N                                       # note that dif[i] may already be shared to some dif[j]
#         isa(r.op[i], Par) && return
#         size(r.dif[i]) == size(r.out0[n]) || continue
#         issparse(r.dif[i]) == nsparse || continue

#     end
#     return free
# end

# initback(r::Net, dy) called at the beginning of a sequence or a
# stand-alone item, never between elements of a sequence.

# function initback1(r::Net, dy; seq=false, a...)
#     fill!(r.dif, nothing)                           # why? (TODO)
#     for n=1:length(r.dif0)
#         y = (n==nops(r) && dy!=nothing ? dy : r.out0[n])
#         initarray(r.dif0, n, y; dense=true) # x and dw may be sparse, dx and w always dense
#         if r.toincr[n]
#             initarray(r.tmp, n, r.dif0[n])
#         end
#     end
#     for n=1:length(r.dif0)
#         r.toincr[n] && fill!(r.dif0[n], 0)           # zeroed by back at every item in sequence
#     end
#     if seq
#         for w in params(r)       # zeroed only once at the beginning of the sequence
#             isdefined(w,:diff) && isdefined(w,:inc) && fill!(w.diff,0)
#         end
#     end
# end


text:
using GLPlot, GLVisualize, GLAbstraction, Colors, GeometryTypes, Plots, FileIO
using Reactive, GLWindow, NIfTI
import Quaternions: qrotation
# load a volume
vol = niread(joinpath(homedir(), "Desktop", "lesson", "volume.nii")).raw;
vol = Float32.(vol ./ maximum(vol));

# prepare the slices
ranges = ntuple(i-> linspace(-2f0, 2f0, size(vol, i)), 3)
slices = ntuple(3) do i
    range_s = play_widget(1:size(vol, i))
    slice = map(range_s) do j
        s = vol[ifelse(i == 1, j, :), ifelse(i == 2, j, :) , ifelse(i == 3, j, :)]
        reinterpret(Intensity{1, Float32}, s)
    end
    slice, range_s
end

rotations = (
    Mat4f0(qrotation(unit(Vec3f0, 1), Float32(pi/2))),
    rotationmatrix_x(Float32(pi/2)) * Mat4f0(qrotation(unit(Vec3f0, 2), Float32(pi/2))),
    rotationmatrix_z(Float32(pi/2))
)

plots = map(1:3) do i
    rot = rotations[i]
    slice, slider = slices[i]
    axis = unit(Vec3f0, 1)
    model = map(slider) do r
        trans = unit(Vec3f0, mod1(3-i, 3)) * ranges[i][r]
        translationmatrix(trans) * rot
    end
    rs = ranges[mod1(i + 1, 3)], ranges[mod1(i + 2, 3)]
    glplot(
        slice,
        model = model,
        ranges = rs,
        preferred_camera = :perspective,
        stroke_width = 0.0f0,
        stroke_color = RGBA{Float32}(0,0,0,0.4),
        color_norm = Vec2f0(0, 1)
    )
end

# low level hack to share attributes... Gotta figure out a nice API for this
attributes = GLPlot.to_edit_dict(plots[1].children[])
for p in plots[2:end]
    for (k, v) in attributes
        k == :intensity && continue # don't share the intensity values
        p.children[].uniforms[k] = v
    end
end


text:
module GLLists

using ModernGL

export GL_QUERY_BUFFER
export GL_SHADER_BINARY_FORMAT_SPIR_V_ARB

export DEFAULT_SHADER_VERSION
export DEFAULT_SHADER_CODE

export CLIP_DISTANCE
export LIST_STATUS
export LIST_BUFFER
export LIST_OPTIONS
export LIST_COLOR_MASK
export LIST_COMPARE
export LIST_INFO
export LIST_SHADER
export LIST_DRAW_MODE
export LIST_ERROR
export LIST_TYPE
export LIST_TYPE_ELEMENTS
export LIST_TYPE_STRING

# missing constants
const GL_QUERY_BUFFER = convert(GLenum, 0x9192)
const GL_SHADER_BINARY_FORMAT_SPIR_V_ARB = convert(GLenum, 0x9551)

""" TODO """
const LIST_STATUS = Dict{Symbol,Dict{Symbol,Function}}(
	:SHADER		=> Dict{Symbol,Function}(:STATE => glGetShaderiv, :INFO => glGetShaderInfoLog),
	:PROGRAM	=> Dict{Symbol,Function}(:STATE => glGetProgramiv, :INFO => glGetProgramInfoLog),
)

const LIST_BUFFER = Dict{Symbol,GLenum}(
	:ARRAY_BUFFER								=> GL_ARRAY_BUFFER,
	:ATOMIC_COUNTER_BUFFER			=> GL_ATOMIC_COUNTER_BUFFER,
	:COPY_READ_BUFFER						=> GL_COPY_READ_BUFFER,
	:COPY_WRITE_BUFFER					=> GL_COPY_WRITE_BUFFER,
	:DISPATCH_INDIRECT_BUFFER		=> GL_DISPATCH_INDIRECT_BUFFER,
	:DRAW_INDIRECT_BUFFER				=> GL_DRAW_INDIRECT_BUFFER,
	:ELEMENT_ARRAY_BUFFER				=> GL_ELEMENT_ARRAY_BUFFER,
	:PIXEL_PACK_BUFFER					=> GL_PIXEL_PACK_BUFFER,
	:PIXEL_UNPACK_BUFFER				=> GL_PIXEL_UNPACK_BUFFER,
	:QUERY_BUFFER								=> GL_QUERY_BUFFER,
	:SHADER_STORAGE_BUFFER			=> GL_SHADER_STORAGE_BUFFER,
	:TEXTURE_BUFFER							=> GL_TEXTURE_BUFFER,
	:TRANSFORM_FEEDBACK_BUFFER	=> GL_TRANSFORM_FEEDBACK_BUFFER,
	:UNIFORM_BUFFER							=> GL_UNIFORM_BUFFER,
)

const CLIP_DISTANCE = [
	GL_CLIP_DISTANCE0,
	GL_CLIP_DISTANCE1,
	GL_CLIP_DISTANCE2,
	GL_CLIP_DISTANCE3,
	GL_CLIP_DISTANCE4,
	GL_CLIP_DISTANCE5,
	GL_CLIP_DISTANCE6,
	GL_CLIP_DISTANCE7,
	#GL_CLIP_DISTANCE8,
]

#	GL_CLIP_ORIGIN,
#	GL_CLIP_DEPTH_MODE
#	GL_CLIPPING_INPUT_PRIMITIVES_ARB
#	GL_CLIPPING_OUTPUT_PRIMITIVES_ARB

const LIST_OPTIONS = Dict{Symbol,GLuint}(
	:ALPHA_TEST                     => 0,                                # GL_ALPHA_TEST
	:AUTO_NORMAL                    => GL_AUTO_GENERATE_MIPMAP,          # GL_AUTO_NORMAL
	:BLEND                          => GL_BLEND,                         # GL_BLEND
	:CLIP_PLANEi                    => 0,                                # GL_CLIP_PLANEi
	:CLIP_DISTANCEi                 => 0,                                # GL_CLIP_DISTANCEi
	:COLOR_LOGIC_OP                 => GL_COLOR_LOGIC_OP,                # GL_COLOR_LOGIC_OP
	:COLOR_MATERIAL                 => 0,                                # GL_COLOR_MATERIAL
	:COLOR_SUM                      => 0,                                # GL_COLOR_SUM
	:COLOR_TABLE                    => 0,                                # GL_COLOR_TABLE
	:CONVOLUTION_1D                 => 0,                                # GL_CONVOLUTION_1D
	:CONVOLUTION_2D                 => 0,                                # GL_CONVOLUTION_2D
	:CULL_FACE                      => GL_CULL_FACE,                     # GL_CULL_FACE
	:DEBUG_OUTPUT                   => GL_DEBUG_OUTPUT,                  # GL_DEBUG_OUTPUT
	:DEBUG_OUTPUT_SYNCHRONOUS       => GL_DEBUG_OUTPUT_SYNCHRONOUS,      # GL_DEBUG_OUTPUT_SYNCHRONOUS
	:DEPTH_CLAMP                    => GL_DEPTH_CLAMP,                   # GL_DEPTH_CLAMP
	:DEPTH_TEST                     => GL_DEPTH_TEST,                    # GL_DEPTH_TEST
	:DITHER                         => GL_DITHER,                        # GL_DITHER
	:FOG                            => 0,                                # GL_FOG
	:FRAMEBUFFER_SRGB               => GL_FRAMEBUFFER_SRGB,              # GL_FRAMEBUFFER_SRGB
	:HISTOGRAM                      => 0,                                # GL_HISTOGRAM
	:INDEX_LOGIC_OP                 => 0,                                # GL_INDEX_LOGIC_OP
	:LIGHTi                         => 0,                                # GL_LIGHTi
	:LIGHTING                       => 0,                                # GL_LIGHTING
	:LINE_SMOOTH                    => GL_LINE_SMOOTH,                   # GL_LINE_SMOOTH
	:LINE_STIPPLE                   => 0,                                # GL_LINE_STIPPLE
	:MAP1_COLOR_4                   => 0,                                # GL_MAP1_COLOR_4
	:MAP1_INDEX                     => 0,                                # GL_MAP1_INDEX
	:MAP1_NORMAL                    => 0,                                # GL_MAP1_NORMAL
	:MAP1_TEXTURE_COORD_1           => 0,                                # GL_MAP1_TEXTURE_COORD_1
	:MAP1_TEXTURE_COORD_2           => 0,                                # GL_MAP1_TEXTURE_COORD_2
	:MAP1_TEXTURE_COORD_3           => 0,                                # GL_MAP1_TEXTURE_COORD_3
	:MAP1_TEXTURE_COORD_4           => 0,                                # GL_MAP1_TEXTURE_COORD_4
	:MAP1_VERTEX_3                  => 0,                                # GL_MAP1_VERTEX_3
	:MAP1_VERTEX_4                  => 0,                                # GL_MAP1_VERTEX_4
	:MAP2_COLOR_4                   => 0,                                # GL_MAP2_COLOR_4
	:MAP2_INDEX                     => 0,                                # GL_MAP2_INDEX
	:MAP2_NORMAL                    => 0,                                # GL_MAP2_NORMAL
	:MAP2_TEXTURE_COORD_1           => 0,                                # GL_MAP2_TEXTURE_COORD_1
	:MAP2_TEXTURE_COORD_2           => 0,                                # GL_MAP2_TEXTURE_COORD_2
	:MAP2_TEXTURE_COORD_3           => 0,                                # GL_MAP2_TEXTURE_COORD_3
	:MAP2_TEXTURE_COORD_4           => 0,                                # GL_MAP2_TEXTURE_COORD_4
	:MAP2_VERTEX_3                  => 0,                                # GL_MAP2_VERTEX_3
	:MAP2_VERTEX_4                  => 0,                                # GL_MAP2_VERTEX_4
	:MINMAX                         => 0,                                # GL_MINMAX
	:MULTISAMPLE                    => GL_MULTISAMPLE,                   # GL_MULTISAMPLE
	:NORMALIZE                      => 0,                                # GL_NORMALIZE
	:POINT_SMOOTH                   => 0,                                # GL_POINT_SMOOTH
	:POINT_SPRITE                   => 0,                                # GL_POINT_SPRITE
	:POLYGON_OFFSET_FILL            => GL_POLYGON_OFFSET_FILL,           # GL_POLYGON_OFFSET_FILL
	:POLYGON_OFFSET_LINE            => GL_POLYGON_OFFSET_LINE,           # GL_POLYGON_OFFSET_LINE
	:POLYGON_OFFSET_POINT           => GL_POLYGON_OFFSET_POINT,          # GL_POLYGON_OFFSET_POINT
	:POLYGON_SMOOTH                 => GL_POLYGON_SMOOTH,                # GL_POLYGON_SMOOTH
	:POLYGON_STIPPLE                => 0,                                # GL_POLYGON_STIPPLE
	:POST_COLOR_MATRIX_COLOR_TABLE  => 0,                                # GL_POST_COLOR_MATRIX_COLOR_TABLE
	:POST_CONVOLUTION_COLOR_TABLE   => 0,                                # GL_POST_CONVOLUTION_COLOR_TABLE
	:PRIMITIVE_RESTART              => GL_PRIMITIVE_RESTART,             # GL_PRIMITIVE_RESTART
	:PRIMITIVE_RESTART_FIXED_INDEX  => GL_PRIMITIVE_RESTART_FIXED_INDEX, # GL_PRIMITIVE_RESTART_FIXED_INDEX
	:PROGRAM_POINT_SIZE             => GL_PROGRAM_POINT_SIZE,            # GL_PROGRAM_POINT_SIZE
	:RASTERIZER_DISCARD             => GL_RASTERIZER_DISCARD,            # GL_RASTERIZER_DISCARD
	:RESCALE_NORMAL                 => 0,                                # GL_RESCALE_NORMAL
	:SAMPLE_ALPHA_TO_COVERAGE       => GL_SAMPLE_ALPHA_TO_COVERAGE,      # GL_SAMPLE_ALPHA_TO_COVERAGE
	:SAMPLE_ALPHA_TO_ONE            => GL_SAMPLE_ALPHA_TO_ONE,           # GL_SAMPLE_ALPHA_TO_ONE
	:SAMPLE_COVERAGE                => GL_SAMPLE_COVERAGE,               # GL_SAMPLE_COVERAGE
	:SAMPLE_SHADING                 => GL_SAMPLE_SHADING,                # GL_SAMPLE_SHADING
	:SAMPLE_MASK                    => GL_SAMPLE_MASK,                   # GL_SAMPLE_MASK
	:SEPARABLE_2D                   => 0,                                # GL_SEPARABLE_2D
	:SCISSOR_TEST                   => GL_SCISSOR_TEST,                  # GL_SCISSOR_TEST
	:STENCIL_TEST                   => GL_STENCIL_TEST,                  # GL_STENCIL_TEST
	:TEXTURE_1D                     => GL_TEXTURE_1D,                    # GL_TEXTURE_1D
	:TEXTURE_2D                     => GL_TEXTURE_2D,                    # GL_TEXTURE_2D
	:TEXTURE_3D                     => GL_TEXTURE_3D,                    # GL_TEXTURE_3D
	:TEXTURE_CUBE_MAP               => GL_TEXTURE_CUBE_MAP,              # GL_TEXTURE_CUBE_MAP
	:TEXTURE_CUBE_MAP_SEAMLESS      => GL_TEXTURE_CUBE_MAP_SEAMLESS,     # GL_TEXTURE_CUBE_MAP_SEAMLESS
	:TEXTURE_GEN_Q                  => 0,                                # GL_TEXTURE_GEN_Q
	:TEXTURE_GEN_R                  => 0,                                # GL_TEXTURE_GEN_R
	:TEXTURE_GEN_S                  => 0,                                # GL_TEXTURE_GEN_S
	:TEXTURE_GEN_T                  => 0,                                # GL_TEXTURE_GEN_T
	:VERTEX_PROGRAM_POINT_SIZE      => GL_VERTEX_PROGRAM_POINT_SIZE,     # GL_VERTEX_PROGRAM_POINT_SIZE
	:VERTEX_PROGRAM_TWO_SIDE        => 0,                                # GL_VERTEX_PROGRAM_TWO_SIDE
)

const LIST_COLOR_MASK = Dict{Symbol,GLuint}(
	:SRC_COLOR                 => GL_SRC_COLOR,
	:ONE_MINUS_SRC_COLOR       => GL_ONE_MINUS_SRC_COLOR,
	:DST_COLOR                 => GL_DST_COLOR,
	:ONE_MINUS_DST_COLOR       => GL_ONE_MINUS_DST_COLOR,
	:SRC_ALPHA                 => GL_SRC_ALPHA,
	:ONE_MINUS_SRC_ALPHA       => GL_ONE_MINUS_SRC_ALPHA,
	:DST_ALPHA                 => GL_DST_ALPHA,
	:ONE_MINUS_DST_ALPHA       => GL_ONE_MINUS_DST_ALPHA,
	:CONSTANT_COLOR            => GL_CONSTANT_COLOR,
	:ONE_MINUS_CONSTANT_COLOR  => GL_ONE_MINUS_CONSTANT_COLOR,
	:CONSTANT_ALPHA            => GL_CONSTANT_ALPHA,
	:ONE_MINUS_CONSTANT_ALPHA  => GL_ONE_MINUS_CONSTANT_ALPHA,
	:SRC_ALPHA_SATURATE        => GL_SRC_ALPHA_SATURATE,
	:SRC1_COLOR                => GL_SRC1_COLOR,
	:ONE_MINUS_SRC1_COLOR      => GL_ONE_MINUS_SRC1_COLOR,
	:SRC1_ALPHA                => GL_SRC1_ALPHA,
	:ONE_MINUS_SRC1_ALPHA      => GL_ONE_MINUS_SRC1_ALPHA,
)

const LIST_COMPARE = Dict{Symbol,GLuint}(
	:NEVER     => GL_NEVER,
	:LESS      => GL_LESS,
	:EQUAL     => GL_EQUAL,
	:LEQUAL    => GL_LEQUAL,
	:GREATER   => GL_GREATER,
	:NOTEQUAL  => GL_NOTEQUAL,
	:GEQUAL    => GL_GEQUAL,
	:ALWAYS    => GL_ALWAYS,
)

const LIST_INFO = Dict{Symbol,GLuint}(
	:RENDERER =>                 GL_RENDERER,
	:VENDOR =>                   GL_VENDOR,
	:VERSION =>                  GL_VERSION,
	:SHADING_LANGUAGE_VERSION => GL_SHADING_LANGUAGE_VERSION,
	:EXTENSIONS =>               GL_EXTENSIONS,
)

const DEFAULT_SHADER_VERSION = 130

function DEFAULT_SHADER_CODE(typ::Symbol)
	code = LIST_SHADER[typ][:CODE]
	string("#version ",DEFAULT_SHADER_VERSION,"\n",code[1],"\n",code[2],"\nvoid main(){\n",code[3],"\n}")
end

const LIST_SHADER = Dict(
	:VERTEX           => Dict(:TYPE => GLuint(GL_VERTEX_SHADER),					:EXT => "vert", :CODE => ("in vec3 iVertex; uniform mat4 iMVP = mat4(vec4(1,0,0,0),vec4(0,1,0,0),vec4(0,0,1,0),vec4(0,0,0,1));", "", "gl_Position = iMVP * vec4(iVertex,1);")),
	:FRAGMENT         => Dict(:TYPE => GLuint(GL_FRAGMENT_SHADER),				:EXT => "frag", :CODE => ("", "out vec4 oFragColor;", "oFragColor=vec4(vec3(0.5),1.0);")),
	:TESS_CONTROL			=> Dict(:TYPE => GLuint(GL_TESS_CONTROL_SHADER),		:EXT => "tesc", :CODE => ("uniform uint iTessLevel = 1;", "layout (vertices = 3) out;", "#define ID gl_InvocationID\n if (ID == 0) for(int i=0; i<3; ++i){	if(i<2) gl_TessLevelInner[i] = iTessLevel; gl_TessLevelOuter[i] = iTessLevel; }; gl_out[ID].gl_Position = gl_in[ID].gl_Position;")),
	:TESS_EVALUATION	=> Dict(:TYPE => GLuint(GL_TESS_EVALUATION_SHADER),	:EXT => "tese", :CODE => ("layout (triangles, equal_spacing, ccw) in;", "", "gl_Position = vec4(gl_in[0].gl_Position.xyz * gl_TessCoord.x + gl_in[1].gl_Position.xyz * gl_TessCoord.y + gl_in[2].gl_Position.xyz * gl_TessCoord.z,1.0);")),
	:GEOMETRY         => Dict(:TYPE => GLuint(GL_GEOMETRY_SHADER),				:EXT => "geom",	:CODE => ("layout(triangles) in;", "layout(triangle_strip, max_vertices=3) out;", "for(int i=0; i<gl_in.length(); ++i){ gl_Position = gl_in[i].gl_Position; EmitVertex(); }; EndPrimitive();")),
	:COMPUTE          => Dict(:TYPE => GLuint(GL_COMPUTE_SHADER),					:EXT => "comp",	:CODE => ("", "", "")),
)

const LIST_DRAW_MODE = Dict{Symbol,GLuint}(
	:POINTS										=> GL_POINTS,										#0x0000
	:LINES										=> GL_LINE,											#0x0001
	:LINE_LOOP								=> GL_LINE_LOOP,								#0x0002
	:LINE_STRIP     					=> GL_LINE_STRIP,								#0x0003
	:TRIANGLES 								=> GL_TRIANGLES,								#0x0004
	:TRIANGLE_STRIP  					=> GL_TRIANGLE_STRIP,						#0x0005
	:TRIANGLE_FAN							=> GL_TRIANGLE_FAN,							#0x0006
	:QUADS										=> GL_QUADS,										#0x0007
	:QUAD_STRIP								=> GL_QUAD_STRIP,								#0x0008 (deprecated since OpenGL3)
	:POLYGON									=> GL_POLYGON,									#0x0009 (deprecated since OpenGL3)
	:LINES_ADJACENCY					=> GL_LINES_ADJACENCY,					#0x000A
	:LINE_STRIP_ADJACENCY			=> GL_LINE_STRIP_ADJACENCY,			#0x000B
	:TRIANGLES_ADJACENCY			=> GL_TRIANGLES_ADJACENCY,			#0x000C
	:TRIANGLE_STRIP_ADJACENCY	=> GL_TRIANGLE_STRIP_ADJACENCY,	#0x000D
	:PATCHES 									=> GL_PATCHES,									#0x000E
)

const LIST_ERROR = Dict{GLuint,AbstractString}(
	GL_NO_ERROR											 =>	"",
	GL_INVALID_ENUM									 =>	"INVALID_ENUM: An unacceptable value is specified for an enumerated argument. The offending command is ignored and has no other side effect than to set the error flag.",
	GL_INVALID_VALUE								 =>	"INVALID_VALUE: A numeric argument is out of range. The offending command is ignored and has no other side effect than to set the error flag.",
	GL_INVALID_OPERATION						 =>	"INVALID_OPERATION: The specified operation is not allowed in the current state. The offending command is ignored and has no other side effect than to set the error flag.",
	GL_INVALID_FRAMEBUFFER_OPERATION => "INVALID_FRAMEBUFFER_OPERATION: The framebuffer object is not complete. The offending command is ignored and has no other side effect than to set the error flag.",
	GL_OUT_OF_MEMORY								 =>	"OUT_OF_MEMORY: There is not enough memory left to execute the command. The state of the GL is undefined, except for the state of the error flags, after this error is recorded.",
	GL_STACK_UNDERFLOW							 =>	"STACK_UNDERFLOW",
	GL_STACK_OVERFLOW								 =>	"STACK_OVERFLOW",
)

const LIST_TYPE = Dict{DataType,GLuint}(
	Bool => GL_BOOL,
	Int8 => GL_BYTE,
	Int16 => GL_SHORT,
	Int32 => GL_INT, #gl_low_int, gl_medium_int, gl_high_int
	Int64 => GL_INT, # no 64 bit support
	Int128 => GL_INT, # no 128 bit support
	UInt8 => GL_UNSIGNED_BYTE,
	UInt16 => GL_UNSIGNED_SHORT,
	UInt32 => GL_UNSIGNED_INT,
	UInt64 => GL_UNSIGNED_INT, # no 64 bit support
	UInt128 => GL_UNSIGNED_INT, # no 128 bit support
	Float16 => GL_HALF_FLOAT,
	Float32 => GL_FLOAT,
	Float64 => GL_DOUBLE,
)

const LIST_TYPE_ELEMENTS = Dict{GLuint,UInt32}(
	# Float
	GL_FLOAT =>      1,
	GL_FLOAT_VEC2 => 2,
	GL_FLOAT_VEC3 => 3,
	GL_FLOAT_VEC4 => 4,

	# Doubles
	GL_DOUBLE =>      1,
	GL_DOUBLE_VEC2 => 2,
	GL_DOUBLE_VEC3 => 3,
	GL_DOUBLE_VEC4 => 4,

	# Int
	GL_INT =>      1,
	GL_INT_VEC2 => 2,
	GL_INT_VEC3 => 3,
	GL_INT_VEC4 => 4,

	# Unsigned Int
	GL_UNSIGNED_INT =>      1,
	GL_UNSIGNED_INT_VEC2 => 2,
	GL_UNSIGNED_INT_VEC3 => 3,
	GL_UNSIGNED_INT_VEC4 => 4,

	# Bool
	GL_BOOL =>      1,
	GL_BOOL_VEC2 => 2,
	GL_BOOL_VEC3 => 3,
	GL_BOOL_VEC4 => 4,

	# Float Matrix
	GL_FLOAT_MAT2 =>   4,
	GL_FLOAT_MAT3 =>   9,
	GL_FLOAT_MAT4 =>   16,
	GL_FLOAT_MAT2x3 => 6,
	GL_FLOAT_MAT2x4 => 8,
	GL_FLOAT_MAT3x2 => 6,
	GL_FLOAT_MAT3x4 => 12,
	GL_FLOAT_MAT4x2 => 8,
	GL_FLOAT_MAT4x3 => 12,

	# Double Matrix
	GL_DOUBLE_MAT2 =>   4,
	GL_DOUBLE_MAT3 =>   9,
	GL_DOUBLE_MAT4 =>   16,
	GL_DOUBLE_MAT2x3 => 6,
	GL_DOUBLE_MAT2x4 => 8,
	GL_DOUBLE_MAT3x2 => 6,
	GL_DOUBLE_MAT3x4 => 12,
	GL_DOUBLE_MAT4x2 => 8,
	GL_DOUBLE_MAT4x3 => 12,
)

const LIST_TYPE_STRING = Dict{GLuint,AbstractString}(
	# Float
	GL_FLOAT =>      "float",
	GL_FLOAT_VEC2 => "vec2f",
	GL_FLOAT_VEC3 => "vec3f",
	GL_FLOAT_VEC4 => "vec4f",

	# Doubles
	GL_DOUBLE =>      "double",
	GL_DOUBLE_VEC2 => "vec2d",
	GL_DOUBLE_VEC3 => "vec3d",
	GL_DOUBLE_VEC4 => "vec4d",

	# Int
	GL_INT =>      "int",
	GL_INT_VEC2 => "vec2i",
	GL_INT_VEC3 => "vec3i",
	GL_INT_VEC4 => "vec4i",

	# Unsigned Int
	GL_UNSIGNED_INT =>      "uint",
	GL_UNSIGNED_INT_VEC2 => "vec2u",
	GL_UNSIGNED_INT_VEC3 => "vec3u",
	GL_UNSIGNED_INT_VEC4 => "vec4u",

	# Bool
	GL_BOOL =>      "bool",
	GL_BOOL_VEC2 => "vec2b",
	GL_BOOL_VEC3 => "vec3b",
	GL_BOOL_VEC4 => "vec4b",

	# Float Matrix
	GL_FLOAT_MAT2 =>   "mat2f",
	GL_FLOAT_MAT3 =>   "mat3f",
	GL_FLOAT_MAT4 =>   "mat4f",
	GL_FLOAT_MAT2x3 => "mat2x3f",
	GL_FLOAT_MAT2x4 => "mat2x4f",
	GL_FLOAT_MAT3x2 => "mat3x2f",
	GL_FLOAT_MAT3x4 => "mat3x4f",
	GL_FLOAT_MAT4x2 => "mat4x2f",
	GL_FLOAT_MAT4x3 => "mat4x3f",

	# Double Matrix
	GL_DOUBLE_MAT2 =>   "mat2d",
	GL_DOUBLE_MAT3 =>   "mat3d",
	GL_DOUBLE_MAT4 =>   "mat4d",
	GL_DOUBLE_MAT2x3 => "mat2x3d",
	GL_DOUBLE_MAT2x4 => "mat2x4d",
	GL_DOUBLE_MAT3x2 => "mat3x2d",
	GL_DOUBLE_MAT3x4 => "mat3x4d",
	GL_DOUBLE_MAT4x2 => "mat4x2d",
	GL_DOUBLE_MAT4x3 => "mat4x3d",

	# Sampler
	GL_SAMPLER_1D =>                   "sampler1D",
	GL_SAMPLER_2D =>                   "sampler2D",
	GL_SAMPLER_3D =>                   "sampler3D",
	GL_SAMPLER_CUBE =>                 "sampler_Cube",
	GL_SAMPLER_CUBE_SHADOW =>          "sampler_Cube_Shadow",
	GL_SAMPLER_1D_SHADOW =>            "sampler1D_Shadow",
	GL_SAMPLER_2D_SHADOW =>            "sampler2D_Shadow",
	GL_SAMPLER_1D_ARRAY =>             "sampler1D_Array",
	GL_SAMPLER_2D_ARRAY =>             "sampler2D_Array",
	GL_SAMPLER_1D_ARRAY_SHADOW =>      "sampler1D_Array_Shadow",
	GL_SAMPLER_2D_ARRAY_SHADOW =>      "sampler2D_Array_Shadow",
	GL_SAMPLER_2D_MULTISAMPLE =>       "sampler2D_Multisample",
	GL_SAMPLER_2D_MULTISAMPLE_ARRAY => "sampler2D_Multisample_Array",
	GL_SAMPLER_BUFFER =>               "sampler_Buffer",
	GL_SAMPLER_2D_RECT =>              "sampler2D_Rect",
	GL_SAMPLER_2D_RECT_SHADOW =>       "sampler2D_Rect_Shadow",

	# Sampler Int
	GL_INT_SAMPLER_1D =>                   "isampler1D",
	GL_INT_SAMPLER_2D =>                   "isampler2D",
	GL_INT_SAMPLER_3D =>                   "isampler3D",
	GL_INT_SAMPLER_CUBE =>                 "isampler_Cube",
	GL_INT_SAMPLER_1D_ARRAY =>             "isampler1D_Array",
	GL_INT_SAMPLER_2D_ARRAY =>             "isampler2D_Array",
	GL_INT_SAMPLER_2D_MULTISAMPLE =>       "isampler2D_Multisample",
	GL_INT_SAMPLER_2D_MULTISAMPLE_ARRAY => "isampler2D_Multisample_Array",
	GL_INT_SAMPLER_BUFFER =>               "isampler_Buffer",
	GL_INT_SAMPLER_2D_RECT =>              "isampler2D_Rect",

	# Sampler Unsigned Int
	GL_UNSIGNED_INT_SAMPLER_1D =>                   "usampler1D",
	GL_UNSIGNED_INT_SAMPLER_2D =>                   "usampler2D",
	GL_UNSIGNED_INT_SAMPLER_3D =>                   "usampler3D",
	GL_UNSIGNED_INT_SAMPLER_CUBE =>                 "usampler_Cube",
	GL_UNSIGNED_INT_SAMPLER_1D_ARRAY =>             "usampler1D_Array",
	GL_UNSIGNED_INT_SAMPLER_2D_ARRAY =>             "usampler2D_Array",
	GL_UNSIGNED_INT_SAMPLER_2D_MULTISAMPLE =>       "usampler2D_Multisample",
	GL_UNSIGNED_INT_SAMPLER_2D_MULTISAMPLE_ARRAY => "usampler2D_Multisample_Array",
	GL_UNSIGNED_INT_SAMPLER_BUFFER =>               "usampler_Buffer",
	GL_UNSIGNED_INT_SAMPLER_2D_RECT =>              "usampler2D_Rect",
)

const BASE_TYPES = Dict{DataType, NamedTuple{(:id, :name),Tuple{GLenum,String}}}(
	Float32 => (id=GL_FLOAT, name="float"),
	Float64 => (id=GL_DOUBLE, name="double"),
	Int32 => (id=GL_INT, name="int"),
	UInt32 => (id=GL_UNSIGNED_INT, name="uint"),
	Bool => (id=GL_BOOL, name="bool"),
)

const UNIFORMS = Dict{Integer,Dict{Integer,Dict{DataType,Function}}}(
	0 => Dict(
		0 => Dict(
			UInt32 => glUniform1ui,
			UInt64 => glUniform1ui,
			Int32 => glUniform1i,
			Int64 => glUniform1i,
			Float32 => glUniform1f,
			Float64 => glUniform1d,
		),
	),
	1 => Dict(
		1 => Dict(
			UInt32 => glUniform1uiv,
			UInt64 => glUniform1uiv,
			Int32 => glUniform1iv,
			Int64 => glUniform1iv,
			Float32 => glUniform1fv,
			Float64 => glUniform1dv,
		),
		2 => Dict(
			UInt32 => glUniform2uiv,
			UInt64 => glUniform2uiv,
			Int32 => glUniform2iv,
			Int64 => glUniform2iv,
			Float32 => glUniform2fv,
			Float64 => glUniform2dv,
		),
		3 => Dict(
			UInt32 => glUniform3uiv,
			UInt64 => glUniform3uiv,
			Int32 => glUniform3iv,
			Int64 => glUniform3iv,
			Float32 => glUniform3fv,
			Float64 => glUniform3dv,
		),
		4 => Dict(
			UInt32 => glUniform4uiv,
			UInt64 => glUniform4uiv,
			Int32 => glUniform4iv,
			Int64 => glUniform4iv,
			Float32 => glUniform4fv,
			Float64 => glUniform4dv,
		),
	),
	2 => Dict(
		4 => Dict(
			Float32 => glUniformMatrix2fv,
			Float64 => glUniformMatrix2dv,
		),
		9 => Dict(
			Float32 => glUniformMatrix3fv,
			Float64 => glUniformMatrix3dv,
		),
		16 => Dict(
			Float32 => glUniformMatrix4fv,
			Float64 => glUniformMatrix4dv,
		),
	),
)

end #GLLists

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# functions from github.com/mcreel/Econometrics, extracted to not depend on the whole thing.

using PrettyTables, Statistics, Printf, Random

function prettyprint(a, cnames="", rnames="")
if rnames !=""
    rnames = rnames[:]
    a = [rnames a]
    if cnames != ""
        cnames = cnames[:]
        cnames = vcat("", cnames)
    end
end
if cnames !=""
    pretty_table(a, cnames; formatters=ft_printf("%12.5f"))
else
    pretty_table(a; formatters=ft_printf("%12.5f"))
end
end

function dstats(x, rnames="";short=false, silent=false)
    k = size(x,2)
    if rnames==""
        rnames = 1:k
        rnames = rnames'
    end
    m = mean(x,dims=1)
    mm = median(x,dims=1)
    s = std(x,dims=1)
    sk = m-m
    kt = m-m
    mn = minimum(x,dims=1)
    mx = maximum(x,dims=1)
    q05 = fill(0.0,k)
    q25 = fill(0.0,k)
    q75 = fill(0.0,k)
    q95 = fill(0.0,k)
    if short == false
        for i = 1:size(x,2) q05[i], q25[i], q75[i],q95[i] = quantile(x[:,i], [0.05,0.25,0.75,0.95]) end
        cnames = ["  mean", " median","  std", "IQR", "min", "max", "q05", "q95"]
        stats = [m' mm' s' (q75-q25) mn' mx' q05 q95]
        if !silent prettyprint(stats, cnames, rnames) end
    else
        cnames = ["  mean", " median", "  std", "min", "max"]
        stats = [m' mm' s' mn' mx']
        if !silent prettyprint(stats, cnames, rnames) end
    end
    return stats
end

# Newey-West covariance estimator
function NeweyWest(Z,nlags=0)
#=
    Returns the Newey-West estimator of the asymptotic variance matrix
    INPUTS: Z, a nxk matrix with rows the vector zt'
            nlags, the number of lags
    OUTPUTS: omegahat, the Newey-West estimator of the covariance matrix
=#
    n,k = size(Z)
    # de-mean the variables
    Z = Z .- mean(Z,dims=1)
    omegahat = Z'*Z/n # sample variance
    # automatic lags?
    if nlags == 0
        nlags = max(0, Int(round(n^0.25)))
    end
    # sample autocovariances
    for i = 1:nlags
       Zlag = Z[1:n-i,:]
       ZZ = Z[i+1:n,:]
       gamma = (ZZ'*Zlag)/n
       weight = 1.0 - (i/(nlags+1.0))
       omegahat += weight*(gamma + gamma')
    end
    return omegahat
end

# method using threads and symmetric proposal
@views function mcmc(θ, reps::Int64, burnin::Int64, Prior::Function, lnL::Function, Proposal::Function, report::Bool, nthreads::Int64)
    perthread = reps ÷ nthreads
    chain = zeros(reps, size(θ,1)+1)
    Threads.@threads for t = 1:nthreads # collect the results from the threads
        chain[t*perthread-perthread+1:t*perthread,:] = mcmc(θ, perthread, burnin, Prior, lnL, Proposal, report)
    end
    return chain
end

# MH method for symmetric proposal (as is the case with SNM methods)
# the main loop
@views function mcmc(θ, reps::Int64, burnin::Int64, Prior::Function, lnL::Function, Proposal::Function, report::Bool=true)
    reportevery = Int((reps+burnin)/10)
    lnLθ = lnL(θ)
    chain = zeros(reps, size(θ,1)+1)
    naccept = zeros(size(θ))
    for rep = 1:reps+burnin
        θᵗ = Proposal(θ) # new trial value
        if report
            changed = Int.(.!(θᵗ .== θ)) # find which changed
        end
        # MH accept/reject: only evaluate logL if proposal is in support of prior (avoid crashes)
        pt = Prior(θᵗ)
        accept = false
        if pt > 0.0
            lnLθᵗ = lnL(θᵗ)
            accept = rand() < exp(lnLθᵗ-lnLθ) * pt/Prior(θ)
            if accept
                θ = θᵗ
                lnLθ = lnLθᵗ
            end
        end
        if report
            naccept = naccept .+ changed .* Int.(accept)
        end
        if (mod(rep,reportevery)==0 && report)
            println("current parameters: ", round.(θ,digits=3))
            println("  acceptance rates: ", round.(naccept/reportevery,digits=3))
            naccept = naccept - naccept
        end
        if rep > burnin
            chain[rep-burnin,:] = vcat(θ, accept)
        end
    end
    return chain
end

#


text:
function get_num_trees(input_file, num_right_steps, num_down_steps, tree_symbol)
    num_trees = open(input_file) do file
        line_width = -1
        current_position = 1 # we start on the top left
        num_trees = 0
        lines_to_skip = 0

        for line in eachline(file)
            if lines_to_skip > 0
                lines_to_skip -= 1
                continue
            end

            cleaned_line = collect(strip(line))
            if line_width == -1
                line_width = length(cleaned_line)
            end
            equality = string(cleaned_line[current_position]) == tree_symbol

            if equality
                num_trees += 1
            end

            current_position = current_position + num_right_steps

            # we have periodic boundary positions
            if current_position > line_width
                current_position = current_position % line_width
            end

            # this is just a way of coding skipping lines
            lines_to_skip += (num_down_steps - 1)
        end
        num_trees
    end
    num_trees
end

input_steps = [
    (input_file = "day3_input.txt", num_right_steps = 1, num_down_steps = 1, tree_symbol = "#"),
    (input_file = "day3_input.txt", num_right_steps = 3, num_down_steps = 1, tree_symbol = "#"),
    (input_file = "day3_input.txt", num_right_steps = 5, num_down_steps = 1, tree_symbol = "#"),
    (input_file = "day3_input.txt", num_right_steps = 7, num_down_steps = 1, tree_symbol = "#"),
    (input_file = "day3_input.txt", num_right_steps = 1, num_down_steps = 2, tree_symbol = "#")]

all_trees = []
for input_settings in input_steps
    num_trees = get_num_trees(input_settings...)
    println("$input_settings num_trees: $num_trees")
    append!(all_trees, num_trees)
end

prod(all_trees)


text:
export MtlEvent, MtlSharedEvent, MtlSharedEventHandle

abstract type MtlAbstractEvent end

const MTLEvent = Ptr{MtEvent}
const MTLSharedEvent = Ptr{MtSharedEvent}
const MTLSharedEventHandle = Ptr{MtSharedEventHandle}

mutable struct MtlEvent <: MtlAbstractEvent
	handle::MTLEvent
	device::MtlDevice
end

mutable struct MtlSharedEvent <: MtlAbstractEvent
	handle::MTLSharedEvent
	device::MtlDevice
end

Base.unsafe_convert(::Type{MTLEvent}, ev::MtlAbstractEvent) = convert(MTLEvent, ev.handle)
Base.unsafe_convert(::Type{MTLSharedEvent}, ev::MtlSharedEvent) = ev.handle

Base.:(==)(a::MtlAbstractEvent, b::MtlAbstractEvent) = a.handle == b.handle
Base.hash(ev::MtlAbstractEvent, h::UInt) = hash(ev.handle, h)

function unsafe_destroy!(fun::MtlAbstractEvent)
	mtRelease(fun.handle)
end

function MtlEvent(dev::MtlDevice)
	handle = mtDeviceNewEvent(dev)
	obj = MtlEvent(handle, dev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

function MtlSharedEvent(dev::MtlDevice)
	handle = mtDeviceNewSharedEvent(dev)
	obj = MtlSharedEvent(handle, dev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

## properties

Base.propertynames(::MtlAbstractEvent) = (:device, :label, :signaledValue)

function Base.getproperty(ev::MtlAbstractEvent, f::Symbol)
    if f == :label
        ptr = mtEventLabel(ev)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    elseif ev isa MtlSharedEvent && f == :signaledValue
        mtSharedEventSignaledValue(ev)
    else
        getfield(ev, f)
    end
end

## shared event handle

mutable struct MtlSharedEventHandle
	handle::MTLSharedEventHandle
	event::MtlSharedEvent
end

function MtlSharedEventHandle(ev::MtlSharedEvent)
	handle = mtSharedEventNewHandle(ev)
	obj = MtlSharedEventHandle(handle, ev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

function unsafe_destroy!(evh::MtlSharedEventHandle)
	mtRelease(evh.handle)
end

Base.unsafe_convert(::Type{MTLSharedEventHandle}, evh::MtlSharedEventHandle) = evh.handle

Base.:(==)(a::MtlSharedEventHandle, b::MtlSharedEventHandle) = a.handle == b.handle
Base.hash(evh::MtlSharedEventHandle, h::UInt) = hash(evh.handle, h)


text:
module ChainRulesCore
using Base.Broadcast: broadcasted, Broadcasted, broadcastable, materialize, materialize!
using Base.Meta
using LinearAlgebra
using SparseArrays: SparseVector, SparseMatrixCSC
using Compat: hasfield, hasproperty

export frule, rrule  # core function
# rule configurations
export RuleConfig, HasReverseMode, NoReverseMode, HasForwardsMode, NoForwardsMode
export frule_via_ad, rrule_via_ad
# definition helper macros
export @non_differentiable, @opt_out, @scalar_rule, @thunk, @not_implemented
export ProjectTo, canonicalize, unthunk  # tangent operations
export add!!  # gradient accumulation operations
export ignore_derivatives, @ignore_derivatives
# tangents
export Tangent, NoTangent, InplaceableThunk, Thunk, ZeroTangent, AbstractZero, AbstractThunk

include("compat.jl")
include("debug_mode.jl")

include("tangent_types/abstract_tangent.jl")
include("tangent_types/abstract_zero.jl")
include("tangent_types/thunks.jl")
include("tangent_types/tangent.jl")
include("tangent_types/notimplemented.jl")

include("tangent_arithmetic.jl")
include("accumulation.jl")
include("projection.jl")

include("config.jl")
include("rules.jl")
include("rule_definition_tools.jl")
include("ignore_derivatives.jl")

include("deprecated.jl")

end # module


text:
function recurring_decimal_length(n)
    d = 1
    while true
        if powermod(10, d, n) == 1
            return d
        end
        d += 1
    end
end

function longest_recurring_decimal(upper_limit)
    longest = 0; arg = 0

    for i in 2:upper_limit
        if i % 2 == 0 || i % 5 == 0
            continue
        end

        len = recurring_decimal_length(i)
        if len > longest
            longest = len; arg = i
        end
    end

    return arg
end

println(longest_recurring_decimal(1000))

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
# # Add time series to `System`

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# An example of how to parse add time series data to a `System` using [PowerSystems.jl](github.com/NREL-SIIP/PowerSystems.jl)
#
# For example, a `System` created by [parsing a MATPOWER file](../../notebook/2_PowerSystems_examples/parse_matpower.ipynb)
# doesn't contain any time series data. So a user may want to add time series to the `System`
# ### Dependencies
# Let's use the 5-bus dataset we parsed in the MATPOWER example
using SIIPExamples
using PowerSystems
using JSON3

pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "02_parse_matpower.jl"))

# ### Define pointers to time series files
# For example, if we want to add a bunch of time series files, say one for each load and
# one for each renewable generator, we need to define pointers to each .csv file containing
# the time series in the following format (PowerSystems.jl also supports a CSV format for this file)

FORECASTS_DIR = joinpath(base_dir, "forecasts", "5bus_ts")
fname = joinpath(FORECASTS_DIR, "timeseries_pointers_da.json")
open(fname, "r") do f
    JSON3.@pretty JSON3.read(f)
end

# ### Read and assign time series to `System` using these parameters.
add_time_series!(sys, fname)
sys


text:
using
    OceanTurb,
    LinearAlgebra,
    Test

Timer(t -> println(" "), 0, interval=240)

#
# Run tests
#

steppers = (:ForwardEuler, :BackwardEuler)

include("test_utils.jl")
include("test_solvers.jl")
include("test_grids.jl")
include("test_fields.jl")

# Models
include("test_diffusion.jl")
include("test_pp.jl")
include("test_kpp.jl")
include("test_modularkpp.jl")
include("test_tkemassflux.jl")


text:
using StatsBase, StaticArrays, BitBasis, Random
export measure, measure!, select, select!

function _measure(rng::AbstractRNG, pl::AbstractVector, nshots::Int)
    N = log2i(length(pl))
    sample(rng, basis(BitStr64{N}), Weights(pl), nshots)
end

function _measure(rng::AbstractRNG, pl::AbstractMatrix, nshots::Int)
    B = size(pl, 2)
    res = Matrix{BitStr64{log2i(size(pl, 1))}}(undef, nshots, B)
    for ib = 1:B
        @inbounds res[:, ib] = _measure(rng, view(pl, :, ib), nshots)
    end
    return res
end

YaoBase.measure(
    ::ComputationalBasis,
    reg::ArrayReg{1},
    ::AllLocs;
    nshots::Int = 1,
    rng::AbstractRNG = Random.GLOBAL_RNG,
) = _measure(rng, reg |> probs, nshots)

function YaoBase.measure(
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs;
    nshots::Int = 1,
    rng::AbstractRNG = Random.GLOBAL_RNG,
) where {B}
    pl = dropdims(sum(reg |> rank3 .|> abs2, dims = 2), dims = 2)
    return _measure(rng, pl, nshots)
end

function YaoBase.measure!(
    ::YaoBase.RemoveMeasured,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs;
    rng::AbstractRNG = Random.GLOBAL_RNG,
) where {B}
    state = reg |> rank3
    nstate = similar(reg.state, 1 << nremain(reg), B)
    pl = dropdims(sum(state .|> abs2, dims = 2), dims = 2)
    res = Vector{BitStr64{nactive(reg)}}(undef, B)
    @inbounds for ib = 1:B
        ires = _measure(rng, view(pl, :, ib), 1)[]
        # notice ires is `BitStr` type, can be use as indices directly.
        nstate[:, ib] = view(state, Int64(ires) + 1, :, ib) ./ sqrt(pl[Int64(ires)+1, ib])
        res[ib] = ires
    end
    reg.state = reshape(nstate, 1, :)
    return B == 1 ? res[] : res
end

function YaoBase.measure!(
    ::YaoBase.NoPostProcess,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs;
    rng::AbstractRNG = Random.GLOBAL_RNG,
) where {B}
    state = reg |> rank3
    nstate = zero(state)
    res = measure!(RemoveMeasured(), reg; rng = rng)
    _nstate = reshape(reg.state, :, B)
    indices = Int64.(res) .+ 1
    for ib = 1:B
        @inbounds nstate[indices[ib], :, ib] .= view(_nstate, :, ib)
    end
    reg.state = reshape(nstate, size(state, 1), :)
    return res
end

function YaoBase.measure!(
    rst::YaoBase.ResetTo,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs;
    rng::AbstractRNG = Random.GLOBAL_RNG,
) where {B}
    state = rank3(reg)
    M, N, B1 = size(state)
    nstate = zero(state)
    res = measure!(YaoBase.RemoveMeasured(), reg; rng = rng)
    nstate[Int(rst.x)+1, :, :] = reshape(reg.state, :, B)
    reg.state = reshape(nstate, M, N * B)
    return res
end

import YaoBase: select, select!
select(r::ArrayReg{B}, bits::AbstractVector{T}) where {B,T<:Integer} =
    ArrayReg{B}(r.state[Int64.(bits).+1, :])
select(r::ArrayReg{B}, bit::Integer) where {B} = select(r, [bit])

function select!(r::ArrayReg, bits::AbstractVector{T}) where {T<:Integer}
    r.state = r.state[Int64.(bits).+1, :]
    return r
end

select!(r::ArrayReg, bit::Integer) = select!(r, [bit])


text:
mutable struct Expo <: AbstractTweenEquation
    name::String
    # A function that accepts a time argument and returns a float
    compute::Function # f(t::Float64) -> Float64
    type::Int64

    function Expo_EaseIn(type::Int64)
        o = new()

        if type == EASE_IN
            o.type = EASE_IN
            o.name = "Expo.easeIn"
            o.compute = _exp_ease_in
        elseif type == EASE_OUT
            o.type = EASE_OUT
            o.name = "Expo.easeOut"
            o.compute = _exp_ease_out
        else
            o.type = EASE_INOUT
            o.name = "Expo.easeInOut"
            o.compute = _exp_ease_inout
        end

        o
    end
end

function _exp_ease_in(t::Float64)
    if t == 0.0
        0.0
    else
        2^(10.0 * (t - 1.0))
    end
end

function _exp_ease_out(t::Float64)
    if t == 1.0
        1.0
    else
        -(2.0^(-10.0 * t)) + 1.0
    end
end

function _exp_ease_inout(t::Float64)
    if t == 0.0
        return 0.0
    end

    if t == 1.0
        return 1.0
    end

    a = t * 2.0
    if a < 1.0
        return 0.5 * (2.0^(10.0 * (a - 1.0)))
    end

    0.5 * (-(2.0^(-10.0 * (t - 1.0))) + 2.0)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    StratSim(var₁=>param₁, var₂=>param₂, ...)

Stratigraphy simulation with Markov-Poisson sampling.

## Parameters

* `environment` - geological environment
* `state`       - initial geological state
* `stack`       - stacking scheme (:erosional or :depositional)
* `nepochs`     - number of epochs (default to 10)
* `fillbase`    - fill value for the bottom layer (default to `NaN`)
* `filltop`     - fill value for the top layer (default to `NaN`)

### References

Hoffimann 2018. *Morphodynamic analysis and statistical
synthesis of geormorphic data.*
"""
@simsolver StratSim begin
  @param environment
  @param state = nothing
  @param stack = :erosional
  @param nepochs = 10
  @param fillbase = NaN
  @param filltop = NaN
end

function preprocess(problem::SimulationProblem, solver::StratSim)
  # retrieve problem info
  pdomain = domain(problem)

  @assert ndims(pdomain) == 3 "solver implemented for 3D domain only"

  # result of preprocessing
  preproc = Dict{Symbol,NamedTuple}()

  for (var, V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = StratSimParam()
    end

    # determine environment
    environment = varparams.environment

    # determine initial state
    if varparams.state ≠ nothing
      state = varparams.state
    else
      nx, ny, nz = size(pdomain)
      land = zeros(nx, ny)
      state = LandState(land)
    end

    # determine stacking scheme
    stack = varparams.stack

    # determine number of epochs
    nepochs = varparams.nepochs

    # determine fill values
    fillbase = varparams.fillbase
    filltop  = varparams.filltop

    # save preprocessed input
    preproc[var] = (environment=environment, state=state,
                    stack=stack, nepochs=nepochs,
                    fillbase=fillbase, filltop=filltop)
  end

  preproc
end

function solve_single(problem::SimulationProblem, var::Symbol,
                      solver::StratSim, preproc)
  # retrieve problem info
  pdomain = domain(problem)
  _, __, nz = size(pdomain)

  # get parameters for the variable
  environment, state, stack, nepochs, fillbase, filltop = preproc[var]

  # simulate the environment
  record = simulate(environment, state, nepochs)

  # build stratigraphy
  strata = Strata(record, stack)

  # return voxel model
  model = voxelize(strata, nz, fillbase=fillbase, filltop=filltop)

  # flatten result
  vec(model)
end


text:
# provide some most commonly used options, leave rest as low level
# common options: Orthogonilization type, KSP type, PC type
# use PC context created as part of KSP

export KSP, petscview

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

# solver context
mutable struct KSP{T}
  p::C.KSP{T}
  function KSP(p::C.KSP{T}) where {T}
    o = new{T}(p)
    if (!deactivate_finalizers)
      finalizer(PetscDestroy,o)
    end
    return o
  end
end

comm(a::KSP{T}) where {T} = C.PetscObjectComm(T, a.p.pobj)

function PetscDestroy(o::KSP{T}) where {T}
  PetscFinalized(T) || C.KSPDestroy(Ref(o.p))
end

function petscview(o::KSP{T}) where {T}
  viewer = C.PetscViewer{T}(C_NULL)
  chk(C.KSPView(o.p, viewer))
end

"""
    KSP(A::Mat, PA=A; kws...)
    KSP(pc::PC; kws...)

Create a KSP solver object that can be used to solve equations `Ax=b` with
the matrix `A`, where `PA` (defaults to `A`) is used to construct the
default preconditioner.  Alternatively, you can supply a preconditioner
object (`PC`).

The keyword options are zero or more of the following:

These control the solver and preconditioner characteristics:
* `ksp_type="a"`: use KSP algorithm `a`
* `ksp_pc_side=n`: set preconditioner side to `GridapDistributedPETScWrappers.C.PC_LEFT`, `GridapDistributedPETScWrappers.C.PC_RIGHT`, or `GridapDistributedPETScWrappers.C.PC_SYMMETRIC`
* `ksp_reuse_preconditioner=true`: use initial preconditioner and don't ever compute a new one
* `ksp_diagonal_scale=true`: symmetrically diagonally scale `A` before solving (note that this *changes* `A` and the right-hand side in a solve, unless you also set `ksp_diagonal_scale_fix=true`)
* `ksp_diagonal_scale_fix=true`: undo diagonal scaling after solve
* `ksp_knoll=true`: use preconditioner applied to `b` for initial guess
* `ksp_constant_null_space=true`: add constant null space to Krylov solver matrix
* `ksp_initial_guess_nonzero=true`: use the contents of initial `x` instead of zero for initial guess
* `ksp_fischer_guess="model,size"`: use Fischer initial guess generator (`model=1` or `2`) for repeated linear solves with subspace of dimension `size`

The following keyword options control the stopping criteria for
iterative solvers:
* `ksp_rtol=x`: `x` is relative decrease in residual norm
* `ksp_atol=x`: `x` is absolute decrease in residual norm
* `ksp_divtol=x`: `x` is amount residual can increase before method is considered to be diverging
* `ksp_max_it=n`: `n` is the max number of iterations
* `ksp_converged_use_initial_residual_norm=true`: use initial residual norm for computing relative convergence
* `ksp_converged_use_min_initial_residual_norm=true`: use min of initial residual norm and `b` for computing relative convergence
* `ksp_error_if_not_converged=true`: generate error if solver does not converge
* `ksp_convergence_test=:default` or `:skip`: use the default convergence test (tolerances and `max_it`) or skip convergence tests and run until `max_it` is reached
* `ksp_norm_type=n`: in residual tests, use norm type `n`, one of default (`GridapDistributedPETScWrappers.C.KSP_NORM_DEFAULT`), none (`GridapDistributedPETScWrappers.C.KSP_NORM_NONE`), of the preconditioned residual (`GridapDistributedPETScWrappers.C.KSP_NORM_PRECONDITIONED`), the true residual (`GridapDistributedPETScWrappers.C.KSP_NORM_UNPRECONDITIONED`), or the "natural" norm (`GridapDistributedPETScWrappers.C.KSP_NORM_NATURAL`)
* `ksp_check_norm_iteration=n`: compute residual norm starting on iteration `n`
* `ksp_lag_norm=true`: lag the calculation of the residual norm by one iteration (trades off reduced communication for an additional iteration)

The following options control output that monitors the progress of the
solver (default none).
* `ksp_monitor=filename`: print the residual norm at each iteration to `filename` (`""` for `STDOUT`)
* `ksp_monitor_short=filename`: print preconditioned residual norm with fewer digits
* `ksp_monitor_range=filename`: prints the percentage of residual elements that are more then 10% of the maximum value
* `ksp_monitor_true_residual=filename`: print true residual norm
* `ksp_monitor_singular_value=filename`: print extreme singular values (via Lanczos or Arnoldi process as the linear system is solved)
* `ksp_monitor_solution=true`: plot solution graphically
* `ksp_monitor_lg_residualnorm=true`: plot preconditioned residual norm graphically
* `ksp_monitor_lg_true_residualnorm=true`: plot preconditioned and true residual norm graphically
* `ksp_monitor_lg_range=true`: plot preconditioned residual norm and range of residual values
* `ksp_monitor_cancel=true`: remove any hardwired monitor routines
* `ksp_compute_singularvalues=true`: print extreme singular values (via Lanczos or Arnoldi process as the linear system is solved)

In addition, if default preconditioner is being used,
then any of the preconditioner options (see `PC`) can be specified to control
this preconditioner (e.g. `pc_type`).
"""
function KSP(pc::PC{T}; kws...) where {T}
  ksp_c = Ref{C.KSP{T}}()
  chk(C.KSPCreate(comm(pc), ksp_c))
  ksp = ksp_c[]
  chk(C.KSPSetPC(ksp, pc.p))
  withoptions(T, kws) do
    chk(C.KSPSetFromOptions(ksp))
  end
  return KSP(ksp)
end

KSP(A::Mat{T}, PA::Mat{T}=A; kws...) where {T} = KSP(PC(A, PA; kws...); kws...)

export KSPSetUp!
function KSPSetUp!(ksp::KSP{T}) where {T}
  chk(C.KSPSetUp(ksp.p))
end

# Retrieve a reference to the matrix in the KSP object, as a raw C.Mat
# pointer.  Note that we should not wrap this in a Mat object, or
# call MatDestroy on it, without incrementing its reference count first!
function _ksp_A(ksp::KSP{T}) where {T}
  a = Ref{C.Mat{T}}()
  pa = Ref{C.Mat{T}}()
  chk(C.KSPGetOperators(ksp.p, a, pa))
  return a[]
end

# x = A \ b
function LinearAlgebra.ldiv!(ksp::KSP{T}, b::Vec{T}, x::Vec{T}) where {T}

 assemble(_ksp_A(ksp))
 assemble(b)
 assemble(x)

 chk(C.KSPSolve(ksp.p, b.p, x.p))

 reason = Ref{Cint}()
 chk(C.KSPGetConvergedReason(ksp.p, reason))
 if (reason[] < 0)
   @warn "KSP solve did not converge"
 end

 return x
end

# x = A.' \ b
# Base.At_ldiv_B! does not exist? the non bang version does
function KSPSolveTranspose!(ksp::KSP{T}, b::Vec{T}, x::Vec{T}) where {T}
  assemble(_ksp_A(ksp))
  assemble(b)
  assemble(x)

  chk(C.KSPSolveTranspose(ksp.p, b.p, x.p))

  reason = Ref{Cint}()
  chk(C.KSPGetConvergedReason(ksp.p, reason))
  if (reason[] < 0)
    @warn "KSP solve did not converge"
  end

  return x
end

# is there nice syntax for this?
export KSPSolveTranspose
function KSPSolveTranspose(ksp::KSP{T}, b::Vec{T}) where {T}

  x = similar(b, size(ksp, 1))
  KSPSolveTranspose!(ksp, b, x)
  return x
end

function Base.size(ksp::KSP)
  m = Ref{PetscInt}()
  n = Ref{PetscInt}()
  chk(C.MatGetSize(_ksp_A(ksp), m, n))
  (Int(m[]), Int(n[]))
end
Base.size(ksp::KSP{T}, dim::Integer)  where {T} = dim > 2 ? 1 : size(ksp)[dim]

import Base: \
(\)(ksp::KSP{T}, b::Vec{T}) where {T} = ldiv!(ksp, b, similar(b, size(ksp, 2)))

# Mat fallbacks
(\)(A::Mat{T}, b::Vec{T}) where {T} = KSP(A)\b
LinearAlgebra.ldiv!(A::Mat{T}, b::Vec{T}, x::Vec{T}) where {T} = ldiv!(KSP(A), b, x)

KSPSolveTranspose(A::Mat{T}, b::Vec{T}) where {T} = KSPSolveTranspose(KSP(A), b)
KSPSolveTranspose(A::Mat{T}, b::Vec{T}, x::Vec{T}) where {T} = KSPSOlveTranspose(KSP(A), x, b)


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using MathOptInterface
import GLPK
import Clp
const MOI = MathOptInterface

function example_diet(optimizer, bridge)
    category_data = [
        1800.0 2200.0
          91.0    Inf
           0.0   65.0
           0.0 1779.0
    ]
    cost = [2.49, 2.89, 1.50, 1.89, 2.09, 1.99, 2.49, 0.89, 1.59]
    food_data = [
        410 24 26 730
        420 32 10 1190
        560 20 32 1800
        380  4 19 270
        320 12 10 930
        320 15 12 820
        320 31 12 1230
        100  8 2.5 125
        330  8 10 180
    ]
    bridge_model = if bridge
        MOI.instantiate(optimizer; with_bridge_type=Float64)
    else
        MOI.instantiate(optimizer)
    end
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.Utilities.AUTOMATIC,
    )
    MOI.Utilities.reset_optimizer(model, bridge_model)
    MOI.set(model, MOI.Silent(), true)
    nutrition = MOI.add_variables(model, size(category_data, 1))
    for (i, v) in enumerate(nutrition)
        MOI.add_constraint(model, v, MOI.GreaterThan(category_data[i, 1]))
        MOI.add_constraint(model, v, MOI.LessThan(category_data[i, 2]))
    end
    buy = MOI.add_variables(model, size(food_data, 1))
    MOI.add_constraint.(model, buy, MOI.GreaterThan(0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(cost, buy), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    for (j, n) in enumerate(nutrition)
        f = MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(food_data[:, j], buy),
            0.0,
        )
        push!(f.terms, MOI.ScalarAffineTerm(-1.0, n))
        MOI.add_constraint(model, f, MOI.EqualTo(0.0))
    end
    MOI.optimize!(model)
    term_status = MOI.get(model, MOI.TerminationStatus())
    @assert term_status == MOI.OPTIMAL
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(1.0, [buy[end-1], buy[end]]),
            0.0,
        ),
        MOI.LessThan(6.0),
    )
    MOI.optimize!(model)
    @assert MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    return
end

if length(ARGS) > 0
    bridge = get(ARGS, 2, "") != "--no-bridge"
    println("Running: $(ARGS[1]) $(get(ARGS, 2, ""))")
    if ARGS[1] == "clp"
        @time example_diet(Clp.Optimizer, bridge)
        @time example_diet(Clp.Optimizer, bridge)
    else
        @assert ARGS[1] == "glpk"
        @time example_diet(GLPK.Optimizer, bridge)
        @time example_diet(GLPK.Optimizer, bridge)
    end
    exit(0)
end

# using SnoopCompileCore
# tinf = @snoopi_deep example_diet(Clp.Optimizer, true)
# using SnoopCompile, ProfileView
# fg = flamegraph(tinf)
# ProfileView.view(fg)

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using Chaos
using Documenter

DocMeta.setdocmeta!(Chaos, :DocTestSetup, :(using Chaos); recursive=true)

makedocs(;
    modules=[Chaos],
    authors="ErikQQY",
    repo="https://github.com/ErikQQY/Chaos/blob/{commit}{path}#{line}",
    sitename="Chaos.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ErikQQY.github.io/Chaos",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ErikQQY/Chaos",
    devbranch="master",
)


text:
using Newton
using Test
using ForwardDiff
using LinearAlgebra
using StaticArrays

multiinput_rf!(r::Vector, x::Vector, A::Matrix, b::Vector) = (r .= b .+ A*x)

function setup_cache(x, A, b)
    rf!(r, x) = multiinput_rf!(r, x, A, b)
    return NewtonCache(x, rf!)
end

function determine_solution(A, b)
    x = zero(b)
    rf!(r, x) = multiinput_rf!(r, x, A, b)
    cache=NewtonCache(x, rf!)
    newtonsolve!(x, get_drdx(cache), rf!, cache)
    return x
end

@testset "linsolve!" begin
    A = 2*I + rand(10,10)
    b = rand(10)
    A1, A2 = [copy(A) for _ in 1:2]
    b1, b2 = [copy(b) for _ in 1:2]
    rf!(r, x) = (r .= x)
    @test A1\b1 ≈ Newton.linsolve!(A2, b2, NewtonCache(b, rf!))
end

@testset "newtonsolve!" begin
    nsize = 4
    (a,b,x0) = [rand(nsize) for _ in 1:3]
    tol = 1.e-10

    # Basic functionality
    function rf_solution!(r, x)
        r .= - a + b.*x + exp.(x)
    end
    x = copy(x0)
    cache = NewtonCache(x, rf_solution!)
    drdx = get_drdx(cache)
    @test size(drdx,1) == nsize
    converged = newtonsolve!(x, drdx, rf_solution!; tol=tol)
    r_check = similar(x0)
    @test converged
    @test isapprox(rf_solution!(r_check, x), zero(r_check); atol=tol)
    @test drdx ≈ ForwardDiff.jacobian(rf_solution!, r_check, x)

    function rf_nosolution!(r, x)
        r .= a .+ b.*x.^2
    end
    x = copy(x0)
    converged = newtonsolve!(x, drdx, rf_nosolution!)
    @test !converged

    # Test function with different anynomous functions for cache and solving
    cache = setup_cache(x0, zeros(nsize,nsize), zeros(nsize))   # If these A and b inputs are used, solution is zero
    A = nsize*I+rand(nsize,nsize)
    A += transpose(A)   # Symmetrize to ensure invertible
    b = rand(nsize)

    rf_solve!(r, x) = multiinput_rf!(r, x, A, b)
    drdx = get_drdx(cache)
    @test newtonsolve!(x, drdx, rf_solve!, cache)
    @test isapprox(rf_solve!(r_check, x), zero(r_check); atol=tol)
    @test .!(isapprox(x, zero(x); atol=tol))

    # Test that error is thrown if we try to use automatic differentiation of the newtonsolve!
    diff_fun(y) = determine_solution(A, y)
    failed = false
    try
        df = ForwardDiff.jacobian(diff_fun, rand(nsize))
    catch err
        failed = true
    end
    @test failed

end

@testset "linsolve" begin
    dim = 5
    A = rand(SMatrix{dim,dim})
    x = rand(SVector{dim})
    b = A*x
    @test x ≈ Newton.linsolve(A, b)
end

@testset "newtonsolve" begin
    dim = 6
    (a,b,x0) = [rand(SVector{dim}) for _ in 1:3]

    rf_solution(x) = - a + b.*x + exp.(x)
    rf_nosolution(x) = a .+ b.*x.^2

    converged, x, drdx = newtonsolve(x0, rf_solution; tol=1.e-6)
    @test converged
    @test isapprox(norm(rf_solution(x)), 0.0, atol=1.e-6)
    @test drdx ≈ ForwardDiff.jacobian(rf_solution, x)

    converged, x, drdx = newtonsolve(x0, rf_nosolution; tol=1.e-6, maxiter=4)
    @test ~converged
    @test all(isnan.(x))
    @test all(isnan.(drdx))

end


text:
import Random.shuffle
import Random.shuffle!

function swap_degree!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1, i2 = sample(iy, 2, replace=false)
   n1, n2 = (from=i1.from, to=i2.to), (from=i2.from, to=i1.to)

   while (n2 ∈ iy)|(n1 ∈ iy)
      i1, i2 = sample(iy, 2, replace=false)
      n1, n2 = (from=i1.from, to=i2.to), (from=i2.from, to=i1.to)
   end

   for old_i in [i1, i2, n1, n2]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

function swap_fill!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1 = sample(iy)
   n1 = (from=sample(species(Y; dims=1)), to=sample(species(Y; dims=2)))

   while (n1 ∈ iy)
      n1 = (from=sample(species(Y; dims=1)), to=sample(species(Y; dims=2)))
   end

   for old_i in [i1, n1]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

function swap_vulnerability!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1 = sample(iy)
   n1 = (from=sample(species(Y; dims=1)), to=i1.to)

   while (n1 ∈ iy)
      n1 = (from=sample(species(Y; dims=1)), to=i1.to)
   end

   for old_i in [i1, n1]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

function swap_generality!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1 = sample(iy)
   n1 = (from=i1.from, to=sample(species(Y; dims=2)))

   while (n1 ∈ iy)
      n1 = (from=i1.from, to=sample(species(Y; dims=2)))
   end

   for old_i in [i1, n1]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

"""
    shuffle(N::BinaryNetwork; constraint::Symbol=:degree)

Return a shuffled copy of the network (the original network is not modified).
See `shuffle!` for a documentation of the keyword arguments.
"""
function shuffle(N::BinaryNetwork; constraint::Symbol=:degree)
   Y = copy(N)
   shuffle!(Y; constraint=constraint)
   return Y
end

"""
    shuffle!(N::BinaryNetwork; constraint::Symbol=:degree)

Shuffles interactions inside a network (the network is *modified*), under the
following `constraint`:

- :degree, which keeps the degree distribution intact
- :generality, which keeps the out-degree distribution intact
- :vulnerability, which keeps the in-degree distribution intact
- :fill, which moves interactions around freely

The function will take two interactions, and swap the species establishing them.
By repeating the process a large enough number of times, the resulting network
should be relatively random. Note that this function will conserve the degree
(when appropriate under the selected constraint) of *every* species. Calling the
function will perform **a single** shuffle. If you want to repeat the shuffling
a large enough number of times, you can use something like:

    [shuffle!(n) for i in 1:10_000]

If the keyword arguments are invalid, the function will throw an
`ArgumentError`.
"""
function shuffle!(N::BinaryNetwork; constraint::Symbol=:degree)
   constraint ∈ [:degree, :generality, :vulnerability, :fill] || throw(ArgumentError("The constraint argument you specificied ($(constraint)) is invalid -- see ?shuffle! for a list."))

   f = EcologicalNetworks.swap_degree!
   constraint == :generality && (f = EcologicalNetworks.swap_generality!)
   constraint == :vulnerability && (f = EcologicalNetworks.swap_vulnerability!)
   constraint == :fill && (f = EcologicalNetworks.swap_fill!)

   f(N)
end

function quanti_shuffle_valid(a, i1, i2)
    u1 = length(unique(i1)) == 2
    u2 = length(unique(i2)) == 2
    mi = minimum(a[i1,i2]) != zero(eltype(a))
    return u1 & u2 & mi
end

"""
TODO
"""
function shuffle!(N::BipartiteQuantitativeNetwork{Int64,K}; constraint::Symbol=:degree, lims=(0, Inf)) where {K <: AllowedSpeciesTypes}
    constraint ∈ [:degree, :generality, :vulnerability, :fill] || throw(ArgumentError("The constraint argument you specificied ($(constraint)) is invalid -- see ?shuffle! for a list."))
    first(lims) < last(lims) || throw(ArgumentError("The interval you used for limits, $(lims), is invalid"))
    all(first(lims) .<= N.A .<= last(lims)) || throw(ArgumentError("The interval you used for limits, $(lims), is smaller than the range of network values"))

    # List of compatible swap schemes under different constraints
    same_row = [[-1 1; -1 1], [-1 1; 0 0], [-1 1; 1 -1], [0 0; -1 1], [0 0; 1 -1], [1 -1; -1 1], [1 -1; 0 0], [1 -1; 1 -1],]
    same_col = [[-1 -1; 1 1], [-1 0; 1 0], [-1 1; 1 -1], [0 -1; 0 1], [0 1; 0 -1], [1 -1; -1 1], [1 0; -1 0], [1 1; -1 -1],]
    same_all = filter(x -> x ∈ same_row, same_col)
    same_lnk = unique(vcat(same_row, same_col))

    # Pick the correcty matrix series
    m = same_lnk
    constraint == :degree && (m = same_all)
    constraint == :vulnerability && (m = same_col)
    constraint == :generality && (m = same_row)

    permuted = false
    while !permuted
        s1 = rand(1:richness(N; dims=1), 2)
        s2 = rand(1:richness(N; dims=2), 2)
        while !quanti_shuffle_valid(N.A, s1, s2)
            s1 = rand(1:richness(N; dims=1), 2)
            s2 = rand(1:richness(N; dims=2), 2)
        end
        st1 = filter(a -> all(first(lims).<=a.<=last(lims)), [N[s1,s2] .+ x for x in m])
        st2 = filter(a -> all((N[s1,s2].!=0).==(a.!=0)), st1)
        if length(st2) > 0
            permuted = true
            k = rand(st2)
            for i in 1:2, j in 1:2
                N[s1[i],s2[j]] = k[i,j]
            end
        end
    end
end

"""
TODO
"""
function shuffle!(N::BipartiteQuantitativeNetwork{Float64,K}; ε::Float64=1e-2, constraint::Symbol=:degree, lims=(0.0, Inf)) where {K <: AllowedSpeciesTypes}
    constraint ∈ [:degree, :generality, :vulnerability, :fill] || throw(ArgumentError("The constraint argument you specificied ($(constraint)) is invalid -- see ?shuffle! for a list."))
    first(lims) < last(lims) || throw(ArgumentError("The interval you used for limits, $(lims), is invalid"))
    all(first(lims) .<= N.A .<= last(lims)) || throw(ArgumentError("The interval you used for limits, $(lims), is smaller than the range of network values"))

    # List of compatible swap schemes under different constraints
    same_row = [[-1 1; -1 1], [-1 1; 0 0], [-1 1; 1 -1], [0 0; -1 1], [0 0; 1 -1], [1 -1; -1 1], [1 -1; 0 0], [1 -1; 1 -1],]
    same_col = [[-1 -1; 1 1], [-1 0; 1 0], [-1 1; 1 -1], [0 -1; 0 1], [0 1; 0 -1], [1 -1; -1 1], [1 0; -1 0], [1 1; -1 -1],]
    same_all = filter(x -> x ∈ same_row, same_col)
    same_lnk = unique(vcat(same_row, same_col))

    # Pick the correcty matrix series
    m = same_lnk
    constraint == :degree && (m = same_all)
    constraint == :vulnerability && (m = same_col)
    constraint == :generality && (m = same_row)
    m = m.*ε

    permuted = false
    while !permuted
        s1 = rand(1:richness(N; dims=1), 2)
        s2 = rand(1:richness(N; dims=2), 2)
        while !quanti_shuffle_valid(N.A, s1, s2)
            s1 = rand(1:richness(N; dims=1), 2)
            s2 = rand(1:richness(N; dims=2), 2)
        end
        st1 = filter(a -> all(first(lims).<=a.<=last(lims)), [N[s1,s2] .+ x for x in m])
        st2 = filter(a -> all((N[s1,s2].!=0).==(a.!=0)), st1)
        if length(st2) > 0
            permuted = true
            k = rand(st2)
            for i in 1:2, j in 1:2
                N[s1[i],s2[j]] = k[i,j]
            end
        end
    end
end

"""
TODO
"""
function shuffle!(N::UnipartiteQuantitativeNetwork{Int64,K}; constraint::Symbol=:degree, lims=(0, Inf)) where {K <: AllowedSpeciesTypes}
    constraint ∈ [:degree, :generality, :vulnerability, :fill] || throw(ArgumentError("The constraint argument you specificied ($(constraint)) is invalid -- see ?shuffle! for a list."))
    first(lims) < last(lims) || throw(ArgumentError("The interval you used for limits, $(lims), is invalid"))
    all(first(lims) .<= N.A .<= last(lims)) || throw(ArgumentError("The interval you used for limits, $(lims), is smaller than the range of network values"))

    # List of compatible swap schemes under different constraints
    same_row = [[-1 1; -1 1], [-1 1; 0 0], [-1 1; 1 -1], [0 0; -1 1], [0 0; 1 -1], [1 -1; -1 1], [1 -1; 0 0], [1 -1; 1 -1],]
    same_col = [[-1 -1; 1 1], [-1 0; 1 0], [-1 1; 1 -1], [0 -1; 0 1], [0 1; 0 -1], [1 -1; -1 1], [1 0; -1 0], [1 1; -1 -1],]
    same_all = filter(x -> x ∈ same_row, same_col)
    same_lnk = unique(vcat(same_row, same_col))

    # Pick the correcty matrix series
    m = same_lnk
    constraint == :degree && (m = same_all)
    constraint == :vulnerability && (m = same_col)
    constraint == :generality && (m = same_row)

    permuted = false
    while !permuted
        s1 = rand(1:richness(N; dims=1), 2)
        while !quanti_shuffle_valid(N.A, s1, s1)
            s1 = rand(1:richness(N; dims=1), 2)
        end
        st1 = filter(a -> all(first(lims).<=a.<=last(lims)), [N[s1,s1] .+ x for x in m])
        st2 = filter(a -> all((N[s1,s1].!=0).==(a.!=0)), st1)
        if length(st2) > 0
            permuted = true
            k = rand(st2)
            for i in 1:2, j in 1:2
                N[s1[i],s1[j]] = k[i,j]
            end
        end
    end
end

"""
TODO
"""
function shuffle!(N::UnipartiteQuantitativeNetwork{Float64,K}; ε::Float64=1e-2, constraint::Symbol=:degree, lims=(0.0, Inf)) where {K <: AllowedSpeciesTypes}
    constraint ∈ [:degree, :generality, :vulnerability, :fill] || throw(ArgumentError("The constraint argument you specificied ($(constraint)) is invalid -- see ?shuffle! for a list."))
    first(lims) < last(lims) || throw(ArgumentError("The interval you used for limits, $(lims), is invalid"))
    all(first(lims) .<= N.A .<= last(lims)) || throw(ArgumentError("The interval you used for limits, $(lims), is smaller than the range of network values"))

    # List of compatible swap schemes under different constraints
    same_row = [[-1 1; -1 1], [-1 1; 0 0], [-1 1; 1 -1], [0 0; -1 1], [0 0; 1 -1], [1 -1; -1 1], [1 -1; 0 0], [1 -1; 1 -1],]
    same_col = [[-1 -1; 1 1], [-1 0; 1 0], [-1 1; 1 -1], [0 -1; 0 1], [0 1; 0 -1], [1 -1; -1 1], [1 0; -1 0], [1 1; -1 -1],]
    same_all = filter(x -> x ∈ same_row, same_col)
    same_lnk = unique(vcat(same_row, same_col))

    # Pick the correcty matrix series
    m = same_lnk
    constraint == :degree && (m = same_all)
    constraint == :vulnerability && (m = same_col)
    constraint == :generality && (m = same_row)
    m = m.*ε

    permuted = false
    while !permuted
        s1 = rand(1:richness(N; dims=1), 2)
        while !quanti_shuffle_valid(N.A, s1, s1)
            s1 = rand(1:richness(N; dims=1), 2)
        end
        st1 = filter(a -> all(first(lims).<=a.<=last(lims)), [N[s1,s1] .+ x for x in m])
        st2 = filter(a -> all((N[s1,s1].!=0).==(a.!=0)), st1)
        if length(st2) > 0
            permuted = true
            k = rand(st2)
            for i in 1:2, j in 1:2
                N[s1[i],s1[j]] = k[i,j]
            end
        end
    end
end

"""
    shuffle(N::T; kwargs...) where {T <: QuantitativeNetwork}

Shuffle a quantitative network -- see shuffle!
"""
function shuffle(N::T; kwargs...) where {T <: QuantitativeNetwork}
    Y = copy(N)
    @info kwargs
    shuffle!(Y; kwargs...)
    return Y
end