text:
# shock capturing using any the SBP-SAT discretization of the second derivative
# term

"""
  Applies [`SBPParabolic`](@ref) shock capturing.
"""
function calcShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             capture::SBPParabolicSC{Tsol, Tres},
                             shockmesh::ShockedElements) where {Tsol, Tres}

  computeGradW(mesh, sbp, eqn, opts, capture, shockmesh,
               capture.entropy_vars, capture.diffusion)

  computeVolumeTerm(mesh, sbp, eqn, opts, capture, shockmesh)

  computeFaceTerm(mesh, sbp, eqn, opts, capture, shockmesh, capture.diffusion,
                  capture.penalty)

  #println("after face term, residual norm = ", calcNorm(eqn, eqn.res))

  if shockmesh.isNeumann
    computeNeumannBoundaryTerm(mesh, sbp, eqn, opts, capture, shockmesh)
  else
    computeDirichletBoundaryTerm(mesh, sbp, eqn, opts, capture, shockmesh)
  end

  computeSharedFaceTerm(mesh, sbp, eqn, opts, capture, shockmesh,
                              capture.diffusion, capture.penalty)

  return nothing
end

"""
  Computes:

  [ grad_x q  =  [ lambda_xx lambda_xy  [ Dx * w
    grad_y q]      lambda_yx lambda_yy]   Dy * w]

  and stores it in `capture.grad_w`.  Note that lambda = 0 for elements that
  do not have shocks in them, so the `grad_q` is set to zero there.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * capture: [`SBPParabolicSC`](@ref)
   * shockmesh: the `ShockedElements`
   * entropy_vars: an [`AbstractVariables`](@ref) object
   * diffusion: an [`AbstractDiffusion`](@ref)
"""
function computeGradW(mesh, sbp, eqn, opts, capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements,
                      entropy_vars::AbstractVariables,
                      diffusion::AbstractDiffusion,
                     ) where {Tsol, Tres}

  wxi = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.dim)
  grad_w = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.dim)

  # do local elements
  @simd for i in shockmesh.local_els
    i_full = shockmesh.elnums_all[i]
    @simd for j=1:mesh.numNodesPerElement
      q_j = ro_sview(eqn.q, :, j, i_full)
      w_j = sview(capture.w_el, :, j, i)

      # convert to entropy variables
      convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
    end

    # apply D operator
    q_i = ro_sview(eqn.q, :, :, i_full)
    w_i = ro_sview(capture.w_el, :, :, i)
    coords = ro_sview(mesh.coords, :, :, i_full)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i_full)
    jac_i = ro_sview(mesh.jac, :, i_full)
    fill!(grad_w, 0)
    applyDx(sbp, w_i, dxidx_i, jac_i, wxi, grad_w)

    # apply diffusion tensor
    lambda_gradq_i = sview(capture.grad_w, :, :, :, i)

    applyDiffusionTensor(diffusion, sbp, eqn.params, q_i, w_i, coords, dxidx_i,
                         jac_i, i, grad_w, lambda_gradq_i)
  end

  # the diffusion is zero in the neighboring elements, so convert to entropy
  # but zero out grad_w
  @simd for i in shockmesh.neighbor_els
    i_full = shockmesh.elnums_all[i]
    @simd for j=1:mesh.numNodesPerElement
      q_j = ro_sview(eqn.q, :, j, i_full)
      w_j = sview(capture.w_el, :, j, i)
      convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
    end

    gradw_i = sview(capture.grad_w, :, :, :, i)
    fill!(gradw_i, 0)
  end

  # do shared elements
  for peer=1:shockmesh.npeers
    peer_full = shockmesh.peer_indices[peer]
    data = eqn.shared_data[peer_full]
    metrics = mesh.remote_metrics[peer_full]

    for i in shockmesh.shared_els[peer]
      i_full = getSharedElementIndex(shockmesh, mesh, peer, i)

      # compute entropy variables
      for j=1:mesh.numNodesPerElement
        q_j = ro_sview(data.q_recv, :, j, i_full)
        w_j = sview(capture.w_el, :, j, i)
        convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
      end

      lambda_gradq_i = sview(capture.grad_w, :, :, :, i)

      # some of these have ee = 0, but we don't have the sensor so we can't
      # check
      q_i = ro_sview(data.q_recv, :, :, i_full)
      w_i = ro_sview(capture.w_el, :, :, i)
      coords_i = ro_sview(metrics.coords, :, :, i_full)
      dxidx_i = ro_sview(metrics.dxidx, :, :, :, i_full)
      jac_i = ro_sview(metrics.jac, :, i_full)
      fill!(grad_w, 0)

      applyDx(sbp, w_i, dxidx_i, jac_i, wxi, grad_w)
      applyDiffusionTensor(diffusion, sbp, eqn.params, q_i,  w_i, coords_i,
                           dxidx_i, jac_i, i, grad_w, lambda_gradq_i)

    end  # end i
  end  # end peer

  return nothing
end

"""
  Computes the volume terms, using the intermediate variable calcualted by
  [`computeGradW`](@ref).  The residual is updated with

  [Qx Qy] * [Lambda] * [Dx * w
                        Dy * w]

  Note that this used Qx, not Qx^T, so the term has not been integrated by
  parts.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * capture: [`SBPParbaolicSC`](@ref)
   * shockmesh
"""
@noinline function computeVolumeTerm(mesh, sbp, eqn, opts,
                           capture::SBPParabolicSC{Tsol, Tres},
                           shockmesh::ShockedElements) where {Tsol, Tres}

  # computeGradW computes Lambda * D * w, so all that remains to do is
  # compute Qx * grad_q_x
  # Note that this term is not entropy stable by itself, because Qx was
  # not replaced by -Qx^T + Ex.  The entire discretization should be
  # entropy-stable however.
  work = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.dim)
  op = SummationByParts.Subtract()
  for i=1:shockmesh.numShock
    i_full = shockmesh.elnums_all[i]

    gradq_i = ro_sview(capture.grad_w, :, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i_full)
    res_i = sview(eqn.res, :, :, i_full)
    applyQx(sbp, gradq_i, dxidx_i, work, res_i)

    #applyQxTransposed(sbp, gradq_i, dxidx_i, work, res_i, op)
  end

  return nothing
end

"""
  Computes the interior face term for any [`SBPParabolic`](@ref) shock
  capturing scheme.

  **Inputs**

   * mesh
   * sbp
   * eqn: `eqn.res` is updated with the result
   * opts
   * shockmesh
   * diffusion: the [`AbstractDiffusion`](@ref) to use (must be consistent
                with the one passed to [`computeGradW`](@ref).
   * penalty: [`AbstractDiffusionPenalty`](@ref) specifying which scheme to use

"""
@noinline function computeFaceTerm(mesh, sbp, eqn, opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements, diffusion::AbstractDiffusion,
                      penalty::AbstractDiffusionPenalty) where {Tsol, Tres}

  #TODO: avoid recalculating the diffusion tensor 6 time: add an API for
  #      caching it (called in getFaceData) and then apply it.  Add a clear
  #      cache function at end so future calls to applyDiffusionTensor work
  delta_w = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  theta = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  t1L = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t1R = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t2L = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t2R = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  op = SummationByParts.Subtract()

  for i=1:shockmesh.numInterfaces
    iface_red = shockmesh.ifaces[i].iface
    iface_idx = shockmesh.ifaces[i].idx_orig
    elnumL = shockmesh.elnums_all[iface_red.elementL]
    elnumR = shockmesh.elnums_all[iface_red.elementR]

    # get data needed for next steps
    qL = ro_sview(eqn.q, :, :, elnumL)
    qR = ro_sview(eqn.q, :, :, elnumR)
    wL = ro_sview(capture.w_el, :, :, iface_red.elementL)
    wR = ro_sview(capture.w_el, :, :, iface_red.elementR)
    gradwL = ro_sview(capture.grad_w, :, :, :, iface_red.elementL)
    gradwR = ro_sview(capture.grad_w, :, :, :, iface_red.elementR)

    coordsL = ro_sview(mesh.coords, :, :, elnumL)
    coordsR = ro_sview(mesh.coords, :, :, elnumR)
    nrm_face = ro_sview(mesh.nrm_face, :, :, iface_idx)
    dxidxL = ro_sview(mesh.dxidx, :, :, :, elnumL)
    dxidxR = ro_sview(mesh.dxidx, :, :, :, elnumR)
    jacL = ro_sview(mesh.jac, :, elnumL)
    jacR = ro_sview(mesh.jac, :, elnumR)
    alphas = ro_sview(capture.alpha, :, i)
    resL = sview(eqn.res, :, :, elnumL)
    resR = sview(eqn.res, :, :, elnumR)

    # compute delta w tilde and theta_bar = Dgk w_k + Dgn w_n
    getFaceVariables(capture, mesh.sbpface, iface_red, wL, wR, gradwL, gradwR,
                     nrm_face, delta_w, theta)

    # apply the penalty coefficient matrix
    applyPenalty(penalty, sbp, eqn.params, mesh.sbpface, diffusion,
                 iface_red, delta_w, theta,
                 qL, qR, wL, wR, coordsL, coordsR, nrm_face, alphas, dxidxL,
                 dxidxR, jacL, jacR, t1L, t1R, t2L, t2R)

    # apply Rgk^T, Rgn^T, Dgk^T, Dgn^T
    # need to apply R^T * t1, not R^T * B * t1, so
    # interiorFaceIntegrate won't work.  Use the reverse mode instead
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.numDofPerNode
        t1L[k, j] = -t1L[k, j]  # the SAT has a minus sign in front of it
        t1R[k, j] = -t1R[k, j]
      end
    end
    interiorFaceInterpolate_rev!(mesh.sbpface, iface_red, resL, resR, t1L, t1R)

    # apply Dgk^T and Dgn^T
    # TODO: there is an allocation here caused by the number of arguments
    #       to the function, however the allocation is context dependent.
    #       Pulling this function out of PDESolver removes the allocation.
    #       The allocation can be avoided by packing dxidxL, jacL, and resL
    #       into a tuple (and similarly for the corresponding R arguments),
    #       however the allocation does not appear to cause a significant
    #       performance problem.
    applyDgkTranspose(capture, sbp, eqn.params, mesh.sbpface, iface_red,
                      diffusion, t2L, t2R,
                      qL, qR, wL, wR, coordsL, coordsR, nrm_face, dxidxL,
                      dxidxR, jacL, jacR, resL, resR, op)

  end  # end loop i

  return nothing
end

"""
  Does the same thing as `compuateFaceTerm`, but for the shared faces, updating
  the residual on the local element only
"""
@noinline function computeSharedFaceTerm(mesh, sbp, eqn, opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements, diffusion::AbstractDiffusion,
                      penalty::AbstractDiffusionPenalty) where {Tsol, Tres}

  delta_w = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  theta = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  t1L = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t1R = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t2L = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t2R = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  op = SummationByParts.Subtract()

  # don't care about resR for shared faces
  resR = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)
  for _peer=1:shockmesh.npeers
    # finish parallel communication for this set of faces
    peer = receiveAlpha(capture, shockmesh)

    peer_full = shockmesh.peer_indices[peer]
    metrics = mesh.remote_metrics[peer_full]
    data = eqn.shared_data[peer_full]

    for i=1:shockmesh.numSharedInterfaces[peer]
      iface_red = shockmesh.shared_interfaces[peer][i].iface
      iface_idx = shockmesh.shared_interfaces[peer][i].idx_orig
      elnumL = shockmesh.elnums_all[iface_red.elementL]
      elnumR = getSharedElementIndex(shockmesh, mesh, peer, iface_red.elementR)

      # get data needed for next steps
      qL = ro_sview(eqn.q, :, :, elnumL)
      qR = ro_sview(data.q_recv, :, :, elnumR)
      wL = ro_sview(capture.w_el, :, :, iface_red.elementL)
      wR = ro_sview(capture.w_el, :, :, iface_red.elementR)
      gradwL = ro_sview(capture.grad_w, :, :, :, iface_red.elementL)
      gradwR = ro_sview(capture.grad_w, :, :, :, iface_red.elementR)

      coordsL = ro_sview(mesh.coords   , :, :, elnumL)
      coordsR = ro_sview(metrics.coords, :, :, elnumR)
      nrm_face = ro_sview(mesh.nrm_sharedface[peer_full], :, :, iface_idx)
      dxidxL = ro_sview(mesh.dxidx,    :, :, :, elnumL)
      dxidxR = ro_sview(metrics.dxidx, :, :, :, elnumR)
      jacL = ro_sview(mesh.jac,    :, elnumL)
      jacR = ro_sview(metrics.jac, :, elnumR)
      alphas = ro_sview(capture.alpha_parallel[peer], :, i)
      resL = sview(eqn.res, :, :, elnumL)

      # compute delta w tilde and theta_bar = Dgk w_k + Dgn w_n
      getFaceVariables(capture, mesh.sbpface, iface_red, wL, wR, gradwL, gradwR,
                       nrm_face, delta_w, theta)

      # apply the penalty coefficient matrix
      applyPenalty(penalty, sbp, eqn.params, mesh.sbpface, diffusion,
                   iface_red, delta_w,
                   theta, qL, qR,  wL, wR, coordsL, coordsR, nrm_face, alphas,
                   dxidxL, dxidxR, jacL, jacR, t1L, t1R, t2L, t2R)

      # apply Rgk^T, Rgn^T, Dgk^T, Dgn^T
      # need to apply R^T * t1, not R^T * B * t1, so
      # interiorFaceIntegrate won't work.  Use the reverse mode instead
      for j=1:mesh.numNodesPerFace
        for k=1:mesh.numDofPerNode
          t1L[k, j] = -t1L[k, j]  # the SAT has a minus sign in front of it
          t1R[k, j] = -t1R[k, j]
        end
      end
      #TODO: this could be a boundaryIntegrate now that resR doesn't matter
      interiorFaceInterpolate_rev!(mesh.sbpface, iface_red, resL, resR,
                                   t1L, t1R)

      # apply Dgk^T and Dgn^T
      applyDgkTranspose(capture, sbp, eqn.params, mesh.sbpface, iface_red,
                        diffusion,
                        t2L, qL, wL, coordsL, nrm_face, dxidxL, jacL, resL, op)
    end  # end i
  end  # end peer

  return nothing
end

"""
  Computes a Neumann boundary condition Lambda * grad w = 0 for all faces
  on the boundary.  This term is required to make the entire shock capturing
  scheme entropy stable.

  This term does not depend on what penalty scheme is used, so it does not
  take an `AbstractPenalty` object.  It uses `capture.grad_w`, so it does
  not need an `AbstractDiffusion` object either.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * capture: an [`SBPParabolic`](@ref) object
   * shockmesh
"""
@noinline function computeNeumannBoundaryTerm(mesh::AbstractMesh{Tmsh}, sbp, eqn, opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements,
                      ) where {Tsol, Tres, Tmsh}

  # for shock capturing, apply the Neumann boundary condition
  # Lambda * grad_w = 0.  The resulting term is -R^T * B * Dgk * u

  #@assert mesh.coord_order == 1

  temp_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  temp2_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  nrm_i = zeros(Tmsh, mesh.dim, mesh.numNodesPerFace)
  op = SummationByParts.Subtract()
  for i=1:shockmesh.numBoundaryFaces
    bndry_i = shockmesh.bndryfaces[i].bndry
    idx_orig = shockmesh.bndryfaces[i].idx_orig
    elnum_orig = shockmesh.elnums_all[bndry_i.element]

    if i < shockmesh.bndry_offsets[end]
      #nrm_i = ro_sview(mesh.nrm_bndry, :, :, idx_orig)
      for j=1:mesh.numNodesPerFace
        for d=1:mesh.dim
          nrm_i[d, j] = mesh.nrm_bndry[d, j, idx_orig]
        end
      end
    else
      error("this doesn't work yet")
      # what needs to happen is all boundaries of the shockmesh that are
      # not part of the boundary of the domain need to be included in
      # shockmesh.bndryfaces.  This is because the volume term is
      # not integrated by parts, so we need to subtract off the boundary
      # term
      # Also, need to get the correct normal vector (nbrperm, factor of -1)
      #nrm_i = ro_sview(mesh.nrm_face, :, :, idx_orig)
      #TODO: nbrperm
      fac = shockmesh.bndryfaces[i].fac
      for j=1:mesh.numNodesPerFace
        for d=1:mesh.dim
          nrm_i[d, j] = fac*mesh.nrm_face[d, j, idx_orig]
        end
      end
    end
    res_i = sview(eqn.res, :, :, elnum_orig)

    # grad_w already has Lambda * [Dx; Dy]*u, now apply R and N
    fill!(temp2_face, 0)
    for d1=1:mesh.dim
      gradw_d = ro_sview(capture.grad_w, :, :, d1, bndry_i.element)
      fill!(temp_face, 0)
      boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, gradw_d, temp_face)

      @simd for j=1:mesh.numNodesPerFace
        @simd for k=1:mesh.numDofPerNode
          temp2_face[k, j] += nrm_i[d1, j]*temp_face[k, j]
        end
      end
    end  # end d1

    # apply R^T B
    boundaryFaceIntegrate!(mesh.sbpface, bndry_i.face, temp2_face, res_i, op)
  end  # end i

  return nothing
end

"""
  Computes a Dirichlet BC that is consistent with the inviscid one
"""
@noinline function computeDirichletBoundaryTerm(mesh, sbp, eqn, opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements,
                      ) where {Tsol, Tres}

  for i=1:mesh.numBC
    bc_range = (shockmesh.bndry_offsets[i]:(shockmesh.bndry_offsets[i+1]-1))
    bc_func = mesh.bndry_funcs[i]

    calcBoundaryFlux(mesh, sbp, eqn, opts, shockmesh, capture,
                     capture.penalty, capture.diffusion,
                     capture.entropy_vars, bc_func, bc_range)
  end

  return nothing
end

@noinline function calcBoundaryFlux(mesh::AbstractMesh, sbp, eqn::EulerData, opts,
                      shockmesh::ShockedElements,
                      capture::SBPParabolicSC{Tsol, Tres},
                      penalty::AbstractDiffusionPenalty,
                      diffusion::AbstractDiffusion,
                      entropy_vars::AbstractVariables, bc_func::BCType,
                      bc_range::AbstractVector,
                      ) where {Tsol, Tres}

  delta_w = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  res1 = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  op = SummationByParts.Subtract()

  for i in bc_range
    bndry_i = shockmesh.bndryfaces[i].bndry
    idx_orig = shockmesh.bndryfaces[i].idx_orig
    elnum_orig = shockmesh.elnums_all[bndry_i.element]

    w_i = ro_sview(capture.w_el, :, :, bndry_i.element)
    q_i = ro_sview(eqn.q, :, :, elnum_orig)
    coords_i = ro_sview(mesh.coords_bndry, :, :, idx_orig)
    nrm_i = ro_sview(mesh.nrm_bndry, :, :, idx_orig)
    alpha = capture.alpha_b[i]
    dxidxL = ro_sview(mesh.dxidx, :, :, :, elnum_orig)
    jacL = ro_sview(mesh.jac, :, elnum_orig)
    res_i = sview(eqn.res, :, :, elnum_orig)

    computeDirichletDelta(capture, eqn.params, mesh.sbpface, bndry_i,
                          bc_func, entropy_vars, w_i, q_i, coords_i,
                          nrm_i, delta_w)

    applyDirichletPenalty(penalty, sbp, eqn.params, mesh.sbpface, diffusion,
                          bndry_i, delta_w, q_i, w_i, coords_i, nrm_i, alpha,
                          dxidxL, jacL, res1)

#    println("after pnealty, res1 = ", res1)
    # apply R^T to T_D * delta_w
    scale!(res1, -1)  # SAT has a - sign in front of it

    boundaryFaceInterpolate_rev!(mesh.sbpface, bndry_i.face, res_i, res1)

    # apply B and then Dgk^T to delta_w
    for j=1:mesh.numNodesPerFace
      for i=1:mesh.numDofPerNode
        delta_w[i, j] *= -mesh.sbpface.wface[j]
      end
    end

    applyDgkTranspose(capture, sbp, eqn.params, mesh.sbpface, bndry_i,
                      diffusion, delta_w, q_i, w_i, coords_i, nrm_i, dxidxL,
                      jacL, res_i, op)

#    println("after Dgk^T, res = ", res_i)
  end

  return nothing
end

"""
  Compute delta_w = Rgk * wL - Rgn - wR and theta = Dgk * wL + Dgn * wR
  at the face.

  **Inputs**

   * capture: an [`SBPParabolic`](@ref) object
   * sbpface
   * iface_red: an `Interface` object from the `shockmesh`
   * wL: the entropy variables for the left element, `numDofPerNode` x
         `numNodesPerElement`
   * wR: the entropy variables for the right element, same size as `wL`
   * gradwL: Lambda * [Dx, Dy]*wL (as computed by [`computeGradW`](@ref),
             `numDofPerNode` x `numNodesPerElement` x `dim`
   * gradwR: same as `gradwR`, but for the right element
   * nrm_face: the (scaled) normal vectors for the face, `dim` x
               `numNodesPerFace`

  **Inputs/Outputs**

   * delta_w: `numDofPerNode` x `numNodesPerFace`, overwritten
   * theta: `numDofPerNode` x `numNodesPerFace`, overwritten
"""
function getFaceVariables(capture::SBPParabolicSC{Tsol, Tres},
                          sbpface::AbstractFace, iface_red::Interface,
                          wL::AbstractMatrix, wR::AbstractMatrix,
                          gradwL::Abstract3DArray, gradwR::Abstract3DArray,
                          nrm_face::AbstractMatrix,
                          delta_w::AbstractMatrix, theta::AbstractArray
                         ) where {Tsol, Tres}

  numDofPerNode, numNodesPerElement = size(wL)
  numNodesPerFace = size(delta_w, 2)
  dim = size(gradwL, 3)

  @unpack capture w_faceL w_faceR grad_faceL grad_faceR

  interiorFaceInterpolate!(sbpface, iface_red, wL, wR, w_faceL, w_faceR)

  @simd for j=1:numNodesPerFace
    @simd for k=1:numDofPerNode
      delta_w[k, j] = w_faceL[k, j] - w_faceR[k, j]
    end
  end

  fill!(theta, 0)
  for d=1:dim
    gradwL_d = ro_sview(gradwL, :, :, d)
    gradwR_d = ro_sview(gradwR, :, :, d)
    interiorFaceInterpolate!(sbpface, iface_red, gradwL_d, gradwR_d,
                             grad_faceL, grad_faceR)

    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        theta[k, j] += nrm_face[d, j]*(grad_faceL[k, j] - grad_faceR[k, j])
      end
    end
  end  # end d

  return nothing
end

"""
  Computes the difference between the solution at the face and the
  prescribed Dirichlet boundary condition.  Specifically, computes the
  difference between the entropy variables interpolated to the face
  and the Dirichlet state computed from the interpolated conservative
  variables (for consistency with the inviscid BC).

  This only works for boundary conditions that define `getDirichletState`.

  **Inputs**

   * capture
   * params
   * sbpface
   * bndry: `Boundary` object
   * bndry_node: `BoundaryNode` object
   * func: `BCType` for this boundary condition
   * entropy_vars: an [`AbstractVariables`](@ref)
   * wL: entropy variables for the element, `numDofPerNode` x
         `numNodesPerElement`
   * qL: conservative variables for the element, same size as `wL`
   * coords: `dim` x `numNodesPerFace` array containing the xyz coordinates
             of the face nodes
   * nrm_face: `dim` x `numNodesPerFace` containing the normal vector at the
               face nodes (can be scaled or unscaled, `getDirichletState`
               should not care).

  **Inputs/Outputs**

   * delta_w: array to be overwritten with the result. `numDofPerNode` x
              `numNodesPerFace`
"""
function computeDirichletDelta(capture::SBPParabolicSC{Tsol, Tres},
                               params::ParamType,
                               sbpface::AbstractFace, bndry::Boundary,
                               func::BCType, entropy_vars::AbstractVariables,
                               wL::AbstractMatrix, qL::AbstractMatrix,
                               coords::AbstractMatrix, nrm_face::AbstractMatrix,
                               delta_w::AbstractMatrix
                              ) where {Tsol, Tres}

  numDofPerNode, numNodesPerFace = size(delta_w)
  numNodesPerElement = size(wL, 2)

  #TODO: maybe interpolating the conservative variables and then converting
  #      would be more consistent, but Eq. 31 requires the same interpolation
  #      Rgk * w for both the interface and the boundary term.
  #      For the boundary state qg we have more leeway because it is zero
  #      for the energy stability analysis (-> entropy stability when u = w)
  wface = zeros(Tsol, numDofPerNode, numNodesPerFace)
  boundaryFaceInterpolate!(sbpface, bndry.face, wL, wface)

  qface = zeros(Tsol, numDofPerNode, numNodesPerFace)
  qgface = zeros(Tres, numDofPerNode)
  wgface = zeros(Tres, numDofPerNode)
  aux_vars = Tres[]

  # the inviscid boundary conditions interpolate q to the face and then
  # compute qg, so do that here as well for consistency
  boundaryFaceInterpolate!(sbpface, bndry.face, qL, qface)
  for i=1:numNodesPerFace
    q_i = sview(qface, :, i)
    coords_i = sview(coords, :, i)
    nrm_i = sview(nrm_face, :, i)

    getDirichletState(func, params, q_i, aux_vars, coords_i, nrm_i, qgface)

    convertToEntropy(entropy_vars, params, qgface, wgface)
#    fill!(wgface, 0)  #TODO: undo this
#    println("wgface = ", wgface)

    # compute the delta
    for j=1:numDofPerNode
      delta_w[j, i] = wface[j, i] - wgface[j]
    end
  end

  return nothing
end

"""
  This function applies Dgk^T * t2L and Dgn^T * t2R.

  **Inputs**

   * capture: an [`SBPParabolic`](@ref) object
   * sbp
   * sbpface
   * iface: an `Interface` from `shockmesh`
   * diffusion: an [`AbstractDiffusion`](@ref)
   * t2L: array to multiply Dgk^T with, `numDofPerNode` x `numNodesPerFace`
   * t2R: array to multiply Dgn^T with, `numDofPerNode` x `numNodesPerFace`
   * wL: entropy variables for left element, `numDofPerNode` x
         `numNodesPerElement`
   * wR: entropy variables for the right element, same size as `wL`
   * nrm_face: (scaled) normal vectors at the face, `dim` x `numNodesPerFace`
   * dxidxL: (scaled) mapping jacobian at left element, `dim` x `dim` x
             `numNodesPerElement`
   * dxidxR: (scaled) mapping jacobian at the right element, same size as
             `dxidxL`
   * jacL: mapping jacobian determinant for left element, `numNodesPerElement`
   * jacR: mapping jacobian determinant for the right element,
           `numNodesPerElement`
   * op: a `SummationByParts.UnaryFunctor`, determins if the output arrays are
         added to or subtracted from

  **Inputs/Outputs**

   * resL: array to update with Dgk^T * t2L, `numDofPerNode` x
           `numNodesPerElement`
   * resR: array to update with Dgn^T * t2R, same size as `resL`

"""
function applyDgkTranspose(capture::SBPParabolicSC{Tsol, Tres}, sbp,
                           params::ParamType,
                           sbpface, iface::Interface,
                           diffusion::AbstractDiffusion,
                           t2L::AbstractMatrix, t2R::AbstractMatrix,
                           qL::AbstractMatrix, qR::AbstractMatrix,
                           wL::AbstractMatrix, wR::AbstractMatrix,
                           coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                           nrm_face::AbstractMatrix,
                           dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                           jacL::AbstractVector, jacR::AbstractVector,
                           resL::AbstractMatrix, resR::AbstractMatrix,
                           op::SummationByParts.UnaryFunctor=SummationByParts.Add()) where {Tsol, Tres}

  dim, numNodesPerFace = size(nrm_face)
  numNodesPerElement = size(resL, 2)
  numDofPerNode = size(wL, 1)

  @unpack capture temp1L temp1R temp2L temp2R temp3L temp3R work
  fill!(temp2L, 0); fill!(temp2R, 0)

  # apply N and R^T
  @simd for d=1:dim
    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        temp1L[k, j] =  nrm_face[d, j]*t2L[k, j]
        temp1R[k, j] = -nrm_face[d, j]*t2R[k, j]
      end
    end

    tmp2L = sview(temp2L, :, :, d); tmp2R = sview(temp2R, :, :, d)
    interiorFaceInterpolate_rev!(sbpface, iface, tmp2L, tmp2R, temp1L, temp1R)
  end

  # multiply by D^T Lambda
  applyDiffusionTensor(diffusion, sbp, params, qL, wL, coordsL, dxidxL, jacL,
                       iface.elementL, temp2L, temp3L)
  applyDiffusionTensor(diffusion, sbp, params, qR, wR,  coordsR, dxidxR, jacR,
                       iface.elementR, temp2R, temp3R)

  # saving temp3 to a mesh-wide array and then applying Dx^T would save
  # a lot of flops.
  applyDxTransposed(sbp, temp3L, dxidxL, jacL, work, resL, op)
  applyDxTransposed(sbp, temp3R, dxidxR, jacR, work, resR, op)

  return nothing
end

"""
  Method to apply Dgk^T * t2L only.
"""
function applyDgkTranspose(capture::SBPParabolicSC{Tsol, Tres}, sbp,
                           params::ParamType,
                           sbpface, bndry::Union{Interface, Boundary},
                           diffusion::AbstractDiffusion,
                           t2L::AbstractMatrix,
                           qL::AbstractMatrix,
                           wL::AbstractMatrix,
                           coordsL::AbstractMatrix,
                           nrm_face::AbstractMatrix,
                           dxidxL::Abstract3DArray,
                           jacL::AbstractVector,
                           resL::AbstractMatrix,
                           op::SummationByParts.UnaryFunctor=SummationByParts.Add()) where {Tsol, Tres}

  dim, numNodesPerFace = size(nrm_face)
  numNodesPerElement = size(resL, 2)
  numDofPerNode = size(wL, 1)

  @unpack capture temp1L temp2L temp3L work
  fill!(temp2L, 0)

  # apply N and R^T
  @simd for d=1:dim
    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        temp1L[k, j] =  nrm_face[d, j]*t2L[k, j]
      end
    end

    tmp2L = sview(temp2L, :, :, d);
    boundaryFaceInterpolate_rev!(sbpface, getFaceL(bndry), tmp2L, temp1L)
  end

  # multiply by D^T Lambda
  applyDiffusionTensor(diffusion, sbp, params, qL, wL, coordsL, dxidxL, jacL,
                       getElementL(bndry), temp2L, temp3L)

  applyDxTransposed(sbp, temp3L, dxidxL, jacL, work, resL, op)

  return nothing
end

#------------------------------------------------------------------------------
# helper functions

"""
  Computes the alpha parameter.  Starts parallel communication so alpha
  for remote elements can be computed later.

  **Inputs**

   * caputure: `capture.alpha` is overwritten.
   * mesh
   * shockmesh
"""
function computeAlpha(capture::SBPParabolicSC, mesh::AbstractMesh,
                      shockmesh::ShockedElements)

  # I think the alpha_gk parameter should be 0 for faces on the Neumann boundary
  # Count the number of faces each element has in the mesh to compute alpha_gk
  # such that it sums to 1.
  if size(capture.alpha, 2) < shockmesh.numInterfaces
    capture.alpha = zeros(2, shockmesh.numInterfaces)
  else
    fill!(capture.alpha, 0)
  end

  if !shockmesh.isNeumann
    if length(capture.alpha_b) < shockmesh.numBoundaryFaces
      capture.alpha_b = zeros(shockmesh.numBoundaryFaces)
    else
      fill!(capture.alpha_b, 0)
    end
  end

  oldsize = length(capture.alpha_parallel)
  resize!(capture.alpha_parallel, shockmesh.npeers)
  for peer=1:shockmesh.npeers
    len_i = shockmesh.numSharedInterfaces[peer]
    # don't try to get the size of alpha_parallel[peer] if it is undefined
    if peer > oldsize || size(capture.alpha_parallel[peer], 2) < len_i
      capture.alpha_parallel[peer] = zeros(2, len_i)
    else
      fill!(capture.alpha_parallel[peer], 0)
    end
  end

  # count number of faces of each element not on boundary
  el_counts = zeros(UInt8, shockmesh.numEl)
  for i=1:shockmesh.numInterfaces
    iface_i = shockmesh.ifaces[i].iface
    elnumL = shockmesh.elnums_all[iface_i.elementL]
    elnumR = shockmesh.elnums_all[iface_i.elementR]
    el_counts[iface_i.elementL] += 1
    el_counts[iface_i.elementR] += 1
  end

  if !shockmesh.isNeumann
    for i=1:shockmesh.numBoundaryFaces
      bndry_i = shockmesh.bndryfaces[i].bndry
      el_counts[bndry_i.element] += 1
    end
  end

  for peer=1:shockmesh.npeers
    for i=1:shockmesh.numSharedInterfaces[peer]
      iface_i = shockmesh.shared_interfaces[peer][i].iface
      el_counts[iface_i.elementL] += 1
      # only do local elements, because we can't see all faces of remote
      # elements
    end
  end

  # start parallel communications
  sendAlphas(capture.alpha_comm, shockmesh, el_counts)

  # compute alpha from face counts
  # for neighbor elements alpha doesn't sum to 1, but thats ok because
  # lambda = 0 there, so alpha multiplies zero.
  for i=1:shockmesh.numInterfaces
    iface_i = shockmesh.ifaces[i].iface
    capture.alpha[1, i] = 1/el_counts[iface_i.elementL]
    capture.alpha[2, i] = 1/el_counts[iface_i.elementR]
  end

  if !shockmesh.isNeumann
    for i=1:shockmesh.numBoundaryFaces
      bndry_i = shockmesh.bndryfaces[i].bndry
      capture.alpha_b[i] = 1/el_counts[bndry_i.element]
    end
  end

  for peer=1:shockmesh.npeers
    for i=1:shockmesh.numSharedInterfaces[peer]
      iface_i = shockmesh.shared_interfaces[peer][i].iface
      capture.alpha_parallel[peer][1, i] = 1/el_counts[iface_i.elementL]
    end
  end

  return nothing
end

"""
  Send the number of faces not on a Neumann boundary each element has to
  other processes.

  **Inputs**

   * comm: `AlphaComm` object
   * shockmesh: `ShockedElements`
   * el_counts: array containing the number of faces each local element has
                that are not on a Neumann boundary

"""
function sendAlphas(comm::AlphaComm, shockmesh::ShockedElements,
                    el_counts::AbstractVector{I}) where {I <: Integer}

  # extract the required values into the send buffers

  # allocateArrays checks that previous communications are complete

  # We only need to send number of faces per element, but that would require
  # uniquely identifying the elements on either side of the partition boundary.
  # Instead, send the same value for all faces of the element and let the
  # the receiver sort it out

  # post receives
  for peer=1:shockmesh.npeers
    Irecv!(comm, peer)
  end

  # send
  for peer=1:shockmesh.npeers
    println("Sending alphas to process ", comm.peer_parts_red[peer]); flush(STDOUT)
    # pack the buffer
    for j=1:shockmesh.numSharedInterfaces[peer]
      iface_i = shockmesh.shared_interfaces[peer][j].iface
      comm.send_bufs[peer][j] = el_counts[iface_i.elementL]
    end

    Isend(comm, peer)
  end

  return nothing
end

"""
  Receives the parallel communication started by `sendAlphas`.  Specifically,
  it receives the communication from one peer process and returns the index
  of that process.  This function should be called exactly the same number of
  times that there are peer process.

  On exit, `capture.alpha` will be populated with the correct alpha values
  for the element shared with the peer process.

  **Inputs**

   * capture: `SBPParabolicSC` shock capturing object
   * shockmesh: `ShockedElements`

  **Outputs**

   * idx: the index of the peer that communication was received from.
"""
function receiveAlpha(capture::SBPParabolicSC, shockmesh::ShockedElements)

  peer = Waitany!(capture.alpha_comm)
  comm = capture.alpha_comm

  # unpack buffer and calculate alpha
  for j=1:shockmesh.numSharedInterfaces[peer]
    iface_i = shockmesh.shared_interfaces[peer][j]
    capture.alpha_parallel[peer][2, j] = 1/comm.recv_bufs[peer][j]
    # the buffer was packed with el_counts for the parent element of each
    # face, so all thats required it so take 1/the value
  end

  return peer
end


text:
using KCenters
using Documenter

makedocs(;
    modules=[KCenters],
    authors="Eric S. Tellez",
    repo="https://github.com/sadit/KCenters.jl/blob/{commit}{path}#L{line}",
    sitename="KCenters.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sadit.github.io/KCenters.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Center selection" => "centerselection.md",
        "Clustering" => "clustering.md",
        "Density nets" => "dnet.md",
        "Epsilon nets" => "enet.md",
        "Stop criterions" => "criterions.md",
        "Partitioning function" => "utils.md"
    ],
)

deploydocs(;
    repo="github.com/sadit/KCenters.jl",
    devbranch="main",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#"]
)

ext: jl


text:
module TestUpdateOrder
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType

@testset "buildorder!(b a c)" begin
    pool = CategoricalPool(
        [
            "a",
            "b",
            "c"
        ],
        Dict(
            "a" => convert(DefaultRefType, 1),
            "b" => convert(DefaultRefType, 2),
            "c" => convert(DefaultRefType, 3),
        )
    )

    order = Vector{DefaultRefType}(undef, length(pool.index))

    CategoricalArrays.buildorder!(order, pool.invindex, ["b", "a", "c"])

    @test order[1] == convert(DefaultRefType, 2)
    @test order[2] == convert(DefaultRefType, 1)
    @test order[3] == convert(DefaultRefType, 3)
end

@testset "levels are built correctly" begin
    orig_index = [2, 5, 1, 3, 4]
    orig_levels = [1, 2, 3, 4, 5]
    pool = CategoricalPool(orig_index, orig_levels, true)
    @test orig_index == CategoricalArrays.index(pool)
    @test orig_levels == levels(pool)
    @test CategoricalArrays.index(pool) == levels(pool)[CategoricalArrays.order(pool)]
end

end


text:
export phi_dphi
"""
If f : Rⁿ → R then we have h(θ) = f(x + θ*d).
We can use φ(θ) = h(θ) - h(0) - τ₀*h'(0).
Inputs are a line model and a LSAtT structure.
This function returns  φ(θ) and φ'(θ).
"""
function phi_dphi(h :: LineModel,  state :: LSAtT;  τ₀ :: Float64 = 0.01)
    ht, dht = objgrad(h, state.x)
    update!(state, ht = ht, gt = dht)

    if state.x == 0.0
        φt = 0.0                      # known that φ(0) = 0.0
        dφt = (1.0 - τ₀) * state.g₀   # known that φ'(0) = (1.0 - τ₀) * h'(0)
    else
        φt = state.ht - state.h₀ - τ₀ * state.x* state.g₀
        dφt = state.gt - τ₀ * state.g₀
    end

    return φt, dφt
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Print function of the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Low-level function to print the table using the text backend.
function _pt_text(io, pinfo;
                  border_crayon::Crayon = Crayon(),
                  header_crayon::Union{Crayon,Vector{Crayon}} = Crayon(bold = true),
                  subheader_crayon::Union{Crayon,Vector{Crayon}} = Crayon(foreground = :dark_gray),
                  rownum_header_crayon::Crayon = Crayon(bold = true),
                  text_crayon::Crayon = Crayon(),
                  autowrap::Bool = false,
                  body_hlines::Vector{Int} = Int[],
                  body_hlines_format::Union{Nothing,NTuple{4,Char}} = nothing,
                  crop::Symbol = :both,
                  crop_subheader::Bool = false,
                  columns_width::Union{Integer,AbstractVector{Int}} = 0,
                  equal_columns_width::Bool = false,
                  highlighters::Union{Highlighter,Tuple} = (),
                  hlines::Union{Nothing,Symbol,AbstractVector} = nothing,
                  linebreaks::Bool = false,
                  maximum_columns_width::Union{Integer,AbstractVector{Int}} = 0,
                  minimum_columns_width::Union{Integer,AbstractVector{Int}} = 0,
                  overwrite::Bool = false,
                  noheader::Bool = false,
                  nosubheader::Bool = false,
                  row_name_crayon::Crayon = Crayon(bold = true),
                  row_name_header_crayon::Crayon = Crayon(bold = true),
                  row_number_alignment::Symbol = :r,
                  screen_size::Union{Nothing,Tuple{Int,Int}} = nothing,
                  show_row_number::Bool = false,
                  sortkeys::Bool = false,
                  tf::TextFormat = unicode,
                  title_autowrap::Bool = false,
                  title_crayon::Crayon = Crayon(bold = true),
                  title_same_width_as_table::Bool = false,
                  vlines::Union{Nothing,Symbol,AbstractVector} = nothing,
                  # Deprecated
                  formatter = nothing,
                  same_column_size = nothing)

    @unpack_PrintInfo pinfo

    # Input variables verification and initial setup
    # ==========================================================================

    # Let's create a `IOBuffer` to write everything and then transfer to `io`.
    io_has_color = get(io, :color, false)
    buf_io       = IOBuffer()
    buf          = IOContext(buf_io, :color => io_has_color)
    screen       = Screen(has_color = io_has_color)

    # If the user did not specified the screen size, then get the current
    # display size. However, if cropping is not desired, then just do nothing
    # since the size is initialized with -1.
    if crop != :none
        if screen_size == nothing
            # For files, the function `displaysize` returns the value of the
            # environments variables "LINES" and "COLUMNS". Hence, here we set
            # those to `-1`, so that we can use this information to avoid
            # limiting the output.
            withenv("LINES" => -1, "COLUMNS" => -1) do
                screen.size = displaysize(io)
            end
        else
            screen.size = screen_size
        end

        # If the user does not want to crop, then change the size to -1.
        if crop == :vertical
            screen.size = (screen.size[1],-1)
        elseif crop == :horizontal
            screen.size = (-1,screen.size[2])
        end
    end

    !noheader && num_cols != header_num_cols &&
    error("The header length must be equal to the number of columns.")

    # Additional processing necessary if the user wants to print the header.
    if !noheader
        # If the user do not want to print the sub-header but wants to print the
        # header, then just force the number of rows in header to be 1.
        if nosubheader
            # Now, `header` will be a view of the first line of the matrix that
            # has the header.
            header = @view header[1:header_num_rows:end]
            header_num_rows = 1
        end

        # Transform some keywords that are single elements to vectors.
        if typeof(header_crayon) == Crayon
            header_crayon = [header_crayon for i = 1:num_cols]
        else
            length(header_crayon) != num_cols &&
            error("The length of `header_crayon` must be the same as the number of columns.")
        end

        if typeof(subheader_crayon) == Crayon
            subheader_crayon = [subheader_crayon for i = 1:num_cols]
        else
            length(subheader_crayon) != num_cols &&
            error("The length of `subheader_crayon` must be the same as the number of columns.")
        end
    end

    # Make sure that `highlighters` is always a tuple.
    !(highlighters isa Tuple) && (highlighters = (highlighters,))

    # Make sure that `maximum_columns_width` is always a vector.
    typeof(maximum_columns_width) <: Integer &&
        (maximum_columns_width = ones(Int, num_cols)*maximum_columns_width)

    # Make sure that `minimum_columns_width` is always a vector.
    typeof(minimum_columns_width) <: Integer &&
        (minimum_columns_width = ones(Int, num_cols)*minimum_columns_width)

    # Check which columns must have fixed sizes.
    typeof(columns_width) <: Integer && (columns_width = ones(Int, num_cols)*columns_width)
    length(columns_width) != num_cols && error("The length of `columns_width` must be the same as the number of columns.")
    fixed_col_width = map(w->w > 0, columns_width)

    # Assign the minimum size to those columns that does not have a fixed size.
    for i = 1:num_cols
        !fixed_col_width[i] && (columns_width[i] = minimum_columns_width[i])
    end

    # Increase the number of printed columns if the user wants to show
    # additional columns.
    Δc = show_row_names + show_row_number
    num_printed_cols += Δc

    # If the user wants to horizontally crop the printing, then it is not
    # necessary to process all the lines. We will to process, at most, the
    # number of lines in the screen.
    if screen.size[1] > 0
        num_printed_rows = min(num_printed_rows, screen.size[1])
    end

    # If the user wants to vertically crop the printing, then it is not
    # necessary to process all the columns. However, we cannot know at this
    # stage the size of each column. Thus, this initial algorithm uses the fact
    # that each column printed will have at least 4 spaces.
    if screen.size[2] > 0
        num_printed_cols = min(num_printed_cols, ceil(Int, screen.size[2]/4))
    end

    # Create the string matrices that will be printed
    # ==========================================================================

    # Get the string which is printed when `print` is called in each element of
    # the matrix. Notice that we must create only the matrix with the printed
    # rows and columns.
    header_str = Matrix{String}(undef, header_num_rows, num_printed_cols)
    data_str   = Matrix{Vector{AbstractString}}(undef,
                                                num_printed_rows,
                                                num_printed_cols)
    num_lines_in_row = ones(Int, num_printed_rows)

    # The variable `columns_width` is the specification of the user for the
    # columns width. The variable `cols_width` contains the actual size of each
    # column. This is necessary because if the user asks for a width equal or
    # lower than 0 in a column, then the width will be automatically computed to
    # fit the longest field.
    #
    # By default, if the user wants to show the row number of the row names,
    # those columns will have the size computed automatically.
    cols_width = Vector{Int}(undef, num_printed_cols)

    for i = (Δc+1):num_printed_cols
        cols_width[i] = columns_width[id_cols[i-Δc]]
    end

    # This variable stores the predicted table width. If the user wants
    # horizontal cropping, then it can be use to avoid unnecessary processing of
    # columns that will not be displayed.
    pred_tab_width = 0

    # Auxiliary variables to adjust `alignment` based on the new columns added
    # to the table.
    alignment_add = Symbol[]

    # Row numbers
    # --------------------------------------------------------------------------

    @inbounds if show_row_number
        # Set the header of the row column.
        header_str[1,1]      = "Row"
        header_str[2:end,1] .= ""

        # Set the data of the row column.
        for i = 1:num_printed_rows
            data_str[i,1] = [string(id_rows[i])]
        end

        # The row number width depends on how many digits the total number of
        # rows # has and the length of the header "Row". Notice that if
        # `noheader` is set to `true`, then we should not take the word "Row"
        # into account.
        cols_width[1] = max(noheader ? 0 : 3, floor(Int, log10(num_rows)) + 1)

        # Add information about the row number column to the variables.
        push!(alignment_add, row_number_alignment)
    end

    # Row names
    # --------------------------------------------------------------------------

    @inbounds @views if show_row_names
        # Escape the row name column title.
        header_str[1,Δc]      = _str_escaped(row_name_column_title)
        header_str[2:end,Δc] .= ""

        # Convert the row names to string.
        max_size = 0
        for i = 1:num_printed_rows
            data_str[i,Δc] = [_str_escaped(sprint(print, row_names[i];
                                                  context = :compact => compact_printing))]

            len_i = textwidth(data_str[i,Δc][1])
            len_i > max_size && (max_size = len_i)
        end

        # Obtain the size of the row name column.
        cols_width[Δc] = max(length(header_str[1,Δc]), max_size)

        # Add information about the row name column to the variables.
        push!(alignment_add, row_name_alignment)
    end

    !isempty(alignment_add) && (alignment = vcat(alignment_add, alignment))

    # Table data
    # --------------------------------------------------------------------------

    @inbounds for i = (1+Δc):num_printed_cols
        # Index of the i-th printed column in `data`.
        ic = id_cols[i-Δc]

        if !noheader
            for j = 1:header_num_rows
                header_str[j,i] = _str_escaped(sprint(print, header[(ic-1)*header_num_rows + j];
                                                      context = :compact => compact_printing))

                # Compute the minimum column size to print this string.
                cell_width = textwidth(header_str[j,i])

                # If the user does not want a fixed column width, then we must
                # store the information to automatically compute the field size.
                if !fixed_col_width[ic]
                    # Check if we should consider the sub-header size when
                    # computing the column width.
                    if (j == 1) || (!crop_subheader)
                        # Check if we need to increase the columns size.
                        cols_width[i] < cell_width && (cols_width[i] = cell_width)

                        # Make sure that the maximum column width is respected.
                        if maximum_columns_width[ic] > 0 &&
                           maximum_columns_width[ic] < cols_width[i]
                            cols_width[i] = maximum_columns_width[ic]
                        end
                    end
                end
            end
        end

        for j = 1:num_printed_rows
            # Index of the j-th printed row in `data`.
            jr = id_rows[j]

            # Apply the formatters.
            data_ij = isassigned(data,jr,ic) ? data[jr,ic] : undef

            for f in formatters
                data_ij = f(data_ij, jr, ic)
            end

            # Handle `nothing`, `missing`, and `undef`.
            if ismissing(data_ij)
                data_str_ij = "missing"
            elseif data_ij == nothing
                data_str_ij = "nothing"
            elseif data_ij == undef
                data_str_ij = "#undef"
            elseif data_ij isa Markdown.MD
                r = repr(data_ij)
                len = max(0, min(length(r), 32)-1)
                data_str_ij = len < length(r) - 1 ? first(r, len)*"…" : first(r, len)
            else
                data_str_ij = sprint(print, data_ij;
                                     context = :compact => compact_printing)
            end

            if cell_first_line_only
                tokens = _str_line_breaks(data_str_ij,
                                          autowrap && fixed_col_width[ic],
                                          columns_width[ic])
                data_str[j,i] = [tokens[1]]
                cell_width    = textwidth(data_str[j,i][1])
            elseif linebreaks
                tokens = _str_line_breaks(data_str_ij,
                                          autowrap && fixed_col_width[ic],
                                          columns_width[ic])
                data_str[j,i] = tokens
                num_lines_ij  = length(tokens)

                # Check if we must update the number of lines in this row.
                num_lines_in_row[j] < num_lines_ij && (num_lines_in_row[j] = num_lines_ij)

                # Compute the maximum length to compute the column size.
                cell_width = maximum(textwidth.(tokens))
            else
                data_str_ij_esc = _str_escaped(data_str_ij)
                data_str[j,i]   = [data_str_ij_esc]
                cell_width      = textwidth(data_str_ij_esc)
            end

            # If the user does not want a fixed column width, then we must store
            # the information to automatically compute the field size.
            if !fixed_col_width[ic]
                # Check if we need to increase the columns size.
                cols_width[i] < cell_width && (cols_width[i] = cell_width)

                # Make sure that the maximum column width is respected.
                if maximum_columns_width[ic] > 0 && maximum_columns_width[ic] < cols_width[i]
                    cols_width[i] = maximum_columns_width[ic]
                end
            end
        end

        # If the user horizontal cropping, then check if we need to process
        # another column.
        #
        # TODO: Should we take into account the dividers?
        if screen.size[2] > 0
            pred_tab_width += cols_width[i]

            if pred_tab_width > screen.size[2]
                num_printed_cols = i
                break
            end
        end
    end

    # If the user wants all the columns with the same size, then select the
    # larger.
    equal_columns_width && (cols_width = [maximum(cols_width) for i = 1:num_printed_cols])

    # Compute where the horizontal and vertical lines must be drawn
    # --------------------------------------------------------------------------

    # Create the format of the horizontal lines.
    if body_hlines_format == nothing
        body_hlines_format = (tf.left_intersection, tf.middle_intersection,
                              tf.right_intersection, tf.row)
    end

    hlines == nothing && (hlines = tf.hlines)
    hlines = _process_hlines(hlines, body_hlines, num_printed_rows, noheader)

    # Process `vlines`.
    vlines == nothing && (vlines = tf.vlines)
    vlines = _process_vlines(vlines, num_printed_cols)

    #                           Print the table
    # ==========================================================================

    # Top table line
    # ==========================================================================

    0 ∈ hlines && _draw_line!(screen, buf, tf.up_left_corner,
                              tf.up_intersection, tf.up_right_corner, tf.row,
                              border_crayon, cols_width, vlines)

    # Header
    # ==========================================================================

    # Header and sub-header texts
    # --------------------------------------------------------------------------

    if !noheader
        @inbounds @views for i = 1:header_num_rows
            0 ∈ vlines && _p!(screen, buf, border_crayon, tf.column)

            for j = 1:num_printed_cols
                # Get the information about the alignment and the crayon.
                if j ≤ Δc
                    if show_row_number && (j == 1)
                        crayon_ij    = rownum_header_crayon
                        alignment_ij = alignment[1]
                    elseif show_row_names
                        crayon_ij    = row_name_header_crayon
                        alignment_ij = alignment[Δc]
                    end
                else
                    jc = id_cols[j-Δc]
                    crayon_ij    = (i == 1) ? header_crayon[jc] : subheader_crayon[jc]
                    alignment_ij = header_alignment[jc]

                    # Check for cell alignment override.
                    for f in header_cell_alignment
                        aux = f(header, i, jc)

                        if aux ∈ (:l, :c, :r, :L, :C, :R, :s, :S)
                            alignment_ij = aux
                            break
                        end
                    end

                    # If alignment is `:s`, then we must use the column
                    # alignment.
                    alignment_ij ∈ (:s,:S) && (alignment_ij = alignment[jc+Δc])
                end

                # Prepare the text to be printed.
                header_i_str = " " * _str_aligned(header_str[i,j],
                                                  alignment_ij,
                                                  cols_width[j]) * " "

                flp = j == num_printed_cols

                # Print the text.
                _p!(screen, buf, crayon_ij, header_i_str)

                # Check if we need to draw a vertical line here.
                _pc!(j ∈ vlines, screen, buf, border_crayon, tf.column, " " , flp)

                _eol(screen) && break
            end

            i != header_num_rows && _nl!(screen,buf)
        end

        _nl!(screen,buf)

        # Bottom header line
        #-----------------------------------------------------------------------

        1 ∈ hlines && _draw_line!(screen, buf, tf.left_intersection,
                                  tf.middle_intersection,
                                  tf.right_intersection, tf.row,
                                  border_crayon, cols_width, vlines)
    end

    # Data
    # ==========================================================================

    @inbounds @views for i = 1:num_printed_rows
        ir = id_rows[i]

        for l = 1:num_lines_in_row[i]
            0 ∈ vlines && _p!(screen, buf, border_crayon, tf.column)

            for j = 1:num_printed_cols
                # Get the information about the alignment and the crayon.
                if j ≤ Δc
                    if show_row_number && (j == 1)
                        crayon_ij    = text_crayon
                        alignment_ij = alignment[1]
                    elseif show_row_names
                            crayon_ij    = row_name_crayon
                        alignment_ij = alignment[Δc]
                    end
                else
                    jc = id_cols[j-Δc]
                    crayon_ij    = text_crayon
                    alignment_ij = alignment[jc+Δc]

                    # Check for highlighters.
                    for h in highlighters
                        if h.f(_getdata(data), ir, jc)
                            crayon_ij = h.fd(h, _getdata(data), ir, jc)
                            break
                        end
                    end

                    # Check for cell alignment override.
                    for f in cell_alignment
                        aux = f(_getdata(data), ir, jc)

                        if aux ∈ [:l, :c, :r, :L, :C, :R]
                            alignment_ij = aux
                            break
                        end
                    end
                end

                # Align the string to be printed.
                if length(data_str[i,j]) >= l
                    data_ij_str = " " * _str_aligned(data_str[i,j][l], alignment_ij, cols_width[j]) * " "
                else
                    data_ij_str = " " * _str_aligned("", alignment_ij, cols_width[j]) * " "
                end

                # Print.
                _p!(screen, buf, crayon_ij, data_ij_str)

                flp = j == num_printed_cols

                # Check if we need to draw a vertical line here.
                _pc!(j ∈ vlines, screen, buf, border_crayon, tf.column, " " , flp)

                _eol(screen) && break
            end

            _nl!(screen, buf)

        end

        # Check if we must draw a horizontal line here.
        i != num_printed_rows && (i+!noheader) in hlines &&
            _draw_line!(screen, buf, body_hlines_format..., border_crayon,
                        cols_width, vlines)

        # Here we must check if the vertical size of the screen has been
        # reached. Notice that we must add 4 to account for the command line,
        # the continuation line, the bottom table line, and the last blank line.
        if (screen.size[1] > 0) && (screen.row + 4 >= screen.size[1])
            _draw_continuation_row(screen, buf, tf, text_crayon, border_crayon,
                                   cols_width, vlines)
            break
        end
    end

    # Bottom table line
    # ==========================================================================

    (num_printed_rows + !noheader) ∈ hlines &&
        _draw_line!(screen, buf, tf.bottom_left_corner, tf.bottom_intersection,
                    tf.bottom_right_corner, tf.row, border_crayon, cols_width,
                    vlines)

    # Overwrite table
    # ==========================================================================

    # If `overwrite` is `true`, then delete the exact number of lines of the
    # table. This can be used to replace the table in the screen continuously.

    str_overwrite = overwrite ? "\e[1F\e[2K"^(screen.row - 1) : ""

    # Title
    # ==========================================================================

    if length(title) > 0
        title_width = title_same_width_as_table ? screen.max_col : screen.size[2]

        print(io, title_crayon)

        # If the title width is not higher than 0, then we should only print the
        # title.
        if title_width ≤ 0
            println(io,title)
        # Otherwise, we must check for the alignments.
        else
            title_tokens = _str_line_breaks(title, title_autowrap, title_width)
            for token in title_tokens
                println(io, _str_aligned(token, title_alignment, title_width))
            end
        end

        print(io, text_crayon)
    end

    # Print the buffer
    # ==========================================================================
    print(io, str_overwrite * String(take!(buf_io)))

    return nothing
end


text:
module MosekMathProgSolverInterface
import ..Mosek

# Known issues:
#  - SOCP and QP cannot be mixed, but this is not checked (an error from mosek will be produced, though)
#  - Adding a conic quadratic constraint will add an empty constraint to ensure that the number of values
#    in constraint solution is as expected. The actual constraint solution value is bogus.
#  - Adding rotated conic quadratic constraints will result in a constraint being added, but the constraint soloution
#    for this is pointless. Also, a variable is added, but this is filtered out in the results.
#  - Loading an SOCP problem file will cause some funky problems as information on extra variables etc. is lost.
#  - Dual information is currently useless.
#
#  - The concept of dual values is a bit shaky. Specifically; for a variable x there is a dual for the upper bound,
#    one for the lower bound and one for the conic "bound". The dual value reported will be (slx-sux+snx).

import MathProgBase

status(t::Mosek.Task) = status(t,Mosek.MSK_RES_OK)
function status(t::Mosek.Task, r::Mosek.Rescode)
    if  r == Mosek.MSK_RES_TRM_MAX_ITERATIONS ||
        r == Mosek.MSK_RES_TRM_MAX_NUM_SETBACKS ||
        r == Mosek.MSK_RES_TRM_MAX_TIME ||
        r == Mosek.MSK_RES_TRM_MIO_NEAR_ABS_GAP ||
        r == Mosek.MSK_RES_TRM_MIO_NEAR_REL_GAP ||
        r == Mosek.MSK_RES_TRM_MIO_NUM_BRANCHES ||
        r == Mosek.MSK_RES_TRM_MIO_NUM_RELAXS ||
        r == Mosek.MSK_RES_TRM_NUM_MAX_NUM_INT_SOLUTIONS ||
        r == Mosek.MSK_RES_TRM_OBJECTIVE_RANGE
        :UserLimit
    #elseif r == Mosek.MSK_RES_TRM_STALL
    #    # Yes, I know this is not a standard code, but it doesn't fit any other code
    #    :Stall
    elseif r == Mosek.MSK_RES_TRM_USER_CALLBACK
        :UserBreak
    elseif r.value >= 1000 && r.value < 10000 # awful hack! Means: MSK_RES_ERR_*. r < 1000 is MSK_RES_WRN_*, r >= 10000 means MSK_RES_TRM_*
        :Error
    elseif ! ( (Mosek.solutiondef(t,Mosek.MSK_SOL_ITG) && Mosek.getsolsta(t,Mosek.MSK_SOL_ITG) != Mosek.MSK_SOL_STA_UNKNOWN) ||
               (Mosek.solutiondef(t,Mosek.MSK_SOL_BAS) && Mosek.getsolsta(t,Mosek.MSK_SOL_BAS) != Mosek.MSK_SOL_STA_UNKNOWN) ||
               (Mosek.solutiondef(t,Mosek.MSK_SOL_ITR) && Mosek.getsolsta(t,Mosek.MSK_SOL_ITR) != Mosek.MSK_SOL_STA_UNKNOWN) )
        :Unknown
    else
        sol =
            if     Mosek.solutiondef(t,Mosek.MSK_SOL_ITG) && Mosek.getsolsta(t,Mosek.MSK_SOL_ITG) != Mosek.MSK_SOL_STA_UNKNOWN
                Mosek.MSK_SOL_ITG
            elseif Mosek.solutiondef(t,Mosek.MSK_SOL_BAS) && Mosek.getsolsta(t,Mosek.MSK_SOL_BAS) != Mosek.MSK_SOL_STA_UNKNOWN
                Mosek.MSK_SOL_BAS
            else
                Mosek.MSK_SOL_ITR
            end

        prosta = Mosek.getprosta(t,sol)
        solsta = Mosek.getsolsta(t,sol)

        if  solsta == Mosek.MSK_SOL_STA_PRIM_AND_DUAL_FEAS
            #:Suboptimal # Not a standard code
            :Unknown
        elseif solsta == Mosek.MSK_SOL_STA_PRIM_FEAS
            if Mosek.getnumintvar(t) > 0
                #:Suboptimal # Not a standard code
                :Unknown
            else
                :Unknown
            end
        elseif solsta == Mosek.MSK_SOL_STA_DUAL_ILLPOSED_CER ||
            solsta == Mosek.MSK_SOL_STA_PRIM_ILLPOSED_CER
            :DualityFailure
        elseif solsta == Mosek.MSK_SOL_STA_DUAL_FEAS
            :Unknown
        elseif solsta == Mosek.MSK_SOL_STA_DUAL_INFEAS_CER
            :Unbounded
        elseif solsta == Mosek.MSK_SOL_STA_PRIM_INFEAS_CER
            :Infeasible
        elseif solsta == Mosek.MSK_SOL_STA_OPTIMAL ||
            solsta == Mosek.MSK_SOL_STA_INTEGER_OPTIMAL
            :Optimal
        else
            :Unknown
        end
    end
end

function getsense(t::Mosek.Task)
    sense = Mosek.getobjsense(t)
    if sense == Mosek.MSK_OBJECTIVE_SENSE_MAXIMIZE
        :Max
    elseif sense == Mosek.MSK_OBJECTIVE_SENSE_MINIMIZE
        :Min
    else
        :None
    end
end

function setsense!(t::Mosek.Task,sense)
    if     sense == :Max
        Mosek.putobjsense(t,Mosek.MSK_OBJECTIVE_SENSE_MAXIMIZE)
    elseif sense == :Min
        Mosek.putobjsense(t,Mosek.MSK_OBJECTIVE_SENSE_MINIMIZE)
    end
end

function getobjgap(t::Mosek.Task)
    sol = getsoldef(t)
    if sol == MSK_SOL_ITG
        Mosek.getdouinf(m.task,MSK_DINF_MIO_OBJ_REL_GAP)
    else
        0.0
    end
end

getobjval(t::Mosek.Task) = Mosek.getprimalobj(t,getsoldef(t))

function makebounds(bl_ :: Vector{Float64},
                    bu_ :: Vector{Float64}) :: Tuple{Vector{Mosek.Boundkey},Vector{Float64},Vector{Float64}}
    bk = Vector{Mosek.Boundkey}(undef, length(bl_))
    bl = Vector{Float64}(undef, length(bl_))
    bu = Vector{Float64}(undef, length(bl_))

    for i in 1:length(bl_)
        if bl_[i] > -Inf
            if bu_[i] < Inf
                if abs(bu_[i]-bl_[i]) < 1e-8
                    bk[i] = Mosek.MSK_BK_FX
                    bl[i] = bl_[i]
                    bu[i] = bl_[i]
                else
                    bk[i] = Mosek.MSK_BK_RA
                    bl[i] = bl_[i]
                    bu[i] = bu_[i]
                end
            else # bu_[i] == Inf
                bk[i] = Mosek.MSK_BK_LO
                bl[i] = bl_[i]
                bu[i] = Inf
            end
        else # bl_[i] == -Inf
            if bu_[i] < Inf
                bk[i] = Mosek.MSK_BK_UP
                bl[i] = -Inf
                bu[i] = bu_[i]
            else
                bk[i] = Mosek.MSK_BK_FR
                bl[i] = -Inf
                bu[i] = Inf
            end
        end
    end

    bk,bl,bu
end

#mutable struct MosekMathProgModel <: MathProgBase.SolverInterface.AbstractMathProgModel
mutable struct MosekMathProgModel <: MathProgBase.AbstractMathProgModel
  task :: Mosek.Task
  probtype :: Int

  # numvar
  #   Number of elements used in varmap,barvarij
  numvar   :: Int
  # varmap
  #   Maps model variables to MOSEK variables.

  #   - varmap[i] > 0: it refers to MOSEK linear variable index (varmap[i]),
  #     barvarij[i] is 0.
  #   - varmap[i] < 0: it refers to MOSEK SDP variable index (-varmap[i]),
  #     and barvarij[i] is the linear linear index into the
  #     column-oriented lower triangular part of the SDP variable.
  varmap   :: Vector{Int32}
  barvarij :: Vector{Int64}

  # numbarvar
  #   Number of used elements in barvarmap.
  numbarvar :: Int32
  # barvarmap
  #   Maps Model PSD variable indexes into MOSEK barvar indexes. When
  #   using the Semidefinite interface, these are the semidefinite
  #   variables thar can be accessed.  Semidefinite variables added
  #   thorugh loadconicproblem!() are not mapped here.
  barvarmap :: Vector{Int32}

  # binvarflag
  #   Defines per variable if it is binary
  binvarflag :: Vector{Bool}

  # numcon
  #   Number of elements used in conmap
  numcon   :: Int32
  # conmap
  #   Maps Model constraints to native MOSEK constraints. Auxiliary
  #   constraints are not mapped through. Positive values map to the
  #   corresponding constraint in MOSEK, while 0 indicated a
  #   'NULL-constraint', i.e. A placeholder constraint that exists
  #   from the user's point of view, but doesn't map to a MOSEK
  #   constraint (used specifically as a place-holder constraint for
  #   conic-quadratic constraints).
  conmap   :: Vector{Int32}
  # conslack
  #   Maps Model constraint to corresponding slack variable
  #
  #   - conslack[i] == 0: No slack (for linear constraints)
  #   - conslack[i] > 0: Constraint is conic quadratic and corresponds to MOSEK variable conslack[i]
  #   - conslack[i] < 0: Constraint is PSD conic and corresponds to
  #
  #   MOSEK variable -coneslack[i]. barconij[i] is the linear index
  #   into the column-oriented lower triangular part of the SDP
  #   variable.
  conslack :: Vector{Int32}
  barconij :: Vector{Int64}

  # quadratic constraints
  numqcon :: Int32
  qconmap :: Vector{Int32}

  # options Options from MosekSolver
  options
end

mutable struct MosekMathProgModelError <: Exception
  msg :: AbstractString
end

function loadoptions_internal!(t::Mosek.Task, options)
    # write to console by default
    printstream(msg::AbstractString) = print(msg)

    be_quiet = false
    for (option,val) in options
        parname = string(option)
        if parname == "QUIET"
            be_quiet = convert(Bool,val)
        elseif startswith(parname, "MSK_IPAR_")
            Mosek.putnaintparam(t, parname, convert(Integer, val))
        elseif startswith(parname, "MSK_DPAR_")
            Mosek.putnadouparam(t, parname, convert(AbstractFloat, val))
        elseif startswith(parname, "MSK_SPAR_")
            Mosek.putnastrparam(t, parname, convert(AbstractString, val))
        elseif isa(val, Integer)
            parname = "MSK_IPAR_$parname"
            Mosek.putnaintparam(t, parname, val)
        elseif isa(val, AbstractFloat)
            parname = "MSK_DPAR_$parname"
            Mosek.putnadouparam(t, parname, val)
        elseif isa(val, AbstractString)
            parname = "MSK_SPAR_$parname"
            Mosek.putnastrparam(t, parname, val)
        else
            error("Value $val for parameter $option has unrecognized type")
        end
    end
    if ! be_quiet
        Mosek.putstreamfunc(t,Mosek.MSK_STREAM_LOG,printstream)
    end
end

function getBoundsKey(lb, ub)
    ret = convert(Int32,0)
    if lb == -Inf && ub == Inf
        ret = MSK_BK_FR
    elseif lb == ub
        ret = MSK_BK_FX
    elseif ub  == Inf
        ret = MSK_BK_LO
    elseif lb == -Inf
        ret = MSK_BK_UP
    else
        ret = MSK_BK_RA
    end
    return convert(Int32, ret) #just to be safe
end

#internal
function complbk(bk,bl)
  if bl > -Inf
    if bk in [ MSK_BK_UP, MSK_BK_RA, MSK_BK_FX ]
      MSK_BK_RA
    else
      MSK_BK_LO
    end
  else
    if bk in [ MSK_BK_UP, MSK_BK_RA, MSK_BK_FX ]
      MSK_BK_UP
    else
      MSK_BK_FR
    end
  end
end

#internal
function compubk(bk,bu)
  if bu < Inf
    if bk in [ MSK_BK_LO, MSK_BK_RA, MSK_BK_FX ]
      MSK_BK_RA
    else
      MSK_BK_UP
    end
  else
    if bk in [ MSK_BK_LO, MSK_BK_RA, MSK_BK_FX ]
      MSK_BK_LO
    else
      MSK_BK_FR
    end
  end
end

function getsoldef(t::Mosek.Task)
    for sol in (Mosek.MSK_SOL_ITG, Mosek.MSK_SOL_BAS, Mosek.MSK_SOL_ITR)
        if Mosek.solutiondef(t,sol) && Mosek.getsolsta(t,sol) != Mosek.MSK_SOL_STA_UNKNOWN
            return sol
        end
    end
    throw(MosekMathProgModelError("No solution available"))
end

#include("MosekNLPSolverInterface.jl")
include("MosekLPQCQPInterface.jl")
include("MosekConicInterface.jl")

end


text:
"""
    Nucleotide_Essentials.FastqRecord
# Components
* ID: The unique sequence identifier associated with that entry
* sequence: The nucleotide sequence of that entry
* quality: The quality scores of that entry
* filename: The original file name
"""
struct FastqRecord
    ID::Vector{SubString{String}}
    sequence::Vector{SubString{String}}
    quality::Vector{SubString{String}}
    filename::String
end

"""
    Nucleotide_Essentials.readFastq
    readFastq(Path::String)
.fastq file => readFastq(Path) => FastqRecord(ID, sequence, quality, filename)

supported keyword arguments include:

* ```Path::String```: The full or relative path to a .fastq file

# Example:
```julia
# Supply the path to a .Fastq file that you would like to import
myfastq = readFastq("myfastq.fastq")
```
"""
function readFastq(Path::String)
    try
        global file = rsplit(Path, "/", limit = 2)[2]
    catch err
        if isa(err, BoundsError)
            global file = rsplit(Path, "/", limit = 2)[1]
        end
    end
    input = read(Path, String)
    input2 = replace(input, "\n@" => "Z")
    input3 = replace(input2, "@" => "Z", count = 1)
    input4 = replace(input3, "Z" => "\n")
    split_fastq = split(input4, "\n")
    out = FastqRecord(split_fastq[2:4:end], split_fastq[3:4:end], split_fastq[5:4:end], file)
    return out
end # end function readFastq


text:
module parallel

export Pipe, send!, get!, brain, initialize, distribute!, gather!, quit!

# Use this code in a separate module by doing
# @everywhere using parallel

# Pipe is designed to act like Python's two-ended PIPE
# which means master and child can both put and take things off the queue
type Pipe
    m::Int #Process ID of master
    c::Int #Process ID of child
    mtc::RemoteRef # Used to send message from Master to Child
    ctm::RemoteRef # Used to send message from Child to Master
end

# Automatically create the remote refs by simply specifying the child process ID
function Pipe(child::Int, master::Int = 1)
    return Pipe(master, child, RemoteRef(master), RemoteRef(child))
end

# Send a message to the pipe. Automatically figures out which RemoteRef to use
# based on whether the master or child calls it. Message can be any object,
# but is typically either a Vector{Float64} or a composite type with parameters
function send!(pipe::Pipe, msg)
    id = myid()
    if id == pipe.m
        put!(pipe.mtc, msg)
    elseif id == pipe.c
        put!(pipe.ctm, msg)
    else
        println("Wrong process communication. ID's current: ", id, " master: ", pipe.m,
        " child: ", pipe.c)
    end
end

# Take a message from the pipe. Automatically figures out which RemoteRef to use
# based on whether the master or child calls it
function get!(pipe::Pipe)
    id = myid()
    if id == pipe.m
        return take!(pipe.ctm)
    elseif id == pipe.c
        return take!(pipe.mtc)
    else
        println("Wrong process communication. ID's current: ", id, " master: ", pipe.m,
        " child: ", pipe.c)
    end
end

# This function is started on each child process, looping constantly. It waits
# until the pipe has put any new information on (ie, the proposed parameters).
# Once information appears, the child takes the value, does its work and then
# puts the answer back on the pipe for the main process to consume. Then it goes
# back to waiting for more information.
# initfunc is designed to return the dataset, or any object that gets passed
# to f, specific to this process
# Key can be either a single integer, or a list of integers
function brain(pipe::Pipe, key, initfunc::Function, f::Function)
    id = myid()
    println("Initialized process (worker): $id with key: $key")

    # Load the dataset according to this key
    dset = initfunc(key)

    # Give the master process something to gather on, so that it can wait until
    # all child processes are initialized
    send!(pipe, 0.0)

    while true
        p = get!(pipe)
        #println("Process $id received parameters $p")

        #Check to see if the message is actually the quit signal
        if p == "QUIT"
            println("Process $id completing!")
            return

        # otherwise proceed with the likelihood calculation
        else
            # Likelihood function is declared elsewhere and takes as input the dataset,
            # the parameters, and the key
            ans = f(dset, key, p)

            # Return the answer to the main process
            send!(pipe, ans)
        end
    end
end

# All of the following functions are run ONLY the master process

# Set up nchild child processes and pipes and return an array of the pipes
# function takes three arguments: the dataset, keylist, and vector of parameters
function initialize{T}(nchild::Int, keylist::Vector{T}, initfunc::Function, f::Function)

    pipes = Array(Pipe, nchild)

    # The master process is always labelled 1, so nprocessors = nchild + 1
    # yet the pipes array is indexed from 1 to nchild, hence the offset
    for i=1:nchild
        child_id = i + 1
        pipes[i] = Pipe(child_id) # Create a pipe from the master process to this child ID

        # Initialize the process with the key, init function, and the function
        remotecall(child_id, brain, pipes[i], keylist[i], initfunc, f)
    end

    return pipes
end

# Distribute parameters to all sub processes for likelihood evaluation
# p can be any object, but is traditionally either a Vector{Float64} or a specific
# composite type (e.g., model.jl/Parameters) that contains as named fields all of the
# parameters
function distribute!(pipes::Vector{Pipe}, p)
    for pipe in pipes
        send!(pipe, p)
    end
end

# Collect the results from all of the sub processes and sum together
function gather!(pipes::Vector{Pipe})
    return sum([get!(pipe) for pipe in pipes])
end

# Call this at the very end of the run to close all the pipes
function quit!(pipes::Vector{Pipe})
    for pipe in pipes
        send!(pipe, "QUIT")
    end
    #Kill all the processes
    rmprocs(workers())
end

# might be worthwhile to investigate the driver code mentioned in the manual. Seems like more of an MPI startup approach (everything is loaded everywhere, and we can check to see if we are master or not).

end #module

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using PyJuliaTester
using Test

@testset "PyJuliaTester.runtests" begin
    PyJuliaTester.runtests()
end


text:
using Distributions: Uniform, Normal, MultivariateNormal
using LazySets: DefaultUniform

for N in [Float64]
    P1 = BallInf([0.6, 0.1, -1.3, -0.4, 0.2], 0.6)
    A = [1.0  0.0;
         0.0  1.0;
        -1.0  0.0;
         0.0 -1.0]
    b = [1.0, 2.0, 3.0, 4.0]
    P2 = HPolyhedron(A, b)
    P3 = VPolygon([N[0, -1], N[2, 1]])

    # Test rand samples are contained in the set
    p1 = LazySets.sample(P1)
    @test p1 ∈ P1
    p1_samples = LazySets.sample(P1, 100)
    @test sum(p1_samples .∈ Ref(P1)) == length(p1_samples)

    # default distribution
    @test LazySets.RejectionSampler(P2).distribution == [DefaultUniform(-3.0,1.0), DefaultUniform(-4.0,2.0)]

    # random-walk sampler
    for P in [P1, P2, P3]
        for sampler in [LazySets.RandomWalkSampler(true), LazySets.RandomWalkSampler(false)]
            p_samples = LazySets.sample(P, 100, sampler=sampler)
            @test all(v ∈ P for v in p_samples)
        end
    end

    # random-walk sampler
    for P in [P1, P2, P3]
        p_samples = LazySets.sample(P, 100, sampler=LazySets.CombinedSampler())
        @test all(v ∈ P for v in p_samples)
    end

    # specifying a distribution from Distributions.jl
    @test LazySets.RejectionSampler(P2, Uniform).distribution == [Uniform(-3.0, 1.0), Uniform(-4.0,2.0)]
    @test LazySets.RejectionSampler(P2, Normal).distribution == [Normal(-3.0, 1.0), Normal(-4.0, 2.0)]

    # test univariate distributions
    X = Interval(1, 2.)
    v = sample(X, 5)
    @test all(vi ∈ X for vi in v)
    d = Normal(1.5, 0.5)

    v = sample(X, 5, sampler=LazySets.RejectionSampler(d))
    @test all(vi ∈ X for vi in v)

    v = sample(X, 5, sampler=LazySets.RejectionSampler([d]))
    @test all(vi ∈ X for vi in v)

    v = sample(X, 5, sampler=LazySets.RejectionSampler([DefaultUniform(1, 2)]))
    @test all(vi ∈ X for vi in v)

    v = sample(X, 5, sampler=LazySets.RejectionSampler(DefaultUniform(1, 2)))
    @test all(vi ∈ X for vi in v)

    # test multivariate distribution
    H = rand(Hyperrectangle, dim=2)
    sampler = LazySets.RejectionSampler(MultivariateNormal(N[1 0; 0 1.]), tight=true)
    v = sample(H, 5, sampler=sampler)

    # including vertices
    for k in 0:4
        p1 = sample(P1, 10; include_vertices=k)
        @test length(p1) == 10 + k
    end
    @test length(sample(P1, 10; include_vertices=false)) == 10
    @test length(sample(P1, 10; include_vertices=true)) == 42

    # face sampling
    H = Hyperrectangle(N[1, 2], N[3, 4])
    v0 = sample(H, 10; sampler=LazySets.FaceSampler(0))
    v1 = sample(H, 10; sampler=LazySets.FaceSampler(1))
    v2 = sample(H, 10; sampler=LazySets.FaceSampler(2))
    @test all(all(vi ∈ H for vi in v) for v in [v0, v1, v2])
    vs = vertices_list(H)
    @test all(vi ∈ vs for vi in v0)  # 0-faces are the vertices
    @test all(any(vi[i:i] ∈ project(H, [i]) for i in 1:2) for vi in v1)  # 1-faces are at the border
end


text:
export grid_graph, Partition, props, Singleton, @singleton

mutable struct Singleton <: AbstractSemialgebraicSet
    x::AbstractVector
end

Singleton(x::Real) = Singleton([x])
Singleton(x::Tuple) = Singleton([x...])

macro singleton(x)
    return :(Singleton($(esc(x))))
end

struct Partition
    graph
    get_vertex
end

function grid_graph(x, lb, ub, n; inf_top = zeros(Int64, length(ub)), inf_floor = zeros(Int64, length(lb)))
    @assert all(lb .<= ub)

    x_ranges = []
    for k in 1:length(n)
        n_eff = n[k] - inf_top[k] - inf_floor[k]
        if n_eff < 0
            push!(x_ranges, [inf_floor[k] == 0 ? lb[k] : -Inf, (inf_top[k] == 0 ? ub[k] : Inf)])
        elseif n_eff == 0
            push!(x_ranges, [inf_floor[k] == 0 ? lb[k] : -Inf, (lb[k] + ub[k])/2, (inf_top[k] == 0 ? ub[k] : Inf)])
        else
            push!(x_ranges, [inf_floor[k] == 0 ? lb[k] : -Inf, range(lb[k], stop=ub[k], length=n_eff)..., (inf_top[k] == 0 ? ub[k] : Inf)])
        end
    end
    return grid_graph(x, x_ranges)
end

function grid_graph(x, x_ranges)
    n = length.(x_ranges) .- 1
    nv = prod(n)

    mg = MetaDiGraph(nv)
    for i in 1:nv
        idx = invert_index(i, n)
        subset = FullSpace()
        for k in 1:length(idx)
            if !isinf(x_ranges[k][idx[k]])
                subset = intersect(subset, @set(x[k] >= x_ranges[k][idx[k]]))
            end
            if !isinf(x_ranges[k][idx[k]+1])
                subset = intersect(subset, @set(x[k] <= x_ranges[k][idx[k]+1]))
            end
        end
        set_prop!(mg, i, :cell, subset)
        for k in 1:length(idx)
            if idx[k] < n[k]
                idx[k] += 1
                j = linearize_index(idx, n)
                add_edge!(mg, i, j)
                set_prop!(mg, Edge(i,j), :interface, [@set(x[k] == x_ranges[k][idx[k]])])
                idx[k] -= 1
            end
        end
    end
    get_vertex = function (x)
                    idx = zeros(Int64, length(n))
                    for k in 1:length(n)
                        j = findfirst(m -> m >= x[k], x_ranges[k]) - 1
                        if isnothing(j)
                            error("state outside domain")
                        end
                        idx[k] = max(j, 1)
                    end
                    return linearize_index(idx, n)
                 end
    return mg, get_vertex
end

function check_membership(X::AbstractSemialgebraicSet, x_var, x)
    is_member = check_membership(X.V)
    if is_member
        for ineq in inequalities(X)
            if ineq(x_var => x) < 0
                is_member = false
                break
            end
        end
    end
    return is_member
end

function check_membership(X::AbstractAlgebraicSet, x_var, x)
    is_member = true
    for eq in equalities(X)
        if eq(x_var => x) != 0
            is_member = false
            break
        end
    end
    return is_member
end

check_membership(X::FullSpace, x) = true

# not elegant but at least correct!
function complement(X::BasicSemialgebraicSet, H = FullSpace())
    Ys = BasicSemialgebraicSet[]
    @assert isempty(equalities(X)) "Equality constraints are not allowed when computing the complement"
    ineqs = inequalities(X)
    m = length(ineqs)
    for i in 1:m
        push!(Ys, intersect(H, BasicSemialgebraicSet(algebraicset(Polynomial{true, Float64}[]),
                                                    vcat(ineqs[1:end-1], -ineqs[end]))))
        pop!(ineqs)
    end
    return Ys
end

# in general this is conservative but holds


text:
using AWSInfinidash
using Test

@testset "AWSInfinidash.jl" begin
    @test i8h isa Function
    @test I8h isa Module
end


text:
@enum OraResult::Int32 begin
    ORA_SUCCESS = 0
    ORA_FAILURE = -1
end

"""
This enumeration identifies the mode to use when authorizing connections to the database.
"""
@enum OraAuthMode::UInt32 begin
    ORA_MODE_AUTH_DEFAULT = 0
    ORA_MODE_AUTH_SYSDBA  = 2
    ORA_MODE_AUTH_SYSOPER = 4
    ORA_MODE_AUTH_PRELIM  = 8
    ORA_MODE_AUTH_SYSASM  = 32768
    ORA_MODE_AUTH_SYSBKP  = 131072
    ORA_MODE_AUTH_SYSDGD  = 262144
    ORA_MODE_AUTH_SYSKMT  = 524288
    ORA_MODE_AUTH_SYSRAC  = 1048576
end

@enum OraConnCloseMode::UInt32 begin
    ORA_MODE_CONN_CLOSE_DEFAULT = 0
    ORA_MODE_CONN_CLOSE_DROP    = 1
    ORA_MODE_CONN_CLOSE_RETAG   = 2
end

"""
This enumeration identifies the mode to use when creating connections to the database. Note that the OCI objects mode is always enabled.
"""
@enum OraCreateMode::UInt32 begin
    ORA_MODE_CREATE_DEFAULT  = 0
    ORA_MODE_CREATE_THREADED = 1
    ORA_MODE_CREATE_EVENTS   = 4
end

"""
This enumeration identifies the mode to use when getting sessions from a session pool.
"""
@enum OraPoolGetMode::UInt8 begin
    ORA_MODE_POOL_GET_WAIT      = 0
    ORA_MODE_POOL_GET_NOWAIT    = 1
    ORA_MODE_POOL_GET_FORCEGET  = 2
    ORA_MODE_POOL_GET_TIMEDWAIT = 3
end

"""
This enumeration identifies the purity of the sessions that are acquired when using connection classes during connection creation.
"""
@enum OraPurity::UInt32 begin
    ORA_PURITY_DEFAULT = 0
    ORA_PURITY_NEW     = 1
    ORA_PURITY_SELF    = 2
end

@enum OraShutdownMode::UInt32 begin
    ORA_MODE_SHUTDOWN_DEFAULT             = 0
    ORA_MODE_SHUTDOWN_TRANSACTIONAL       = 1
    ORA_MODE_SHUTDOWN_TRANSACTIONAL_LOCAL = 2
    ORA_MODE_SHUTDOWN_IMMEDIATE           = 3
    ORA_MODE_SHUTDOWN_ABORT               = 4
    ORA_MODE_SHUTDOWN_FINAL               = 5
end

@enum OraStartupMode::UInt32 begin
    ORA_MODE_STARTUP_DEFAULT  = 0
    ORA_MODE_STARTUP_FORCE    = 1
    ORA_MODE_STARTUP_RESTRICT = 2
end

@enum OraOracleTypeNum::UInt32 begin
    ORA_ORACLE_TYPE_NONE          = 2000
    ORA_ORACLE_TYPE_VARCHAR       = 2001
    ORA_ORACLE_TYPE_NVARCHAR      = 2002
    ORA_ORACLE_TYPE_CHAR          = 2003
    ORA_ORACLE_TYPE_NCHAR         = 2004
    ORA_ORACLE_TYPE_ROWID         = 2005
    ORA_ORACLE_TYPE_RAW           = 2006
    ORA_ORACLE_TYPE_NATIVE_FLOAT  = 2007
    ORA_ORACLE_TYPE_NATIVE_DOUBLE = 2008
    ORA_ORACLE_TYPE_NATIVE_INT    = 2009
    ORA_ORACLE_TYPE_NUMBER        = 2010
    ORA_ORACLE_TYPE_DATE          = 2011
    ORA_ORACLE_TYPE_TIMESTAMP     = 2012
    ORA_ORACLE_TYPE_TIMESTAMP_TZ  = 2013
    ORA_ORACLE_TYPE_TIMESTAMP_LTZ = 2014
    ORA_ORACLE_TYPE_INTERVAL_DS   = 2015
    ORA_ORACLE_TYPE_INTERVAL_YM   = 2016
    ORA_ORACLE_TYPE_CLOB          = 2017
    ORA_ORACLE_TYPE_NCLOB         = 2018
    ORA_ORACLE_TYPE_BLOB          = 2019
    ORA_ORACLE_TYPE_BFILE         = 2020
    ORA_ORACLE_TYPE_STMT          = 2021
    ORA_ORACLE_TYPE_BOOLEAN       = 2022
    ORA_ORACLE_TYPE_OBJECT        = 2023
    ORA_ORACLE_TYPE_LONG_VARCHAR  = 2024
    ORA_ORACLE_TYPE_LONG_RAW      = 2025
    ORA_ORACLE_TYPE_NATIVE_UINT   = 2026
    ORA_ORACLE_TYPE_JSON          = 2027
    ORA_ORACLE_TYPE_JSON_OBJECT   = 2028
    ORA_ORACLE_TYPE_JSON_ARRAY    = 2029
    ORA_ORACLE_TYPE_MAX           = 2030
end

@enum OraNativeTypeNum::UInt32 begin
    ORA_NATIVE_TYPE_INT64       = 3000
    ORA_NATIVE_TYPE_UINT64      = 3001
    ORA_NATIVE_TYPE_FLOAT       = 3002
    ORA_NATIVE_TYPE_DOUBLE      = 3003
    ORA_NATIVE_TYPE_BYTES       = 3004
    ORA_NATIVE_TYPE_TIMESTAMP   = 3005
    ORA_NATIVE_TYPE_INTERVAL_DS = 3006
    ORA_NATIVE_TYPE_INTERVAL_YM = 3007
    ORA_NATIVE_TYPE_LOB         = 3008
    ORA_NATIVE_TYPE_OBJECT      = 3009
    ORA_NATIVE_TYPE_STMT        = 3010
    ORA_NATIVE_TYPE_BOOLEAN     = 3011
    ORA_NATIVE_TYPE_ROWID       = 3012
    ORA_NATIVE_TYPE_JSON        = 3013
    ORA_NATIVE_TYPE_JSON_OBJECT = 3014
    ORA_NATIVE_TYPE_JSON_ARRAY  = 3015
    ORA_NATIVE_TYPE_NULL        = 3016
    #ORA_NATIVE_TYPE_NUMBER      = xxxx # see issue #21
end

@enum OraStatementType::UInt16 begin
    ORA_STMT_TYPE_UNKNOWN      = 0
    ORA_STMT_TYPE_SELECT       = 1
    ORA_STMT_TYPE_UPDATE       = 2
    ORA_STMT_TYPE_DELETE       = 3
    ORA_STMT_TYPE_INSERT       = 4
    ORA_STMT_TYPE_CREATE       = 5
    ORA_STMT_TYPE_DROP         = 6
    ORA_STMT_TYPE_ALTER        = 7
    ORA_STMT_TYPE_BEGIN        = 8
    ORA_STMT_TYPE_DECLARE      = 9
    ORA_STMT_TYPE_CALL         = 10
    ORA_STMT_TYPE_EXPLAIN_PLAN = 15
    ORA_STMT_TYPE_MERGE        = 16
    ORA_STMT_TYPE_ROLLBACK     = 17
    ORA_STMT_TYPE_COMMIT       = 21
end

@enum OraPoolCloseMode::UInt32 begin
    ORA_MODE_POOL_CLOSE_DEFAULT = 0x0000
    ORA_MODE_POOL_CLOSE_FORCE   = 0x0001
end


text:
# p2.jl - convergence of periodic spectral method (compare p1.jl)

# For various N (even), set up grid as before:
clf(); axes([.1,.4,.8,.5]);
for N = 2:2:100
    h = 2*pi/N;  x = [ -pi + i*h for i = 1:N ];
    u = @. exp(sin(x)); uprime = @. cos(x)*u;

    # Construct spectral differentiation matrix:
    column = [ 0; @. .5*(-1)^(1:N-1)*cot((1:N-1)*h/2) ];
    D = toeplitz(column,column[[1;N:-1:2]]);

    # Plot max(abs(D*u-uprime)):
    error = norm(D*u-uprime,Inf);
    loglog(N,error,"k.",markersize=6);
end
grid(true); xlabel("N"); ylabel("error");
title("Convergence of spectral differentiation")


text:
module VQE

#export min_expec_val, expec_val, ansatz_gen, run_circuit_with_qiskit

include("GateOps.jl")
include("Utils.jl")

using PyCall
using Optim
using UUIDs

###############################################################################
#               Module setup
###############################################################################

#Track the number of circuits generated
cct_runs = 0

#Create temporary directory hold intermdiate QASM files
cct_dir = "/tmp/cct_dir_" * string(UUIDs.uuid4()) * "/"

function __init__()
    init_env()
    mkdir(cct_dir)
end

"""
    init_env()

Installs necessary packages to run the VQE method. Qiskit is installed into the
local Python environment packaged weith Julia's conda.
"""
function init_env()
    run(`$(PyCall.python) -m pip install qiskit`)
end

"""
Cleanup the generated files.
"""
function refresh_dir()
    rm(cct_dir, recursive=true)
    mkdir(cct_dir)
end

###############################################################################
#               Circuit evaulation functions
###############################################################################

"""
    run_circuit_with_qiskit(cct::String, num_samples::Int64)

Generate an ansatz using qubit register `q_reg` on indices `qubit_indices` with angle `theta`.
Sample ansatz uses R_x(theta) gate.
"""
function run_circuit_with_qiskit(cct::String, num_samples::Int64)
    global cct_runs += 1
    qiskit = PyCall.pyimport("qiskit");
    qiskit_execute = PyCall.pyimport("qiskit.execute");

    qubit_label_q = "qr"
    qubit_label_c = "cr"
    q_count = 1
    c_count = q_count

    qubit_indices = collect(0:q_count-1)

    qasm_file = GateOps.gen_qasm_header(q_count, c_count, qubit_label_q, qubit_label_c )
    qasm_file *= cct

    cct_file = open(cct_dir * "cct_$cct_runs.qasm", "w")
    for i in Utils.format_string_nl(qasm_file)
        write(cct_file, i)
    end
    close(cct_file)

    circ = qiskit.QuantumCircuit.from_qasm_file(cct_dir * "cct_$cct_runs.qasm")

    #Run the circuit
    simulator = qiskit.Aer.get_backend("qasm_simulator")
    job = qiskit_execute.execute(circ, simulator, shots=num_samples)
    result = job.result()
    counts = result.get_counts(circ)
    print("\nTotal count for 00 and 11 are:", counts)

    return counts;
end

"""
    ansatz_gen(q_reg::String, qubit_indices::Vector{Int}, theta::Float64)

Generate an ansatz using qubit register `q_reg` on indices `qubit_indices` with angle `theta`.
Sample ansatz uses R_x(theta) gate.
"""
function ansatz_gen(q_reg::String, qubit_indices::Vector{Int}, theta::Float64)
    ansatz = GateOps.reset_reg(qubit_indices, q_reg)
    for idx in qubit_indices
        ansatz *= GateOps.apply_gate_rx(q_reg, theta, idx)
    end
    return ansatz
end

"""
    expec_val(theta::Real)

Evaluate the expectation value of the given test-case Pauli-Z.
"""
function expec_val(theta::Real, num_samples::Int64)
    q_reg="qr"
    c_reg="cr"
    qubit_indices=[0]
    circ = ""
    circ *= ansatz_gen(q_reg, qubit_indices, theta)
    circ *= GateOps.measure_qubits_all(qubit_indices, qubit_indices, q_reg, c_reg)

    c = run_circuit_with_qiskit(circ, num_samples, )

    p0 = [val for (key, val) in c if key=="0"];
    p1 = [val for (key, val) in c if key=="1"];

    if length(p0) == 0
        p0 = [0]
    elseif length(p1) == 0
        p1 = [0]
    end
    return (p0[1] - p1[1]) / num_samples
end

"""
    min_expec_val(theta::Real)

Find the minimum expectation value of the given Hamiltonian (Pauli-Z).
Uses Optim.jl package to perform classical optimisation.
"""
function min_expec_val(start_val::Real, num_samples::Int64=1000)
    #Clear previous data before performing new experiment
    global cct_runs = 0
    refresh_dir()

    f(x) = expec_val(x[1], num_samples)
    min_val = Optim.optimize(f, [start_val], Optim.BFGS())
    return min_val
end

end