text:
using ABCDistances
using Distributions
import Distributions.length, Distributions._rand!, Distributions._pdf ##So that these can be extended
Libdl.dlopen("/usr/lib/liblapack.so.3", Libdl.RTLD_GLOBAL); ##Needed to avoid PyPlot problems on my work machine
##I sometimes also need run "export LD_LIBRARY_PATH=$HOME/.julia/v0.6/Conda/deps/usr/lib/:$LD_LIBRARY_PATH" - see https://github.com/JuliaPy/Conda.jl/issues/105
using PyPlot
using StatsBase
using ProgressMeter

#########################
##DEFINE MODELS AND PRIOR
#########################
##Define prior: uniform on [0,10]^4
type GKPrior <: ContinuousMultivariateDistribution
end

function length(d::GKPrior)
    4
end

function _rand!{T<:Real}(d::GKPrior, x::AbstractVector{T})
    x = 10.0*rand(4)
end

function _pdf{T<:Real}(d::GKPrior, x::AbstractVector{T})
    if all(0.0 .<= x .<= 10.0)
        return 0.0001
    else
        return 0.0
    end
end

##Define model and summary statistics
quantiles = Int[1250*i for i in 1:7];
ndataset = 10000;

function sample_sumstats(pars::Array{Float64,1})
    success = true
    stats = rgk_os(pars, quantiles, ndataset)
    (success, stats)
end

################################################
##DETAILED ANALYSIS OF A SINGLE OBSERVED DATASET
################################################
theta0 = [3.0,1.0,1.5,0.5]
srand(1);
(success, sobs) = sample_sumstats(theta0)

abcinput = ABCInput();
abcinput.prior = GKPrior();
abcinput.sample_sumstats = sample_sumstats;
abcinput.abcdist = WeightedEuclidean(sobs);
abcinput.nsumstats = length(quantiles);

##Perform ABC-PMC
srand(2);
pmcoutput_alg3V = abcPMC3V(abcinput, 1000, 1/2, 1000000);
srand(2);
pmcoutput_alg5 = abcPMC5(abcinput, 1000, 1/2, 1000000);
srand(2);
pmcoutput_alg3 = abcPMC3(abcinput, 1000, 1/2, 1000000);
abcinput.abcdist = MahalanobisEmp(sobs);
srand(2);
pmcoutput_alg5_Mahalanobis = abcPMC5(abcinput, 1000, 1/2, 1000000);
abcinput.abcdist = WeightedEuclidean(sobs);
srand(2);
pmcoutput_alg4 = abcPMC4(abcinput, 1000, 1/2, 1000000);

##Plot MSEs (and also bias^2, variance)
b1 = parameter_means(pmcoutput_alg3V);
b2 = parameter_means(pmcoutput_alg5);
b3 = parameter_means(pmcoutput_alg3);
b4 = parameter_means(pmcoutput_alg5_Mahalanobis);
b5 = parameter_means(pmcoutput_alg4);
v1 = parameter_vars(pmcoutput_alg3V);
v2 = parameter_vars(pmcoutput_alg5);
v3 = parameter_vars(pmcoutput_alg3);
v4 = parameter_vars(pmcoutput_alg5_Mahalanobis);
v5 = parameter_vars(pmcoutput_alg4);
c1 = pmcoutput_alg3V.cusims ./ 1000;
c2 = pmcoutput_alg5.cusims ./ 1000;
c3 = pmcoutput_alg3.cusims ./ 1000;
c4 = pmcoutput_alg5_Mahalanobis.cusims ./ 1000;
c5 = pmcoutput_alg4.cusims ./ 1000;
PyPlot.figure(figsize=(12,8))
pnames = ("A", "B", "g", "k")
for i in 1:4
    PyPlot.subplot(220+i)
    PyPlot.plot(c3, log10.(v3[i,:] .+ (b3[i,:]-theta0[i]).^2), "b-o")
    PyPlot.plot(c5, log10.(v5[i,:] .+ (b5[i,:]-theta0[i]).^2), "g-^")
    PyPlot.plot(c2, log10.(v2[i,:] .+ (b2[i,:]-theta0[i]).^2), "y-*")
    PyPlot.title(pnames[i])
    PyPlot.xlabel("Number of simulations (000s)")
    PyPlot.ylabel("log₁₀(MSE)")
    PyPlot.legend(["Algorithm 3","Algorithm 4", "Algorithm 5"])
end
PyPlot.tight_layout();
PyPlot.savefig("gk_mse.pdf")

PyPlot.figure()
for i in 1:4
    PyPlot.subplot(220+i)
    PyPlot.plot(c1, log10.(v1[i,:]), "r-x")
    PyPlot.plot(c2, log10.(v2[i,:]), "g-^")
    PyPlot.plot(c3, log10.(v3[i,:]), "b-o")
    PyPlot.plot(c4, log10.(v4[i,:]), "k-|")
    PyPlot.plot(c5, log10.(v5[i,:]), "y-*")
    PyPlot.axis([0,maximum(vcat(c1,c2,c3,c4,c5)),-4,1]);
    PyPlot.title(pnames[i])
    PyPlot.xlabel("Number of simulations (000s)")
    PyPlot.ylabel("log₁₀(estimated variance)")
    PyPlot.legend(["Non-adaptive (alg 4)","Adaptive (alg 5)","Non-adaptive (alg 3)", "Mahalanobis", "Alg 4"])
end

PyPlot.figure()
for i in 1:4
    PyPlot.subplot(220+i)
    PyPlot.plot(c1, log10.((b1[i,:]-theta0[i]).^2), "b-o")
    PyPlot.plot(c2, log10.((b2[i,:]-theta0[i]).^2), "g-^")
    PyPlot.plot(c3, log10.((b3[i,:]-theta0[i]).^2), "r-x")
    PyPlot.plot(c4, log10.((b4[i,:]-theta0[i]).^2), "k-|")
    PyPlot.plot(c5, log10.((b5[i,:]-theta0[i]).^2), "y-*")
    PyPlot.title(pnames[i])
    PyPlot.xlabel("Number of simulations (000s)")
    PyPlot.ylabel("log₁₀(bias squared)")
    PyPlot.legend(["Non-adaptive (alg 4)","Adaptive (alg 5)","Non-adaptive (alg 3)", "Mahalanobis", "Alg 4"])
end

##Posterior summaries
hcat(b3[:,end], b5[:,end], b2[:,end])' ##Rows are algs 3,4,5 of paper
sqrt.(hcat(v3[:,end], v5[:,end], v2[:,end]))'

##Compute weights
w1 = pmcoutput_alg3V.abcdists[1].w;
w2 = Array{Float64}(pmcoutput_alg5.niterations, length(quantiles));
for i in 1:pmcoutput_alg5.niterations
    w2[i,:] = pmcoutput_alg5.abcdists[i].w
end
w3 = pmcoutput_alg3.abcdists[1].w;
w5 = Array{Float64}(pmcoutput_alg4.niterations-1, length(quantiles));
for i in 2:pmcoutput_alg4.niterations
    w5[i-1,:] = pmcoutput_alg4.abcdists[i].w
end

##Plot weights
PyPlot.figure(figsize=(12,4))
PyPlot.plot(quantiles, w3/sum(w3), "b-o")
wlast = w5[pmcoutput_alg4.niterations-1, :] ##IS -1 NEEDED?!?!
PyPlot.plot(quantiles, wlast/sum(wlast), "g-^")
wlast = w2[pmcoutput_alg5.niterations, :]
PyPlot.plot(quantiles, wlast/sum(wlast), "y-*")
##PyPlot.axis([1.0,9.0,0.0,0.35]) ##Sometimes needed to fit legend in
PyPlot.legend(["Algorithm 3", "Algorithm 4\n(last iteration)", "Algorithm 5\n(last iteration)"])
PyPlot.xlabel("Order statistic")
PyPlot.ylabel("Relative weight")
PyPlot.tight_layout();
PyPlot.savefig("gk_weights.pdf")

##Marginal posterior plots
PyPlot.figure();
samplesABC = (pmcoutput_alg3.parameters[:,:,pmcoutput_alg3.niterations], pmcoutput_alg4.parameters[:,:,pmcoutput_alg4.niterations], pmcoutput_alg5.parameters[:,:,pmcoutput_alg5.niterations]);
weightsABC = (pmcoutput_alg3.weights[:,pmcoutput_alg3.niterations], pmcoutput_alg4.weights[:,pmcoutput_alg4.niterations], pmcoutput_alg5.weights[:,pmcoutput_alg5.niterations]);
for i in 1:3 ##Loop over algorithms
    ww = weightsABC[i]
    for j in 1:4 ##Loop over parameters
        pp = samplesABC[i][j,:]
        PyPlot.subplot(140+j)
        PyPlot.plt[:hist](pp, weights=ww, normed=true, alpha=0.5)
    end
end

##Investigate best choice of alpha (as requested by reviewers)
function get_mse(pars::Array{Float64, 2}, w::Array{Float64, 1})
    mse = 0.0
    for i in 1:size(pars)[2]
        mse += sum((pars[:, i] .- theta0).^2) * w[i]
    end
    mse
end

alphas = 0.05:0.1:0.95;
MSEs_alg5 = zeros(alphas);
MSEs_alg4 = zeros(alphas);
MSEs_alg3 = zeros(alphas);
abcinput.abcdist = WeightedEuclidean(sobs);
srand(1);
for i in 1:length(alphas)
    x = abcPMC5(abcinput, 1000, alphas[i], 10^6)
    nits = x.niterations
    MSEs_alg5[i] = (nits == 0) ? Inf : get_mse(x.parameters[:,:,nits], x.weights[:,nits])
    x = abcPMC3(abcinput, 1000, alphas[i], 10^6)
    nits = x.niterations
    MSEs_alg3[i] = (nits == 0) ? Inf : get_mse(x.parameters[:,:,nits], x.weights[:,nits])
    x = abcPMC4(abcinput, 1000, alphas[i], 10^6)
    nits = x.niterations
    MSEs_alg4[i] = (nits == 0) ? Inf : get_mse(x.parameters[:,:,nits], x.weights[:,nits])
end

PyPlot.figure();
plot(alphas, log10.(MSEs_alg5), "r-x");
plot(alphas, log10.(MSEs_alg3), "b-o");
plot(alphas, log10.(MSEs_alg4), "g-*");
##MSE has large flat minimum around alpha in [0.3, 0.7]

###############################
##ANALYSIS OF MULTIPLE DATASETS
###############################
ndatasets = 100;
trueθs = zeros((4, ndatasets));
RMSEs = zeros((4, 5, ndatasets)); ##Indices are: parameters, method, dataset
vars =  zeros((4, 5, ndatasets,));
squaredbiases = zeros((4, 5, ndatasets));

##Returns squared bias, variance and RMSE of weighted posterior sample wrt true parameters
function getError(s::ABCPMCOutput, pobs::Array{Float64, 1})
    n = s.niterations
    p = s.parameters[:,:,n]
    wv = Weights(s.weights[:,n])
    bias2 = (mean(p, wv, 2) - pobs).^2
    bias2 = vec(bias2)
    v = var(p, wv, 2, corrected=false)
    v = vec(v)
    rmse = sqrt.(bias2 + v)
    (bias2, v, rmse)
end

abcinput = ABCInput();
abcinput.prior = GKPrior();
abcinput.sample_sumstats = sample_sumstats;
abcinput.nsumstats = length(quantiles);

srand(1);
prog = Progress(5*ndatasets, 1); ##Progress meter
for i in 1:ndatasets
    theta0 = rand(GKPrior())
    (success, sobs) = sample_sumstats(theta0)
    abcinput.abcdist = WeightedEuclidean(sobs)
    pmcoutput_alg3V = abcPMC3V(abcinput, 1000, 1/2, 1000000, silent=true)
    next!(prog)
    pmcoutput_alg5 = abcPMC5(abcinput, 1000, 1/2, 1000000, silent=true)
    next!(prog)
    pmcoutput_alg3 = abcPMC3(abcinput, 1000, 1/2, 1000000, silent=true)
    next!(prog)
    abcinput.abcdist = MahalanobisEmp(sobs)
    pmcoutput_alg5_Mahalanobis = abcPMC5(abcinput, 1000, 1/2, 1000000, silent=true)
    next!(prog)
    abcinput.abcdist = WeightedEuclidean(sobs)
    pmcoutput_alg4 = abcPMC4(abcinput, 1000, 1/2, 1000000, silent=true)
    next!(prog)

    trueθs[:,i] = theta0
    (squaredbiases[:,1,i], vars[:,1,i], RMSEs[:,1,i]) = getError(pmcoutput_alg3V, theta0)
    (squaredbiases[:,2,i], vars[:,2,i], RMSEs[:,2,i]) = getError(pmcoutput_alg5, theta0)
    (squaredbiases[:,3,i], vars[:,3,i], RMSEs[:,3,i]) = getError(pmcoutput_alg3, theta0)
    (squaredbiases[:,4,i], vars[:,4,i], RMSEs[:,4,i]) = getError(pmcoutput_alg5_Mahalanobis, theta0)
    (squaredbiases[:,5,i], vars[:,5,i], RMSEs[:,5,i]) = getError(pmcoutput_alg4, theta0)
end

##Summarise output
mean(RMSEs, 3)
mean(vars, 3)
mean(squaredbiases, 3)

##Save output for further analysis without lengthy rerun
writedlm("gk_RMSE.txt", RMSEs)
writedlm("gk_vars.txt", vars)
writedlm("gk_bias2.txt", squaredbiases)
writedlm("gk_thetas.txt", trueθs)


text:
quad1DW = [0.2777777777777778, 0.4444444444444444, 0.2777777777777778];

quad1DX = [0.1127016653792583, 0.5, 0.8872983346207417];

#quadX = [[0.659027622374092,  0.231933368553031],
#         [0.659027622374092,  0.109039009072877],
#         [0.231933368553031,  0.659027622374092],
#         [0.231933368553031,  0.109039009072877],
#         [0.109039009072877,  0.659027622374092],
#         [0.109039009072877,  0.231933368553031]];

#quadW = 1.0/12.0*ones(6,1);

quadX = [[1.88409405952072e-01,   7.87659461760847e-01],
         [5.23979067720101e-01,   4.09466864440735e-01],
         [8.08694385677670e-01,   8.85879595127039e-02],
         [1.06170269119576e-01,   7.87659461760847e-01],
         [2.95266567779633e-01,   4.09466864440735e-01],
         [4.55706020243648e-01,   8.85879595127039e-02],
         [2.39311322870805e-02,   7.87659461760847e-01],
         [6.65540678391645e-02,   4.09466864440735e-01],
         [1.02717654809626e-01,   8.85879595127039e-02]];

quadW = [1.93963833059595e-02,
         6.36780850998851e-02,
         5.58144204830443e-02,
         3.10342132895352e-02,
         1.01884936159816e-01,
         8.93030727728709e-02,
         1.93963833059595e-02,
         6.36780850998851e-02,
         5.58144204830443e-02];

PDESystem(;A=spzeros(0,0), b=[], bc=[], DI=Set{Int64}(), vec_ind=Set{Int64}(), qdim=1, Factors=[],
          SystemMatrix=spzeros(0,0), B=spzeros(0,0), state=[], rhs=[]) =
            PDESystem(A, b, bc, DI, vec_ind, qdim, Factors, SystemMatrix, B, state, rhs)

function assembly(S::PDESystem)
  if S.Factors == []
    if S.DI != Set{Int64}()
      S.vec_ind = Set{Int64}()
      for d=1:S.qdim
        union!(S.vec_ind, S.qdim*(collect(S.DI).-1).+d)
      end
      ii = length(S.vec_ind)
      m,n = size(S.A)
      S.B = getDirichletProjection(m, S.DI, qdim=S.qdim)

      S.SystemMatrix = [S.A            S.B';
                        S.B  spzeros(ii,ii)]
    else
      S.SystemMatrix = S.A
    end
    S.Factors = lu(S.SystemMatrix)
  end
end

function refresh(S::PDESystem)
  S.Factors = []
  assembly(S)
end

"""
  `solve(S::PDESystem)`

First tries to set up the system matrix with multipliers for Dirichlet conditions.
If the system has already been used before, this step is skipped.
This is determined depending on an existing factorization of the system matrix.
If the stiffness matrix or Dirichlet conditions have changes, one should invole refresh() first.

Finally the system is solved via matrix factorization.
"""
function solve(S::PDESystem)
  assembly(S)
  S.rhs = [S.b; S.B*S.bc]
  S.state = (S.Factors\S.rhs)[1:length(S.b)]
end

"""
  `getDirichletProjection(jj::Int64, DI::Set{Int64};qdim=1)`

Build the projection onto the Dirichlet nodes where the input jj is understood as qdim*nnodes.
"""

function getDirichletProjection(jj::Int64, DI::Set{Int64}; qdim=1)
  vec_ind = Set{Int64}()
  for d=1:qdim
    union!(vec_ind, qdim*(collect(DI).-1).+d)
  end
  ii = length(vec_ind)
  return sparse(1:ii, collect(vec_ind), ones(ii), ii, jj)
end

"""
  `evaluateMeshFunction(mesh::Mesh, f::Function; region=[], qdim=1)`

Evaluate a given function object f on all nodes of the mesh. Or, if region is specified,
it is assumed to be a set or vector of node indices, where f should be evaluated.
"""
function evaluateMeshFunction(mesh::Mesh, f::Function; region=[], qdim=1)
  v = zeros(Float64, qdim*mesh.nnodes)
  nodes = 1:mesh.nnodes
  if region != []
    nodes = collect(region)
  end
  if(qdim==1)
    for i in nodes
      v[i] = f(mesh.Nodes[i])
    end
  else
    for i in nodes
      v[qdim*(i-1)+1:qdim*i] = f(mesh.Nodes[i])
    end
  end
  return v
end

function Phi(ii::Int64, x::AbstractVector)
  if ii==1
    return 1.0-x[1]-x[2]
  elseif ii==2
    return x[1]
  elseif ii==3
    return x[2]
  end
end

function PhiEdge(ii::Int64, x::Float64)
  if ii==1
    return 1.0-x
  elseif ii==2
    return x
  end
end

function gradPhi(ii::Int64)
  grad = zeros(2,1)
  if ii==1
    grad = [-1.0; -1.0]
  elseif ii==2
    grad = [1.0; 0.0]
  elseif ii==3
    grad = [0.0; 1.0]
  end
  return grad
end

"""
  `function asmDirichletCondition(SM, DI::Set{Int64}, rhs=[], bc=[]; qdim=1, insert=1.0)`

Modify a stiffness matrix and a right hand side according to the given Dirichlet conditions.
DI has to be the set of node indices for which the condition should be active.
For vector valued states either DI can be set to each component that should have a
Dirichlet condtion or qdim is set, if all components should have the condition.
The value insert is put as diagonal element. Usually you want a 1.0 here.
"""
function asmDirichletCondition(SM, DI::Set{Int64}; rhs=[], bc=[], qdim=1, insert=1.0)
  if rhs != [] && bc != []
    for i in DI
      for d=1:qdim
        ii = qdim*(i-1)+d
        bcind = SM[:, ii].nzind
        rhs[bcind] -= SM[bcind, ii]*bc[ii]
      end
    end
    for i in DI
      for d=1:qdim
        ii = qdim*(i-1)+d
        rhs[ii] = bc[ii]
      end
    end
  end
  for i in DI
    for d=1:qdim
      ii = qdim*(i-1)+d
      SM[ii, SM[ii,:].nzind] *= 0.0;
      SM[:, ii] *= 0.0;
      SM[ii, ii] = insert;
    end
  end
end

function asmSparseMatrix(D::Dict{Tuple{Int64,Int64}, Float64}; nrows=0, ncols=0)
  nn = length(D)
  II = zeros(Int64, nn)
  JJ = zeros(Int64, nn)
  AA = zeros(Float64, nn)
  for (i, pair) in enumerate(D)
    II[i] = pair[1][1]
    JJ[i] = pair[1][2]
    AA[i] = pair[2]
  end
  return sparse(II,JJ,AA, maximum([nrows,maximum(II)]), maximum([ncols,maximum(JJ)]))
end

"""
  `asmLaplacian(mesh::Mesh)`

Assemble the Laplacian stiffness matrix for all elements in the mesh.
"""
function asmLaplacian(mesh::Mesh)

  AA = zeros(Float64, mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(3,3)
    for i=1:3
      for j=1:3
        elemMat[i,j] += (J*gradPhi(i))'*(J*gradPhi(j))
      end
    end
    elemMat *= sum(quadW)*detJ

    for i=1:3
      for j=1:3
        n = n+1
        II[n] = nodes[i]
        JJ[n] = nodes[j]
        AA[n] = elemMat[i,j]
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `asmMassMatrix(mesh::Mesh; qdim=1)`

Assemble a mass matrix for all elements of the given mesh.
"""
function asmMassMatrix(mesh::Mesh; qdim=1)

  AA = zeros(Float64, qdim^2 * mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(3,3)
    for i=1:3
      for j=1:3
        for (q, x) in enumerate(quadX)
          elemMat[i,j] += Phi(i, x) * Phi(j, x) * quadW[q] * detJ
        end
      end
    end

    for d=1:qdim
      for i=1:3
        for j=1:3
          n = n+1
          II[n] = qdim*(nodes[i]-1)+d
          JJ[n] = qdim*(nodes[j]-1)+d
          AA[n] = elemMat[i,j]
        end
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `asmBoundaryMassMatrix(mesh::Mesh; Edges=Set{Int64}(), qdim=1)`

Assemble a mass matrix for the given set of boundary edges.
"""
function asmBoundaryMassMatrix(mesh::Mesh; Edges=Set{Int64}(), qdim=1)

  AA = zeros(Float64, qdim^2 * mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  if(isempty(Edges))
    Edges = 1:mesh.nedges
  end

  for el in Edges
    edge = mesh.Edges[el]
    detJ = EdgeJacobian(mesh, el)
    elemMat = zeros(2,2)
    for i=1:2
      for j=1:2
        for (q, x) in enumerate(quad1DX)
          elemMat[i,j] += PhiEdge(i, x) * PhiEdge(j, x) * quad1DW[q] * detJ
        end
      end
    end

    for d=1:qdim
      for i=1:2
        for j=1:2
          n = n+1
          II[n] = qdim*(edge[i]-1)+d
          JJ[n] = qdim*(edge[j]-1)+d
          AA[n] = elemMat[i,j]
        end
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n], qdim*mesh.nnodes, qdim*mesh.nnodes)
end

"""
  `asmBoundarySource(mesh::Mesh, S::Array{Float64,1}, BoundaryEdges=Set{Int64}(-1); qdim=1)`

Assemble the source S on the given set of boundary edges.
"""
function asmBoundarySource(mesh::Mesh, S::Array{Float64,1}, BoundaryEdges=Set{Int64}(-1); qdim=1)
  V = zeros(mesh.nnodes)

  if(in(-1, BoundaryEdges))
    BoundaryEdges = 1:mesh.nedges
  end

  for el in BoundaryEdges
    edge = mesh.Edges[el]
    detJ = EdgeJacobian(mesh, el)
    for i=1:2
      for j=1:2
        for (q, x) in enumerate(quad1DX)
          V[edge[i]] += S[edge[j]]*PhiEdge(j, x) * PhiEdge(i, x) * quad1DW[q] * detJ
        end
      end
    end
  end

  return V
end

"""
  `L2norm(Mass::SparseMatrixCSC{Float64,Int64}, v::AbstractVector; qdim=1)`

Computes the L2 norm of a vector v according to the mass matrix Mass.
"""
function L2norm(Mass::SparseMatrixCSC{Float64,Int64}, v::AbstractVector; qdim=1)
  return sqrt(v'*Mass*v)
end

"""
  `asmCubicTerm(mesh::Mesh, y::AbstractVector)`

The cubic term of the standard semilinear parabolic equation.
"""
function asmCubicTerm(mesh::Mesh, y::AbstractVector)
  V = zeros(mesh.nnodes)
  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    y_cubic = zeros(length(quadW))
    for (q, x) in enumerate(quadX)
      y_quad = 0
      for i=1:3
        y_quad += y[nodes[i]]*Phi(i, x)
      end
      y_cubic[q] = y_quad^3
    end

    for i=1:3
      for (q, x) in enumerate(quadX)
        V[nodes[i]] += y_cubic[q] * Phi(i, x) * quadW[q] * detJ
      end
    end

  end

  return V
end

"""
  `asmCubicDerivativeMatrix(mesh::Mesh, y::AbstractVector)`

Assembly of the linearization of the cubic term of the standard semilinear elliptic equation.
"""
function asmCubicDerivativeMatrix(mesh::Mesh, y::AbstractVector)

  AA = zeros(Float64, mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    y_quadratic = zeros(length(quadW))
    for (q, x) in enumerate(quadX)
      y_quad = 0
      for i=1:3
        y_quad += y[nodes[i]]*Phi(i, x)
      end
      y_quadratic[q] = y_quad^2
    end

    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(3,3)
    for i=1:3
      for j=1:3
        for (q, x) in enumerate(quadX)
          elemMat[i,j] += 3.0* y_quadratic[q] * Phi(i, x) * Phi(j, x) * quadW[q] * detJ
        end
      end
    end

    for i=1:3
      for j=1:3
        n = n+1
        II[n] = nodes[i]
        JJ[n] = nodes[j]
        AA[n] = elemMat[i,j]
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `asmCubicSecondDerivativeMatrix(mesh::Mesh, y::AbstractVector, p::AbstractVector)`

Assembly of the second derivative of the cubic term of the standard semilinear elliptic equation
around the state y.
"""
function asmCubicSecondDerivativeMatrix(mesh::Mesh, y::AbstractVector, p::AbstractVector)

  AA = zeros(Float64, mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    y_quad = zeros(length(quadW))
    p_quad = zeros(length(quadW))
    for i=1:3
      for (q, x) in enumerate(quadX)
        y_quad[q] += y[nodes[i]]*Phi(i, x)
        p_quad[q] += p[nodes[i]]*Phi(i, x)
      end
    end

    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(3,3)
    for i=1:3
      for j=1:3
        for (q, x) in enumerate(quadX)
          elemMat[i,j] += 6.0 * y_quad[q] * p_quad[q] * Phi(i, x) * Phi(j, x) * quadW[q] * detJ
        end
      end
    end

    for i=1:3
      for j=1:3
        n = n+1
        II[n] = nodes[i]
        JJ[n] = nodes[j]
        AA[n] = elemMat[i,j]
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `strainTensor(grad::AbstractMatrix)`

Strain tensor for linear elasticity with gradient as argument.
"""
function strainTensor(grad::AbstractMatrix)
  return 0.5*(grad + grad')
end

"""
  `stressTensor(grad::AbstractMatrix, lambda::Float64, mu::Float64)`

Strain tensor for linear elasticity with gradient as argument.
"""
function stressTensor(grad::AbstractMatrix, lambda::Float64, mu::Float64)
  epsilon = strainTensor(grad)
  return lambda * tr(epsilon) * [1.0 0.0;0.0 1.0] + 2.0 * mu * epsilon
end

"""
  `asmElasticity(mesh::Mesh, lambda::Float64, mu::Float64)`

Assembly of the linear elasticity stiffness matrix for constant coefficients lambda and mu.
"""
function asmElasticity(mesh::Mesh, lambda::Float64, mu::Float64)
  qdim = 2

  AA = zeros(Float64, qdim^2 * mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(qdim*3,qdim*3)
    for i=1:3
      for j=1:3
        for ic=1:qdim
          for jc=1:qdim
            grad_i = zeros(2,2)
            grad_i[:,ic] = J*gradPhi(i)
            grad_j = zeros(2,2)
            grad_j[:,jc] = J*gradPhi(j)
            elemMat[qdim*(i-1)+ic,qdim*(j-1)+jc] += dot(stressTensor(grad_j,lambda,mu), strainTensor(grad_i))
          end
        end
      end
    end
    elemMat *= sum(quadW)*detJ

    for i=1:3
      for j=1:3
        for ic=1:qdim
          for jc=1:qdim
            n = n+1
            II[n] = qdim*(nodes[i]-1)+ic
            JJ[n] = qdim*(nodes[j]-1)+jc
            AA[n] = elemMat[qdim*(i-1)+ic,qdim*(j-1)+jc]
          end
        end
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `computeGradient(mesh::Mesh, y::AbstractVector; qdim=1)`

Compute the gradient of a state y on a given mesh.
"""
function computeGradient(mesh::Mesh, y::AbstractVector; qdim=1)
  grad = zeros(2*qdim*mesh.nelems)
  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)

    for q=1:qdim
      g = view(grad, (1:2).+(2*qdim*(el-1) + 2*(q-1)))
      for i=1:3
        g[:] += y[qdim*(nodes[i]-1)+q]*(J*gradPhi(i))
      end
    end
  end
  return grad
end

"""
  `asmGradient(mesh::Mesh; qdim=1)`

Assembles the linear mapping from a state on the given mesh to the gradient.
"""
function asmGradient(mesh::Mesh; qdim=1)

  AA = zeros(Float64, qdim^2 * mesh.nelems * 3 * 2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    G = zeros(2*qdim, 3*qdim)

    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)

    for q=1:qdim
      for i=1:3
        g = view(G, [2*(q-1)+1, 2*(q-1)+2], qdim*(i-1)+q)
        g[:] += J*gradPhi(i)
      end
    end

    for i=1:2
      for j=1:3
        for ic=1:qdim
          for jc=1:qdim
            n=n+1
            II[n] = 2*qdim*(el-1)+qdim*(i-1)+ic
            JJ[n] = qdim*(nodes[j]-1)+jc
            AA[n] = G[qdim*(i-1)+ic,qdim*(j-1)+jc]
          end
        end
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n], 2*qdim*mesh.nelems, qdim*mesh.nnodes)
end

"""
  `asmIncompressibility(mesh::Mesh)`

Assembles the matrix for a linear incompressibility block in a saddle point system.
"""
function asmIncompressibility(mesh::Mesh)

  AA = zeros(Float64, mesh.nelems * 3^2 * 2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(2*3,3)
    for i=1:3
      for j=1:3
        for (q, x) in enumerate(quadX)
          elemMat[2*(i-1).+[1;2],j] += Phi(j, x) * (J*gradPhi(i)) * quadW[q] * detJ
        end
      end
    end

    for d=1:2
      for i=1:3
        for j=1:3
          n = n+1
          II[n] = 2*(nodes[i]-1)+d
          JJ[n] = nodes[j]
          AA[n] = elemMat[2*(i-1)+d,j]
        end
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `asmVectorLaplacian(mesh::Mesh; qdim=1)`

Assemble the vector-valued Laplacian stiffness matrix for all elements in the mesh.
"""
function asmVectorLaplacian(mesh::Mesh; qdim=1)

  AA = zeros(Float64, qdim^2 * mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(3,3)
    for i=1:3
      for j=1:3
        elemMat[i,j] += (J*gradPhi(i))'*(J*gradPhi(j))
      end
    end
    elemMat *= sum(quadW)*detJ

    for d=1:qdim
      for i=1:3
        for j=1:3
          n = n+1
          II[n] = qdim*(nodes[i]-1)+d
          JJ[n] = qdim*(nodes[j]-1)+d
          AA[n] = elemMat[i,j]
        end
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

"""
  `function asmStokesStabilization(mesh::Mesh)`

Assemble the vector-valued Laplacian stiffness matrix for all elements in the mesh.
"""
function asmStokesStabilization(mesh::Mesh)

  AA = zeros(Float64, mesh.nelems * 3^2)
  II = zeros(Int64, length(AA))
  JJ = zeros(Int64, length(AA))
  n = 0
  diam = MinFEM.getCellDiameter(mesh)

  for el=1:mesh.nelems
    nodes = mesh.Triangles[el]
    (detJ, J) = Jacobian(mesh, el)
    elemMat = zeros(3,3)
    d = diam[el]
    for i=1:3
      for j=1:3
        elemMat[i,j] += d^2 * (J*gradPhi(i))'*(J*gradPhi(j))
      end
    end
    elemMat *= sum(quadW)*detJ

    for i=1:3
      for j=1:3
        n = n+1
        II[n] = nodes[i]
        JJ[n] = nodes[j]
        AA[n] = elemMat[i,j]
      end
    end
  end

  return sparse(II[1:n],JJ[1:n],AA[1:n])
end

ext: jl


text:
using Pkg.Artifacts, ghr_jll
using Pkg.BinaryPlatforms
import ZipFile

release_tag = "v0.3.1"

artifact_toml = joinpath(@__DIR__, "Artifacts.toml")

function unzip(zf_path, dest_path)
    r = ZipFile.Reader(zf_path)
    k_path = joinpath(dest_path, "kaleido")
    rm(k_path , force=true, recursive=true)
    mkdir(k_path)
    for f in r.files
        @show f
        out_path = joinpath(dest_path, f.name)
        if (endswith(f.name, "/") || endswith(f.name, "\\"))
            !isdir(out_path) && mkdir(out_path)
        else
            parent_dir = joinpath(dest_path, dirname(f.name))
            if !isdir(parent_dir)
                mkpath(parent_dir)
            end
            v = Vector{UInt8}(undef, f.uncompressedsize)
            read!(f, v)
            write(out_path, v)
        end
    end
    close(r)
end

_download_os(::Linux) = "linux"
_download_os(::MacOS) = "mac"
_download_os(::Windows) = "win"

function _download_os()
    if Sys.isunix()
        return Sys.isapple() ? "mac" : "linux"
    end
    if Sys.iswindows()
        return "win"
    end
    error("Only mac, linux, and windows are supported")
end

function download_kaleido(k_dir, os=_download_os())
    url = "https://github.com/plotly/Kaleido/releases/download/v0.0.3/kaleido_$os-0.0.3.zip"
    dest = joinpath(k_dir, "kaleido.zip")
    download(url, dest)
    unzip(dest, k_dir)
    rm(dest)  # remove zip file
    if Sys.isunix() && (os ∈ ["mac", "linux"])
        ex1_path = joinpath(k_dir, "kaleido", "kaleido")
        run(`chmod +x $ex1_path`)

        ex2_path = joinpath(k_dir, "kaleido", "bin", "kaleido")
        run(`chmod +x $ex2_path`)
    end
    return
end

for platform in [Linux(:x86_64), MacOS(:x86_64), Windows(:x86_64)]
# for platform in [MacOS(:x86_64)]
    @show platform
    kaleido_hash = artifact_hash("kaleido", artifact_toml, platform=platform)
    # If the name was not bound, or the hash it was bound to does not exist, create it!
    if kaleido_hash === nothing || !artifact_exists(kaleido_hash)
        # create_artifact() returns the content-hash of the artifact directory once we're finished creating it
        kaleido_hash = create_artifact() do k_dir
            download_kaleido(k_dir, _download_os(platform))
        end
    end

    local tarball_hash
    gz_name = "kaleido_$(_download_os(platform)).tar.gz"
    mktempdir() do dir
        gz_path = joinpath(dir, gz_name)
        tarball_hash = archive_artifact(kaleido_hash, gz_path)
        ghr() do ghr
            run(`$(ghr) -u sglyon -r PlotlyBase.jl -replace $(release_tag) $(gz_path)`)
        end
    end
    bind_artifact!(
        artifact_toml,
        "kaleido",
        kaleido_hash;
        platform=platform,
        force=true,
        download_info = [
            ("https://github.com/sglyon/PlotlyBase.jl/releases/download/$(release_tag)/$(gz_name)", tarball_hash),
        ],
    )
end


text:
# TODO: Add tests for came configurations

using REPL.TerminalMenus
export configure_game

function cofigure_n_players()
    n_players = nothing
    while true
        println("Enter number of players")
        n_players = readline()
        try
            n_players = parse(Int, n_players)
            2 ≤ n_players ≤ 10 && break
            println("Number of players must be between 2 and 10 (inclusive)")
        catch
            println("Number of players must be an integer between 2 and 10 (inclusive)")
        end
    end
    @assert n_players ≠ nothing
    @info "$n_players playing"
    return n_players
end

function cofigure_bank_roll(blinds)
    bank_roll = nothing
    while true
        println("Enter bank roll:")
        bank_roll = readline()
        try
            bank_roll = parse(Float64, bank_roll)
            bank_roll ≥ blinds.big && break
            println("Bank roll must be larger than or equal to the big blind")
        catch
            println("Bank roll must be a Float64 larger than or equal to the big blind")
        end
    end
    @assert bank_roll ≠ nothing
    @info "Players starting with a bank roll of \$$bank_roll"
    return bank_roll
end

function cofigure_blinds()
    big_blinds = nothing
    while true
        println("Enter big blind amount:")
        big_blinds = readline()
        try
            big_blinds = parse(Float64, big_blinds)
            big_blinds > sqrt(eps(Float64)) && break
            println("Blinds must be larger than `sqrt(eps(Float64))`")
        catch
            println("Blinds must be a float larger than `sqrt(eps(Float64))`")
        end
    end
    @assert big_blinds ≠ nothing
    @info "Blinds: $(big_blinds/2), $(big_blinds)"
    return Blinds(big_blinds/2, big_blinds)
end

function configure_basic_heads_up_game()
    bank_roll = 100
    blinds = Blinds(1, 2)
    players = (
        Player(Human(), 1; bank_roll=bank_roll),
        Player(ai_to_use(), 2; bank_roll=bank_roll)
    )
    return Game(players; blinds=blinds)
end

function configure_basic_1v4_game()
    bank_roll = 100
    blinds = Blinds(1, 2)
    players = (
        Player(Human(), 1; bank_roll=bank_roll),
        Player(ai_to_use(), 2; bank_roll=bank_roll),
        Player(ai_to_use(), 3; bank_roll=bank_roll),
        Player(ai_to_use(), 4; bank_roll=bank_roll),
        Player(ai_to_use(), 5; bank_roll=bank_roll),
    )
    return Game(players; blinds=blinds)
end

function configure_basic_2_bots_game()
    bank_roll = 100
    blinds = Blinds(1, 2)
    players = (
        Player(ai_to_use(), 1; bank_roll=bank_roll),
        Player(ai_to_use(), 2; bank_roll=bank_roll),
    )
    return Game(players; blinds=blinds)
end

function configure_basic_4_bots_game()
    bank_roll = 100
    blinds = Blinds(1, 2)
    players = (
        Player(ai_to_use(), 1; bank_roll=bank_roll),
        Player(ai_to_use(), 2; bank_roll=bank_roll),
        Player(ai_to_use(), 3; bank_roll=bank_roll),
        Player(ai_to_use(), 4; bank_roll=bank_roll),
    )
    return Game(players; blinds=blinds)
end

function configure_human_players(n_players)
    options = string.(1:n_players)
    menu = MultiSelectMenu(options) # charset=:unicode is not supported in earlier Julia versions
    choices = request("Select which players are human:", menu)
    println("$(length(choices)) human players ($(join(sort(collect(choices)), ", ")))")
    length(choices) > 0 || println("No human players.")
    human_seat_numbers = collect(choices)
    return human_seat_numbers
end

function configure_custom_game()

    n_players = cofigure_n_players()
    blinds = cofigure_blinds()
    human_seat_numbers = configure_human_players(n_players)
    bank_roll = cofigure_bank_roll(blinds)

    players = ntuple(n_players) do i
        if i in human_seat_numbers
            Player(Human(), i; bank_roll=bank_roll)
        else
            Player(ai_to_use(), i; bank_roll=bank_roll)
        end
    end

    return Game(players; blinds=blinds)
end

function configure_game()
    options = [
        "1) Heads up-- you vs AI, \$1-\$2 blinds, \$100 bank roll",
        "2) 1v4-- you vs 4 AI, \$1-\$2 blinds, \$100 bank roll",
        "3) 2 bots!, \$1-\$2 blinds, \$100 bank roll",
        "4) 4 bots!, \$1-\$2 blinds, \$100 bank roll",
        "5) custom game"
    ]
    menu = RadioMenu(options, pagesize=4)
    choice = request("Configure game:", menu)
    if choice != -1
        choice == 1 && return configure_basic_heads_up_game()
        choice == 2 && return configure_basic_1v4_game()
        choice == 3 && return configure_basic_2_bots_game()
        choice == 4 && return configure_basic_4_bots_game()
        choice == 5 && return configure_custom_game()
        error("Uncaught case")
    else
        println("Menu canceled.")
    end
end


text:
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

function code_div(b1::Int, b2::Int, res::Int)
    code = @dice begin
        function uniform(b::Int, w::Int) # b is the bits for uniform, w is the bitwidth
            x = Vector(undef, b)
            for i = b:-1:1
                x[i] = flip(0.5)
            end
            return add_bits(DistInt(x), w - b)
        end
        a = (uniform(b1, b1+1) + 1)[1]
        b = (uniform(b2, b2))
        y = (a/b)
        # println(y)
        # y[1]
        prob_equals(y[1], res) & !y[2]
    end
    code
end

# BDD analysis
code = code_div(1, 1, 1)
bdd = compile(code)
# println(infer(code, :bdd))
@assert infer(code, :bdd) ≈ 0.25

code = code_div(2, 1, 2)
bdd = compile(code)
@assert infer(code, :bdd) ≈ 0.125

code = code_div(2, 1, 4)
bdd = compile(code)
@assert infer(code, :bdd) ≈ 0.125

code = code_div(1, 2, 0)
bdd = compile(code)
@assert infer(code, :bdd) ≈ 0.375

code = code_div(1, 2, 1)
bdd = compile(code)
@assert infer(code, :bdd) ≈ 0.25

code = code_div(1, 2, 2)
bdd = compile(code)
@assert infer(code, :bdd) ≈ 0.125
# #         bdd = compile(code)
# # num_flips(bdd)
# # num_nodes(bdd)
# # @assert infer(code, :bdd) == 0.0625

# IR analysis
# to_dice_ir(code)
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)


text:
@testset "Strings" begin
    @testset "Strings: Is Palindrome" begin
        s = "AA"     # A palindrome
        @test is_palindrome(s) == true
        s = "ABC"     # Not a palindrome
        @test is_palindrome(s) == false
        s = "Statistics"    # Not a palindrome
        @test is_palindrome(s) == false
        s = "Stats"    # A palindrome
        @test is_palindrome(s) == true
        s = "Racecar"      # A palindrome
        @test is_palindrome(s) == true
        s = "Hello"      # Not a palindrome
        @test is_palindrome(s) == false
    end

    @testset "Strings: Detect Anagrams" begin
        s = "anagram"    # No matches
        samples = ["grammar", "tutor", "kilogram"]
        @test detect_anagrams(s, samples) == []
        s = "nose"  # Matches words with different case
        samples = ["Eons", "ONES"]
        @test detect_anagrams(s, samples) == ["Eons", "ONES"]
        s = "listen"  # Does not include the same word
        samples = ["LiSten", "enlist"]
        @test detect_anagrams(s, samples) == ["enlist"]
    end

    @testset "Strings: KMP SubString Search" begin
        # contain sub string
        string = "seems like kmp search works"
        sub_string = "seems like kmp search"
        @test contain_substring_with_kmp(string, sub_string, false) == true
        string = "ABABDABACDABABCABAB"
        sub_string = "ABABCABAB"
        @test contain_substring_with_kmp(string, sub_string, false) == true
        string = "abcxabcdabcdabcy"
        sub_string = "abcdabcy"
        @test contain_substring_with_kmp(string, sub_string, false) == true
        # does not contain sub string
        string = "ABABDABACDABABCABAB"
        sub_string = "X"
        @test contain_substring_with_kmp(string, sub_string, false) == false
        string = "abcxabcdabcdabcy"
        sub_string = "so this"
        @test contain_substring_with_kmp(string, sub_string, false) == false
        # ignore cases
        string = "SeEms Like IgNOrE CaSe Work"
        sub_string = "seems like ignore case work"
        @test contain_substring_with_kmp(string, sub_string, true) == true
        string = "123 456 AbC"
        sub_string = "123 456 ABC"
        @test contain_substring_with_kmp(string, sub_string, true) == true
        string = "abcdefg"
        sub_string = "ABCDEFG"
        @test contain_substring_with_kmp(string, sub_string, true) == true
    end

    @testset "Strings: Pangrams" begin
        input = "Pack my box with five dozen liquor jugs"    # Pangram
        @test ispangram(input) == true
        input = "Hello World"    # Not Pangram
        @test ispangram(input) == false
        input = "The quick brown fox jumps over the lazy dog"    # Pangram
        @test ispangram(input) == true
    end

    @testset "Strings: word_count" begin
        sentence = "The quick brown fox jumps over the lazy dog"    # 8 different words
        @test word_count(sentence) == Dict{Any,Any}(
            "jumps" => 1,
            "the" => 2,
            "brown" => 1,
            "over" => 1,
            "quick" => 1,
            "lazy" => 1,
            "dog" => 1,
            "fox" => 1,
        )
        sentence = "Hello World"    # 2 different words
        @test word_count(sentence) == Dict{Any,Any}("hello" => 1, "world" => 1)
        sentence = "the sky is blue and beautiful"    # 6 different words
        @test word_count(sentence) == Dict{Any,Any}(
            "and" => 1,
            "the" => 1,
            "sky" => 1,
            "blue" => 1,
            "is" => 1,
            "beautiful" => 1,
        )
    end

    @testset "Strings: hamming_distance" begin
        input1 = "karolin"
        input2 = "kathrin"
        @test hamming_distance(input1, input2) == 3
        input1 = "0000"
        input2 = "1111"
        @test hamming_distance(input1, input2) == 4
        input1 = "ABCDEFG"
        input2 = "AZCDWGM"
        @test hamming_distance(input1, input2) == 4
    end

    @testset "Strings: Rabin-Karp algorithm" begin
        text = "a"
        pattern = "airplane"
        @test rabin_karp(text, pattern) == []
        text = "Hello, World!!"
        pattern = "World"
        @test rabin_karp(text, pattern) == [8]
        text = "Banana"
        pattern = "an"
        @test rabin_karp(text, pattern) == [2, 4]
        text = "aaaaaa"
        pattern = "aaa"
        @test rabin_karp(text, pattern) == [1, 2, 3, 4]
    end
    @testset "Strings: Longest common subsequence" begin
        @test LCS("ABCD", "EFGHIJ") == ""

        @test LCS("AAAAA", "AAA") == "AAA"

        @test LCS("ABABABA", "ABBABBA") == "ABBABA"

        @test LCS("AGCDEG", "BGFDG") == "GDG"
    end
end

max_stars_repo_name: Nikola-Mircic/Julia

max_issues_repo_name: Nikola-Mircic/Julia

max_forks_repo_name: Nikola-Mircic/Julia


text:
rand_lat() = 180 * (rand()-0.5)
rand_lon() = 360 * (rand()-0.5)
rand_lon(m::Int64, n::Int64) = 360.0.*rand(Float64, m, n).-0.5
rand_datum() = rand() > 0.2 ? "WGS-84" : rand(geodetic_datum)

# Fill :misc with garbage
function rand_misc(N::Integer)
  D = Dict{String, Any}()
  for n = 1:N
    t = code2typ(rand(OK))
    k = randstring(rand(2:12))
    if Bool(t <: Real) == true
      D[k] = rand(t)
    elseif Bool(t <: Complex) == true
      D[k] = rand(Complex{real(t)})
    elseif Bool(t <: Array) == true
      y = eltype(t)
      if Bool(y <: Number) == true
        D[k] = rand(y, rand(1:1000))
      end
    end
  end
  (haskey(D, "hc")) && (delete!(D, "hc"))
  return D
end

function repop_id!(S::GphysData; s::Bool=true)
  while length(unique(S.id)) < length(S.id)
    for i = 1:S.n
      fs = S.fs[i]
      if fs == 0.0
        S.id[i] = rand_irr_id()
      else
        S.id[i], S.units[i] = rand_reg_id(fs, s)
      end
    end
  end
  return nothing
end


text:
# We will use the logistic resource type II consumer model - obviously the most interesting bifurcation diagrams

using Distributed
using PyPlot

#parallel set up
addprocs(length(Sys.cpu_info())-1)

@everywhere using Parameters
@everywhere using DifferentialEquations
@everywhere using ForwardDiff
@everywhere using LinearAlgebra

@everywhere include("/home/chrisgg/julia/TheoryClub/src/julia/CR/cr_model.jl")

#non parallel set up
include("cr_model.jl")

# Min max version
## Parallel version
@everywhere function calc_minmax(eval, tend)
    u0 = [0.5, 0.5]
    t_span = (0.0, tend)
    remove_transient = 300.0:1.0:tend
    p = CRPar(e = eval)
    prob = ODEProblem(cr_logII_model, u0, t_span, p)
    sol = DifferentialEquations.solve(prob, reltol = 1e-8)
    solend = sol(remove_transient)
    minC = minimum(solend[2, 1:end])
    maxC = maximum(solend[2, 1:end])

    return [minC, maxC]
end

function parallel_minmax(tend)
    effrange = 0.43:0.01:0.9
    data = pmap(eval -> calc_minmax(eval, tend), effrange)
    minmaxC = zeros(length(effrange),2)
    for i in 1:length(effrange)
        minmaxC[i,1] = data[i][1]
        minmaxC[i,2] = data[i][2]
    end
    return hcat(collect(effrange), minmaxC)
end

let
    data = parallel_minmax(400.0)
    minmaxplot = figure()
    plot(data[:,1], data[:,2])
    plot(data[:,1], data[:,3])
    xlabel("Efficiency")
    ylabel("Consumer Max/Min")
    return minmaxplot
end

## non-parallel version
function minmaxbifurc(tend)
    effrange = 0.43:0.01:0.9
    minC = zeros(length(effrange))
    maxC = zeros(length(effrange))

    u0 = [0.5, 0.5]
    t_span = (0.0, tend)
    remove_transient = 300.0:1.0:tend

    for (ei, evals) in enumerate(effrange)
        p = CRPar(e = evals)
        prob = ODEProblem(cr_logII_model, u0, t_span, p)
        sol = DifferentialEquations.solve(prob, reltol = 1e-8)
        solend = sol(remove_transient)
        solend[2, 1:end]

        minC[ei] = minimum(solend[2, 1:end])
        maxC[ei] = maximum(solend[2, 1:end])
    end

    return hcat(effrange, minC, maxC)
end

let
    data = minmaxbifurc(400.0)
    minmaxplot = figure()
    plot(data[:,1], data[:,2])
    plot(data[:,1], data[:,3])
    xlabel("Efficiency")
    ylabel("Consumer Max/Min")
    return minmaxplot
end

# Check mark diagram
function eff_maxeigen_data()
    evals = 0.4:0.0001:0.9
    max_eig = zeros(length(evals))

    for (ei, eval) in enumerate(evals)
        p = CRPar(e = eval)
        equ = logII_eq(p)
        max_eig[ei] = λ_stability(jac(equ, cr_logII_model, p))
    end
    return hcat(collect(evals), max_eig)
end

let
    data = eff_maxeigen_data()
    maxeigen_plot = figure()
    plot(data[:,1], data[:,2], color = "black")
    ylabel("Re(λₘₐₓ)", fontsize = 15)
    xlim(0.4, 0.8)
    ylim(-0.35, 0.1)
    xlabel("Efficiency", fontsize = 15)
    hlines(0.0, 0.4,0.8, linestyles = "dashed", linewidth = 0.5)
    vlines([0.441, 0.710], ymin = -0.35, ymax = 0.1, linestyles = "dashed")
    return maxeigen_plot
end

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
using NNlib: DenseConvDims
import NNlib: conv!, ∇conv_filter!, ∇conv_data!, conv_bias_act!

using CUDA.CUDNN: scalingParameter, CUDNN_CONVOLUTION, convdims,
                  cudnnConvolutionDescriptor, cudnnConvolutionBwdDataAlgoPerf,
                  cudnnConvolutionForward!, cudnnConvolutionBwdFilterAlgoPerf,
                  cudnnConvolutionBackwardData, cudnnConvolutionBackwardFilter,
                  cudnnConvolutionBackwardBias

const CUDNNFloat = Union{Float16,Float32,Float64}

function cudnnConvolutionDescriptor(cdims::DenseConvDims, x::DenseCuArray{T}) where T
    mode=(NNlib.flipkernel(cdims) ? CUDNN_CROSS_CORRELATION : CUDNN_CONVOLUTION)
    cudnnConvolutionDescriptor(convdims(nnlibPadding(cdims),size(x),0),
                               convdims(NNlib.stride(cdims),size(x),1),
                               convdims(NNlib.dilation(cdims),size(x),1),
                               mode,
                               cudnnDataType(T),
                               math_mode(),
                               CUDNN_DEFAULT_REORDER,
                               Cint(NNlib.groupcount(cdims)))
end

function conv!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}, cdims::DenseConvDims;
               alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    d = cudnnConvolutionDescriptor(cdims, x)
    cudnnConvolutionForward!(y, w, x, d; alpha, beta, z=y)
end

function conv_bias_act!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T},
                        cdims::DenseConvDims, bias::DenseCuArray{T}, σ=identity;
                        z::DenseCuArray{T}=y, alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    d = cudnnConvolutionDescriptor(cdims, x)
    # only relu and identity are supported by cudnnConvolutionForward!
    activation = (σ == NNlib.relu ? CUDNN_ACTIVATION_RELU : CUDNN_ACTIVATION_IDENTITY)
    cudnnConvolutionForward!(y, w, x, d; z, bias, activation, alpha, beta)
    if activation === CUDNN_ACTIVATION_IDENTITY && σ ∉ (nothing, identity)
        y = σ.(y)
    end
    return y
end

function ∇conv_data!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, w::DenseCuArray{T},
                     cdims::DenseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(dx), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(w)
    convDesc = cudnnConvolutionDescriptor(cdims, dx)
    p = cudnnConvolutionBwdDataAlgoPerf(wDesc, w, yDesc, dy, convDesc, xDesc, dx)
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardData(handle(), alpha, wDesc, w, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, xDesc, dx)
    end
    return dx
end

function ∇conv_filter!(dw::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                       cdims::DenseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(x), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(dw)
    convDesc = cudnnConvolutionDescriptor(cdims, x)
    p = cudnnConvolutionBwdFilterAlgoPerf(xDesc, x, yDesc, dy, convDesc, wDesc, dw);
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardFilter(handle(), alpha, xDesc, x, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, wDesc, dw);
    end
    return dw
end

function ∇conv_bias!(db::DenseCuArray{T}, dy::DenseCuArray{T}; alpha=1, beta=0) where T<:CUDNNFloat
    alpha,beta = scalingParameter(T,alpha), scalingParameter(T,beta)
    bDesc, yDesc = cudnnTensorDescriptor.((db,dy))
    cudnnConvolutionBackwardBias(handle(), alpha, yDesc, dy, beta, bDesc, db)
    return db
end


text:
# TODO
"""
Vector-based `SpatialIndex`.
While insertion is `O(1)`, the search is `O(N)`.

Generally should not be used, except for performance comparisons
or when the number of stored elements is expected to be very small (<100).
"""
struct SimpleSpatialIndex{T,N,V} <: SpatialIndex{T,N,V}
    elems::Vector{V}
end

Base.length(si::SimpleSpatialIndex) = length(si.elems)
Base.isempty(si::SimpleSpatialIndex) = isempty(si.elems)


text:
println("Processing Ecocyc...")
include("get_gene_information_ecocyc.jl")
println("Done!\n")
println("Processing RegulonDB...")
include("get_gene_information_regulonDB.jl")
println("Done!\n")
println("Processing search results...")
include("database_processing.jl")
println("Done!\n")
println("Designing sequences...")
include("design_sequences.jl")
println("Done!\n")


text:
# test_sequences.jl

function test_sequences()

end


text:
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Debug)
global_logger(logger)

ode = @ODEmodel(
    N'(t) = - N(t) * mu_N - N(t) * P(t) * delta_NE,
    E'(t) = N(t) * P(t) * delta_NE - E(t)^2 * mu_EE - E(t) * delta_EL + E(t) * P(t) * rho_E,
    S'(t) = S(t) * delta_EL - S(t) * delta_LM - S(t)^2 * mu_LL - E(t) * S(t) * mu_LE,
    M'(t) = S(t) * delta_LM - mu_M * M(t),
    P'(t) = P(t)^2 * rho_P - P(t) * mu_P - E(t) * P(t) * mu_PE - S(t) * P(t) * mu_PL,
    y1(t) = N(t),
    y2(t) = E(t) + S(t),
    y3(t) = M(t)
)

@time println(assess_identifiability(ode))


text:
## ---  Test equivalence to Base/Stdlib for other summary statistics

    # Test vsum
    for nd = 1:5
        @info "Testing vsum: $nd-dimensional arrays"
        # Generate random array
        A = rand((1 .+ (1:nd))...)

        # Test equivlalence when reducing over all dims
        Σ = sum(A)
        @test vsum(A, multithreaded=false) ≈ Σ
        @test vsum(A, multithreaded=true) ≈ Σ

        # Test equivalence when reducing over a single dimension
        for i = 1:nd
            @info "Testing vsum: reduction over dimension $i"
            Σ = sum(A, dims=i)
            @test vsum(A, dims=i, multithreaded=false) ≈ Σ
            @test vsum(A, dims=i, multithreaded=true) ≈ Σ
        end

        # Test equivalence when reducing over two dimensions
        if nd > 1
            for i = 2:nd
                for j = 1:i-1
                    @info "Testing vsum: reduction over dimensions $((j,i))"
                    Σ = sum(A, dims=(j,i))
                    @test vsum(A, dims=(j,i), multithreaded=false) ≈ Σ
                    @test vsum(A, dims=(j,i), multithreaded=true) ≈ Σ
                end
            end
        end

        # Test equivalence when reducing over three dimensions
        if nd > 2
            for i = 3:nd
                for j = 2:i-1
                    for k = 1:j-1
                        @info "Testing vsum: reduction over dimensions $((k,j,i))"
                        Σ = sum(A, dims=(k,j,i))
                        @test vsum(A, dims=(k,j,i), multithreaded=false) ≈ Σ
                        @test vsum(A, dims=(k,j,i), multithreaded=true) ≈ Σ
                    end
                end
            end
        end
    end

    # Test vmean
    for nd = 1:5
        @info "Testing vmean: $nd-dimensional arrays"
        # Generate random array
        A = rand((2 .+ (1:nd))...)

        # Test equivlalence when reducing over all dims
        μ = mean(A)
        @test vmean(A, multithreaded=false) ≈ μ
        @test vmean(A, multithreaded=true) ≈ μ

        # Test equivalence when reducing over a single dimension
        for i = 1:nd
            @info "Testing vmean: reduction over dimension $i"
            μ = mean(A, dims=i)
            @test vmean(A, dims=i, multithreaded=false) ≈ μ
            @test vmean(A, dims=i, multithreaded=true) ≈ μ
        end

        # Test equivalence when reducing over two dimensions
        if nd > 1
            for i = 2:nd
                for j = 1:i-1
                    @info "Testing vmean: reduction over dimensions $((j,i))"
                    μ = mean(A, dims=(j,i))
                    @test vmean(A, dims=(j,i), multithreaded=false) ≈ μ
                    @test vmean(A, dims=(j,i), multithreaded=true) ≈ μ
                end
            end
        end

        # Test equivalence when reducing over three dimensions
        if nd > 2
            for i = 3:nd
                for j = 2:i-1
                    for k = 1:j-1
                        @info "Testing vmean: reduction over dimensions $((k,j,i))"
                        μ = mean(A, dims=(k,j,i))
                        @test vmean(A, dims=(k,j,i), multithreaded=false) ≈ μ
                        @test vmean(A, dims=(k,j,i), multithreaded=true) ≈ μ
                    end
                end
            end
        end

        # Test equivalence when reducing over four dimensions
        if nd > 3
            for i = 4:nd
                for j = 2:i-1
                    for k = 1:j-1
                        for l = 1:k-1
                            @info "Testing vmean: reduction over dimensions $((l,k,j,i))"
                            μ = mean(A, dims=(l,k,j,i))
                            @test vmean(A, dims=(l,k,j,i), multithreaded=false) ≈ μ
                            @test vmean(A, dims=(l,k,j,i), multithreaded=true) ≈ μ
                        end
                    end
                end
            end
        end
    end

    # Test vvar
    for nd = 1:5
        @info "Testing vvar: $nd-dimensional arrays"
        # Generate random array
        A = randn((2 .+ (1:nd))...)

        # Test equivlalence when reducing over all dims
        σ² = var(A, corrected=false)
        @test vvar(A, corrected=false, multithreaded=false) ≈ σ²
        @test vvar(A, corrected=false, multithreaded=true) ≈ σ²

        # Test equivalence when reducing over a single dimension
        for i = 1:nd
            @info "Testing vvar: reduction over dimension $i"
            σ² = var(A, dims=i, corrected=false)
            @test vvar(A, dims=i, corrected=false, multithreaded=false) ≈ σ²
            @test vvar(A, dims=i, corrected=false, multithreaded=true) ≈ σ²
        end

        # Test equivalence when reducing over two dimensions
        if nd > 1
            for i = 2:nd
                for j = 1:i-1
                    @info "Testing vvar: reduction over dimensions $((j,i))"
                    σ² = var(A, dims=(j,i), corrected=false)
                    @test vvar(A, dims=(j,i), corrected=false, multithreaded=false) ≈ σ²
                    @test vvar(A, dims=(j,i), corrected=false, multithreaded=true) ≈ σ²
                end
            end
        end

        # Test equivalence when reducing over three dimensions
        if nd > 2
            for i = 3:nd
                for j = 2:i-1
                    for k = 1:j-1
                        @info "Testing vvar: reduction over dimensions $((k,j,i))"
                        σ² = var(A, dims=(k,j,i))
                        @test vvar(A, dims=(k,j,i), multithreaded=false) ≈ σ²
                        @test vvar(A, dims=(k,j,i), multithreaded=true) ≈ σ²
                    end
                end
            end
        end
    end

    # Test vstd
    for nd = 1:5
        @info "Testing vstd: $nd-dimensional arrays"
        # Generate random array
        A = randn((2 .+ (1:nd))...)

        # Test equivlalence when reducing over all dims
        @test vstd(A, multithreaded=true) ≈ vstd(A, multithreaded=false) ≈ std(A)
        @test vstd(A, corrected=false, multithreaded=true) ≈ vstd(A, corrected=false, multithreaded=false) ≈ std(A, corrected=false)

        # Test equivalence when reducing over a single dimension
        for i = 1:nd
            @info "Testing vstd: reduction over dimension $i"
            @test vstd(A, dims=i, multithreaded=true) ≈ vstd(A, dims=i, multithreaded=false) ≈ std(A, dims=i)
        end

        # Test equivalence when reducing over two dimensions
        if nd > 1
            for i = 2:nd
                for j = 1:i-1
                    @info "Testing vstd: reduction over dimensions $((j,i))"
                    @test vstd(A, dims=(j,i)) ≈ std(A, dims=(j,i))
                end
            end
        end

        # Test equivalence when reducing over three dimensions
        if nd > 2
            for i = 3:nd
                for j = 2:i-1
                    for k = 1:j-1
                        @info "Testing vstd: reduction over dimensions $((k,j,i))"
                        @test vstd(A, dims=(k,j,i)) ≈ std(A, dims=(k,j,i))
                    end
                end
            end
        end

        # Test equivalence when reducing over four dimensions
        if nd > 3
            for i = 4:nd
                for j = 2:i-1
                    for k = 1:j-1
                        for l = 1:k-1
                            @info "Testing vstd: reduction over dimensions $((l,k,j,i))"
                            @test vstd(A, dims=(l,k,j,i)) ≈ std(A, dims=(l,k,j,i))
                        end
                    end
                end
            end
        end
    end

    # Test fallbacks for complex reductions
    A = randn((2 .+ (1:6))...);
    @test vmean(A, dims=(4,5,6)) ≈ mean(A, dims=(4,5,6))
    @test vstd(A, dims=(4,5,6)) ≈ std(A, dims=(4,5,6))
    @test vstd(A, dims=(4,5,6)) ≈ vstd(A, dims=(4,5,6), mean=vmean(A, dims=(4,5,6)))

    # Test results when diminsions are out of range
    A = rand(10,10)
    @test vminimum(A, dims=3) == A
    @test vminimum(A, dims=(1,3)) == minimum(A, dims=(1,3))
    @test vsum(A, dims=3, multithreaded=false) == vmean(A, dims=(3,4), multithreaded=false) == A
    @test vsum(A, dims=3, multithreaded=true) == vmean(A, dims=(3,4), multithreaded=true) == A
    @test vmean(A, dims=3, multithreaded=false) == vmean(A, dims=(3,4), multithreaded=false) == A
    @test vmean(A, dims=3, multithreaded=true) == vmean(A, dims=(3,4), multithreaded=true) == A
    @test isequal(vstd(A, dims=(3,4), multithreaded=false), fill(NaN, 10, 10))
    @test isequal(vstd(A, dims=(3,4), multithreaded=true), fill(NaN, 10, 10))
    @test vstd(A, dims=(1,3)) == vstd(A, dims=1)

## -- End of File


text:
"""
    plotresults(;flux=flux, posterior=posterior, [σ2=1e-4*ones(size(flux, 1))])

"""
function plotresults(; flux=flux, posterior=posterior, g=g, randx = randx, σ2=1e-4*ones(size(flux, 1)), showstd=false)

  linemeansigma = marginaline(posterior, σ2)

  sampleinter   = [randx() for i=1:20]

  # project flux obs on line to find minimum, maximum t for plotting purposes

  D = round(Int, length(posterior)/2)

  ā, b̄ = mean(posterior)[1:D], mean(posterior)[D+1:2*D]

  proj(x) = dot(ā, x-b̄) / dot(ā, ā)

  tproj = [proj(flux[:,i]) for i=1:size(flux, 2)]

  tinter = proj(mean(sampleinter))

  tmin, tmax = min(tinter,minimum(tproj)), max(tinter, maximum(tproj))

  tmin, tmax = sign(tmin) * abs(tmin * 1.1), sign(tmax) * abs(tmax * 1.1) # make them a bit longer

  trange = collect(LinRange(tmin, tmax, 400))

  predictions     = map(linemeansigma, trange)

  meanpredictions = reduce(hcat, [p[1] for p in predictions])

  covpredictions  = reduce(hcat, [diag(p[2]) for p in predictions])

  figure()

  plot3D(flux[1,:], flux[2,:], flux[3,:], "mo", label="observed total flux", alpha=0.75)

  # plot3D(vec(meanpredictions[1,:]), vec(meanpredictions[2,:]), vec(meanpredictions[3,:]), "m-", label="mean fit")

  tmp = reduce(hcat, sampleinter)
  plot3D(tmp[1,:],tmp[2,:],tmp[3,:], "yo", label = "intersection samples", markersize=6)

  plotvector(g/norm(g) * norm(ā*tinter .+ b̄)*1.2, "g", "galaxy vector")

  if showstd

      circle = [cos.(LinRange(0.0,2π, 100)) sin.(LinRange(0.0,2π, 100))]'

      for t in 1:length(trange)

            μ = predictions[t][1]
            c = predictions[t][2]

            # principal components
            E = eigen(c)

            vaxis  = [E.vectors[:,1]*(E.values[1])*20 E.vectors[:,2]*(E.values[2])*20]
            v     = eigen(c).vectors[:, 3]

            pnts = vaxis*circle
            plot3D(pnts[1,:].+μ[1],pnts[2,:].+μ[2],pnts[3,:].+μ[3],"-c", alpha=0.1)

      end

      # plot3D(vec(meanpredictions[1,:]) .+ 2*vec(sqrt.(covpredictions[1,:])),
      #   vec(meanpredictions[2,:]),
      #   vec(meanpredictions[3,:]), "m-", label="2std", alpha=0.5)
      #
      # plot3D(vec(meanpredictions[1,:]),
      #   vec(meanpredictions[2,:] .+ 2*vec(sqrt.(covpredictions[2,:]))),
      #   vec(meanpredictions[3,:]), "m-", label="2std", alpha=0.5)
      #
      # plot3D(vec(meanpredictions[1,:]),
      #   vec(meanpredictions[2,:]),
      #   vec(meanpredictions[3,:]) .+ 2*vec(sqrt.(covpredictions[3,:])), "m-", label="2std", alpha=0.5)

  end

  # plot sample lines

  legend()

    xlim(0.0, (xlim()[2]))
    ylim(0.0, (ylim()[2]))
    zlim(0.0, (zlim()[2]))

    xlabel("band at 3670")
    ylabel("band at 4826")
    zlabel("band at 6223")

  nothing

end


text:
using Distributed
@everywhere using DelimitedFiles, JuMP, Gurobi, Statistics, StatsBase, DataFrames, Gadfly
include("./init.jl")
include("../input_data/read_data.jl")
@everywhere include("../input_data/data_manipulation.jl")
@everywhere include("./basic_functions.jl")
@everywhere include("./validations.jl")
@everywhere include("../heuristics/indicators.jl")
@everywhere include("../heuristics/construction_algorithms.jl")
@everywhere include("../heuristics/local_search.jl")
include("../JuMP/JuMP_functions.jl")
include("../JuMP/JuMP_models.jl")
@everywhere include("../solution_methods/GRASP_filter_multi.jl")
# @everywhere include("../solution_methods/GRASP_local_search_multi.jl")
# @everywhere include("../solution_methods/GRASP_reactive_multi.jl")
# @everywhere include("../solution_methods/GRASP_algorithm_multi.jl")
# include("../solution_methods/GRASP_algorithm.jl")
include("../output_data/results.jl")

#Read Data
    #ATOM
    #doc = "./data/Benchmark/60C_10Type_LessDense_2QC.txt"
    #doc = "./data/Benchmark/240C_20Type_LessDense_2QC.txt"
    #doc = "./data/Benchmark/500C_100Type_LessDense_4QC.txt"
    #doc = "./data/Benchmark/1000C_100Type_UniformDense_4QC.txt"

results = Dict{String, Array{Tuple{Number, Number}, 1}}()
# for i = 1:10
for file in reverse(readdir("../../data/Benchmark/"))
    println("-------------")
    println(file)
    if file != ".ipynb_checkpoints" && file != "00_DataDescription.txt"
    # if file == "500C_100Type_Scattered_4QC.txt"
        if file[1:3] == "60C"
            if file[5:6] == "10"
                beta = -0.5
            else
                beta = 0.25
            end
        elseif file[1:3] == "240"
            if file[6:7] == "20"
                beta = -0.25
            else
                beta = 0.25
            end
        elseif file[1:3] == "500"
            if file[6:7] == "20"
                beta = 0.75
            else
                beta = -0.25
            end
        else
            if file[6:7] == "20"
                beta = 0.75
            elseif file[6:7] == "60"
                beta = 0.5
            else
                beta = 0.25
            end
        end

        doc = "../../data/Benchmark/"*file
        C,P,CP,Q,J,tt,d,ci,pj,bj,prejj,cpij = read_data(doc)

        possible_tasks = Array{LTask, 1}()
        tasks_by_position = Dict{Int, Array{LTask, 1}}()
        prec = Dict{Int, Array}()
        task_times = zeros(Int, P, C)

        for p = 1:P
            l=Array{LTask, 1}()
            for t=1:CP
                target_bay = bj[cpij[t,1]+1]
                if p==P
                    #create array with all tasks
                    push!(possible_tasks, LTask(cpij[t,1]+1, target_bay, cpij[t,2]+1, cpij[t,3]))
                    #create matrix with all task times
                    task_times[cpij[t,1]+1, cpij[t,2]+1] = cpij[t,3]
                end
                #create dict with the tasks sorted by position
                if cpij[t,1]+1==p
                    push!(l,LTask(cpij[t,1]+1, target_bay, cpij[t,2]+1, cpij[t,3]))
                end
            end
            tasks_by_position[p]=sort!(l, by=x->x.t)

            #create dict with the position precedences
            s=Array{Int,1}()
            for j = 1:P
                if prejj[p,j]==1
                    push!(s,j)
                end
            end
            prec[p]=s

        end

        H = horizon_plan(C,P,J,Q,tt,d,bj,tasks_by_position)
        global CTS = Constants(C,P,J,Q,H,tt,d)
        global it_results, best_makespan, best_LS, exec_time
        start = time()
        makespan_1, best_LS_1 = GRASP_multi_thread(0 , prec, task_times, bj, CTS)
        makespan_2, best_LS_2 = GRASP_multi_thread(beta , prec, task_times, bj, CTS)
        if makespan_1 < makespan_2
            best_makespan = makespan_1
            best_LS = best_LS_1
        else
            best_makespan = makespan_2
            best_LS = best_LS_2
        end
        exec_time = time() - start
        #write_results(doc, it_results)
        # write_time_results(doc, exec_time)

        println(best_makespan)
        println(exec_time)
        #println(best_LS)
        if haskey(results, file)
            push!(results[file], (best_makespan, exec_time))
        else
            results[file] = [(best_makespan, exec_time)]
        end
    end
end
# end

# plot_solution(best_LS, best_makespan, CTS)

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
@testset "perft1" begin
    b = Bobby.set_board()
    pt = Bobby.perft(b, 4, b.player_color)
    @test pt.nodes == [20, 400, 8902, 197281]
end

@testset "perft2" begin
    @test test_fen(
        "r3k2r/p1ppqpb1/bn2pnp1/3PN3/1p2P3/2N2Q1p/PPPBBPPP/R3K2R w KQkq - 1 1",
         2, [48, 2039])
    # [48, 2039, 97862, 4085603]
end

@testset "perft3" begin
    @test test_fen("8/2p5/3p4/KP5r/1R3p1k/8/4P1P1/8 w - - 1 1", 2, [14, 191])
    # [14, 191, 2812, 43238]
end

@testset "perft4" begin
    @test test_fen(
        "r3k2r/Pppp1ppp/1b3nbN/nP6/BBP1P3/q4N2/Pp1P2PP/R2Q1RK1 w kq - 0 1",
        2, [6, 264])
    # [6, 264, 9467, 422333]
end

@testset "perft4.1" begin
    @test test_fen(
        "r2q1rk1/pP1p2pp/Q4n2/bbp1p3/Np6/1B3NBn/pPPP1PPP/R3K2R b KQ - 0 1",
        2, [6, 264])
    # [6, 264, 9467, 422333]
end

@testset "perft5" begin
    @test test_fen(
        "rnbq1k1r/pp1Pbppp/2p5/8/2B5/8/PPP1NnPP/RNBQK2R w KQ - 1 8",
        3, [44, 1486, 62379])
    # [44, 1486, 62379, 2103487]
end

@testset "perft6" begin
    @test test_fen(
        "r4rk1/1pp1qppp/p1np1n2/2b1p1B1/2B1P1b1/P1NP1N2/1PP1QPPP/R4RK1 w - - 0 10",
        2, [46, 2079])
    # [46, 2079, 89890, 3894594]
end

@testset "perft7" begin
    @test test_fen(
        "r3k2r/1bp2pP1/5n2/1P1Q4/1pPq4/5N2/1B1P2p1/R3K2R b KQkq c3 0 1",
        2, [60, 2608])
    # [60, 2608, 113742, 4812099]
end


text:
import ArgParse
import ArgParse: @project_version, @add_arg_table!, parse_args, ArgParseSettings

function parse_inputs(args; error_handler = ArgParse.default_handler)
    s = ArgParseSettings(
        description = "Reproducible Input Ouput (I/O) Pattern Application (RIOPA)",
        add_version = true,
        version = @project_version,
        commands_are_required = false,
        exc_handler = error_handler,
    )
    @add_arg_table! s begin
        "hello"
        help = """Run in "hello" mode (minimal functionality)"""
        action = :command
        "generate-config", "gencfg"
        help = """Create default configuration file
        (as "$(default_config_filename())" unless --config option is used)"""
        action = :command
        "--config", "-c"
        help = "Specify name of (YAML) config file to generate I/O"
        arg_type = String
    end

    inputs = parse_args(args, s)

    return inputs
end


text:
abstract type Extinction{R,W} <: CellRule{R,W} end

"""
    RandomExtinction
"""
struct RandomExtinction{R,W,C} <: Extinction{R,W}
    probability::C
end
function RandomExtinction{R,W}(; probability::PT = 0.1) where {R,W,PT}
    RandomExtinction{R,W,PT}(probability)
end

function RandomExtinction(names::T; probability = 0.1) where {T<:Vector{Symbol}}
    rules = Ruleset()
    for n in names
        rules += RandomExtinction{n}(probability = probability)
    end
    return rules
end

function DynamicGrids.applyrule(data, rule::RandomExtinction, N, I)
    N > zero(N) || return zero(N)
    rand() < rule.probability || return (N)
    return zero(N)
end

"""
    AbioticExtinction
"""
struct AbioticExtinction{R,W,L,C,F} <: Extinction{R,W}
    envlayer::L
    baseprobability::C
    mapping::F
end
function AbioticExtinction{R,W}(
    envlayer;
    baseprobability = 0.1,
    mapping = ExponentialFitness(),
) where {R,W}
    AbioticExtinction{R,W}(envlayer, baseprobability, mapping)
end
function DynamicGrids.applyrule(data, rule::AbioticExtinction, O, I)
    E = data[rule.envlayer][I...]
    O > zero(O) || return zero(O)
    extinctprob = rule.baseprobability * rule.mapping(E)  # goes to 0 as E goes to infinity with default Expontential Fitness
    rand() < extinctprob && return one(O)
    return zero(O)
end


text:
# This file is a part of SimilaritySearch.jl

export PivotedSearch

"""
    PivotedSeach(index::PivotTable, db::AbstractVector, dist::PreMetric)

Defines a search index for Pivot tables
"""
struct PivotedSearch{DataType<:AbstractVector, DistanceType<:PreMetric} <: AbstractSearchContext
    dist::DistanceType
    db::DataType
    pivots::DataType
    table::Vector{Vector{Float32}} # pivot table
    dqp::Vector{Float32} # query mapped to the pivot space
end

Base.copy(index::PivotedSearch;
        dist=index.dist,
        db=index.db,
        pivots=index.pivots,
        table=index.table,
        dqp=zeros(Float32, length(pivots))
    ) = PivotedSearch(dist, db, pivots, table, dqp)

Base.string(p::PivotedSearch) = "{PivotedSearch: dist=$(p.dist), n=$(length(p.db)), pivs=$(length(p.pivots))}"

PivotedSearch(dist::PreMetric, db, pivots, table) =
    PivotedSearch(dist, db, pivots, table, zeros(Float32, length(pivots)))

"""
    PivotedSearch(::Type{PivotTable}, dist::PreMetric, db::AbstractVector{T}, pivots::Vector{T})
    PivotedSearch(::Type{PivotTable}, dist::PreMetric, db::AbstractVector{T}, numpivots::Integer)

Creates a `PivotTable` index with the given pivots. If the number of pivots is specified,
then they will be randomly selected from the dataset.
"""
function PivotedSearch(dist::PreMetric, db::AbstractVector{T}, pivots::AbstractVector{T}) where T
    @info "Creating a pivot table with $(length(pivots)) pivots and distance=$(dist)"
    table = Vector{Vector{Float32}}(undef, length(db))
    for i in 1:length(db)
        u = db[i]
        table[i] = [evaluate(dist, u, pivots[j]) for j in 1:length(pivots)]
    end

    PivotedSearch(dist, db, pivots, table)
end

function PivotedSearch(dist::PreMetric, db::AbstractVector{T}, numpivots::Integer) where T
    pivots = rand(db, numpivots)
    PivotedSearch(dist, db, pivots)
end

"""
    search(index::PivotedSeach, q, res::KnnResult)

Solves a query with the pivot index.
"""
function search(index::PivotedSearch, q, res::KnnResult)
    @inbounds for i in eachindex(index.pivots)
        index.dqp[i] = evaluate(index.dist, q, index.pivots[i])
    end

    for i in eachindex(index.db)
        dpu = index.table[i]
        need_eval = true

        for pivID in 1:length(index.pivots)
            if abs(index.dqp[pivID] - dpu[pivID]) > covrad(res)
                need_eval = false
                break
            end
        end

        if need_eval
            d = evaluate(index.dist, q, index[i])
            push!(res, i, d)
        end
    end

    res
end


text:
# see also https://github.com/tlienart/Franklin.jl/issues/330
@testset "locvar" begin
    s = raw"""
        @def va = 5
        @def vb = 7
        ```julia:ex
        #hideall
        println(locvar("va")+locvar("vb"))
        ```
        \output{ex}
        """ |> fd2html_td
    @test isapproxstr(s, """

12

        """)
end


text:
function hlp_JuMP_inflow_discharge(model::JuMP.Model, r::Reservoir, s::Int, t::Int;
                                   imposeDiverted::Union{Bool, JuMP.JuMPArray}=true)
  if imposeDiverted === false # Comparison on value and type.
    error("When not imposing the value of the diverted rivers, must provide a set of variables.")
  end

  return sum([hlp_JuMP_inflow_discharge(model, r, riv, s, t, imposeDiverted=imposeDiverted) for riv in eachriver(r)])
end

hlp_JuMP_inflow_discharge(model::JuMP.Model, r::Reservoir, nr::NaturalRiver, s::Int, t::Int; kwargs...) =
  getScenario(nr, s, t)

function hlp_JuMP_inflow_discharge(model::JuMP.Model, r::Reservoir, dr::DivertedRiver, s::Int, t::Int;
                                   imposeDiverted::Union{Bool, JuMP.JuMPArray, JuMP.Variable}=true)
  if imposeDiverted === false # Comparison on value and type.
    error("When not imposing the value of the diverted rivers, must provide a set of variables in the argument imposeDiverted.")
  end

  if imposeDiverted === true # Comparison on value and type.
    return maxAllowableFlow(dr, s, t)
  else
    if isa(imposeDiverted, JuMP.JuMPArray)
      imposeDiverted = imposeDiverted[s, t, getName(dr)]
    end

    @constraint(model, imposeDiverted <= maxAllowableFlow(dr, s, t))
    return imposeDiverted
  end
end


text:
using Primes

ν(n::Int) = sum(digits(n, base=2))
λ(n::Int) = floor(Int, log2(n))
function lowerbound(n::Int)
    if n ≤ 327_678 || ν(n) ≤ 16
        # We add 1 because Julia indexs from 1
        λ(n) + ceil(Int, log2(ν(n))) + 1
    else
        # We add 1 because Julia indexs from 1
        ceil(Int, log2(n) + log2(ν(n)) - 2.13)
    end
end

abstract type BoundingSequence end
abstract type BoundingSequenceA <: BoundingSequence end
abstract type BoundingSequenceC <: BoundingSequence end

function bounds(::Type{BoundingSequenceA}, n::Int, lb::Int)
    seq = zeros(Int, lb + 1)
    seq[end] = n
    for i in 1:lb
        seq[end-i] = ceil(Int, seq[end-i+1] / 2)
    end
    seq
end

function bounds(::Type{BoundingSequenceC}, n::Int, lb::Int)
    t = iseven(n) ? factor(n)[2] : zero(n)

    seq = zeros(Int, lb + 1)
    seq[end] = n
    for i in 1:lb-t-1
        seq[i] = ceil(Int, n / (3 * 2^(lb - (i + 1))))
    end
    for i in lb-t:lb
        seq[i] = ceil(Int, n / 2^(lb - (i - 1)))
    end
    seq
end

function bounds(n::Int, lb::Int)
    if n % 5 != 0
        vertical = bounds(BoundingSequenceC, n, lb)
        vertical, vertical
    else
        vertical = bounds(BoundingSequenceC, n, lb)
        slant = bounds(BoundingSequenceA, n, lb)
        vertical, slant
    end
end

retain(bound::Int, aᵢ::Int) = aᵢ ≥ bound
retain(bound::Int, aᵢ₋₁::Int, aᵢ::Int) = aᵢ + aᵢ₋₁ ≥ bound
retain(n::Int, lb::Int, i::Int, aᵢ::Int) = n != 2^(lb - i + 1) * aᵢ

function retain(n::Int, lb::Int, v::Int, s::Int, i::Int, aᵢ₋₁::Int, aᵢ::Int)
    retain(v, aᵢ) && retain(n, lb, i, aᵢ) && retain(s, aᵢ₋₁, aᵢ)
end

function stackchildren!(n::Int, stack::Vector{Vector{Int}}; verbose=false)
    aᵢ = stack[end][end]
    segment = Int[]
    for i in 1:length(stack), j in i:length(stack)
        aᵢ₊₁ = stack[i][end] + stack[j][end]
        if aᵢ < aᵢ₊₁ ≤ n
            push!(segment, aᵢ₊₁)
        end
    end
    unique!(sort!(segment))
    push!(stack, segment)
end

function backup(stack::Vector{Vector{Int}}; verbose=false)
    while length(stack) > 2
        pop!(stack[end])
        if isempty(stack[end])
            pop!(stack)
        else
            break
        end
    end
    length(stack) == 2
end

function thurber(n::Int; verbose=false)
    if n < one(n)
        error("no chains defined for integers less than 1")
    elseif isone(n)
        return 0
    end

    stack = Vector{Int}[[1],[2]]
    lb = lowerbound(n)

    loop = 1
    while true
        vertical, slant = bounds(n, lb)
        verbose && @info "Outer Loop" lb vertical slant
        while true
            i = length(stack)
            verbose && @info "State" n stack lb loop
            if i ≤ lb
                aᵢ₋₁, aᵢ = stack[i-1][end], stack[i][end]

                if retain(n, lb, vertical[i], slant[i+1], i, aᵢ₋₁, aᵢ)
                    verbose && @info "Retained" aᵢ
                    if aᵢ == n
                        return length(stack) - 1
                    end
                    stackchildren!(n, stack; verbose)
                else
                    verbose && @info "Did not retain" aᵢ loop
                    if backup(stack; verbose)
                        loop += 1
                        break
                    end
                end
            else
                verbose && @info "Reached lower bound" loop
                if backup(stack; verbose)
                    loop += 1
                    break
                end
            end
            loop += 1
        end
        lb += 1
    end
end


text:
#Import Julia Libs
using LinearAlgebra
using LinearAlgebra.BLAS
using DifferentialEquations
using SparseArrays

using Parameters
using HDF5
using DelimitedFiles

using ProgressMeter

#Import structures
include("./structs/Operators.jl")
include("./structs/Parameters.jl")

#Operators
include("./operators/Diff_Opps.jl")
include("./operators/Pos_Opps.jl")
include("./operators/absorption.jl")
include("./operators/boundary_conditions.jl")
include("./operators/potential.jl")

#Simulation helpers
include("./sim_utils/initial_conditions.jl")
include("./sim_utils/create_structures.jl")
include("./sim_utils/setup_folders.jl")
include("./sim_utils/setup_routine.jl")

#Solvers
include("./solvers/KG_Solver.jl")

#Include other utils
include("./utils/print_param_info.jl")
include("./utils/data_read.jl")
include("./utils/data_read_filename.jl")
include("./utils/get_energy.jl")


text:
@show LongMultiplication.longmult(big(2) ^ 64, big(2) ^ 64)
@show LongMultiplication.longmult("18446744073709551616", "18446744073709551616")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# References:
#
# - https://lammps.sandia.gov/doc/Howto_library.htm
# - https://lammps.sandia.gov/doc/Python_library.html
#
# Some wrapper methods were generated with help of Clang.jl

"""
Julia wrapper around `liblammps`.
"""
module LAMMPSWrapper

# Dependencies
using Base: cconvert

using Reexport
@reexport using FilePathsBase
@reexport using MPI

# JLL dependencies
using LAMMPS_jll

const LAMMPSPtr = Ptr{Cvoid}

"List of available LAMMPS packages"
const PACKAGES = String[]

# TODO: Design a way to check overlap of MPI communicators for different
# LAMMPS instances.
"""
    LAMMPSObject

Wrapper around LAMMPS objects.

# Fields
- `args::Vector{String}`: LAMMPS command-line arguments.
- `comm::MPI.Comm`: MPI communicator.
- `ptr::LAMMPSPtr`: Pointer to the running instance.
"""
mutable struct LAMMPSObject
    args::Vector{String}
    comm::MPI.Comm
    ptr::LAMMPSPtr

    @doc """
        LAMMPSObject(; args = String[], comm = MPI.COMM_WORLD)

    Creates a LAMMPS instance.

    # Keywords
    - `args::Vector{String} = String[]`: LAMMPS command-line arguments for this instance.
    - `comm::MPI.Comm = MPI.COMM_WORLD`: MPI communicator.
    """
    function LAMMPSObject(; args = String[], comm = MPI.COMM_WORLD)
        ptrref = Ref{LAMMPSPtr}(C_NULL)

        ccall(
            (:lammps_open, liblammps),
            Cvoid,
            (Cint, Ref{Cstring}, MPI.Comm, Ptr{LAMMPSPtr}),
            length(args), args, comm, ptrref
        )

        o = new(args, comm, ptrref[])
        finalizer(unset!, o)
        return o
    end
end

"""
    initialized(o::LAMMPSObject) :: Bool

Returns true if `o` has been instantiated from `liblammps`.
"""
initialized(o::LAMMPSObject) = o.ptr != C_NULL

"""
    check_state(o::LAMMPSObject)

Throws is instance `o` is currently inactive.
"""
function check_state(o::LAMMPSObject)
    if !initialized(o)
        error("LAMMPS instance is inactive")
    end
    return nothing
end

"""
    reset!(o::LAMMPSObject; args = o.args, comm = o.comm)

Restarts an inactive LAMMPS instance `o`.

# Arguments
- `o::LAMMPSObject`: Previously defined, but currently inactive, LAMMPS instance.

# Keywords
- `args::Vector{String} = o.args`: LAMMPS command-line arguments for this instance.
- `comm::MPI.Comm = o.comm`: MPI communicator.
"""
function reset!(o::LAMMPSObject; args = o.args, comm = o.comm)
    if initialized(o)
        error("cannot reset active LAMMPS objects")
    end

    ptrref = Ref{LAMMPSPtr}(C_NULL)

    ccall(
        (:lammps_open, liblammps),
        Cvoid,
        (Cint, Ref{Cstring}, MPI.Comm, Ptr{LAMMPSPtr}),
        length(args), args, comm, ptrref
    )

    o.args = args
    o.comm = comm
    o.ptr = ptrref[]
    return o
end

"""
    unset!(o::LAMMPSObject)

Closes LAMMPS instance `o` from `liblammps` and incativates it.
"""
function unset!(o::LAMMPSObject)
    if initialized(o)
        ccall((:lammps_close, liblammps), Cvoid, (LAMMPSPtr,), o.ptr)
        # Clean the pointer address so we can check safely before running any
        # command on a given instance
        o.ptr = C_NULL
    end
    return nothing
end

"""
    version(o::LAMMPSObject) :: Int

Returns the version of the underlying LAMMPS library as an integer with format
`YYYYMMDD`.
"""
function version(o::LAMMPSObject)
    check_state(o)
    return Int(ccall((:lammps_version, liblammps), Cint, (LAMMPSPtr,), o.ptr))
end

"""
    load!(o::LAMMPSObject, filename::AbstractPath)

Runs the LAMMPS input file `filename` on instance `o`.
"""
function load!(o::LAMMPSObject, path::AbstractPath)
    check_state(o)
    if !isfile(path)
        throw(ArgumentError("$filename does not exist or is not a file"))
    end
    ccall(
        (:lammps_file, liblammps),
        Cvoid,
        (LAMMPSPtr, Cstring),
        o.ptr, last(path.segments)
    )
    return nothing
end
"""
    load!(o::LAMMPSObject, cmd::String)

Runs LAMMPS command `cmd` on instance `o`.
"""
function load!(o::LAMMPSObject, cmd::String)
    check_state(o)
    if ispath(cmd)
        throw(ArgumentError("for running input files call `run(o, p\"$cmd\")` instead"))
    elseif '\n' in cmd
        # There is no point in supporting `:lammps_commands_string` from
        # `liblammps` if we already have support for `:lammps_commands_list`
        throw(ArgumentError(
            "for running multiple commands call `run(o, split($cmd, '\n'))` instead"
        ))
    end
    try
        ccall((:lammps_command, liblammps), Cstring, (LAMMPSPtr, Cstring), o.ptr, cmd)
    catch e
        throw(e)
    end
    return nothing
end
"""
    load!(o::LAMMPSObject, cmds::Vector{String})

Runs a list of LAMMPS commands `cmds` on instance `o`.
"""
function run(o::LAMMPSObject, cmds::Vector{String})
    check_state(o)
    ccall(
        (:lammps_commands_list, liblammps),
        Cvoid,
        (LAMMPSPtr, Cint, Ref{Cstring}),
        o.ptr, length(cmds), cmds
    )
    return nothing
end

"""
    sizeof(o::LAMMPSObject, typename::String)

Size of the LAMMPS data type `typename`.
"""
function sizeof(o::LAMMPSObject, typename::String)
    check_state(o)
    size = ccall(
        (:lammps_extract_setting, liblammps),
        Cint,
        (LAMMPSPtr, Cstring),
        o.ptr, typename
    )
    return Int(size)
end

"""
    get_box(o::LAMMPSObject) :: NamedTuple

Box information of LAMMPS instance `o`.

# Returns
NamedTuple with entries
- `o::Vector{Float64}`: min-coordinates corner.
- `hi::Vector{Float64}`: max-coordinates corner.
- `xy::Float64`: tilt factor for XY face.
- `yz::Float64`: tilt factor for YZ face.
- `xz::Float64`: tilt factor for XZ face.
- `periodicity::Vector{Int}`
- `box_change::Int`
"""
function get_box(o::LAMMPSObject)
    check_state(o)

    # See https://sourceforge.net/p/lammps/mailman/message/36726313 for
    # reference on the output structure to `lammps_extract_box`
    lo = Vector{Cdouble}(undef, 3)
    hi = Vector{Cdouble}(undef, 3)
    xy = Ref{Cdouble}(0)
    yz = Ref{Cdouble}(0)
    xz = Ref{Cdouble}(0)
    periodicity = Vector{Cint}(undef, 3)
    box_change = Ref{Cint}(0)

    ccall(
        (:lammps_extract_box, liblammps),
        Cvoid,
        (
            LAMMPSPtr, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
            Ref{Cdouble}, Ref{Cdouble}, Ref{Cint}, Ref{Cint}
        ),
        o.ptr, lo, hi, xy, yz, xz, periodicity, box_change
    )

    return (
        lo = lo,
        hi = hi,
        xy = xy[],
        yz = yz[],
        xz = xz[],
        periodicity = Vector{Int}(periodicity),
        box_change = Int(box_change[]),
    )
end

# TODO: Find out allowed `property` values and their respective data types
function get_global_ptr(o::LAMMPSObject, property::String)
    check_state(o)
    return ccall(
        (:lammps_extract_global, liblammps),
        Ptr{Cvoid},
        (LAMMPSPtr, Cstring),
        o.ptr, property
    )
end

# TODO: Find out allowed `property` values and their respective data types
function get_atom_ptr(o::LAMMPSObject, property::String)
    check_state(o)
    return ccall(
        (:lammps_extract_atom, liblammps),
        Ptr{Cvoid},
        (LAMMPSPtr, Cstring),
        o.ptr, property
    )
end

# TODO: Find out allowed `compute` values and their respective data types
function get_compute_ptr(o::LAMMPSObject, name, i, j)
    check_state(o)
    return ccall(
        (:lammps_extract_compute, liblammps),
        Ptr{Cvoid},
        (LAMMPSPtr, Cstring, Cint, Cint),
        o.ptr, name, i, j
    )
end

# TODO: Find out allowed `fix` values and their respective data types
function get_fix_ptr(o::LAMMPSObject, name, id, style, i, j)
    check_state(o)
    return ccall(
        (:lammps_extract_fix, liblammps),
        Ptr{Cvoid},
        (LAMMPSPtr, Cstring, Cint, Cint, Cint, Cint),
        o.ptr, name, id, style, i, j
    )
end

# TODO: Find out allowed `variable` values and their respective data types
function get_variable_ptr(o::LAMMPSObject, name, group)
    check_state(o)
    return ccall(
        (:lammps_extract_variable, liblammps),
        Ptr{Cvoid},
        (LAMMPSPtr, Cstring, Cstring),
        o.ptr, name, group
    )
end

const thermo_keywords = [
    "step", "elapsed", "elaplong", "dt", "time", "cpu", "tpcpu", "spcpu",
    "cpuremain", "part", "timeremain", "atoms", "temp", "press", "pe", "ke",
    "etotal", "enthalpy", "evdwl", "ecoul", "epair", "ebond", "eangle",
    "edihed", "eimp", "emol", "elong", "etail", "vol", "density", "lx", "ly",
    "lz", "xlo", "xhi", "ylo", "yhi", "zlo", "zhi", "xy", "xz", "yz", "xlat",
    "ylat", "zlat", "bonds", "angles", "dihedrals", "impropers", "pxx", "pyy",
    "pzz", "pxy", "pxz", "pyz", "fmax", "fnorm", "nbuild", "ndanger", "cella",
    "cellb", "cellc", "cellalpha",
]

const thermo_regexps = [
    r"^(c|f)_\w+(\[(\d*\*\d*|\d+)\])?$",
    r"^(c|f)_\w+\[\d+\]\[\d+\]",
    r"^v_\w+(\[\d+\])$",
]

"""
    get_thermo(o::LAMMPSObject, keyword::String) :: Float64

Returns the value of a thermo keyword. A non-exhaustive list of allowed
`keyword` values is stored in the `LAMMPS.thermo_keywords` variable. For a
detailed description of each valid `keyword`, refer to the [`thermo_style`]
section of the LAMMPS manual.

[`thermo_style`]: https://lammps.sandia.gov/doc/thermo_style.html
"""
function get_thermo(o::LAMMPSObject, keyword::String)
    check_state(o)
    is_valid_keyword = (
        keyword in thermo_keywords ||
        any(occursin(r, keyword) for r in thermo_regexps)
    )
    if !is_valid_keyword
        throw("invalid keyword")
    end
    value = ccall(
        (:lammps_get_thermo, liblammps),
        Cdouble,
        (LAMMPSPtr, Cstring),
        o.ptr,
        keyword
    )
    return Float64(value)
end

"""
    get_natoms(o::LAMMPSObject) :: Int

Returns the total number of atoms in the system.
"""
function get_natoms(o::LAMMPSObject)
    check_state(o)
    return Int(ccall((:lammps_get_natoms, liblammps), Cint, (LAMMPSPtr,), o.ptr))
end

"""
    set_variable!(o::LAMMPSObject, old::String, new::String)

Sets an existing string-style variable to a new string value.
"""
function set_variable!(o::LAMMPSObject, variable::String, name::String)
    check_state(o)
    flag = ccall(
        (:lammps_set_variable, liblammps),
        Cint,
        (LAMMPSPtr, Cstring, Cstring),
        o.ptr, variabl, name
    )
    if flag != 0
        throw("could not set variable $variable to $name")
    end
    return nothing
end

"""
    reset_box!(o::LAMMPSObject, lo::Vector{<:Real}, hi::Vector{<:Real},
               xy::Real, yz::Real, xz::Real)

Resets the size and shape of the simulation box.

# Arguments
- `o::LAMMPSObject`: LAMMPS instance.
- `lo::Vector{<:Real}`: min-coordinates corner.
- `hi::Vector{<:Real}`: max-coordinates corner.
- `xy::Real`: XY-plane tilt value.
- `yz::Real`: YZ-plane tilt value.
- `xz::Real`: XZ-plane tilt value.
"""
function reset_box!(
    o::LAMMPSObject, lo::Vector{<:Real}, hi::Vector{<:Real}, xy::Real, yz::Real, xz::Real
)
    check_state(o)
    @assert length(lo) == length(hi) == 3
    new_lo = cconvert(Vector{Cdouble}, lo)
    new_hi = cconvert(Vector{Cdouble}, hi)
    ccall(
        (:lammps_reset_box, liblammps),
        Cvoid,
        (LAMMPSPtr, Ref{Cdouble}, Ref{Cdouble}, Cdouble, Cdouble, Cdouble),
        o.ptr, new_lo, new_hi, xy, yz, xz
    )
    return nothing
end

#function gather_atoms(o::LAMMPSObject, arg2, arg3, arg4, arg5)
#    check_state(o)
#    ccall(
#        (:lammps_gather_atoms, liblammps),
#        Cvoid,
#        (LAMMPSPtr, Cstring, Cint, Cint, Ptr{Cvoid}),
#        o.ptr, arg2, arg3, arg4, arg5
#    )
#end

#function gather_atoms_concat(o::LAMMPSObject, arg2, arg3, arg4, arg5)
#    check_state(o)
#    ccall(
#        (:lammps_gather_atoms_concat, liblammps),
#        Cvoid,
#        (LAMMPSPtr, Cstring, Cint, Cint, Ptr{Cvoid}),
#        o.ptr, arg2, arg3, arg4, arg5
#    )
#end

#function gather_atoms_subset(o::LAMMPSObject, arg2, arg3, arg4, arg5, arg6, arg7)
#    check_state(o)
#    ccall(
#        (:lammps_gather_atoms_subset, liblammps),
#        Cvoid,
#        (LAMMPSPtr, Cstring, Cint, Cint, Cint, Ptr{Cint}, Ptr{Cvoid}),
#        o.ptr, arg2, arg3, arg4, arg5, arg6, arg7
#    )
#end

#function scatter_atoms(o::LAMMPSObject, arg2, arg3, arg4, arg5)
#    check_state(o)
#    ccall(
#        (:lammps_scatter_atoms, liblammps),
#        Cvoid,
#        (LAMMPSPtr, Cstring, Cint, Cint, Ptr{Cvoid}),
#        o.ptr, arg2, arg3, arg4, arg5
#    )
#end

#function scatter_atoms_subset(o::LAMMPSObject, arg2, arg3, arg4, arg5, arg6, arg7)
#    check_state(o)
#    ccall(
#        (:lammps_scatter_atoms_subset, liblammps),
#        Cvoid,
#        (LAMMPSPtr, Cstring, Cint, Cint, Cint, Ptr{Cint}, Ptr{Cvoid}),
#        o.ptr, arg2, arg3, arg4, arg5, arg6, arg7
#    )
#end

function set_callback!(o::LAMMPSObject, arg2, callback, arg4)
    check_state(o)
    ccall(
        (:lammps_set_fix_external_callback, liblammps),
        Cvoid,
        (LAMMPSPtr, Cstring, Ptr{Cvoid}, Ptr{Cvoid}),
        o.ptr, arg2, callback, arg4
    )
end

has_package(name::String) = uppercase(name) in PACKAGES

function has_gzip_support()
    return Bool(ccall((:lammps_config_has_gzip_support, liblammps), Cint, ()))
end

function has_png_support()
    return Bool(ccall((:lammps_config_has_png_support, liblammps), Cint, ()))
end

function has_jpeg_support()
    return Bool(ccall((:lammps_config_has_jpeg_support, liblammps), Cint, ()))
end

function has_ffmpeg_support()
    return Bool(ccall((:lammps_config_has_ffmpeg_support, liblammps), Cint, ()))
end

function has_exceptions()
    return Bool(ccall((:lammps_config_has_exceptions, liblammps), Cint, ()))
end

function pair_neighbor_list_index(
    o::LAMMPSObject, style::String; exact = true, nsub = 0, request = 0
)
    index = ccall(
        (:lammps_find_pair_neighlist, liblammps),
        Cint,
        (LAMMPSPtr, Cstring, Cint, Cint, Cint),
        o.ptr, style, exact, nsub, request
    )
    if index == -1
        throw("neighbor list index not found")
    end
    return Int(index)
end

function fix_neighbor_list_index(o::LAMMPSObject, id; request = 0)
    index = ccall(
        (:lammps_find_fix_neighlist, liblammps),
        Cint,
        (LAMMPSPtr, Cstring, Cint),
        o.ptr, id, request
    )
    if index == -1
        throw("neighbor list index not found")
    end
    return Int(index)
end

function compute_neighbor_list_index(o::LAMMPSObject, id; request = 0)
    index = ccall(
        (:lammps_find_compute_neighlist, liblammps),
        Cint,
        (LAMMPSPtr, Cstring, Cint),
        o.ptr, id, request
    )
    if index == -1
        throw("neighbor list index not found")
    end
    return Int(index)
end

function neighbor_list_size(o::LAMMPSObject, index)
    n = ccall(
        (:lammps_neighlist_num_elements, liblammps),
        Cint,
        (LAMMPSPtr, Cint),
        o.ptr, index
    )
    return Int(n)
end

function neighbor_list(o::LAMMPSObject, idx, element)
    check_state(o)
    nref = Ref{Cint}()
    idref = Ref{Cint}()
    ptrref = Ref{Ptr{Cint}}(C_NULL)
    ccall(
        (:lammps_neighlist_element_neighbors, liblammps),
        Cvoid,
        (LAMMPSPtr, Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{Ptr{Cint}}),
        o.ptr, idx, element, idref, nref, ptrref
    )
    return idref[], unsafe_wrap(Array, ptrref[], nref[])
end

function create_atoms!(o::LAMMPSObject, ids, types, xs, vs, image, shrink)
    check_state(o)
    c_ids = cconvert(Vector{Cint}, ids)
    c_types = cconvert(Vector{Cint}, types)
    c_xs = cconvert(Vector{Cdouble}, xs)
    c_vs = cconvert(Vector{Cdouble}, vs)
    c_image = cconvert(Vector{Cdouble}, image)
    ccall(
        (:lammps_create_atoms, liblammps),
        Cvoid,
        (
            LAMMPSPtr, Cint, Ref{Cint}, Ref{Cint},
            Ref{Cdouble}, Ref{Cdouble}, Ref{Cint}, Cint
        ),
        o.ptr, length(c_ids), c_ids, c_types, c_xs, c_vs, c_image, shrink
    )
    return nothing
end

function __init__()
    if !MPI.Initialized()
        MPI.Init()
    end

    # Write the installed packages list to global variable `PACKAGES`
    package_count = ccall((:lammps_config_package_count, liblammps), Cint, ())
    buffer = Vector{UInt8}(undef, 32)
    buffer_ptr = pointer(buffer)
    empty!(PACKAGES)
    for i in 0:(package_count - 1)
        ccall(
            (:lammps_config_package_name, liblammps),
            Cint,
            (Cint, Cstring, Cint),
            i, buffer_ptr, length(buffer)
        )
        push!(PACKAGES, unsafe_string(buffer_ptr))
    end
end

end # module


text:
# built-in Autumn types: Int, Bool, Cell, Position, Click

# RETURN TYPE: OBJECT
"""
function objClicked(click::Union{Click, Nothing}, objects::AbstractArray)::Object
function removeObj(obj::Object)::Object
function updateObj(obj::Object, field::String, value)::Object
function rotate(object::Object)::Object
function rotateNoCollision(object::Object)::Object
function move(object::Object, position::Position)::Object
function move(object::Object, x::Int, y::Int)::Object
function moveLeft(object::Object)::Object
function moveRight(object::Object)::Object
function moveUp(object::Object)::Object
function moveDown(object::Object)::Object
function moveNoCollision(object::Object, position::Position)::Object
function moveNoCollision(object::Object, x::Int, y::Int)
function moveLeftNoCollision(object::Object)::Object
function moveRightNoCollision(object::Object)::Object
function moveUpNoCollision(object::Object)::Object
function moveDownNoCollision(object::Object)::Object
function moveWrap(object::Object, position::Position)::Object
function moveLeftWrap(object::Object)::Object
function moveRightWrap(object::Object)::Object
function moveUpWrap(object::Object)::Object
function moveDownWrap(object::Object)::Object
function nextLiquid(object::Object)::Object
function nextSolid(object::Object)::Object
"""

# RETURN TYPE: INT
"""
function distance(position1::Position, position2::Position)::Int
function distance(object1::Object, object2::Object)::Int
function distance(object::Object, position::Position)::Int
"""

# RETURN TYPE: BOOL
"""
function isWithinBounds(obj::Object)::Bool
function clicked(click::Union{Click, Nothing}, object::Object)::Bool
function clicked(click::Union{Click, Nothing}, objects::AbstractArray)::Bool
function clicked(click::Union{Click, Nothing}, x::Int, y::Int)::Bool
function clicked(click::Union{Click, Nothing}, pos::Position)::Bool
function intersects(obj1::Object, obj2::Object)::Bool
function intersects(obj1::Object, obj2::Array{<:Object})::Bool
function intersects(list1, list2)::Bool
function intersects(object::Object)::Bool
function isWithinBounds(position::Position)::Bool
function isFree(position::Position)::Bool
function isFree(click::Union{Click, Nothing})::Bool
function adjacent(position1::Position, position2::Position)::Bool
function adjacent(cell1::Cell, cell2::Cell)::Bool

"""


text:
using ChemistryFeaturization
using Test

@testset "SMILES_functions" begin
    include("SMILES_tests.jl")
end

@testset "Featurization" begin
    include("atomfeat_tests.jl")
end

@testset "Atomic Graphs" begin
    include("atomgraph_tests.jl")
end


text:
@test ("hi", "b") ⊕ ("!", ".") == ("hi!", "b.")
@test neutral(Tuple{String, String, String}) == ("", "", "")