text:
@doc """
  This function calculates the potential flux psi from entropy stability
  theory.  Methods are available for both 2 and 3D
"""

function getPsi(params::ParamType{2, :conservative}, q_vals::AbstractVector, nrm::AbstractVector)
  return nrm[1]*q_vals[2] + nrm[2]*q_vals[3]
end

function getPsi(params::ParamType{3, :conservative}, q_vals::AbstractVector, nrm::AbstractVector)
  return nrm[1]*q_vals[2] + nrm[2]*q_vals[3] +nrm[3]*q_vals[4]
end

function getPsi(params, qL::AbstractMatrix{Tsol}, qR::AbstractMatrix{Tsol}, nrm::AbstractVector) where Tsol

  numDofPerNode, numNodesPerElement = size(qL)
  psiL = zeros(Tsol, numNodesPerElement)
  psiR = zeros(psiL)
  dim = length(nrm)
  # calculate psi at volume nodes
  for i=1:numNodesPerElement
    qL_i = qL[:, i]
    qR_i = qR[:, i]

    psi_vecL = psi_vec(params, qL_i)
    psi_vecR = psi_vec(params, qR_i)
    for d=1:dim
      psiL[i] += nrm[d]*psi_vecL[d]
      psiR[i] += nrm[d]*psi_vecR[d]
    end
  end

  return psiL, psiR
end

@doc """
  This function computes the integral of the potential flux over an interface

  Inputs:
    params: a ParamType
    iface: the Interface
    sbpface: an AbstractFace
    nrm_scaled: the scaled normal vector at the face nodes in x-y space
    qL: the solution at the volume nodes of elementL
    qR: the solution at the volume nodes of elementR

  Outputs:
    rhs: the value of the integral
"""
function computeInterfacePotentialFlux(
                params::ParamType{Tdim, :conservative, Tsol, Tres},
                iface::Interface, sbpface, nrm_scaled::AbstractMatrix,
                qL::AbstractMatrix, qR::AbstractMatrix) where {Tdim, Tsol, Tres}
# compute the potential flux then compute the reduction with Eface

  rhs = reduceEface(params, iface, sbpface, nrm_scaled, dim, qL, qR)
  return rhs
end

@doc """
  Computes 1^T(E_faceL * psiL + E_faceR * psiR), where E_face performs the
  face integral over the specified face and psi is the potential flux
  (from entropy stability theory).  The dir argument indicates whether
  the normal vector used for the integral is x, y, or z, ie. in

  integral psi dot n dGamma

  it determines the normal vector n

  TODO: this is currently broken for SparseFace SBP operators
"""
function reduceEface(
params::ParamType{Tdim, :conservative, Tsol, Tres, Tmsh},
iface::Interface, sbpface, nrm_scaled::AbstractMatrix,
qL::AbstractMatrix, qR::AbstractMatrix) where {Tdim, Tsol, Tres, Tmsh}
  # compute Ex_gamma kappa * psiL + Ex_gamma_nu * psiR, where x is one
  # of either x or y, as specified by dir

  # interpolate q to the faces
  qfaceL = zeros(Tsol, size(qL, 1), sbpface.numnodes)
  qfaceR = zeros(qfaceL)
  interiorFaceInterpolate!(sbpface, iface, qL, qR, qfaceL, qfaceR)

  psiL = zeros(Tres, sbpface.numnodes)
  psiR = zeros(Tres, sbpface.numnodes)
  # compute psi dot n
  for i=1:sbpface.numnodes
    nrm_i = sview(nrm_scaled, :, i)
    psiL[i] = getPsi(params, sview(qfaceL, :, i), nrm_i)
    psiR[i] = -getPsi(params, sview(qfaceR, :, i), nrm_i)  # negative sign for
                                                           # direction
  end

  val = integrateBoundaryFunctional!(sbpface, iface.faceL, psiL)
  val += integrateBoundaryFunctional!(sbpface, iface.faceR, psiR)

  return val
end

#TODO: this can be made more efficient once SBP stores E
function computeVolumePotentialFlux(params::ParamType{Tdim, :conservative, Tsol, Tres}, sbp, q_i::AbstractMatrix, dxidx) where {Tdim, Tsol, Tres}

  numDofPerNode, numNodesPerElement = size(q_i)

  nrm = params.nrm
  # calculate psi vector
  psi = zeros(numNodesPerElement, 2)
  for j=1:numNodesPerElement
    q_j = ro_sview(q_i, :, j)
    for d=1:Tdim
      fill!(nrm, 0.0)
      nrm[d] = 1
      psi[j, d] = getPsi(params, q_j, nrm)
    end
  end

  rhs_reduced = zero(Tres)
  for d=1:Tdim
#    println("dimension = ", d)
    E_d = (sbp.Q[:, :, d] + sbp.Q[:, :, d].')
    psi_nrm = dxidx[d, 1, 1]*psi[:, 1] + dxidx[d, 2, 1]*psi[:, 2]
#    println("psi_nrm = ", psi_nrm)
    val = sum(E_d*psi_nrm)
#    println("rhs_reduced = ", val)

    rhs_reduced += val
  end

  return -rhs_reduced
end

@doc """
  Calculate the integral of entropy over the entire domain for the given
  solution vector q_vec.
  This performs an MPI blocking collective operation, so all processes must
  call this function at the same time.

  The (scalar) value is returned.
"""
function calcEntropyIntegral(mesh::AbstractMesh{Tmsh},
             sbp, eqn::EulerData{Tsol, Tres, Tdim}, opts, q_vec::AbstractVector) where {Tsol, Tres, Tmsh, Tdim}

  val = zero(Tsol)
  for i=1:mesh.numDofPerNode:mesh.numDof
    q_vals = ro_sview(eqn.q_vec, i:(i+(mesh.numDofPerNode - 1)))
#      s = calcEntropy(eqn.params, q_vals)
    s = calcEntropyIR(eqn.params, q_vals)
#      val += real(q_vals[1]*s)*eqn.M[i]
    val += s*eqn.M[i]
  end

  val2 = MPI.Allreduce(val, MPI.SUM, eqn.comm)
  return val2
end

@doc """
  Compute w^T * res_vec, where w is the vector of entropy variables.
  This performs an MPI blocking collective operation, so all processes must
  call this funciton at the same time.

  The (scalar) value is returned.
"""
function contractResEntropyVars(
             mesh::AbstractDGMesh{Tmsh},
             sbp, eqn::EulerData{Tsol, Tres, Tdim}, opts, q_vec::AbstractVector,
             res_vec::AbstractVector) where {Tsol, Tres, Tmsh, Tdim}

  val = zero(Tres)
  w_vals = eqn.params.contract_res_entropy_vars_data.w_vals
  for i=1:mesh.numDofPerNode:mesh.numDof
    q_vals_i = ro_sview(eqn.q_vec, i:(i+mesh.numDofPerNode - 1))
    convertToEntropy(eqn.params, q_vals_i, w_vals)
    scale!(w_vals, 1./eqn.params.gamma_1)  # the IR entropy variables are
                                           # scaled by 1/gamma compared to
                                           # Hugh's
    res_vals = sview(res_vec, i:(i+mesh.numDofPerNode - 1))
    for p=1:mesh.numDofPerNode
      val += w_vals[p]*res_vals[p]
    end
  end

  val = MPI.Allreduce(val, MPI.SUM, eqn.comm)
  return val
end

"""
  Like contractResEntropyVars, but uses a special summation technique
"""
function contractResEntropyVars2(
             mesh::AbstractDGMesh{Tmsh},
             sbp, eqn::EulerData{Tsol, Tres, Tdim}, opts, q_vec::AbstractVector,
             res_vec::AbstractVector) where {Tsol, Tres, Tmsh, Tdim}

  vals = zeros(Tres, mesh.numNodes)
  w_vals = eqn.params.contract_res_entropy_vars_data.w_vals
  idx = 1
  for i=1:mesh.numDofPerNode:mesh.numDof
    q_vals_i = ro_sview(eqn.q_vec, i:(i+mesh.numDofPerNode - 1))
    convertToEntropy(eqn.params, q_vals_i, w_vals)
    scale!(w_vals, 1./eqn.params.gamma_1)  # the IR entropy variables are
                                           # scaled by 1/gamma compared to
                                           # Hugh's
    res_vals = sview(res_vec, i:(i+mesh.numDofPerNode - 1))
    for p=1:mesh.numDofPerNode
      vals[idx] += w_vals[p]*res_vals[p]
    end
    idx += 1
  end

  # now do the sum in the order of decreasing absolute value
  pvec = sortperm(abs(vals))
  val = zero(Tres)
  for i=mesh.numNodes:-1:1  # decreasing absolute value
    val += vals[pvec[i]]
  end

  val = MPI.Allreduce(val, MPI.SUM, eqn.comm)
  return val
end

@doc """
  Compute the SBP approximation to integral q dOmega, ie. the mass matrix
  times the vector of conservative variables at each node in the mesh.

  This function returns an array of length numDofPerNode containing the
  result for each variable.
"""
function integrateQ( mesh::AbstractDGMesh{Tmsh},
             sbp, eqn::EulerData{Tsol, Tres, Tdim}, opts, q_vec::AbstractVector) where {Tsol, Tres, Tmsh, Tdim}

  vals = zeros(Tsol, mesh.numDofPerNode)
  for i=1:mesh.numDofPerNode:mesh.numDof
    q_vals_i = ro_sview(eqn.q_vec, i:(i+mesh.numDofPerNode - 1))
    w_val = eqn.M[i]
    for j=1:length(q_vals_i)
      vals[j] += w_val*q_vals_i[j]
    end
  end

  vals2 = zeros(vals)
  MPI.Allreduce!(vals, vals2, MPI.SUM, eqn.comm)

  return vals2
end

@doc """
  Computes the net potential flux integral over all interfaces, where the
  potential flux is calculated from q_arr.
  Does not work in parallel
"""
function calcInterfacePotentialFlux(
           mesh::AbstractDGMesh{Tmsh}, sbp,
           eqn::EulerData{Tsol, Tres, Tdim}, opts,
           q_arr::Abstract3DArray) where {Tsol, Tres, Tdim, Tmsh}

  if mesh.commsize != 1
    throw(ErrorException("cannot perform reduction over interfaces in parallel"))
  end

  val = zero(Tres)
  for i=1:mesh.numInterfaces
    iface = mesh.interfaces[i]
    elL = iface.elementL
    elR = iface.elementR
    qL = ro_sview(q_arr, :, :, elL)
    qR = ro_sview(q_arr, :, :, elR)
    nrm_scaled = ro_sview(mesh.nrm_face, :, :, i)

    bndry_potentialflux = -computeInterfacePotentialFlux(eqn.params, iface, mesh.sbpface, nrm_scaled, qL, qR)
    val += bndry_potentialflux
  end

  return val
end

@doc """
  Computes the net potential flux from the volume terms, calulated from
  q_arr.
  This performs a blocking MPI collective operation, so all processes must
  call this function at the same time.
"""
function calcVolumePotentialFlux(mesh::AbstractMesh{Tmsh},
               sbp, eqn::EulerData{Tsol, Tres}, opts,
               q_arr::Abstract3DArray) where {Tsol, Tres, Tmsh}

  val = zero(Tres)
  for i=1:mesh.numEl
    q_i = ro_sview(q_arr, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    volume_potentialflux = -computeVolumePotentialFlux(eqn.params, sbp, q_i, dxidx_i)
    val += volume_potentialflux
  end

  val = MPI.Allreduce(val, MPI.SUM, mesh.comm)

  return val
end

"""
  Calculates ( 1/(2V) )*integral(rho * omega dot omega dV), where V is the
  volume of the mesh omega is the vorticity vector, and rho is the density.
  3D, conservative variables only.  This should work for CG and DG, but
  has only been tested for the latter

  Inputs:
    mesh: an AbstractMesh
    sbp: an SBP operator
    eqn: an EulerData object
    opts: options dictionary
    q_arr: a 3D array of conservative variables

  Outputs:
    val: the value of the integral (over the entire domain, not just the
         part owned by this process)

  Aliasing restrictions: see calcVorticity
"""
function calcEnstrophy(mesh::AbstractMesh{Tmsh}, sbp,
                       eqn::EulerData{Tsol, Tres}, opts,
                       q_arr::Abstract3DArray{Tsol}) where {Tsol, Tres, Tmsh}

  @assert mesh.dim == 3
  Tdim = 3

  val = zero(Tres)
  vorticity = zeros(Tres, Tdim, mesh.numNodesPerElement)
  for i=1:mesh.numEl
    q_i = ro_sview(q_arr, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    jac_i = ro_sview(mesh.jac, :, i)

    calcVorticity(eqn.params, sbp, q_i, dxidx_i, jac_i, vorticity)

    for j=1:mesh.numNodesPerElement

      rho = q_i[1, j]

      # accumulate vorticity dot vorticity
      vorticity_mag = zero(Tres)
      for k=1:Tdim
        vorticity_mag += vorticity[k, j]*vorticity[k, j]
      end

      val += sbp.w[j]*rho*vorticity_mag/jac_i[j]
    end  # end loop j
  end  # end loop i

  val = MPI.Allreduce(val, MPI.SUM, mesh.comm)

  return 0.5*val/mesh.volume
end

function getVorticity(mesh::AbstractMesh{Tmsh}, sbp, eqn::EulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}
# get an array of the vorticity at every node
# also computes enstrophy

  new_field = zeros(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
  vorticity = zeros(Tres, mesh.dim, mesh.numNodesPerElement)
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    jac_i = ro_sview(mesh.jac, :, i)

    calcVorticity(eqn.params, sbp, q_i, dxidx_i, jac_i, vorticity)

    # copy into new_field
    for j=1:mesh.numNodesPerElement
      for k=1:3
        new_field[k, j, i] = real(vorticity[k, j])
      end
      vorticity_mag = real(sqrt(vorticity[1, j]^2 + vorticity[2, j]^2 + vorticity[3, j]^2))
      new_field[4, j, i] = real(vorticity_mag)
      new_field[5, j, i] = real(eqn.q[1, j, i]*vorticity_mag)

    end

  end

  return new_field
end

"""
  This function calculates ( 1/(2*V) )*integral(rho * v dot v dV), where
  V is the volume of the mesh, v is the velocity vector, and rho is the density.
  This is the total kinetic energy normalized by the volume of the domain
  Conservative variables only.

  This function contains an MPI blocking collective operation.  It must be
  called by all processes at the same time.

  This function relies on the sequential numbering of dofs on the same node

  Inputs:
    mesh: a mesh
    sbp: an SBP Operator
    eqn: an EulerData object
    opts: options dictionary
    q_vec: the vector of conservative variables for the entire mesh

  Outputs:
    val: the value of the integral (over the entire domain, not just teh part
         owned by this procss)

  Aliasing restrictions: none
"""
function calcKineticEnergy(mesh::AbstractMesh{Tmsh}, sbp,
   eqn::EulerData{Tsol, Tres, Tdim}, opts,
   q_vec::AbstractVector{Tsol}) where {Tsol, Tres, Tdim, Tmsh}

  val = zero(Tsol)
  for i=1:mesh.numDofPerNode:mesh.numDof
#    println(f, "node ", i)
    rho_i = q_vec[i]

    # compute v dot v
    v_magnitude = zero(Tsol)
    for j=1:Tdim
      v_j = q_vec[i+j]/rho_i
      v_magnitude += v_j*v_j
    end

    val += eqn.M[i]*rho_i*v_magnitude
  end

  val = MPI.Allreduce(val, MPI.SUM, mesh.comm)

  val = 0.5*val/mesh.volume

  return val
end

"""
  This function calclates the time derivative of calcKineticEnergy.

  The idea is to expand the left hand side of d rho*u/dt = res using the
  product rule and solve for du/dt, then use it to compute the integral.
  Inputs:
    mesh: a mesh
    sbp: a SBP operator
    eqn: an EulerData
    opts: options dictionary
    q_vec: vector of conservative variables for the entire mesh
    res_vec: residual vector (dq/dt) of entire mesh

  Aliasing restrictions: none
"""
function calcKineticEnergydt(mesh::AbstractMesh{Tmsh},
      sbp, eqn::EulerData{Tsol, Tres, Tdim}, opts,
      q_vec::AbstractVector{Tsol},
      res_vec::AbstractVector{Tres}) where {Tsol, Tres, Tdim, Tmsh}

  val = zero(Tres)
  for i=1:mesh.numDofPerNode:mesh.numDof
    rho_i = q_vec[i]
    drhodt = res_vec[i]  # time derivative of rho
    # accumulate v dot rho*dv/dt
    term_i = zero(Tres)
    for j=1:Tdim
      v_j = q_vec[i + j]/rho_i
      rho_dvdt = res_vec[i + j] - drhodt*v_j

      term_i += v_j*rho_dvdt
    end

    val += eqn.M[i]*term_i
  end

  val = MPI.Allreduce(val, MPI.SUM, mesh.comm)

  return val/mesh.volume
end


text:
using Base.Threads
using LoopVectorization
import TriangleMesh
using Printf
using CairoMakie
import AbstractPlotting.GeometryBasics

#--------------------------------------------------------------------#

function tplot(Mesh, v)
    p   = [AbstractPlotting.GeometryBasics.Polygon( Point2f0[ (Mesh.point[1, Mesh.cell[j,i]],Mesh.point[2, Mesh.cell[j,i]]) for j=1:3]    )
    for i in 1:Mesh.n_cell]
    display(poly(p, color = v ))
    return
end

function Tanalytic!( mesh, xc , yc, T, nvert_ele, a, b, c, d, alp, bet, vec )
    # Evaluate T analytic on barycentres: vec 0 or 1
    if vec==0 # Loop version vectorized at compilation using @avx
        @avx for i=1:mesh.n_cell
            xc[i] = 0;
            yc[i] = 0;
            for j=1:nvert_ele
                xc[i] += 1.0/nvert_ele*( mesh.point[1,mesh.cell[j,i]] )
                yc[i] += 1.0/nvert_ele*( mesh.point[2,mesh.cell[j,i]] )
            end
        end
        @avx for i=1:mesh.n_cell
            T[i] = exp(alp.sin(a*xc[i] + c*yc[i]) + bet*cos(b*xc[i] + d*yc[i]));
        end
    else # Standard vectorization
        # X    = mesh.point[1,mesh.cell[:,:]]'  # X coordinates of vertices of each element
        # Y    = mesh.point[2,mesh.cell[:,:]]'  # Y coordinates of vertices of each element
        # sumX = sum(X, dims=2)
        # sumY = sum(Y, dims=2)
        # @. xc = 1.0/nvert_ele * sumX[:]
        # @. yc = 1.0/nvert_ele * sumY[:]
    for j=1:nvert_ele # as fast as above and less cryptic
        @. xc = 0.0
        @. yc = 0.0
        @. xc = xc + 1.0/nvert_ele * mesh.point[1,mesh.cell[j,:]]
        @. yc = yc + 1.0/nvert_ele * mesh.point[2,mesh.cell[j,:]]
    end
    @. T = exp(alp*sin(a*xc + c*yc) + bet*cos(b*xc + d*yc));
end
return
end

#--------------------------------------------------------------------#

function main()

# Create sides of mesh
xmin, xmax = 0, 1
ymin, ymax = 0, 1
nx, ny     = 20, 20

# Four corners of the domain
px = [0.0 1.0 1.0 0.0]
py = [0.0 0.0 1.0 1.0]
sx = [ 1 2 3 4 ]
sy = [ 2 3 4 1 ]
st = [ 1 1 1 1 ]          # segment markers
p  = vcat(px, py)         # points
s  = vcat(sx, sy)         # segments

# TriangleMesh.Polygon_pslg:
# n_point           :: Int32
# point             :: Array{Float64,2}
# n_point_marker    :: Int32
# point_marker      :: Array{Int32,2}
# n_point_attribute :: Int32
# point_attribute   :: Array{Float64,2}
# n_segment         :: Int32
# segment           :: Array{Int32,2}
# segment_marker    :: Array{Int32,1}
# n_hole            :: Int32
# hole              :: Array{Float64,2}

# Triangulation
holes    = Array{Float64}(undef,2,0)
domain   = TriangleMesh.Polygon_pslg(size(p,2), p, 0, Array{Int64}(undef,2,0), 0, Array{Float64}(undef,2,0),  size(s,2), s, st[:], 0, holes)
switches = "Dpenq33IAa0.01"#"Dpenq33o2IAa"
mesh     = TriangleMesh.create_mesh(domain, switches)
nvert_el = 3 # vertices per element
println(size(mesh.cell))

# Generate function to be visualised on mesh
alp = 0.1; bet = 0.3; a = 5.1; b = 4.3; c = -6.2; d = 3.4;
xc = zeros(mesh.n_cell)
yc = zeros(mesh.n_cell)
T  = zeros(mesh.n_cell)
ncalls = 4
# A) Loop version with @avx
vec = 0
@printf("Looped, %d times:\n", ncalls)
for icall=1:ncalls
    @time Tanalytic!(mesh, xc ,yc , T, nvert_el, a, b, c, d, alp, bet, vec)
end
# B) Standard vectorisation
vec = 1
@printf("Vectorized, %d times:\n", ncalls)
for icall=1:ncalls
    @time Tanalytic!(mesh, xc ,yc , T, nvert_el, a, b, c, d, alp, bet, vec)
end
# # Visualise
tplot(mesh, T)

return
end

main()

ext: jl


text:
using DDD
using Test, SparseArrays, StaticArrays
cd(@__DIR__)

@testset "Forces" begin
    fileDislocationParameters = "./testData/sampleDislocationParameters.json"
    fileMaterialParameters = "./testData/sampleMaterialParameters.json"
    fileFEMParameters = "./testData/sampleFEMParameters.json"
    fileIntegrationParameters = "./testData/sampleIntegrationParameters.json"
    fileSlipSystem = "./testData/BCC.json"
    fileDislocationLoop = "./testData/samplePrismShear.json"
    dlnParams, matParams, femParams, intParams, slipSystems, dislocationLoop =
        loadParameters(
            fileDislocationParameters,
            fileMaterialParameters,
            fileFEMParameters,
            fileIntegrationParameters,
            fileSlipSystem,
            fileDislocationLoop,
        )
    regularCuboidMesh = buildMesh(matParams, femParams)
    numFEMNode = regularCuboidMesh.numNode
    dx = regularCuboidMesh.dx
    dy = regularCuboidMesh.dy
    dz = regularCuboidMesh.dz
    segLen = dx / 8

    f = spzeros(3 * numFEMNode)
    f[[112, 118, 133, 141, 213, 244, 262, 272, 317]] = [
        0.43048187784858616,
        0.22724536603830137,
        0.4340867899691503,
        0.6660863546953892,
        0.30358515797696106,
        0.2945958951093859,
        0.7278367502911502,
        0.7095924334694701,
        0.1642050526375538,
    ]
    fHat = spzeros(3 * numFEMNode)
    fHat[[32, 48, 55, 88, 138, 148, 191, 230, 253, 335]] = [
        0.09706224225842108,
        0.07773687633248638,
        0.13682398802299178,
        0.4752286167553166,
        0.7423196193496164,
        0.8286077556473421,
        0.7023632196408749,
        0.9813639162461198,
        0.5296701796678411,
        0.5523797553266823,
    ]
    u = spzeros(3 * numFEMNode)
    u[[30, 127, 195, 221, 316, 325, 338, 348, 370]] = [
        0.8792592573507609,
        0.8430664083925272,
        0.4711050560756602,
        0.4860071865093816,
        0.7905698600135145,
        0.39047211692578077,
        0.6545538020629462,
        0.5446700211111557,
        0.8865721648558644,
    ]
    uHat = spzeros(3 * numFEMNode)
    uHat[[91, 126, 130, 195, 217, 226, 229, 256, 281, 293, 309, 342]] = [
        0.5231621885339968,
        0.5771429489788034,
        0.7151190318538345,
        0.7283662326812077,
        0.6314274719472075,
        0.9814688915693632,
        0.5672795171250207,
        0.002712918060655989,
        0.1788941754890383,
        0.188299784057536,
        0.8489027048214433,
        0.029995302953659708,
    ]
    forceDisplacement =
        ForceDisplacement(nothing, uHat * 1000, u * 1000, nothing, fHat * 1000, f * 1000)

    prismSquare = DislocationLoop(;
        loopType = loopPrism(),    # Prismatic loop, all segments are edge segments.
        numSides = 4,   # 5-sided loop.
        nodeSide = 2,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 2,  # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystems,  # Slip plane of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(
            0 + segLen,
            0 + segLen,
            0 + segLen,
            dx - segLen,
            dy - segLen,
            dz - segLen,
        ),  # Distribution range
        dist = Rand(),  # Loop distribution.
    )

    shearSquare = DislocationLoop(;
        loopType = loopShear(),    # Prismatic loop, all segments are edge segments.
        numSides = 4,   # 5-sided loop.
        nodeSide = 2,   # One node per side, if 1 nodes will be in the corners.
        numLoops = 2,  # Number of loops of this type to generate when making a network.
        segLen = segLen * SVector{8}(ones(8)),  # Length of each segment between nodes, equal to the number of nodes.
        slipSystemIdx = 1, # Slip System (assuming slip systems are stored in a file, this is the index).
        slipSystem = slipSystems,  # Slip plane of the segments.
        label = SVector{8, nodeTypeDln}(1, 1, 1, 1, 1, 1, 1, 1),    # Node labels, has to be equal to the number of nodes.
        buffer = 0,   # Buffer to increase the dislocation spread.
        range = SMatrix{3, 2, Float64}(
            0 + segLen,
            0 + segLen,
            0 + segLen,
            dx - segLen,
            dy - segLen,
            dz - segLen,
        ),  # Distribution range
        dist = Rand(),  # Loop distribution.
    )
    network = DislocationNetwork([prismSquare, shearSquare])
    numSeg = network.numSeg[1]
    network.coord[:, 1:numSeg] = [
        1201.5837679702274 1283.233426063 1424.6547823003095 1566.076138537619 1484.4264804448464 1402.7768223520739 1261.3554661147643 1119.9341098774548 304.6214867750916 386.2711448678642 527.6925011051737 669.1138573424832 587.4641992497106 505.81454115693805 364.39318491962854 222.97182868231897 1270.8257124306385 1352.4753705234111 1467.9454243613363 1583.4154781992615 1501.765820106489 1420.1161620137163 1304.6461081757911 1189.176054337866 462.2031932442291 543.8528513370018 659.3229051749269 774.792959012852 693.1433009200794 611.4936428273068 496.02358898938166 380.5535351514565
        696.4782710631935 614.8286129704209 756.2499692077304 897.6713254450399 979.3209835378125 1060.970641630585 919.5492853932756 778.1279291559661 302.5604363129223 220.91077822014972 362.33213445745923 503.75349069476874 585.4031487875413 667.0528068803139 525.6314506430044 384.210094405695 1055.3880794417687 973.7384213489961 858.2683675110709 742.7983136731457 824.4479717659183 906.0976298586909 1021.5676836966161 1137.0377375345413 1195.4446644351433 1113.7950063423707 998.3249525044455 882.8548986665203 964.5045567592929 1046.1542148520655 1161.6242686899907 1277.094322527916
        1219.6759003965913 1382.9752165821365 1382.9752165821365 1382.9752165821365 1219.6759003965913 1056.376584211046 1056.376584211046 1056.376584211046 1604.4473849682922 1767.7467011538374 1767.7467011538374 1767.7467011538374 1604.4473849682922 1441.148068782747 1441.148068782747 1441.148068782747 951.679209138293 1114.9785253238383 999.5084714859131 884.038417647988 720.7391014624428 557.4397852768976 672.9098391148227 788.3798929527478 1431.0547935274708 1594.354109713016 1478.8840558750908 1363.4140020371656 1200.1146858516204 1036.8153696660752 1152.2854235040004 1267.7554773419254
    ]

    fPKTest = [
        0.0 0.0 -0.0
        0.0 0.0 0.0
        0.0 0.0 0.0
        -0.0 0.0 0.0
        0.0298244946801427 -3.935141923601883 -1.982483209141013
        -0.4565798368140385 0.4565798368140385 40.759564106652746
        11.942933454784617 -11.942933454784617 50.64827964686894
        0.0 0.0 -0.0
        0.0 0.0 -0.0
        0.0 0.0 0.0
        0.0 0.0 0.0
        -0.0 0.0 0.0
        -0.0 0.0 0.0
        0.0 -0.0 0.0
        0.0 -0.0 0.0
        0.0 0.0 -0.0
        34.47756746101095 25.45149478892432 -4.513036336043308
        -13.183687913414278 -7.188653151905362 -5.995034761508917
        -7.178181870597174 -0.03305076155172504 -7.1451311090454555
        -0.0 0.0 0.0
        -41.951995497071024 -22.03609089968592 9.957952298692552
        25.62507197894962 -0.6852147795546628 26.310286758504304
        60.72042574006913 -1.1681574909891737 61.88858323105836
        108.2078469749699 72.74142517615954 -17.73321089940518
        0.0 0.0 -0.0
        0.0 0.0 -0.0
        0.0 0.0 -0.0
        -0.0 0.0 0.0
        0.015530976306194433 0.027503198516044464 0.005986111104925015
        0.006498597861991796 0.016097702989461442 -0.00959910512746965
        0.0 -0.0 0.0
        0.0 0.0 -0.0
    ]

    fPK = calcPKForce(regularCuboidMesh, forceDisplacement, network)
    @test isapprox(fPK', fPKTest)

    network.segForce .= 0
    calcPKForce!(regularCuboidMesh, forceDisplacement, network)
    @test isapprox(fPK, network.segForce[:, 1, 1:numSeg] * 2)
    @test isapprox(network.segForce[:, 1, 1:numSeg], network.segForce[:, 2, 1:numSeg])

    idx = 5
    fPK = calcPKForce(regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(fPK, fPKTest[idx, :])

    network.segForce .= 0
    calcPKForce!(regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(fPK, network.segForce[:, 1, idx] * 2)
    @test isapprox(network.segForce[:, 1, idx], network.segForce[:, 2, idx])

    idx = [5; 10; 9; 23]
    fPK = calcPKForce(regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(fPK', fPKTest[idx, :])

    network.segForce .= 0
    calcPKForce!(regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(fPK, network.segForce[:, 1, idx] * 2)
    @test isapprox(network.segForce[:, 1, idx], network.segForce[:, 2, idx])

    idx = 15:32
    fPK = calcPKForce(regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(fPK', fPKTest[idx, :])

    network.segForce .= 0
    calcPKForce!(regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(fPK, network.segForce[:, 1, idx] * 2)
    @test isapprox(network.segForce[:, 1, idx], network.segForce[:, 2, idx])

    dlnParamsPar = DislocationParameters(;
        coreRad = dlnParams.coreRad,
        coreRadMag = dlnParams.coreRadMag,
        coreEnergy = dlnParams.coreEnergy,
        minSegLen = dlnParams.minSegLen,
        maxSegLen = dlnParams.maxSegLen,
        minArea = dlnParams.minArea,
        maxArea = dlnParams.maxArea,
        remesh = dlnParams.remesh,
        collision = dlnParams.collision,
        separation = dlnParams.separation,
        virtualRemesh = dlnParams.virtualRemesh,
        parCPU = true,
        parGPU = dlnParams.parGPU,
        dragCoeffs = dlnParams.dragCoeffs,
        mobility = dlnParams.mobility,
    )

    pentagon = DislocationLoop(;
        loopType = loopPrism(),
        numSides = 5,
        nodeSide = 1,
        numLoops = 2,
        segLen = 10 * ones(5),
        slipSystemIdx = 4,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )
    network = DislocationNetwork(pentagon; memBuffer = 1)
    network.coord[:, 6:end] .+= [10; 10; 10]
    for i in eachindex(network.segForce)
        network.segForce[i] = i
    end
    for i in eachindex(network.nodeVel)
        network.nodeVel[i] = -i
    end

    selfForce = calcSelfForce(dlnParams, matParams, network)

    f1 = [
        -1.109241905058071 0.758809517064033 -0.350432387994039
        -0.118537684538429 1.035981587065609 0.917443902527180
        1.035981587065609 -0.118537684538429 0.917443902527180
        0.758809517064033 -1.109241905058071 -0.350432387994039
        -0.567011514533144 -0.567011514533144 -1.134023029066286
        -1.109241905058071 0.758809517064033 -0.350432387994039
        -0.118537684538429 1.035981587065609 0.917443902527180
        1.035981587065609 -0.118537684538429 0.917443902527180
        0.758809517064033 -1.109241905058071 -0.350432387994039
        -0.567011514533144 -0.567011514533144 -1.134023029066286
    ]
    f2 = [
        1.109241905058071 -0.758809517064033 0.350432387994039
        0.118537684538429 -1.035981587065609 -0.917443902527180
        -1.035981587065609 0.118537684538429 -0.917443902527180
        -0.758809517064033 1.109241905058071 0.350432387994039
        0.567011514533144 0.567011514533144 1.134023029066286
        1.109241905058071 -0.758809517064033 0.350432387994039
        0.118537684538429 -1.035981587065609 -0.917443902527180
        -1.035981587065609 0.118537684538429 -0.917443902527180
        -0.758809517064033 1.109241905058071 0.350432387994039
        0.567011514533144 0.567011514533144 1.134023029066286
    ]
    @test isapprox(selfForce[1], f1')
    @test isapprox(selfForce[2], f2')

    idx = rand(1:(network.numNode[1]), Int(network.numNode[1] / 2))
    selfIdx = calcSelfForce(dlnParams, matParams, network, idx)
    @test isapprox(selfForce[1][:, idx], selfIdx[1])
    @test isapprox(selfForce[2][:, idx], selfIdx[2])

    idx = rand(1:(network.numNode[1]))
    selfIdx = calcSelfForce(dlnParams, matParams, network, idx)
    @test isapprox(selfForce[1][:, idx], selfIdx[1])
    @test isapprox(selfForce[2][:, idx], selfIdx[2])

    pentagon = DislocationLoop(;
        loopType = loopPrism(),
        numSides = 5,
        nodeSide = 1,
        numLoops = 2,
        segLen = 10 * ones(5),
        slipSystemIdx = 4,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )
    network = DislocationNetwork(pentagon, memBuffer = 1)
    remoteForceSer = calcSegSegForce(dlnParams, matParams, network)
    remoteForcePar = calcSegSegForce(dlnParamsPar, matParams, network)
    f1 = [
        0.000330019407456 0.001179598001408 0.001509617408864
        0.001578610251880 -0.000645615383196 0.000932994868684
        0.000645615383195 -0.001578610251879 -0.000932994868684
        -0.001179598001402 -0.000330019407454 -0.001509617408857
        -0.001374647041125 0.001374647041125 -0.000000000000000
        0.000330019407456 0.001179598001408 0.001509617408865
        0.001578610251881 -0.000645615383196 0.000932994868685
        0.000645615383195 -0.001578610251879 -0.000932994868684
        -0.001179598001403 -0.000330019407455 -0.001509617408857
        -0.001374647041125 0.001374647041125 -0.000000000000000
    ]
    f2 = [
        0.000330019407454 0.001179598001403 0.001509617408857
        0.001578610251880 -0.000645615383196 0.000932994868684
        0.000645615383196 -0.001578610251880 -0.000932994868684
        -0.001179598001408 -0.000330019407456 -0.001509617408864
        -0.001374647041126 0.001374647041126 0.000000000000000
        0.000330019407454 0.001179598001402 0.001509617408856
        0.001578610251878 -0.000645615383195 0.000932994868683
        0.000645615383196 -0.001578610251881 -0.000932994868685
        -0.001179598001408 -0.000330019407456 -0.001509617408864
        -0.001374647041126 0.001374647041126 0.000000000000000
    ]
    @test isapprox(remoteForceSer[:, 1, :], f1')
    @test isapprox(remoteForceSer[:, 2, :], f2')
    @test isapprox(remoteForcePar[:, 1, :], f1')
    @test isapprox(remoteForcePar[:, 2, :], f2')

    pentagon = DislocationLoop(;
        loopType = loopPrism(),
        numSides = 5,
        nodeSide = 1,
        numLoops = 2,
        segLen = 10 * ones(5),
        slipSystemIdx = 4,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )
    network = DislocationNetwork(pentagon, memBuffer = 1)
    network.coord[:, 6:end] .+= [20; 20; 20]
    remoteForceSer = calcSegSegForce(dlnParams, matParams, network)
    remoteForcePar = calcSegSegForce(dlnParamsPar, matParams, network)
    f1 = [
        0.000217825304625 0.000839481257576 0.001128278735983
        0.001146987202523 -0.000492397907541 0.000704213490464
        0.000532132124482 -0.001268053797081 -0.000764723861474
        -0.000989719230921 -0.000261682028865 -0.001314774875872
        -0.001103824816939 0.001093373550371 0.000005225633284
        0.000258286008162 0.000976109417854 0.001296054439544
        0.001312903402718 -0.000551465905248 0.000792350410865
        0.000507633513349 -0.001180265566845 -0.000725717309154
        -0.000831056838483 -0.000215895485883 -0.001116145458036
        -0.001051183335307 0.001060779799870 -0.000004798232281
    ]
    f2 = [
        0.000215895485882 0.000831056838483 0.001116145458036
        0.001180265566845 -0.000507633513349 0.000725717309154
        0.000551465905248 -0.001312903402718 -0.000792350410864
        -0.000976109417853 -0.000258286008162 -0.001296054439544
        -0.001060779799870 0.001051183335308 0.000004798232281
        0.000261682028865 0.000989719230921 0.001314774875872
        0.001268053797081 -0.000532132124482 0.000764723861474
        0.000492397907540 -0.001146987202523 -0.000704213490464
        -0.000839481257576 -0.000217825304625 -0.001128278735983
        -0.001093373550372 0.001103824816939 -0.000005225633284
    ]
    @test isapprox(remoteForceSer[:, 1, :], f1')
    @test isapprox(remoteForceSer[:, 2, :], f2')
    @test isapprox(remoteForcePar[:, 1, :], f1')
    @test isapprox(remoteForcePar[:, 2, :], f2')

    idx = rand(1:(network.numSeg[1]), Int(network.numSeg[1] / 2))
    serIdx = calcSegSegForce(dlnParams, matParams, network, idx)
    @test isapprox(remoteForceSer[:, :, idx], serIdx)
    @test isapprox(remoteForcePar[:, :, idx], serIdx)

    idx = rand(1:(network.numSeg[1]))
    serIdx = calcSegSegForce(dlnParams, matParams, network, idx)
    @test isapprox(remoteForceSer[:, :, idx], serIdx)
    @test isapprox(remoteForcePar[:, :, idx], serIdx)

    hexagonPris = DislocationLoop(;
        loopType = loopPrism(),
        numSides = 6,
        nodeSide = 1,
        numLoops = 1,
        segLen = 10 * ones(6),
        slipSystemIdx = 1,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1; 2],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )
    hexagonShear = DislocationLoop(;
        loopType = loopShear(),
        numSides = 6,
        nodeSide = 1,
        numLoops = 1,
        segLen = 10 * ones(6),
        slipSystemIdx = 1,
        slipSystem = slipSystems,
        label = nodeTypeDln[1; 2; 1; 2; 1; 2],
        buffer = 0.0,
        range = Float64[-100 100; -100 100; -100 100],
        dist = Zeros(),
    )
    network = DislocationNetwork([hexagonPris, hexagonShear], memBuffer = 1)
    remoteForceSer = calcSegSegForce(dlnParams, matParams, network)
    remoteForcePar = calcSegSegForce(dlnParamsPar, matParams, network)

    f1 = [
        0.000148330200377 0.000867173203941 -0.001163833604695
        -0.001162438253910 0.000146934849592 -0.000868568554726
        -0.001310620554548 -0.000720386254088 0.000295117150230
        -0.000148330200377 -0.000867173203940 0.001163833604693
        0.001162438253908 -0.000146934849592 0.000868568554724
        0.001310620554548 0.000720386254089 -0.000295117150230
        0.000487619084780 0.000552951524211 -0.000501376468106
        -0.000840679364846 -0.000186252085257 -0.000213189755230
        -0.001337886532722 -0.000747434694617 0.000295225919052
        -0.000487619084782 -0.000552951524212 0.000501376468101
        0.000840679364847 0.000186252085255 0.000213189755231
        0.001337886532723 0.000747434694617 -0.000295225919053
    ]
    f2 = [
        0.000146934849592 0.000868568554725 -0.001162438253909
        -0.001163833604693 0.000148330200376 -0.000867173203941
        -0.001310620554545 -0.000720386254089 0.000295117150228
        -0.000146934849592 -0.000868568554727 0.001162438253910
        0.001163833604695 -0.000148330200376 0.000867173203942
        0.001310620554545 0.000720386254088 -0.000295117150228
        0.000480410434831 0.000546521015274 -0.000507348104804
        -0.000843647543328 -0.000196923065664 -0.000210680448989
        -0.001337886532724 -0.000747434694618 0.000295225919053
        -0.000480410434829 -0.000546521015273 0.000507348104808
        0.000843647543328 0.000196923065666 0.000210680448988
        0.001337886532722 0.000747434694618 -0.000295225919052
    ]
    @test isapprox(remoteForceSer[:, 1, :], f1')
    @test isapprox(remoteForceSer[:, 2, :], f2')
    @test isapprox(remoteForcePar[:, 1, :], f1')
    @test isapprox(remoteForcePar[:, 2, :], f2')

    PKForce = calcPKForce(regularCuboidMesh, forceDisplacement, network)
    selfForce = calcSelfForce(dlnParams, matParams, network)
    remoteForceSer = calcSegSegForce(dlnParams, matParams, network)
    remoteForcePar = calcSegSegForce(dlnParamsPar, matParams, network)
    sumForceSer = (
        selfForce[1] .+ remoteForceSer[:, 1, :] + 0.5 * PKForce,
        selfForce[2] .+ remoteForceSer[:, 2, :] + 0.5 * PKForce,
    )
    sumForcePar = (
        selfForce[1] .+ remoteForcePar[:, 1, :] + 0.5 * PKForce,
        selfForce[2] .+ remoteForcePar[:, 2, :] + 0.5 * PKForce,
    )

    totalForceSer =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    totalForcePar =
        calcSegForce(dlnParamsPar, matParams, regularCuboidMesh, forceDisplacement, network)

    @test isapprox(totalForceSer[:, 1, :], sumForceSer[1])
    @test isapprox(totalForcePar[:, 2, :], sumForcePar[2])
    @test isapprox(totalForceSer[:, 1, :], totalForcePar[:, 1, :])
    @test isapprox(totalForceSer[:, 2, :], totalForcePar[:, 2, :])

    idx = rand(1:(network.numSeg[1]), Int(network.numSeg[1] / 2))
    totalForceIdx = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network,
        idx,
    )
    totalForceSer =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    totalForcePar =
        calcSegForce(dlnParamsPar, matParams, regularCuboidMesh, forceDisplacement, network)
    @test isapprox(totalForceSer[:, :, idx], totalForceIdx)
    @test isapprox(totalForcePar[:, :, idx], totalForceIdx)

    idx = rand(1:(network.numSeg[1]))
    totalForceIdx = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network,
        idx,
    )
    totalForceSer =
        calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    totalForcePar =
        calcSegForce(dlnParamsPar, matParams, regularCuboidMesh, forceDisplacement, network)
    @test isapprox(totalForceSer[:, :, idx], totalForceIdx)
    @test isapprox(totalForcePar[:, :, idx], totalForceIdx)

    # In-place functions.
    numSeg = network.numSeg[1]
    idx = rand(1:numSeg)
    self = calcSelfForce(dlnParams, matParams, network)
    selfIdx = calcSelfForce(dlnParams, matParams, network, idx)
    @test isapprox(self[1][:, idx], selfIdx[1])
    @test isapprox(self[2][:, idx], selfIdx[2])
    network.segForce .= 0
    calcSelfForce!(dlnParams, matParams, network)
    @test isapprox(self[1], network.segForce[:, 1, 1:numSeg])
    isapprox(self[2], network.segForce[:, 2, 1:numSeg])
    network.segForce .= 0
    calcSelfForce!(dlnParams, matParams, network, idx)
    @test isapprox(self[1][:, idx], network.segForce[:, 1, idx])
    @test isapprox(self[2][:, idx], network.segForce[:, 2, idx])

    ser = calcSegSegForce(dlnParams, matParams, network)
    network.segForce .= 0
    calcSegSegForce!(dlnParams, matParams, network)
    @test isapprox(network.segForce[:, :, 1:numSeg], ser)
    network.segForce .= 0
    calcSegSegForce!(dlnParams, matParams, network, idx)
    @test isapprox(network.segForce[:, :, idx], ser[:, :, idx])

    par = calcSegSegForce(dlnParamsPar, matParams, network)
    network.segForce .= 0
    calcSegSegForce!(dlnParamsPar, matParams, network)
    @test isapprox(network.segForce[:, :, 1:numSeg], par)
    network.segForce .= 0
    calcSegSegForce!(dlnParamsPar, matParams, network, idx)
    @test isapprox(network.segForce[:, :, idx], par[:, :, idx])
    @test isapprox(par, ser)

    tot = calcSegForce(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    @test isapprox(tot[:, 1, :], self[1] + ser[:, 1, :] + 0.5 * PKForce)
    @test isapprox(tot[:, 2, :], self[2] + ser[:, 2, :] + 0.5 * PKForce)
    network.segForce .= 0
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network)
    @test isapprox(tot, network.segForce[:, :, 1:numSeg])

    tot = calcSegForce(
        dlnParams,
        matParams,
        regularCuboidMesh,
        forceDisplacement,
        network,
        idx,
    )
    @test isapprox(tot[:, 1, :], self[1][:, idx] + ser[:, 1, idx] + 0.5 * PKForce[:, idx])
    @test isapprox(tot[:, 2, :], self[2][:, idx] + ser[:, 2, idx] + 0.5 * PKForce[:, idx])
    network.segForce .= 0
    calcSegForce!(dlnParams, matParams, regularCuboidMesh, forceDisplacement, network, idx)
    @test isapprox(tot, network.segForce[:, :, idx])

    network.coord .+= dx * dy * dz
    PKForce = calcPKForce(regularCuboidMesh, forceDisplacement, network)
    @test isapprox(PKForce, zeros(3, numSeg))
end

@testset "Calculating sigma_hat" begin
    fileDislocationParameters = "../inputs/simParams/sampleDislocationParameters.json"
    fileMaterialParameters = "../inputs/simParams/sampleMaterialParameters.json"
    fileFEMParameters = "../inputs/simParams/sampleFEMParameters.json"
    fileIntegrationParameters = "../inputs/simParams/sampleIntegrationParameters.json"
    fileSlipSystem = "../data/slipSystems/BCC.json"
    fileDislocationLoop = "../inputs/dln/samplePrismShear.json"
    fileIntVar = "../inputs/simParams/sampleIntegrationTime.json"
    dlnParams, matParams, femParams, intParams, slipSystems, dislocationLoop =
        loadParameters(
            fileDislocationParameters,
            fileMaterialParameters,
            fileFEMParameters,
            fileIntegrationParameters,
            fileSlipSystem,
            fileDislocationLoop,
        )

    regularCuboidMesh = buildMesh(matParams, femParams)
    numFEMNode = regularCuboidMesh.numNode

    f = spzeros(3 * numFEMNode)
    f[[112, 118, 133, 141, 213, 244, 262, 272, 317]] = [
        0.43048187784858616,
        0.22724536603830137,
        0.4340867899691503,
        0.6660863546953892,
        0.30358515797696106,
        0.2945958951093859,
        0.7278367502911502,
        0.7095924334694701,
        0.1642050526375538,
    ]
    fHat = spzeros(3 * numFEMNode)
    fHat[[32, 48, 55, 88, 138, 148, 191, 230, 253, 335]] = [
        0.09706224225842108,
        0.07773687633248638,
        0.13682398802299178,
        0.4752286167553166,
        0.7423196193496164,
        0.8286077556473421,
        0.7023632196408749,
        0.9813639162461198,
        0.5296701796678411,
        0.5523797553266823,
    ]
    u = spzeros(3 * numFEMNode)
    u[[30, 127, 195, 221, 316, 325, 338, 348, 370]] = [
        0.8792592573507609,
        0.8430664083925272,
        0.4711050560756602,
        0.4860071865093816,
        0.7905698600135145,
        0.39047211692578077,
        0.6545538020629462,
        0.5446700211111557,
        0.8865721648558644,
    ]
    uHat = spzeros(3 * numFEMNode)
    uHat[[67, 68, 69, 195, 217, 226, 229, 256, 281, 293, 309, 342]] = [
        0.5231621885339968,
        0.5771429489788034,
        0.7151190318538345,
        0.7283662326812077,
        0.6314274719472075,
        0.9814688915693632,
        0.5672795171250207,
        0.002712918060655989,
        0.1788941754890383,
        0.188299784057536,
        0.8489027048214433,
        0.029995302953659708,
    ]

    forceDisplacement =
        ForceDisplacement(nothing, uHat * 1000, u * 1000, nothing, fHat * 1000, f * 1000)
    regularCuboidMesh.coord[:, 91]
    σHat = calc_σHat(regularCuboidMesh, forceDisplacement, [1800.0, 1400.0, 200.0])
    σHatTest = [
        -1.5186579125763586 -0.6876907109455002 -1.2291546581681494
        -0.6876907109455002 -1.5861338631323667 -0.35243484259464375
        -1.2291546581681494 -0.35243484259464375 -0.8481461758746458
    ]
    @test isapprox(σHat, σHatTest)

    σHat = calc_σHat(regularCuboidMesh, forceDisplacement, [1893.0, 1356.0, 345.0])
    σHatTest = [
        -0.5560725333509386 -0.169126060231533 -0.9758794760722914
        -0.169126060231533 -0.4428115859951268 0.11892368539106829
        -0.9758794760722914 0.11892368539106829 -0.5401850945259058
    ]
    @test isapprox(σHat, σHatTest)
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Represents a Universally Unique Identifier (UUID).
    Can be built from one `UInt128` (all byte values), two `UInt64`, or four `UInt32`.
    Conversion from a string will check the UUID validity.
"""
struct UUID
    value::UInt128
end
UUID(u::UUID) = u
UUID(u::NTuple{2, UInt64}) = UUID((UInt128(u[1]) << 64) | UInt128(u[2]))
UUID(u::NTuple{4, UInt32}) = UUID((UInt128(u[1]) << 96) | (UInt128(u[2]) << 64) |
                                  (UInt128(u[3]) << 32) | UInt128(u[4]))

function convert(::Type{NTuple{2, UInt64}}, uuid::UUID)
    bytes = uuid.value
    hi = UInt64((bytes >> 64) & 0xffffffffffffffff)
    lo = UInt64(bytes & 0xffffffffffffffff)
    return (hi, lo)
end

function convert(::Type{NTuple{4, UInt32}}, uuid::UUID)
    bytes = uuid.value
    hh = UInt32((bytes >> 96) & 0xffffffff)
    hl = UInt32((bytes >> 64) & 0xffffffff)
    lh = UInt32((bytes >> 32) & 0xffffffff)
    ll = UInt32(bytes & 0xffffffff)
    return (hh, hl, lh, ll)
end

UInt128(u::UUID) = u.value

let
@noinline throw_malformed_uuid(s) = throw(ArgumentError("Malformed UUID string: $(repr(s))"))
@inline function uuid_kernel(s, i, u)
    _c = UInt32(@inbounds codeunit(s, i))
    d = __convert_digit(_c, UInt32(16))
    d >= 16 && throw_malformed_uuid(s)
    u <<= 4
    u | d
end

global UUID
function UUID(s::AbstractString)
    u = UInt128(0)
    ncodeunits(s) != 36 && throw_malformed_uuid(s)
    for i in 1:8
        u = uuid_kernel(s, i, u)
    end
    @inbounds codeunit(s, 9) == UInt8('-') || @goto error
    for i in 10:13
        u = uuid_kernel(s, i, u)
    end
    @inbounds codeunit(s, 14) == UInt8('-') || @goto error
    for i in 15:18
        u = uuid_kernel(s, i, u)
    end
    @inbounds codeunit(s, 19) == UInt8('-') || @goto error
    for i in 20:23
        u = uuid_kernel(s, i, u)
    end
    @inbounds codeunit(s, 24) == UInt8('-') || @goto error
    for i in 25:36
        u = uuid_kernel(s, i, u)
    end
    return Base.UUID(u)
    @label error
    throw_malformed_uuid(s)
end
end

parse(::Type{UUID}, s::AbstractString) = UUID(s)
function tryparse(::Type{UUID}, s::AbstractString)
    try
        return parse(UUID, s)
    catch e
        if isa(e, ArgumentError)
            return nothing
        end
        rethrow(e)
    end
end

let groupings = [36:-1:25; 23:-1:20; 18:-1:15; 13:-1:10; 8:-1:1]
    global string
    function string(u::UUID)
        u = u.value
        a = Base.StringVector(36)
        for i in groupings
            @inbounds a[i] = hex_chars[1 + u & 0xf]
            u >>= 4
        end
        @inbounds a[24] = a[19] = a[14] = a[9] = '-'
        return String(a)
    end
end

print(io::IO, u::UUID) = print(io, string(u))
show(io::IO, u::UUID) = print(io, "UUID(\"", u, "\")")

isless(a::UUID, b::UUID) = isless(a.value, b.value)

# give UUID scalar behavior in broadcasting
Base.broadcastable(x::UUID) = Ref(x)

max_stars_repo_name: FelipeLema/julia

max_issues_repo_name: FelipeLema/julia

max_forks_repo_name: FelipeLema/julia


text:
module KissCaches
export cached, @cached, iscached

function cached(storage, f, args...; kw...)
    key = makekey(storage, f, args...; kw...)
    get!(storage, key) do
        @debug """
        Adding to cache:
        f = $f
        args = $args
        kw = $kw
        key = $key
        cache = $cache
        """
        f(args...; kw...)
    end
end
# makekey(storage, f, args...; kw...) = Base.hash((f, args, kw))
makekey(storage, f, args...; kw...) = (f, args, (;kw...))

function iscached(storage, f, args...; kw...)
    k = makekey(storage, f, args...; kw...)
    haskey(storage, k)
end

function deletecached!(storage, f, args...; kw...)
    k = makekey(storage, f, args...; kw...)
    delete!(storage, k)
end

const GLOBAL_CACHE = Dict()

macro cached(ex)
    c = GLOBAL_CACHE
    esc(cachedmacro(c, ex))
end

macro cached(cache, ex)
    esc(cachedmacro(cache, ex))
end

function cachedmacro(cache, ex)
    if Meta.isexpr(ex, :call)
        if length(ex.args) >= 2 && Meta.isexpr(ex.args[2], :parameters)
            Expr(:call, cached, ex.args[2], cache, ex.args[1], ex.args[3:end]...)
        else
            Expr(:call, cached, cache, ex.args...)
        end
    else
        ex
    end
end

# function apply_cache_to_all_calls(cache, expr)
#     MacroTools.postwalk(expr) do ex
#         if Meta.isexpr(ex, :call)
#             Expr(:call, cached, cache, ex.args...)
#         else
#             ex
#         end
#     end
# end Write your package code here.

end


text:
using IncrementalInference
using Test

# setSerializationNamespace!("Main" => Main)

@testset "test solve by saving and loading basic jld..." begin

fgt = initfg()

addVariable!(fgt, :x1, ContinuousScalar)
addFactor!(fgt, [:x1], Prior(Normal()))

addVariable!(fgt, :x2, ContinuousScalar)
addFactor!(fgt, [:x1;:x2], LinearConditional(Normal(10,1)))

savejld(fgt)
# savejld(fgt, file=joinpath(Pkg.dir("IncrementalInference"),"test","testdata","compatibility_test_fg.jld"))
fgl, = loadjld()

fct = solverData(getFactor(fgt, :x1x2f1)).fnc.cpt[1].factormetadata
fcl = solverData(getFactor(fgl, :x1x2f1)).fnc.cpt[1].factormetadata

@test fct.solvefor == fcl.solvefor
@test length(fct.variableuserdata) == length(fcl.variableuserdata)

batchSolve!(fgt)
batchSolve!(fgl)

# test again after solve
@test fct.solvefor == fcl.solvefor
@test length(fct.variableuserdata) == length(fcl.variableuserdata)

# save load and repeat tests
savejld(fgt, file="tempfg.jld2")
fgl, = loadjld(file="tempfg.jld2")

Base.rm("tempfg.jld2")

fct = solverData(getFactor(fgt, :x1x2f1)).fnc.cpt[1].factormetadata
fcl = solverData(getFactor(fgl, :x1x2f1)).fnc.cpt[1].factormetadata

# @test fct.solvefor == fcl.solvefor # defaults to :null during load
@test length(fct.variableuserdata) == length(fcl.variableuserdata)

end

@testset "test backwards compatibility of loadjld from previous versions of IIF (from an existing test data file)..." begin

filename = joinpath(dirname(pathof(IncrementalInference)), "..", "test","testdata","compatibility_test_fg.jld2")
# savejld(fgt, file=filename)  # when changing to a new compat test file (human required)
fgprev, = loadjld(file=filename )

fc = solverData(getFactor(fgprev, :x1x2f1))
@test length(fc.fnc.cpt[1].factormetadata.variableuserdata) == 2
@test fc.fnc.cpt[1].factormetadata.solvefor == :null

# batchSolve!(fgprev)
solveTree!(fgprev)

end

#


text:
@testset "Closures" begin
    closures = Dict(
        "HNC" => HypernettedChain(),
        "MSA" => MeanSpherical(),
        "SMSA" => SoftMeanSpherical(),
        "HMSA" => HMSA(0.1),
        "PY" => PercusYevick(),
        "MV" => ModifiedVerlet(),
    )

    for (k, v) in closures
        @testset "$(k)" begin
            result = oz_sw(v)
            if 1.0 ∈ result.r
                first_contact = result.gr[result.r .== 1.0][1]
            else
                first_contact = maximum(result.gr)
            end
            @show "$(k)" first_contact
            # We don't know exactly the value at contact, so just check that it is a valid number
            @test first_contact isa Real
        end
    end
end


text:
using Test, Printf

import JuliaSyntax.Tokenize

# Takes 10s to run
# include("lex_yourself.jl")
@testset "lexer" begin
include("lexer.jl")
end


text:
include("../../src/AlphaShapes.jl")
using LinearAlgebraicRepresentation, ViewerGL
using TimerOutputs
Lar = LinearAlgebraicRepresentation
GL =  ViewerGL
using BenchmarkTools

"""
	pointsRand(V, VV, n, m)

Generate random points inside and otuside `(V, VV)`.
"""
function pointsRand(
		V::Lar.Points, EV::Lar.Cells, n = 1000, m = 0
	)::Tuple{Lar.Points, Lar.Points, Lar.Cells, Lar.Cells}
	classify = Lar.pointInPolygonClassification(V, EV)
	Vi = [0;0]
	Ve = [0;0]
	k1 = 0
	k2 = 0
	while k1 < n || k2 < m
		queryPoint = [rand();rand()]
		inOut = classify(queryPoint)

		if k1 < n && inOut == "p_in"
			Vi = hcat(Vi, queryPoint)
			k1 = k1 + 1;
		end
		if k2 < m && inOut == "p_out"
			Ve = hcat(Ve, queryPoint)
			k2 = k2 + 1;
		end
	end
	VVi = [[i] for i = 1 : n]
	VVe = [[i] for i = 1 : m]
	return Vi[:,2:end], Ve[:,2:end], VVi, VVe
end

filename = "./examples/examples2D/svg_files/Lar2.svg";
V,EV = Lar.svg2lar(filename);

Vi, Ve, VVi, VVe = pointsRand(V, EV, 1000, 10000);

GL.VIEW([
	GL.GLGrid(Vi, VVi, GL.COLORS[1], 1)
	GL.GLGrid(Ve, VVe, GL.COLORS[12], 1)
])

filtration = AlphaShapes.alphaFilter(Vi);

VV,EV,FV = AlphaShapes.alphaSimplex(Vi, filtration, 0.02)

points = [[p] for p in VV]
faces = [[f] for f in FV]
edges = [[e] for e in EV]
GL.VIEW(
	GL.GLExplode(Vi, [edges; faces], 1.5, 1.5, 1.5, 99, 1)
 );

filter_key = sort(unique(values(filtration)))

granular = 10

reduced_filter = filter_key[sort(abs.(rand(Int, granular).%length(filter_key)))]
reduced_filter = [reduced_filter; max(filter_key...)]

#
# Arlecchino's Lar
#

for α in reduced_filter
	@show α
	VV,EV,FV = AlphaShapes.alphaSimplex(Vi, filtration, α)
	# GL.VIEW(
	# 	GL.GLExplode(
	# 		Vi,
	# 		[[[f] for f in FV]; [[e] for e in EV]],
	# 		1., 1., 1.,	# Explode Ratio
	# 		99, 1		# Colors
	# 	)
	# )
end

#
# Appearing Colors
#

reduced_filter = [
	0.002;	0.003;	0.004;	0.005;  0.006;
	0.007;	0.008;	0.009;	0.010;	0.013;
	0.015;	0.020;	0.050;	1.000
]

for i = 2 : length(reduced_filter)
	VV0, EV0, FV0 = AlphaShapes.alphaSimplex(Vi, filtration, reduced_filter[i-1])
	VV,  EV,  FV  = AlphaShapes.alphaSimplex(Vi, filtration, reduced_filter[i])
	EV0mesh = GL.GLGrid(Vi, EV0)
	FV0mesh = GL.GLGrid(Vi, FV0)
	EVmesh = GL.GLGrid(Vi, setdiff(EV, EV0), GL.COLORS[2], 1)
	FVmesh = GL.GLGrid(Vi, setdiff(FV, FV0), GL.COLORS[7], 1)
	#GL.VIEW([EV0mesh; FV0mesh; EVmesh; FVmesh])
end

print_timer(AlphaShapes.to);


text:
immutable DataFunction{T<:Number,D<:Number,A<:AbstractArray} <: AbstractData
    index::Vector
    values::A
    modifyvalues::Bool
    generator::Function
    args::Tuple
    DataFunction(i::AbstractVector{T},d::AbstractArray{D},w::Bool,g::Function,args...) = new(collect(i),d,w,g,args)
end

function _data(::Type{Val{:function!}},i::AbstractVector,d::AbstractArray,g::Function,args...)
    @assert length(i) == size(d,1)
    DataFunction{eltype(i),eltype(d),typeof(d)}(i,d,true,g,args...)
end

function _data(::Type{Val{:function}},i::AbstractVector,g::Function,args...)
    d = g(args...)::AbstractArray
    @assert length(i) == size(d,1)
    DataFunction{eltype(i),eltype(d),typeof(d)}(i,d,false,g,args...)
end

@inline _generate(::Type{Val{true}},d::DataFunction) = d.generator(d.values,d.args...)
@inline _generate(::Type{Val{false}},d::DataFunction) = copy!(d.values,d.generator(d.args...))

numvalues(d::DataFunction) = size(d.values,1)
numvars(d::DataFunction) = size(d.values,2)
eltype(d::DataFunction) = eltype(d.values)
generate!(d::DataFunction) = (_generate(Val{d.modifyvalues},d) ; d.values)
dataindex(d::DataFunction) = d.index
datavalues(d::DataFunction) = d.values
datatypename(d::DataFunction) = "DataFunction"


text:
## ------------------------------------------------------------------
# Meta

# TODO: Add check to the data being entered.
# e.g: A label can not be its own anthagonist

# RBTagMeta
get_similars(obj::RBTagMeta) = getproperty!(() -> Set{String}(), obj, :similar)
set_similars!(obj::RBTagMeta,  t::String, ts::String...) = (_push_csv!(empty!(get_similars(obj)), t, ts...); obj)

get_antagonists(obj::RBTagMeta) = getproperty!(() -> Set{String}(), obj, :antagonist)
set_antagonists!(obj::RBTagMeta,  t::String, ts::String...) = (_push_csv!(empty!(get_antagonists(obj)), t, ts...); obj)

get_description(t::RBTagMeta) = getproperty(t, :description, "")
set_description!(t::RBTagMeta, text::String) = setproperty!(t, :description, strip(text))

get_supertag(t::RBTagMeta) = getproperty(t, :supertag, nothing)
set_supertag!(t::RBTagMeta, tag::String) = setproperty!(t, :supertag, strip(tag))

get_book(d::RBTagMeta) = get_parent(d)
set_book!(d::RBTagMeta, book::RBook) = set_parent!(d, book)

## ------------------------------------------------------------------
# Show
function Base.show(io::IO, tm::RBTagMeta)
    println(io, "RBTagMeta(\"", get_label(tm), "\")")
    print(io, " - path: '", localrelpath(tm), "'")
    text = get_description(tm)
    !isempty(text) && print(io, "\n - description: \"", _preview(io, text), "\"")
    return nothing
end


text:
# This file was generated, do not modify it. # hide
misclassification_rate(mode.(ŷ), ytrain)


text:
@testset "printing" begin
    angle = 45.0

    str = format_angle(deg2dms(angle))
    @test str == "45:0:0.0"
    strd = format_angle(deg2hms(angle), delim=["d", "m", "s"])
    @test strd == "3d0m0.0s"
    @test_throws BoundsError format_angle(deg2dms(angle), delim=(':', ' '))
end

function test_print_integration(angles, f1, f2)
    @test all(angles) do angle
        parts = f1(angle)
        str = format_angle(parts)
        res = f2(str)
        return all(res .≈ parts)
    end

    @test all(angles) do angle
        parts = f1(angle)
        str = format_angle(parts, delim=" ")
        res = f2(str)
        return all(res .≈ parts)
    end

    @test all(angles) do angle
        parts = f1(angle)
        str = format_angle(parts, delim=(":", "m", ""))
        res = f2(str)
        return all(res .≈ parts)
    end
end

@testset "printing integration" begin
    degrees = randdegree(rng, 100)
    test_print_integration(degrees, deg2dms, parse_dms)
    radians = randrad(rng, 100)
    test_print_integration(radians, rad2dms, parse_dms)
    has = randha(rng, 100)
    test_print_integration(has, ha2dms, parse_dms)
end

@testset "negatives" begin
    @test format_angle(parse_dms("-0:0:1.0")) == "-0:0:1.0"
    @test format_angle(parse_hms("-0:0:1.0")) == "-0:0:1.0"
end


text:
import NLPModels: grad, cons, jac

"""
unconstrained_check: return the infinite norm of the gradient of the objective function

`unconstrained_check( :: AbstractNLPModel, :: NLPAtX; pnorm :: Real = Inf, kwargs...)`

Require `state.gx` (filled if not provided)

See also `unconstrained2nd_check`, `optim_check_bounded`, `KKT`
"""
function unconstrained_check(pb::AbstractNLPModel, state::NLPAtX; pnorm::Real = Inf, kwargs...)
  if state.gx == _init_field(typeof(state.gx)) # should be filled if empty
    update!(state, gx = grad(pb, state.x))
  end

  return norm(state.gx, pnorm)
end

"""
unconstrained2nd_check: check the norm of the gradient and the smallest
                   eigenvalue of the hessian.

`unconstrained2nd_check( :: AbstractNLPModel, :: NLPAtX; pnorm :: Real = Inf, kwargs...)`

Require are `state.gx`, `state.Hx` (filled if not provided).

See also `unconstrained_check`, `optim_check_bounded`, `KKT`
"""
function unconstrained2nd_check(pb::AbstractNLPModel, state::NLPAtX; pnorm::Real = Inf, kwargs...)
  if state.gx == _init_field(typeof(state.gx)) # should be filled if empty
    update!(state, gx = grad(pb, state.x))
  end
  if state.Hx == _init_field(typeof(state.Hx))
    update!(state, Hx = hess(pb, state.x).data)
  end

  res = max(norm(state.gx, pnorm), max(-eigmin(Symmetric(state.Hx, :L)), 0.0))

  return res
end

"""
optim\\_check\\_bounded: gradient of the objective function projected

`optim_check_bounded( :: AbstractNLPModel, :: NLPAtX; pnorm :: Real = Inf, kwargs...)`

Require `state.gx` (filled if not provided).

See also `unconstrained_check`, `unconstrained2nd_check`, `KKT`
"""
function optim_check_bounded(pb::AbstractNLPModel, state::NLPAtX; pnorm::Real = Inf, kwargs...)
  if state.gx == _init_field(typeof(state.gx)) # should be filled if void
    update!(state, gx = grad(pb, state.x))
  end

  proj = max.(min.(state.x - state.gx, pb.meta.uvar), pb.meta.lvar)
  gradproj = state.x - proj

  return norm(gradproj, pnorm)
end

"""
constrained: return the violation of the KKT conditions
length(lambda) > 0
"""
function _grad_lagrangian(pb::AbstractNLPModel, state::NLPAtX)
  if (pb.meta.ncon == 0) & !has_bounds(pb)
    return state.gx
  elseif pb.meta.ncon == 0
    return state.gx + state.mu
  else
    return state.gx + state.mu + state.Jx' * state.lambda
  end
end

function _sign_multipliers_bounds(pb::AbstractNLPModel, state::NLPAtX)
  if has_bounds(pb)
    return vcat(
      min.(max.(state.mu, 0.0), -state.x + pb.meta.uvar),
      min.(max.(-state.mu, 0.0), state.x - pb.meta.lvar),
    )
  else
    return zeros(0)
  end
end

function _sign_multipliers_nonlin(pb::AbstractNLPModel, state::NLPAtX)
  if pb.meta.ncon == 0
    return zeros(0)
  else
    return vcat(
      min.(max.(state.lambda, 0.0), -state.cx + pb.meta.ucon),
      min.(max.(-state.lambda, 0.0), state.cx - pb.meta.lcon),
    )
  end
end

function _feasibility(pb::AbstractNLPModel, state::NLPAtX)
  if pb.meta.ncon == 0
    return vcat(max.(state.x - pb.meta.uvar, 0.0), max.(-state.x + pb.meta.lvar, 0.0))
  else
    return vcat(
      max.(state.cx - pb.meta.ucon, 0.0),
      max.(-state.cx + pb.meta.lcon, 0.0),
      max.(state.x - pb.meta.uvar, 0.0),
      max.(-state.x + pb.meta.lvar, 0.0),
    )
  end
end

"""
KKT: verifies the KKT conditions

`KKT( :: AbstractNLPModel, :: NLPAtX; pnorm :: Real = Inf, kwargs...)`

Note: `state.gx` is mandatory + if bounds `state.mu` + if constraints `state.cx`, `state.Jx`, `state.lambda`.

See also `unconstrained_check`, `unconstrained2nd_check`, `optim_check_bounded`
"""
function KKT(pb::AbstractNLPModel, state::NLPAtX; pnorm::Real = Inf, kwargs...)

  #Check the gradient of the Lagrangian
  gLagx = _grad_lagrangian(pb, state)
  #Check the complementarity condition for the bounds
  dual_res_bounds = _sign_multipliers_bounds(pb, state)
  #Check the complementarity condition for the constraints
  res_nonlin = _sign_multipliers_nonlin(pb, state)
  #Check the feasibility
  feas = _feasibility(pb, state)

  res = vcat(gLagx, feas, dual_res_bounds, res_nonlin)

  return norm(res, pnorm)
end


text:
using SafeTestsets

@safetestset "cholesky tests" begin
    include("cholesky_tests.jl")
end

@safetestset "display tests" begin
    include("display_tests.jl")
end

@safetestset "eval tests" begin
    include("eval_tests.jl")
end

@safetestset "find_path tests" begin
    include("find_path_tests.jl")
end

@safetestset "make_adjacency_matrix tests" begin
    include("make_adjacency_matrix_tests.jl")
end

@safetestset "make_cue_matrix tests" begin
    include("make_cue_matrix_tests.jl")
end

@safetestset "make_semantic_matrix tests" begin
    include("make_semantic_matrix_tests.jl")
end

@safetestset "make_yt_matrix tests" begin
    include("make_yt_matrix_tests.jl")
end

@safetestset "output_matrix tests" begin
    include("output_tests.jl")
end

@safetestset "preprocess tests" begin
    include("preprocess_tests.jl")
end

@safetestset "test_combo tests" begin
    include("test_combo_tests.jl")
end

@safetestset "wh tests" begin
    include("wh_tests.jl")
end


text:
using Test
using LinearAlgebraicRepresentation
using DataStructures
using PyCall
using Plasm

include("./views.jl")
include("./graphic_text.jl")


text:
using SCEDCAutoCorr
using Test

@testset "SCEDCAutoCorr.jl" begin
    # Write your tests here.
end


text:
# InvFourierLayerState
immutable InvFourierLayerState <: AbstractScatteredLayerState
    blobs::Vector{Mocha.Blob}
    layer::InvFourierLayer
end

function InvFourierLayerState(
        backend::Mocha.CPUBackend,
        layer::InvFourierLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    blobs = Vector{Mocha.Blob}(length(inputs))
    for idblob in eachindex(inputs)
        innodes = inputs[idblob].nodes
        outnodes = DataStructures.SortedDict{Path,InvComplexFourierNode,
            Base.Order.ForwardOrdering}()
        for path in keys(innodes)
            outnodes[path] = InvComplexFourierNode(
                innodes[path],
                findin(innodes[path].ranges, layer.pathkeys),
                layer.flags,
                layer.timelimit)
        end
        blobs[idblob] = ScatteredBlob(outnodes)
    end
    # TODO: build diffs
    return InvFourierLayerState(blobs, layer)
end

function Mocha.setup(
        backend::Mocha.Backend,
        layer::InvFourierLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    return InvFourierLayerState(backend, layer, inputs, diffs)
end


text:
n = 10 # rand(1:10)
m = 9  # rand(1:10)
@test matrixdepot("kahan", n) == matrixdepot("kahan", n, 1.2, 25.)
@test matrixdepot("kahan", m ,n, 4, 5) == matrixdepot("kahan", Float64, m, n, 4, 5)

@test matrixdepot("kahan", n, 0.5*pi, 1) ≈ Matrix(1.0I, n,n)
println("'kahan' passed test...")


text:
# Integers
export DistInt, add_bits, max_bits

struct DistInt <: Dist{Int}
    mgr
    # first index is least significant bit
    # most significant bits that are always false are trimmed
    bits::Vector{DistBool}
end

function DistInt(mgr, i::Int)
    @assert i >= 0
    # num_b = num_bits(i)
    num_b = ndigits(i, base = 2)
    bits = Vector{DistBool}(undef, num_b)
    for bit_idx = 1:num_b
        b::Bool = i & 1
        @inbounds bits[bit_idx] = DistBool(mgr, b)
        i = i >> 1
    end
    DistInt(mgr, bits)
end

function DistInt(bits::Vector)
    @assert !isempty(bits) "`DistInt` requires a bit width of at least 1"
    DistInt(bits[1].mgr, bits)
end

function infer(d::DistInt)
    mb = max_bits(d)
    ans = Vector(undef, 2^mb)
    non_zero_index = 0
    for i=0:2^mb - 1
        a = infer(prob_equals(d, i))
        if !(a ≈ 0)
            non_zero_index = i+1
        end
        ans[i + 1] = a
    end
    ans[1:non_zero_index]
end

max_bits(i::DistInt) =
    length(i.bits)

# @inline flip(mgr, ::Type{DistInt}, bits::Int) =
#     DistInt(mgr, [flip(mgr) for i = 1:bits])

function prob_equals(x::DistInt, y::DistInt)
    shared = min(max_bits(x), max_bits(y))
    eq = true
    for i=max_bits(x):-1:shared+1
        eq &= !(x.bits[i])
        !issat(eq) && return eq
    end
    for i=max_bits(y):-1:shared+1
        eq &= !(y.bits[i])
        !issat(eq) && return eq
    end
    for i=shared:-1:1
        # println(i)
        eq &= prob_equals(x.bits[i], y.bits[i])
        !issat(eq) && return eq
    end
    eq
end

prob_equals(x::DistInt, y::Int) =
    prob_equals(x, DistInt(x.mgr, y))

prob_equals(x::Int, y::DistInt) =
    prob_equals(y, x)

function ifelse(cond::DistBool, then::DistInt, elze::DistInt)
    if isvalid(cond)
        then
    elseif !issat(cond)
        elze
    else
        last_sat_bit = 1
        mbthen, mbelze = max_bits(then), max_bits(elze)
        mb = max(mbthen, mbelze)
        z = Vector{DistBool}(undef, mb)
        for i=1:mb
            z[i] = if i > mbthen
                !cond & elze.bits[i]
            elseif i > mbelze
                 cond & then.bits[i]
            else
                ifelse(cond, then.bits[i], elze.bits[i])
            end
            if issat(z[i])
                last_sat_bit = i
            end
        end
        DistInt(cond.mgr, z[1:last_sat_bit])
    end
end

#TODO: Look for optimizations here
function Base.:>(x::DistInt, y::DistInt)
    mx, my = max_bits(x), max_bits(y)
    shared = min(mx, my)
    eq = false
    prop = true
    for i=mx:-1:my+1
        prop &= !x.bits[i]
    end
    eq |= !prop
    for i = my:-1:mx+1
        prop &= !y.bits[i]
    end
    for i=shared:-1:1
        eq |= ifelse(prop, ifelse(x.bits[i] & !y.bits[i], true, false), false)
        prop &= prob_equals(x.bits[i], y.bits[i])
    end
    eq
end

Base.:>(x::DistInt, y::Int) = x > DistInt(x.mgr, y)

Base.:>(x::Int, y::DistInt) = DistInt(y.mgr, x) > y

# No canonical bitwidth
function Base.:+(p1::DistInt, p2::DistInt)
    last_sat_bit = 0
    mb1, mb2 = max_bits(p1), max_bits(p2)
    if (mb1 > mb2)
        t1, t2 = p1, p2
        mab, mib = mb1, mb2
    else
        t1, t2 = p2, p1
        mab, mib = mb2, mb1
    end

    z = Vector{DistBool}(undef, mab)
    carry = false
    for i = 1:mib
        z[i] = (!t1.bits[i] & t2.bits[i]) | (t1.bits[i] & !t2.bits[i])
        z[i] = (!z[i] & carry) | (z[i] & !carry) #combine in one line or maybe add XOR too
        carry = (t1.bits[i] & t2.bits[i]) | (carry & (t1.bits[i] | t2.bits[i]))
        # if issat(z[i])
        #     last_sat_bit = i
        # end
    end
    for i = mib+1:mab
        z[i] = (!t1.bits[i] & carry) | (t1.bits[i] & !carry)
        carry = t1.bits[i] & carry
        # if issat(z[i])
        #     last_sat_bit = i
        # end
    end
    DistInt(t1.mgr, z), carry
end

Base.:+(p1::DistInt, p2::Int) =
    p1 + DistInt(p1.mgr, p2)

Base.:+(p1::Int, p2::DistInt) =
    DistInt(p2.mgr, p1) + p2

function Base.:+(p1::Tuple{DistInt, DistBool}, p2::DistInt)
    a = p1[1] + p2
    a[1], a[2] | p1[2]
end

function Base.:+(p1::DistInt, p2::Tuple{DistInt, DistBool})
    a = p1 + p2[1]
    a[1], a[2] | p2[2]
end

function Base.:+(p2::Tuple{DistInt, DistBool}, p1::Tuple{DistInt, DistBool})
    a = p1[1] + p2[1]
    a[1], a[2] | p1[2] | p2[2]
end

function Base.:-(t1::DistInt, t2::DistInt)
    last_sat_bit = 0
    mb1, mb2 = max_bits(t1), max_bits(t2)
    mab = max(mb1, mb2)
    mib = min(mb1, mb2)

    z = Vector{DistBool}(undef, mab)
    borrow = false
    for i = 1:mib
        z[i] = ifelse(!borrow, (!t1.bits[i] & t2.bits[i]) | (t1.bits[i] & !t2.bits[i]), prob_equals(t1.bits[i], t2.bits[i]))
        borrow = ifelse(!borrow, !t1.bits[i] & t2.bits[i], !t1.bits[i] | t2.bits[i])
    end
    for i = mib+1:mb1
        z[i] = (!t1.bits[i] & borrow) | (t1.bits[i] & !borrow)
        borrow = !t1.bits[i] & borrow
    end

    for i = mib+1:mb2
        z[i] = (!t2.bits[i] & borrow) | (t2.bits[i] & !borrow)
        borrow = t2.bits[i] | borrow
    end

    for i = 1:mab
        z[i] = ifelse(borrow, !z[i], z[i])
    end

    ans = ifelse(borrow, (DistInt(z) + 1)[1], DistInt(z))
    ans, borrow
end

Base.:-(p1::DistInt, p2::Int) =
    p1 - DistInt(p1.mgr, p2)

Base.:-(p1::Int, p2::DistInt) =
    DistInt(p2.mgr, p1) - p2

function Base.:*(p1::DistInt, p2::DistInt)
    mbp1 = max_bits(p1)
    mbp2 = max_bits(p2)
    if (mbp1 > mbp2)
        t1, t2 = p1, p2
        mb1, mb2 = mbp1, mbp2
    else
        t1, t2 = p2, p1
        mb1, mb2 = mbp2, mbp1
    end

    P = DistInt(p1.mgr, 0)
    P = add_bits(P, mb1 - 1)
    shifted_bits = t1.bits
    carry = false
    for i=1:mb2
        if (i != 1)
            carry |= ifelse(t2.bits[i], shifted_bits[mb1], false)
            println(mb1)
            println(shifted_bits)
            shifted_bits = vcat(DistBool(p1.mgr, false), shifted_bits[1:mb1 - 1])
        end
        added = ifelse(t2.bits[i], DistInt(p1.mgr, shifted_bits), DistInt(p1.mgr, 0))
        res = (P + added)
        P = res[1]
        carry |= res[2]
    end
    P, carry
end

Base.:*(p1::DistInt, p2::Int) =
    p1 * DistInt(p1.mgr, p2)

Base.:*(p1::Int, p2::DistInt) =
    DistInt(p2.mgr, p1) * p2

function Base.:*(p1::Tuple{DistInt, DistBool}, p2::DistInt)
    a = p1[1] * p2
    a[1], a[2] | p1[2]
end

function Base.:*(p1::DistInt, p2::Tuple{DistInt, DistBool})
    a = p1 * p2[1]
    a[1], a[2] | p2[2]
end

function Base.:*(p2::Tuple{DistInt, DistBool}, p1::Tuple{DistInt, DistBool})
    a = p1[1] * p2[1]
    a[1], a[2] | p1[2] | p2[2]
end

function Base.:/(p1::DistInt, p2::DistInt) #p1/p2
    mb1 = max_bits(p1)
    mb2 = max_bits(p2)
    p1_bits = Vector(undef, 0)
    ans = Vector(undef, 0)

    is_zero = DistBool(p1.mgr, true)
    for i=1:mb2
        is_zero &= !p2.bits[i]
    end

    for i = mb1:-1:1
        p1_proxy = DistInt(p1.mgr, p1_bits)
        p1_bits = vcat(p1.bits[i:i], ifelse(p2 > p1_proxy, p1_proxy, (p1_proxy - p2)[1]).bits)
        ans = vcat(ifelse(p2 > p1_proxy, DistBool(p1.mgr, false), DistBool(p1.mgr, true)), ans)
    end
    p1_proxy = DistInt(p1.mgr, p1_bits)
    ans = vcat(ifelse(p2 > p1_proxy, DistBool(p1.mgr, false), DistBool(p1.mgr, true)), ans)
    ifelse(is_zero, p2, DistInt(p1.mgr, ans)), is_zero
end

function Base.:%(p1::DistInt, p2::DistInt) #p1%p2
    mb1 = max_bits(p1)
    mb2 = max_bits(p2)
    p1_bits = Vector(undef, 0)

    is_zero = DistBool(p1.mgr, true)
    for i=1:mb2
        is_zero &= !p2.bits[i]
    end

    for i = mb1:-1:1
        p1_proxy = DistInt(p1.mgr, p1_bits)
        p1_bits = vcat(p1.bits[i], ifelse(p2 > p1_proxy, p1_proxy, (p1_proxy - p2)[1]).bits)
    end
    p1_proxy = DistInt(p1.mgr, p1_bits)
    p1_bits = ifelse(p2 > p1_proxy, p1_proxy, (p1_proxy - p2)[1]).bits
    ifelse(is_zero, p2, DistInt(p1.mgr, p1_bits)), is_zero
    # DistInt(p1.mgr, p1_bits), DistBool(p1.mgr, false)
    # end
end

bools(i::DistInt) = i.bits

#No safety checks yet
function add_bits(p::DistInt, w::Int)
    mb = max_bits(p)
    ext = Vector(undef, w)
    for i= w:-1:1
        ext[i] = DistBool(p.mgr, false)
    end
    DistInt(p.mgr, vcat(p.bits, ext))
end


text:
using DataFrames: DataFrames
using PyCall, PyPlot

pandas = ArviZ.pandas

@testset "utils" begin
    @testset "replacemissing" begin
        @test isnan(ArviZ.replacemissing(missing))
        x = 1.0
        @test ArviZ.replacemissing(x) === x
        x = [1.0, 2.0, 3.0]
        @test ArviZ.replacemissing(x) === x
        x = [1.0, missing, 2.0]
        x2 = ArviZ.replacemissing(x)
        @test eltype(x2) <: Real
        @test x2[1] == 1.0
        @test isnan(x2[2])
        @test x2[3] == 2.0
        x = Any[1.0, missing, 2.0]
        x2 = ArviZ.replacemissing(x)
        @test eltype(x2) <: Real
        @test x2[1] == 1.0
        @test isnan(x2[2])
        @test x2[3] == 2.0

        x = [[1.0, missing, 2.0]]
        x2 = ArviZ.replacemissing([[1.0, missing, 2.0]])
        @test x2 isa Vector
        @test x2[1] isa Vector
        @test x2[1][1] == 1.0
        @test isnan(x2[1][2])
        @test x2[1][3] == 2.0
    end

    @testset "frompytype" begin
        x = 1.0
        @test ArviZ.frompytype(x) === x
        x2 = PyObject(x)
        @test ArviZ.frompytype(x2) == x
        @test ArviZ.frompytype([x2]) == [x]
        @test ArviZ.frompytype(Any[x2]) == [x]
        @test eltype(ArviZ.frompytype(Any[x2])) <: Real
        @test ArviZ.frompytype([[x2]]) == [[x]]
    end

    @testset "topandas" begin
        @testset "DataFrames.DataFrame -> pandas.DataFrame" begin
            colnames = [:a, :b, :c]
            index = ["d", "e"]
            rowvals = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
            df = DataFrames.DataFrame([
                :i => ["d", "e"],
                :a => [1.0, 4.0],
                :b => [2.0, 5.0],
                :c => [3.0, 6.0],
            ])
            pdf = ArviZ.topandas(Val(:DataFrame), df; index_name = :i)
            @test pyisinstance(pdf, pandas.DataFrame)
            pdf_exp = pandas.DataFrame(rowvals; columns = colnames, index = index)
            @test py"($(pdf) == $(pdf_exp)).all().all()"
        end

        @testset "DataFrames.DataFrame -> pandas.Series" begin
            df2 = DataFrames.DataFrame([:a => [1.0], :b => [2.0], :c => [3.0]])
            ps = ArviZ.topandas(Val(:Series), df2)
            @test pyisinstance(ps, pandas.Series)
            ps_exp = pandas.Series([1.0, 2.0, 3.0], [:a, :b, :c])
            @test py"($(ps) == $(ps_exp)).all()"
        end

        @testset "DataFrames.DataFrame -> ELPDData" begin
            idata = load_arviz_data("centered_eight")
            df = loo(idata)
            elpd_data = ArviZ.arviz.loo(idata)
            @test all(df == ArviZ.todataframes(elpd_data))
            ps = ArviZ.topandas(Val(:ELPDData), df)
            @test pyisinstance(ps, ArviZ.arviz.stats.ELPDData)
            @test py"($(ps) == $(elpd_data)).all()"
        end
    end

    @testset "todataframes" begin
        @testset "pandas.DataFrame -> DataFrames.DataFrame" begin
            colnames = [:a, :b, :c]
            index = ["d", "e"]
            rowvals = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
            pdf = pandas.DataFrame(rowvals; columns = colnames, index = index)
            df = ArviZ.todataframes(pdf; index_name = :i)
            @test df isa DataFrames.DataFrame
            @test df == DataFrames.DataFrame([
                :i => ["d", "e"],
                :a => [1.0, 4.0],
                :b => [2.0, 5.0],
                :c => [3.0, 6.0],
            ])
            @test df == ArviZ.todataframes(pdf; index_name = :i)
        end

        @testset "pandas.Series -> DataFrames.DataFrame" begin
            ps = pandas.Series([1.0, 2.0, 3.0], [:a, :b, :c])
            df2 = ArviZ.todataframes(ps)
            @test df2 isa DataFrames.DataFrame
            @test df2 == DataFrames.DataFrame([:a => [1.0], :b => [2.0], :c => [3.0]])
            @test df2 == ArviZ.todataframes(ps)
        end
    end

    @testset "styles" begin
        @test ArviZ.styles() isa AbstractArray{String}
        @test "arviz-darkgrid" ∈ ArviZ.styles()
        @test ArviZ.styles() == ArviZ.arviz.style.available
    end

    @testset "use_style" begin
        ArviZ.use_style("arviz-darkgrid")
        ArviZ.use_style("default")
    end
end


text:
struct Spherical <: AbstractPredictorCorrectorAlgorithm{true} end

"""

    predict(alg::Affine, H, J_H, ∂H∂t, x, t, Δt)

Spherical tangent predictor following the formulation of Chen[^1]. Denote by ``\hat{H}`` the lift of the homotopy ``H``
as a projective mapping.
Euler's method yields
```math
\begin{bmatrix}
D_x\hat{H}^(x,t) \\
x^{H}
\end{bmatrix}
\dot{x} =
\begin{bmatrix}
-D_t\hat{H}^(x,t) \\
0
\end{bmatrix}
```
and we can then obtain the spherical projection:
```math
\cos(\norm{\dot{x}}_2 Δt) x + \sin(\norm{\dot{x}}_2 Δt) \frac{\dot{x}}{\norm{\dot{x}}_2}
```

[^1]: Tianran Chen and Tien-Yien Li.
    “Spherical projective path tracking for homotopy continuation methods”.
    Communications in Information and Systems 12(3):195-220 (2012)
"""
function predict(alg::Spherical, H::AbstractHomotopy, J_H, ∂H∂t, x, t, Δt)
    dot_x = \([J_H(x,t); x'], [-∂H∂t(x,t); 0])
    norm_dot_x = norm(dot_x)
    # TODO: conditining? i.e. check singular values
    cos(-norm_dot_x*Δt) .* x .+ (sin(-norm_dot_x*Δt) / norm_dot_x) .* dot_x
end

function correct!(
    u::Vector{T},
    alg::Spherical,
    H::AbstractHomotopy,
    J_H::F,
    x::Vector{T},
    t,
    tol::Float64,
    max_iterations::Int
) where {T,F}
    N = length(x)
    res = zeros(x)
    dx_xt = zeros(T, N, N)
    Δx = zeros(x)
    u .= x

    for k = 1:max_iterations
        res[1:N-1] = -evaluate(H, u, t)
        # println("newton iteration: $(k), res: $(norm(res))")
        if norm(res) < tol
            return true
        end

        dx_xt[1:(N-1), :] .= J_H(u, t)
        dx_xt[N, :] = u'

        Δx .= \(dx_xt, res)
        norm_Δx = norm(Δx)
        u .= cos(norm_Δx) .* u .+ (sin(norm_Δx) / norm_Δx) .* Δx
    end

    norm(res) < tol
end


text:
# tests for addHybridizationUpdate! for the 5 taxon network
# we have in the ipad summary
# Claudia September 2014
##########################################################

warn("to run tests_5taxon.jl, you need to set updateGammaz to return true always. this is because we have bad triangles in here")
# types in "types.jl"
include("../src/types.jl")

# needed modules:
using Base.Collections # for updateInCycle with priority queue

# test functions
include("test_functions_5taxon.jl")

warn("BUG IN CASE C: sometimes it shows errors, but if Julia is closed and reopened, no more error shown")

tests = ["F","G","H","J","I"];
wrong = String[];

for t in tests
    include("../src/functions.jl")
    include("tree_example.jl");
    tp = string("add_hybrid_case","$(t).jl");
    println("running $(tp)-----");
    try
        include(tp)
    catch
        println("error in $(tp)");
        push!(wrong,t);
    end
end

if(!isempty(wrong))
    for t in wrong
        include("../src/functions.jl")
        include("tree_example.jl");
        tp = string("add_hybrid_case","$(t).jl");
        println("running $(tp)");
        include(tp)
    end
else
    println("----------NO ERRORS!----------");
end


text:
begin # Create (Deterministic) Limit Order Generator
    using Base.Iterators: cycle, take, zip, flatten
    orderid_iter = Base.Iterators.countfrom(1)
    sign_iter = cycle([1,-1,-1,1,1,-1])
    side_iter = ( s>0 ? :ASK : :BID for s in sign_iter )
    spread_iter = cycle([3 2 3 2 2 2 3 2 3 4 2 2 3 2 3 2 3 3 2 2 3 2 5 2 2 2 2 2 4 2 3 6 5 6 3 2 3 5 4]*1e-2)
    price_iter = ( Float32(100.0 + sgn*δ) for (δ,sgn) in zip(spread_iter,sign_iter) )
    size_iter = cycle([2,5,3,4,10,15,1,6,13,11,4,1,5])
    # zip them all together
    lmt_order_info_iter = zip(orderid_iter,price_iter,size_iter,side_iter)
end

begin # Create (Deterministic) Market Order Generator
    mkt_size_iter = cycle([10,20,30,15,25,5,7])
    mkt_side_iter = cycle([:BID,:ASK,:ASK,:BID,:ASK,:BID])
    mkt_order_info_iter = zip(mkt_size_iter,mkt_side_iter)
end

@testset "Submit and Cancel 1" begin # Add and delete all orders, verify book is empty, verify account tracking
    ob = LimitOrderBook.OrderBook() #Initialize empty book
    order_info_lst = take(lmt_order_info_iter,50000)
    # Add a bunch of orders
    for (orderid, price, size, side) in order_info_lst
        LimitOrderBook.submit_limit_order!(ob,orderid,price,size,side,acct_id=10101)
    end
    @test length(ob.acct_map[10101]) == 50000 # Check account order tracking
    # Cancel them all
    for (orderid, price, size, side) in order_info_lst
        LimitOrderBook.cancel_limit_order!(ob,orderid,price,side,acct_id=10101)
    end
    # Check emptiness
    @test isempty(ob.bid_orders)
    @test isempty(ob.ask_orders)
    @test isempty(ob.acct_map[10101])
end

@testset "MO Liquidity Wipe" begin # Wipe out book completely, try MOs on empty book
    ob = LimitOrderBook.OrderBook() #Initialize empty book
    # Add a bunch of orders
    for (orderid, price, size, side) in Base.Iterators.take( lmt_order_info_iter, 50 )
        LimitOrderBook.submit_limit_order!(ob,orderid,price,size,:BID)
    end
    mo_matches, mo_flag = LimitOrderBook.submit_market_order!(ob,:BID,100000)

    # Tests
    @test length( mo_matches ) == 50
    @test mo_flag == :INCOMPLETE
    @test isempty(LimitOrderBook.submit_market_order!(ob,:BID,10000)[1] )
    @test isempty(ob.bid_orders)
    @test isempty(ob.ask_orders)
end

@testset "Order match exact - bid" begin # Test correctness in order matching system / Stat calculation (:BID)
    ob = LimitOrderBook.OrderBook() #Initialize empty book
    # record order book info before
    order_lst_tmp = Base.Iterators.take( Base.Iterators.filter( x-> x[4]==:BID, lmt_order_info_iter), 7 ) |> collect

    # Add a bunch of orders
    for (orderid, price, size, side) in order_lst_tmp
        LimitOrderBook.submit_limit_order!(ob,orderid,price,size,:BID)
    end

    # record information from before
    expected_bid_volm_before = sum( x[3] for x in order_lst_tmp )
    expected_bid_n_orders_before =  length(order_lst_tmp)

    # execute MO
    mo_matches, mo_flag = LimitOrderBook.submit_market_order!(ob,:BID,30)
    mo_match_sizes = [o.size for o in mo_matches]

    # record what is expected to be seen
    expected_bid_volm_after = expected_bid_volm_before - 30
    expected_bid_n_orders_after = expected_bid_n_orders_before - 5
    expected_best_bid_after = Float32(99.97)

    # record what is expected of MO result
    expected_mo_match_size = [5,15,6,1,2,1]
    expected_mo_flag = :COMPLETE

    # Compute realized values
    book_info_after = LimitOrderBook.book_depth_info(ob,max_depth=1000)
    realized_bid_volm_after  = sum(book_info_after[:BID][:volume])
    realized_bid_n_orders_after = sum(book_info_after[:BID][:orders])
    realized_best_bid_after = first(book_info_after[:BID][:price])

    # Check all expected vs realized values
    @test realized_bid_volm_after == expected_bid_volm_after
    @test realized_bid_n_orders_after == expected_bid_n_orders_after
    @test realized_best_bid_after == expected_best_bid_after
    @test mo_match_sizes == expected_mo_match_size
    @test mo_flag == expected_mo_flag
end


text:
# Parser functionality for folders adhering to the
# ETSF Nanoquanta file format, for details see http://www.etsf.eu/fileformats/

using NCDatasets
using JSON

struct EtsfFolder
    folder
    gsr
    den
    eig
    pspmap
end

"""
Initialize a EtsfFolder from the path to the folder which contains
the data in the ETSF Nanoquanta format.
"""
function EtsfFolder(folder::AbstractString)
    if !isfile(joinpath(folder, "out_GSR.nc"))
        error("Did not find file $folder/out_GSR.nc")
    end
    gsr = Dataset(joinpath(folder, "out_GSR.nc"))

    den = nothing
    eig = nothing
    pspmap = Dict{Int, Any}()

    if isfile(joinpath(folder, "out_DEN.nc"))
        den = Dataset(joinpath(folder, "out_DEN.nc"))
    end
    if isfile(joinpath(folder, "out_EIG.nc"))
        eig = Dataset(joinpath(folder, "out_EIG.nc"))
    end
    if isfile(joinpath(folder, "pspmap.json"))
        dict = JSON.parsefile(joinpath(folder, "pspmap.json"))
        pspmap = Dict{Int, Any}(parse(Int, k) => dict[k] for k in keys(dict))
    end

    EtsfFolder(folder, gsr, den, eig, pspmap)
end

"""
Load a DFTK-compatible lattice object from the ETSF folder
"""
load_lattice(T, folder::EtsfFolder) = Mat3{T}(folder.gsr["primitive_vectors"][:])
load_lattice(folder; kwargs...) = load_lattice(Float64, folder; kwargs...)

"""
Load a DFTK-compatible atoms object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_atoms(T, folder::EtsfFolder)
    atoms = Dict{ElementPsp, Vector{Vec3{T}}}()
    n_species = length(folder.gsr["atomic_numbers"])
    for ispec in 1:n_species
        atnum = Int(folder.gsr["atomic_numbers"][ispec])
        spec = ElementPsp(atnum, psp=load_psp(folder.pspmap[atnum]))

        mask_species = findall(isequal(ispec), folder.gsr["atom_species"][:])
        positions = folder.gsr["reduced_atom_positions"][:, mask_species]
        atoms[spec] = [Vec3{T}(positions[:, m]) for m in mask_species]
    end
    collect(pairs(atoms))
end
load_atoms(folder; kwargs...) = load_atoms(Float64, folder; kwargs...)

"""
Load a DFTK-compatible model object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_model(T, folder::EtsfFolder)
    # Parse functional ...
    functional = Vector{Symbol}()
    ixc = Int(folder.gsr["ixc"][:])
    if ixc == 0
        # No XC functional
    elseif ixc == 1 || ixc == -20
        push!(functional, :lda_xc_teter93)
    elseif ixc == -1007
        push!(functional, :lda_x, :lda_c_vwn)
    elseif ixc == -101130 || ixc == 11
        push!(functional, :gga_x_pbe, :gga_c_pbe)
    else
        error("Unknown Functional ixc==$ixc")
    end

    # Parse smearing-scheme ...
    smearing = strip(join(folder.gsr["smearing_scheme"][:]))
    Tsmear = 0
    smearing_function = nothing
    if smearing == "none"
        # No smearing
    elseif smearing == "Fermi-Dirac"
        smearing_function = Smearing.FermiDirac()
    elseif smearing == "Methfessel and Paxton"
        smearing_function = Smearing.MethfesselPaxton1()
    elseif smearing == "gaussian"
        smearing_function = Smearing.Gaussian()
    else
        error("Unknown Smearing scheme: $smearing")
    end
    smearing_function !== nothing && (Tsmear = folder.gsr["smearing_width"][:])

    # Build model
    lattice = load_lattice(T, folder)
    atoms = load_atoms(T, folder)
    model_DFT(Array{T}(lattice), atoms, functional;
              smearing=smearing_function, temperature=Tsmear)
end
load_model(folder; kwargs...) = load_model(Float64, folder; kwargs...)

"""
Load a DFTK-compatible basis object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_basis(T, folder::EtsfFolder)
    model = load_model(T, folder)
    atoms = load_atoms(T, folder)

    Ecut = folder.gsr["kinetic_energy_cutoff"][:]

    n_kpoints = size(folder.gsr["reduced_coordinates_of_kpoints"], 2)
    kcoords = Vector{Vec3{T}}(undef, n_kpoints)
    for ik in 1:n_kpoints
        kcoords[ik] = Vec3{T}(folder.gsr["reduced_coordinates_of_kpoints"][:, ik])
    end

    kgrid_size = Vector{Int}(folder.gsr["monkhorst_pack_folding"])
    kcoords_new, ksymops = bzmesh_ir_wedge(kgrid_size, model.lattice, atoms)
    @assert kcoords_new ≈ kcoords

    PlaneWaveBasis(model, Ecut, kcoords, ksymops)
end
load_basis(folder; kwargs...) = load_basis(Float64, folder; kwargs...)

"""
Load a DFTK-compatible density object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_density(T, folder::EtsfFolder)
    # So far this function is untested
    basis = load_basis(T, folder)

    n_comp = size(folder.den["density"], 5)
    n_real_cpx = size(folder.den["density"], 1)
    @assert n_comp == 1
    @assert n_real_cpx == 1

    ρ_real = Array{Complex{T}}(folder.den["density"][1, :, :, :, 1])
    @assert basis.fft_size == size(ρ_real)

    from_real(basis, ρ_real)
end
load_density(folder; kwargs...) = load_density(Float64, folder; kwargs...)


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ d7b19f60-6348-11ec-11b8-abe5c0a6cee7
md"""
# Day 1: Sonar Sweep

**Reference:** [Day 1 (Advent of Code 2021)](https://adventofcode.com/2021/day/1)

**Source:** [ljk233/AdventOfCode2021/days/01/d1.jl (GitHub)](https://github.com/ljk233/AdventOfCode2021/blob/main/days/01/d1.jl)
"""

# ╔═╡ aeba97ec-fd96-4ef4-99b0-726665015217
md"""
## Part 1

**How many measurements are larger than the previous measurement?**
"""

# ╔═╡ 7f03ad5e-bf73-4b16-958b-592ba4b6d923
"""
	main1(f: String) -> Integer
"""
function main1(f)
	nums = parse.(Int, readlines(f))
	return sum(diff(nums) .≥ 1)
end

# ╔═╡ 8cc55754-6815-4c08-8538-b6f4c74d5c28
md"""
**Test #1.**
Expected result = 7.
"""

# ╔═╡ d0d5ecd6-fa56-4fca-a293-fb0a682e4a7b
main1("test.in")

# ╔═╡ 63a48804-1a55-48c6-832f-74f52ac1d583
md"""
**Solution.**
Confirmed result = 1466.
"""

# ╔═╡ 6c58e558-21fe-4221-b324-1efd90c7911a
main1("data.in")

# ╔═╡ 0e681b8a-82bc-425c-a374-e36d360559d2
md"""
## Part 2

Consider sums of a three-measurement sliding window.

**How many sums are larger than the previous sum?**
"""

# ╔═╡ 774f7d4b-9d47-41eb-a6d7-e9d0fc579fbb
"""
	main2(f: String) -> Integer
"""
function main2(f)
	nums = parse.(Int, readlines(f))
	agg = Int[]
	for i ∈ 1:length(nums)-2
    	push!(agg, sum(nums[i:i+2]))
	end
	return sum(diff(agg) .≥ 1)
end

# ╔═╡ 2e90f5d0-e5ee-4abb-932b-5638afb3a146
md"""
**Test #2.**
Expected result = 5.
"""

# ╔═╡ b5db6a00-eb1d-4cf0-bfad-d04b0d7d9968
main2("test.in")

# ╔═╡ 4d9b89d5-7f30-4031-a619-57fdb5936dc5
md"""
**Solution.**
Confirmed result = 1491.
"""

# ╔═╡ 17e05631-6ee0-49f0-b273-c151812220af
main2("data.in")

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─d7b19f60-6348-11ec-11b8-abe5c0a6cee7
# ╟─aeba97ec-fd96-4ef4-99b0-726665015217
# ╠═7f03ad5e-bf73-4b16-958b-592ba4b6d923
# ╟─8cc55754-6815-4c08-8538-b6f4c74d5c28
# ╠═d0d5ecd6-fa56-4fca-a293-fb0a682e4a7b
# ╟─63a48804-1a55-48c6-832f-74f52ac1d583
# ╠═6c58e558-21fe-4221-b324-1efd90c7911a
# ╟─0e681b8a-82bc-425c-a374-e36d360559d2
# ╠═774f7d4b-9d47-41eb-a6d7-e9d0fc579fbb
# ╟─2e90f5d0-e5ee-4abb-932b-5638afb3a146
# ╠═b5db6a00-eb1d-4cf0-bfad-d04b0d7d9968
# ╟─4d9b89d5-7f30-4031-a619-57fdb5936dc5
# ╠═17e05631-6ee0-49f0-b273-c151812220af
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Algorithms to find connected components

###########################################################
#
#   Connected components of undirected graph
#
##########################################################

function connected_components{V}(graph::AbstractGraph{V})
    @graph_requires graph vertex_list vertex_map adjacency_list

    !is_directed(graph) || error("graph must be undirected.")

    cmap = zeros(Int, num_vertices(graph))
    components = Array(Vector{V}, 0)

    for v in vertices(graph)
        if cmap[vertex_index(v, graph)] == 0
            visitor = VertexListVisitor{V}(0)
            traverse_graph(graph, BreadthFirst(), v, visitor, colormap=cmap)
            push!(components, visitor.vertices)
        end
    end

    components
end

###########################################################
#
#   Connected components of directed graph
#
##########################################################

#
# Recursively compute the strongly connected components of a directed graph.
#
# Adapted from
# http://en.wikipedia.org/wiki/Tarjan%27s_strongly_connected_components_algorithm
# on April 29, 2013.
#
function strongly_connected_components_recursive{V}(graph::AbstractGraph{V})
    @graph_requires graph vertex_list vertex_map adjacency_list

    is_directed(graph) || error("graph must be directed.")

    preorder_idx = 1
    stack        = Array(V, 0)
    lowlinks     = zeros(Int, num_vertices(graph))
    indices      = zeros(Int, num_vertices(graph))
    components   = Array(Vector{V}, 0)

    function strongconnect(v)
        v_idx = vertex_index(v, graph)
        indices[v_idx]  = preorder_idx
        lowlinks[v_idx] = preorder_idx
        preorder_idx   += 1
        push!(stack, v)

        # consider successors of v
        for w in out_neighbors(v, graph)
            w_idx = vertex_index(w, graph)
            if indices[v_idx(w)] == 0 # w is un-visited
                strongconnect(w)
                lowlinks[v_idx] = min(lowlinks[v_idx], lowlinks[w_idx])
            elseif (w in stack)
                lowlinks[v_idx] = min(lowlinks[v_idx], indices[w_idx])
            end
        end

        # if v is a root node, pop the stack and generate an SCC
        if lowlinks[v_idx] == indices[v_idx]
            component = Array(typeof(v), 0)
            while !isempty(stack)
                w = pop!(stack)
                push!(component, w)
                if w == v
                    break
                end
            end
            push!(components, component)
        end
    end

    for v in vertices(graph)
        if indices[vertex_index(v, graph)] == 0
            strongconnect(v)
        end
    end
    components
end

type TarjanVisitor{G<:AbstractGraph,V} <: AbstractGraphVisitor
    graph::G
    stack::Vector{V}
    lowlink::Vector{Int}
    index::Vector{Int}
    components::Vector{Vector{V}}
end

TarjanVisitor{V}(graph::AbstractGraph{V}) = TarjanVisitor{typeof(graph),V}(graph,
        V[], Int[], zeros(Int, num_vertices(graph)), Vector{V}[])

function discover_vertex!(vis::TarjanVisitor, v)
    iv = vertex_index(v, vis.graph)
    vis.index[iv] = length(vis.stack) + 1
    push!(vis.lowlink, length(vis.stack) + 1)
    push!(vis.stack, v)
    return true
end

function examine_neighbor!(vis::TarjanVisitor, v, w, w_color::Int, e_color::Int)
    if w_color == 1 # 1 means added seen, but not explored
        while vis.index[vertex_index(w, vis.graph)] < vis.lowlink[end]
            pop!(vis.lowlink)
        end
    end
end

function close_vertex!(vis::TarjanVisitor, v)
    iv = vertex_index(v, vis.graph)
    if vis.index[iv] == vis.lowlink[end]
        component = vis.stack[vis.index[iv]:end]
        splice!(vis.stack, vis.index[iv]:length(vis.stack))
        pop!(vis.lowlink)
        push!(vis.components, component)
    end
end

#
# Computes the strongly connected components of a directed graph.
#
# Adapated from
# http://code.activestate.com/recipes/578507-strongly-connected-components-of-a-directed-graph/
# on April 30, 2013.
#
function strongly_connected_components{V}(graph::AbstractGraph{V})
    @graph_requires graph vertex_list vertex_map adjacency_list

    cmap = zeros(Int, num_vertices(graph))
    components = Array(Vector{V}, 0)

    for v in vertices(graph)
        if cmap[vertex_index(v, graph)] == 0 # 0 means not visited yet
            visitor = TarjanVisitor(graph)
            traverse_graph(graph, DepthFirst(), v, visitor, vertexcolormap=cmap)
            for component in visitor.components
                push!(components, component)
            end
        end
    end
    components
end


text:
using MaxwellFDFD
using Test
using Statistics: mean
using AbbreviatedTypes

Base.isapprox(a::Tuple, b::Tuple; kws...) = all(p -> isapprox(p...; kws...), zip(a,b))

# @testset "MaxwellFDFD" begin

include("source.jl")

# end  # @testset "MaxwellFDFD"