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::SBPParabolicReduced2SC{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)

  # this doesn't seem to help
#  computeDirichletBoundaryTerm(mesh, sbp, eqn, opts, capture, shockmesh)
  #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: [`SBPParabolicReduced2SC`](@ref)
   * shockmesh: the `ShockedElements`
   * entropy_vars: an [`AbstractVariables`](@ref) object
   * diffusion: an [`AbstractDiffusion`](@ref)
"""
function computeGradW(mesh, sbp, eqn, opts, capture::SBPParabolicReduced2SC{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::SBPParabolicReduced2SC{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::SBPParabolicReduced2SC{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)
    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)

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

    # 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::SBPParabolicReduced2SC{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
=#

"""
  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::SBPParabolicReduced2SC{Tsol, Tres},
                          sbpface::AbstractFace, iface_red::Interface,
                          wL::AbstractMatrix, wR::AbstractMatrix,
                          gradwL::Abstract3DArray, gradwR::Abstract3DArray,
                          nrm_face::AbstractMatrix,
                          delta_w::AbstractMatrix,
                         ) 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

  return nothing
end

#------------------------------------------------------------------------------
# Dirichlet boundary condition

const ParabolicReduceds{Tsol, Tres} = Union{SBPParabolicReducedSC{Tsol, Tres},
                                            SBPParabolicReduced2SC{Tsol, Tres}
                                           }

"""
  Computes a Dirichlet BC that is consistent with the inviscid one
"""
@noinline function computeDirichletBoundaryTerm(mesh, sbp, eqn, opts,
                      capture::ParabolicReduceds{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::ParabolicReduceds{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()  # the - sign is applied with B

  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]
    alpha = 1
    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)
#    boundaryFaceInterpolate_rev!(mesh.sbpface, bndry_i.face, res_i, delta_w)

    # 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]  # the sat has a - sign in
                                                 # front of it
      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)  # the op is because the B_gamma has
                                        # a negative sign

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

  return nothing
end

"""
  Method to apply Dgk^T * t2L only.
"""
function applyDgkTranspose(capture::ParabolicReduceds{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

"""
  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::ParabolicReduceds{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

#------------------------------------------------------------------------------
# penalty functions

"""
  Applies the penalty for the SBP-SAT generalization of the modified scheme
  of Bassi and Rebay (BR2).
"""
function applyReduced2Penalty(penalty::BR2Penalty{Tsol, Tres}, sbp, params::AbstractParamType,
                      sbpface,
                      diffusion::AbstractDiffusion, iface::Interface,
                      delta_w::AbstractMatrix{Tsol}, theta::AbstractMatrix{Tres},
                      qL_el::AbstractMatrix, qR_el::AbstractMatrix,
                      wL::AbstractMatrix, wR::AbstractMatrix,
                      coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                      jacL::AbstractVector, jacR::AbstractVector,
                      res1L::AbstractMatrix, res1R::AbstractMatrix,
                     ) where {Tsol, Tres}

  fill!(res1L, 0); fill!(res1R, 0)
  numDofPerNode, numNodesPerFace = size(delta_w)
  numNodesPerElement = length(jacL)
  dim = size(nrm_face, 1)

  @unpack penalty delta_w_n qL qR t1L t1R t2L t2R
  fill!(qL, 0); fill!(qR, 0)

  #--------------------------
  # apply T1
  # multiply by normal vector, then R^T B
  #alpha_g = 1/(dim + 1)  # = 1/number of faces of a simplex
  @simd for d1=1:dim
    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        delta_w_n[k, j] = 0.25*delta_w[k, j]*nrm_face[d1, j]
      end
    end

    qL_d = sview(qL, :, :, d1); qR_d = sview(qR, :, :, d1)
    interiorFaceIntegrate!(sbpface, iface, delta_w_n, qL_d, qR_d)

    # interiorFaceIntegrates -= the second output
    @simd for j=1:numNodesPerElement
      @simd for k=1:numDofPerNode
        qR_d[k, j] = -qR_d[k, j]
      end
    end
  end

  # apply Lambda matrix
  applyDiffusionTensor(diffusion, sbp, params, qL_el, wL, coordsL, dxidxL, jacL, iface.elementL, sbpface, iface.faceL,
                       qL, t1L)
  applyDiffusionTensor(diffusion, sbp, params, qR_el, wR, coordsR, dxidxR, jacR,
                       iface.elementR, sbpface, iface.faceR,
                       qR, t1R)

  # apply inverse mass matrix, then apply B*Nx*R*t2L_x + B*Ny*R*t2L_y
  @simd for d1=1:dim
    @simd for j=1:numNodesPerElement
      facL = jacL[j]/sbp.w[j]
      facR = jacR[j]/sbp.w[j]
      #facL = (3)*jacL[j]/sbp.w[j]
      #facR = (3)*jacR[j]/sbp.w[j]

      @simd for k=1:numDofPerNode
        t1L[k, j, d1] *= facL
        t1R[k, j, d1] *= facR
      end
    end

    #TODO: make t2L 2D rather than 3D
    t1L_d = ro_sview(t1L, :, :, d1); t1R_d = ro_sview(t1R, :, :, d1)
    t2L_d = sview(t2L, :, :, d1);    t2R_d = sview(t2R, :, :, d1)
    interiorFaceInterpolate!(sbpface, iface, t1L_d, t1R_d, t2L_d, t2R_d)

    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        val = sbpface.wface[j]*nrm_face[d1, j]*(t2L_d[k, j] + t2R_d[k, j])
        res1L[k, j] += val
        res1R[k, j] -= val  # delta_w is reversed for elementR
      end
    end

  end  # end d1

#=
  #----------------------------
  # apply T2 and T3
  @simd for j=1:numNodesPerFace
    @simd for k=1:numDofPerNode
      val1 =  0.5*sbpface.wface[j]*theta[k, j]
      val2 = -0.5*sbpface.wface[j]*delta_w[k, j]
      res1L[k, j] += val1
      res1R[k, j] += val1
      res2L[k, j] += val2
      res2R[k, j] -= val2  # delta_w is reversed for elementR
    end
  end
=#
  return nothing
end


text:
#TODO: rename *temperature to more genereral terms
@system LeafLengthTracker(BetaFunction) begin
    pheno: phenology ~ ::Phenology(override)
    T: temperature ~ ::int(override)
    #FIXME: leaves_potential is already max(leaves_generic, leaves_total)?
    n(np=pheno.leaves_potential, ng=pheno.leaves_generic): leaf_count => max(np, ng) ~ track
    Tn: minimum_temperature => 0 ~ preserve
    To(n): optimal_temperature => 0.88n ~ track
    Tx(n): maximum_temperature => 1.64n ~ track
end

@system LeafColdInjury begin
    T: temperature ~ hold

    _a: cold_injury_factor1 => -0.1 ~ preserve(u"K^-1", parameter)
    _b: cold_injury_factor2 => 1.6 ~ preserve(parameter)
    _Tc: cold_injury_critical_temperature => 0 ~ preserve(u"°C", parameter)

    CID(T, _Tc): cold_injury_duration => begin
        T < Tc ? 1 : -1
    end ~ accumulate(u"hr", min=0)

    "preliminary cold injury effect (2019-05-23: KDY)"
    ACIE(_a, _b, _Tc, T, CID): apparent_cold_injury_effect => begin
        a = T < Tc ? log(a * (T - Tc) + b) : 0
        1 - a / exp(1u"hr" / CID)
    end ~ track(min=0, max=1)

    _enable => true ~ flag(parameter)

    CIE(CIE, ACIE, CID): cold_injury_effect => begin
        iszero(CID) ? 1 : min(CIE, ACIE)
    end ~ track(init=1, when=_enable)
end

@system Leaf(Organ, LeafColdInjury) begin
    rank ~ ::int(override) # preserve

    # cm dd-1 Fournier and Andrieu 1998 Pg239.
    # This is the "potential" elongation rate with no water stress Yang
    # elongation_rate => 0.564 ~ preserve(u"cm/d", parameter)

    # max elongation rate (cm per day) at optimal temperature
    # (Topt: 31C with Tbase = 9.8C using 0.564 cm/dd rate from Fournier 1998 paper above
    LER_max: maximum_elongation_rate => 12 ~ preserve(u"cm/d", parameter)

    LM_min: minimum_length_of_longest_leaf => 60 ~ preserve(u"cm", parameter)

    # leaf lamina width to length ratio
    length_to_width_ratio => begin
        #0.106 # for maize
        0.05 # for Radish
    end ~ preserve(parameter)

    # leaf area coeff with respect to L*W (A_LW)
    area_ratio => 0.75 ~ preserve(parameter)

    # staygreen trait of the hybrid
    # stay green for this value times growth period after peaking before senescence begins
    # An analogy for this is that with no other stresses involved,
    # it takes 15 years to grow up, stays active for 60 years,
    # and age the last 15 year if it were for a 90 year life span creature.
    # Once fully grown, the clock works differently so that the hotter it is quicker it ages
    SG: stay_green => 3.5 ~ preserve(parameter)

    maximum_aging_rate(LER_max) => LER_max ~ track(u"cm/d")

    #############
    # Variables #
    #############

    #FIXME
    potential_leaves(pheno.leaves_potential) ~ track

    #FIXME
    extra_leaves(np=potential_leaves, ng=pheno.leaves_generic) => (np - ng) ~ track

    k: maximum_length_of_longest_leaf_adjustment => begin
        # no length adjustment necessary for Radish, unlike MAIZE (KY, 2016-10-12)
        # 24.0
        0
    end ~ preserve(u"cm^2", parameter)

    maximum_length(LM_min, extra_leaves, k): maximum_length_of_longest_leaf => begin
        sqrt(LM_min^2 + k * extra_leaves)
    end ~ track(u"cm")

    maximum_width(l=maximum_length, r=length_to_width_ratio) => begin
        # Fournier and Andrieu(1998) Pg242 YY
        l * r
    end ~ track(u"cm")

    maximum_area(l=maximum_length, w=maximum_width, r=area_ratio) => begin
        # daughtry and hollinger (1984) Fournier and Andrieu(1998) Pg242 YY
        l * w * r
    end ~ track(u"cm^2")

    area_from_length(; L(u"cm")) => begin
        #HACK ensure zero area for zero length
        # for Radish, see JH's thesis
        l = Cropbox.deunitfy(L)
        iszero(l) ? l : 0.639945 + 0.954957l + 0.005920l^2
    end ~ call(u"cm^2")

    area_increase_from_length(length) => begin
        # for Radish, see JH's thesis
        l = Cropbox.deunitfy(length)
        0.954957 + 2*0.005920l
    end ~ track(u"cm^2")

    #TODO better name, shared by growth_duration and potential_area
    #TODO should be a plant parameter not leaf (?)
    rank_effect(rank, n=potential_leaves, weight=1) => begin
        n_m = 5.93 + 0.33n # the rank of the largest leaf. YY
        a = (-10.61 + 0.25n) * weight
        b = (-5.99 + 0.27n) * weight
        # equation 7 in Fournier and Andrieu (1998). YY

        # equa 8(b)(Actually eqn 6? - eqn 8 deals with leaf age - DT)
        # in Fournier and Andrieu(1998). YY
        scale = rank / n_m - 1
        exp(a * scale^2 + b * scale^3)
    end ~ track

    potential_length_tracker(context, pheno, temperature=rank) ~ ::LeafLengthTracker
    potential_length(maximum_length, β=potential_length_tracker.ΔT) => begin
        # for MAIZSIM
        #self.maximum_length * self.rank_effect(weight=0.5)
        # for beta fn calibrated from JH's thesis for SP and KM varieties, 8/10/15, SK
        maximum_length * β
    end ~ track(u"cm")

    # from CLeaf::calc_dimensions()
    # LM_min is a length characteristic of the longest leaf,in Fournier and Andrieu 1998, it was 90 cm
    # LA_max is a fn of leaf no (Birch et al, 1998 fig 4) with largest reported value near 1000cm2. This is implemented as lfno_effect below, SK
    # LM_min of 115cm gives LA of largest leaf 1050cm2 when totalLeaves are 25 and Nt=Ng, SK 1-20-12
    # Without lfno_effect, it can be set to 97cm for the largest leaf area to be at 750 cm2 with Nt ~= Ng (Lmax*Wmax*0.75) based on Muchow, Sinclair, & Bennet (1990), SK 1-18-2012
    # Eventually, this needs to be a cultivar parameter and included in input file, SK 1-18-12
    # the unit of k is cm^2 (Fournier and Andrieu 1998 Pg239). YY
    # L_max is the length of the largest leaf when grown at T_peak. Here we assume LM_min is determined at growing Topt with minmal (generic) leaf no, SK 8/2011
    # If this routine runs before TI, totalLeaves = genericLeafNo, and needs to be run with each update until TI and total leaves are finalized, SK
    GD(potential_length, LER_max): growth_duration => begin
        # shortest possible linear phase duration in physiological time (days instead of GDD) modified
        days = potential_length / LER_max
        # for Radish
        1.5days
    end ~ track(u"d")

    phase1_delay(rank) => begin
        # not used in MAIZSIM because LTAR is used to initiate leaf growth.
        # Fournier's value : -5.16+1.94*rank;equa 11 Fournier and Andrieu(1998) YY, This is in plastochron unit
        -5.16 + 1.94rank
    end ~ track(min=0)

    leaf_number_effect(potential_leaves) => begin
        # Fig 4 of Birch et al. (1998)
        exp(-1.17 + 0.047potential_leaves)
    end ~ track(min=0.5, max=1.0)

    potential_area(potential_length, area_from_length) => begin
        # for MAIZSIM
        # equa 6. Fournier and Andrieu(1998) multiplied by Birch et al. (1998) leaf no effect
        # LA_max the area of the largest leaf
        # PotentialArea potential final area of a leaf with rank "n". YY
        #self.maximum_area * self.leaf_number_effect * self.rank_effect(weight=1)
        # for Radish
        area_from_length(potential_length)
    end ~ track(u"cm^2")

    green_ratio(senescence_ratio) => (1 - senescence_ratio) ~ track

    green_area(green_ratio, area) => (green_ratio * area) ~ track(u"cm^2")

    #TODO implement Parent and Tardieu (2011, 2012) approach for leaf elongation in response to T and VPD, and normalized at 20C, SK, Nov 2012
    # elongAge indicates where it is now along the elongation stage or duration.
    # duration is determined by totallengh/maxElongRate which gives the shortest duration to reach full elongation in the unit of days.
    elongation_age(pheno.BF.ΔT) ~ accumulate(when=growing, u"d")

    #TODO move to common module (i.e. Organ?)
    beta_growth(t_b=0u"d", delta=1; t(u"d"), t_e(u"d"), c_m(u"cm/d")) => begin
        t = clamp(t, zero(t), t_e)
        t_m = t_e / 2 #TODO: allow custom value again?
        t_et = t_e - t
        t_em = t_e - t_m
        t_tb = t - t_b
        t_mb = t_m - t_b
        c_m * ((t_et / t_em) * (t_tb / t_mb)^(t_mb / t_em))^delta
    end ~ call(u"cm/d")

    potential_elongation_rate(beta_growth, elongation_age, LER_max, GD) => begin
        #TODO proper integration with scipy.integrate?
        beta_growth(elongation_age, GD, LER_max)
    end ~ track(when=growing, u"cm/d")

    #TODO: incorporate stress effects
    actual_elongation_rate(potential_elongation_rate, cold_injury_effect) => begin
        potential_elongation_rate * cold_injury_effect
    end ~ track(u"cm/d")

    temperature_effect_func(; T_grow(u"°C"), T_peak(u"°C"), T_base(u"°C")) => begin
        # T_peak is the optimal growth temperature at which the potential leaf size determined in calc_mophology achieved.
        # Similar concept to fig 3 of Fournier and Andreiu (1998)

        # phyllochron corresponds to PHY in Lizaso (2003)
        # phyllochron needed for next leaf appearance in degree days (GDD8) - 08/16/11, SK.
        #phyllochron = (dv->get_T_Opt()- Tb)/(dv->get_Rmax_LTAR());

        T_ratio = (T_grow - T_base) / (T_peak - T_base)
        # final leaf size is adjusted by growth temperature determining cell size during elongation
        max(T_ratio * exp(1 - T_ratio), zero(T_ratio))
    end ~ call

    temperature_effect => begin
        #temperature_effect_func(self.p.pheno.growing_temperature, 18.7, 8.0) # for MAIZSIM
        #FIXME Radish model uses current temperature, not average growing temperature
        #temperature_effect_func(self.p.pheno.temperature, self.p.pheno.optimal_temperature, 0) # for Radish
        #FIXME Radish model does not actually use temperature effect on final leaf size calculation
        1.0 # for Radish
    end ~ track

    potential_expansion_rate(t=elongation_age, t_e=GD, w_max=potential_area) => begin
        # t_e = 1.5 * w_max / c_m
        t = min(t, t_e)
        #FIXME can we introduce new w_max here when w_max in t_e (growth duration) supposed to be potential length?
        # c_m from Eq. 9, r (= dw/dt / c_m) from Eq. 7 of Yin (2003)
        #HACK can be more simplified
        #c_m = 1.5 / t_e * w_max
        #r = 4t * (t_e - t) / t_e^2
        t_m = t_e / 2
        c_m = (2t_e - t_m) / (t_e * (t_e - t_m)) * (t_m / t_e)^(t_m / (t_e - t_m)) * w_max
        r = (t_e - t) / (t_e - t_m) * (t / t_m)^(t_m / (t_e - t_m))
        #FIXME dt here is physiological time, whereas timestep multiplied in potential_area_increase is chronological time
        c_m * r # dw/dt
    end ~ track(u"cm^2/d")

    potential_area_increase(area_from_length, length, actual_length_increase, area) => begin
        ##area = max(0, water_effect * T_effect * self.potential_area * (1 + (t_e - self.elongation_age) / (t_e - t_m)) * (self.elongation_age / t_e)**(t_e / (t_e - t_m)))
        #maximum_expansion_rate = T_effect * self.potential_area * (2*t_e - t_m) / (t_e * (t_e - t_m)) * (t_m / t_e)**(t_m / (t_e - t_m))
        # potential leaf area increase without any limitations
        #max(0, maximum_expansion_rate * max(0, (t_e - self.elongation_age) / (t_e - t_m) * (self.elongation_age / t_m)**(t_m / (t_e - t_m))) * self.timestep)
        # for MAIZSIM
        #self.potential_expansion_rate * self.timestep
        # for Radish
        #TODO need common framework dealing with derivatives
        #area_increase_from_length(actual_length_increase)
        area_from_length(length + actual_length_increase) - area
    end ~ track(when=growing, u"cm^2")

    # create a function which simulates the reducing in leaf expansion rate
    # when predawn leaf water potential decreases. Parameterization of rf_psil
    # and rf_sensitivity are done with the data from Boyer (1970) and Tanguilig et al (1987) YY
    #FIXME: unit of call args kPa?
    water_potential_effect_func(; psi_predawn, psi_th) => begin
        #psi_predawn = self.p.soil.WP_leaf_predawn
        # psi_th: threshold wp below which stress effect shows up

        # DT Oct 10, 2012 changed this so it was not as sensitive to stress near -0.5 lwp
        # SK Sept 16, 2014 recalibrated/rescaled parameter estimates in Yang's paper. The scale of Boyer data wasn't set correctly
        # sensitivity = 1.92, LeafWPhalf = -1.86, the sensitivity parameter may be raised by 0.3 to 0.5 to make it less sensitivity at high LWP, SK
        s_f = 0.4258 # 0.5
        psi_f = -1.4251 # -1.0
        e = (1 + exp(s_f * psi_f)) / (1 + exp(s_f * (psi_f - (psi_predawn - psi_th))))
        min(e, 1)
    end ~ call

    water_potential_effect(; threshold) => begin
        # for MAIZSIM
        #water_potential_effect_func(soil.WP_leaf_predawn, threshold)
        # for Radish
        1.0
    end ~ call

    carbon_effect => 1.0 ~ track

    length(actual_elongation_rate) ~ accumulate(u"cm", time=elongation_age, timeunit=u"d")
    actual_length_increase(actual_elongation_rate) ~ capture(u"cm", time=elongation_age, timeunit=u"d")

    # actual area
    area(water_potential_effect, carbon_effect, temperature_effect, area_from_length, length) => begin
        # See Kim et al. (2012) Agro J. for more information on how this relationship has been determined based on multiple studies and is applicable across environments
        #FIXME: unit of threshold kPa?
        we = water_potential_effect(-0.8657)

        # place holder
        ce = carbon_effect
        te = temperature_effect

        # growth temperature effect is now included here, outside of potential area increase calculation
        #TODO water and carbon effects are not multiplicative?
        min(we, ce) * te * area_from_length(length)
    end ~ track(u"cm^2")

    #TODO remove if unnecessary
    # @property
    # def actual_area_increase(self):
    #     #FIXME area increase tracking should be done by some global state tracking manager
    #     raise NotImplementedError("actual_area_increase")
    #
    # @property
    # def relative_area_increase(self):
    #     #HACK meaning changed from 'relative to other leaves' (spatial) to 'relative to previous state' (temporal)
    #     # adapted from CPlant::calcPerLeafRelativeAreaIncrease()
    #     #self.potential_area_increase / self.nodal_unit.plant.area.potential_leaf_increase
    #     da = self.actual_area_increase
    #     a = self.area - da
    #     if a > 0:
    #         da / a
    #     else:
    #         0

    stay_green_water_stress_duration(water_potential_effect, scale=0.5, threshold=-4.0) => begin
        # One day of cumulative severe water stress (i.e., water_effect = 0.0 around -4MPa) would result in a reduction of leaf lifespan in relation staygreeness and growthDuration, SK
        # if scale is 1.0, one day of severe water stress shortens one day of stayGreenDuration
        #TODO remove WaterStress and use general Accumulator with a lambda function?
        scale * (1 - water_potential_effect(threshold))
    end ~ accumulate(when=mature, u"d")

    stay_green_duration(SG, GD, stay_green_water_stress_duration) => begin
        # SK 8/20/10: as in Sinclair and Horie, 1989 Crop sciences, N availability index scaled between 0 and 1 based on
        #nitrogen_index = max(0, (2 / (1 + exp(-2.9 * (self.g_content - 0.25))) - 1))
        SG * GD - stay_green_water_stress_duration
    end ~ track(u"d", min=0)

    # Assumes physiological time for senescence is the same as that for growth though this may be adjusted by stayGreen trait
    # a peaked fn like beta fn not used here because aging should accelerate with increasing T not slowing down at very high T like growth,
    # instead a q10 fn normalized to be 1 at T_opt is used, this means above Top aging accelerates.
    active_age(pheno.Q10.ΔT) ~ accumulate(when=mature & !aging, u"d")

    senescence_water_stress_duration(water_potential_effect, scale=0.5, threshold=-4.0) => begin
        # if scale is 0.5, one day of severe water stress at predawn shortens one half day of agingDuration
        scale * (1 - water_potential_effect(threshold))
    end ~ accumulate(when=aging, u"d")

    senescence_duration(GD, senescence_water_stress_duration) => begin
        # end of growth period, time to maturity
        GD - senescence_water_stress_duration
    end ~ track(u"d", min=0)

    #TODO active_age and senescence_age could share a tracker with separate intervals
    #TODO support clipping with @rate option or sub-decorator (i.e. @active_age.clip)
    senescence_age(pheno.Q10.ΔT) ~ accumulate(when=aging & !dead, u"d")

    #TODO confirm if it really means the senescence ratio, not rate
    senescence_ratio(maximum_aging_rate, senescence_age, length) => begin
        # for MAIZSIM
        # t = self.senescence_age
        # t_e = self.senescence_duration
        # if t >= t_e:
        #     1
        # else:
        #     t_m = t_e / 2
        #     r = (1 + (t_e - t) / (t_e - t_m)) * (t / t_e)**(t_e / (t_e - t_m))
        #     clip(r, 0., 1.)
        # for Radish
        if iszero(length)
            0
        else
            maximum_aging_rate * senescence_age / length
        end
    end ~ track(min=0, max=1)

    senescent_area(senescence_ratio, area) => begin
        # Leaf senescence accelerates with drought and heat. see http://www.agry.purdue.edu/ext/corn/news/timeless/TopLeafDeath.html
        # rate = self._growth_rate(self.senescence_age, self.senescence_duration)
        # rate * self.timestep * self.area
        senescence_ratio * area
    end ~ track(u"cm^2")

    SLA(area, mass): specific_leaf_area => begin
        # temporary for now - it should vary by age. Value comes from some of Soo's work
        #200.0
        area / mass
    end ~ track(u"cm^2/g")
    LMA(SLA): leaf_mass_per_area => 1 / SLA ~ track(u"g/cm^2")

    # Maturity

    #HACK: tracking should happen after plant emergence (due to implementation of original beginFromEmergence)
    maturity(pheno.GD.r) ~ accumulate(when=pheno.emerged & !mature, u"K")

    # Nitrogen

    #FIXME avoid incorrect cycle detection (nitrogen member vs. module) - ?
    N: nitrogen => begin
        #TODO is this default value needed?
        # no N stress
        3.0
        #TODO remove self.p.* referencing
        #FIXME enable Nitrogen trait
        #self.p.nitrogen.leaf_content
    end ~ track

    ##########
    # States #
    ##########

    initiated => begin
        # no explicit initialize() here
        true
    end ~ flag

    appeared(rank, l=pheno.leaves_appeared) => (rank <= l) ~ flag

    growing(appeared, mature) => (appeared && !mature) ~ flag

    mature(elongation_age, GD, area, potential_area) => begin
        elongation_age >= GD || area >= potential_area
    end ~ flag

    aging(mature, physiological_age, SG, maturity) => begin
        # for MAIZSIM
        #active_age >= stay_green_duration
        # for Radish
        mature && physiological_age > SG * maturity
    end ~ flag

    dead(senescence_ratio) => begin
        #senescent_area >= area
        senescence_ratio >= 1
        #senescence_age >= senescence_duration?
    end ~ flag

    dropped(mature, dead) => (mature && dead) ~ flag
end

ext: jl


text:
## "Abstract" definitions

"""
    AbstractFPGroup

An Abstract type representing finitely presented groups. Every instance `` must implement
 * `KnuthBendix.alphabet(G::MyFPGroup)`
 * `rewriting(G::MyFPGroup)` : return the rewriting object which must implement
 > `KnuthBendix.rewrite_from_left!(u, v, rewriting(G))`.
By default `alphabet(G)` is returned, which amounts to free rewriting in `G`.
 * `relations(G::MyFPGroup)` : return a set of defining relations.

AbstractFPGroup may also override `word_type(::Type{MyFPGroup}) = Word{UInt16}`,
which controls the word type used for group elements. If a group has more than `255` generators you need to define e.g.
> `word_type(::Type{MyFPGroup}) = Word{UInt16}`
"""
abstract type AbstractFPGroup <: GroupsCore.Group end

word_type(G::AbstractFPGroup) = word_type(typeof(G))
# the default:
word_type(::Type{<:AbstractFPGroup}) = Word{UInt8}

"""
    rewriting(G::AbstractFPGroup)
Return a "rewriting object" for elements of `G`. The rewriting object must must implement
    KnuthBendix.rewrite_from_left!(
        u::AbstractWord,
        v::AbstractWord,
        rewriting(G)
    )

For example if `G` is a `FreeGroup` then `alphabet(G)` is returned which results in free rewriting. For `FPGroup` a rewriting system is returned which may (or may not) rewrite word `v` to its normal form.
"""
function rewriting end

Base.@propagate_inbounds function (G::AbstractFPGroup)(
    word::AbstractVector{<:Integer},
)
    @boundscheck @assert all(
        l -> 1 <= l <= length(KnuthBendix.alphabet(G)),
        word,
    )
    return FPGroupElement(word_type(G)(word), G)
end

## Group Interface

Base.one(G::AbstractFPGroup) = FPGroupElement(one(word_type(G)), G)

Base.eltype(::Type{FPG}) where {FPG<:AbstractFPGroup} =
    FPGroupElement{FPG,word_type(FPG)}

include("iteration.jl")

GroupsCore.ngens(G::AbstractFPGroup) = length(G.gens)

function GroupsCore.gens(G::AbstractFPGroup, i::Integer)
    @boundscheck 1 <= i <= GroupsCore.ngens(G)
    l = alphabet(G)[G.gens[i]]
    return FPGroupElement(word_type(G)([l]), G)
end
GroupsCore.gens(G::AbstractFPGroup) =
    [gens(G, i) for i in 1:GroupsCore.ngens(G)]

# TODO: ProductReplacementAlgorithm
function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{<:AbstractFPGroup},
)
    l = rand(10:100)
    G = rs[]
    nletters = length(alphabet(G))
    return FPGroupElement(word_type(G)(rand(1:nletters, l)), G)
end

Base.isfinite(::AbstractFPGroup) = (
    @warn "using generic isfinite(::AbstractFPGroup): the returned `false` might be wrong"; false
)

## FPGroupElement

abstract type AbstractFPGroupElement{Gr} <: GroupElement end

mutable struct FPGroupElement{Gr<:AbstractFPGroup,W<:AbstractWord} <:
               AbstractFPGroupElement{Gr}
    word::W
    savedhash::UInt
    parent::Gr

    FPGroupElement(
        word::W,
        G::AbstractFPGroup,
        hash::UInt = UInt(0),
    ) where {W<:AbstractWord} = new{typeof(G),W}(word, hash, G)

    FPGroupElement{Gr,W}(word::AbstractWord, G::Gr) where {Gr,W} =
        new{Gr,W}(word, UInt(0), G)
end

Base.show(io::IO, ::Type{<:FPGroupElement{Gr}}) where {Gr} =
    print(io, FPGroupElement, "{$Gr, …}")

word(f::AbstractFPGroupElement) = f.word

#convenience
KnuthBendix.alphabet(g::AbstractFPGroupElement) = alphabet(parent(g))

function Base.show(io::IO, f::AbstractFPGroupElement)
    f = normalform!(f)
    return KnuthBendix.print_repr(io, word(f), alphabet(f))
end

## GroupElement Interface for FPGroupElement

Base.parent(f::AbstractFPGroupElement) = f.parent

function Base.:(==)(g::AbstractFPGroupElement, h::AbstractFPGroupElement)
    @boundscheck @assert parent(g) === parent(h)
    normalform!(g)
    normalform!(h)
    hash(g) != hash(h) && return false
    return equality_data(g) == equality_data(h)
end

function Base.deepcopy_internal(g::FPGroupElement, stackdict::IdDict)
    return FPGroupElement(copy(word(g)), parent(g), g.savedhash)
end

function Base.inv(g::GEl) where {GEl<:AbstractFPGroupElement}
    G = parent(g)
    return GEl(inv(alphabet(G), word(g)), G)
end

function Base.:(*)(g::GEl, h::GEl) where {GEl<:AbstractFPGroupElement}
    @boundscheck @assert parent(g) === parent(h)
    return GEl(word(g) * word(h), parent(g))
end

GroupsCore.isfiniteorder(g::AbstractFPGroupElement) =
    isone(g) ? true :
    (
        @warn "using generic isfiniteorder(::AbstractFPGroupElement): the returned `false` might be wrong"; false
    )

# additional methods:
Base.isone(g::AbstractFPGroupElement) = (normalform!(g); isempty(word(g)))

## Free Groups

struct FreeGroup{T,O} <: AbstractFPGroup
    gens::Vector{T}
    ordering::O

    function FreeGroup(gens, ordering::KnuthBendix.WordOrdering)
        @assert length(gens) == length(unique(gens))
        L = KnuthBendix.letters(alphabet(ordering))
        @assert all(l -> l in L, gens)
        return new{eltype(gens),typeof(ordering)}(gens, ordering)
    end
end

FreeGroup(gens, A::Alphabet) = FreeGroup(gens, KnuthBendix.LenLex(A))

function FreeGroup(A::Alphabet)
    @boundscheck @assert all(
        KnuthBendix.hasinverse(l, A) for l in KnuthBendix.letters(A)
    )
    ltrs = KnuthBendix.letters(A)
    gens = Vector{eltype(ltrs)}()
    invs = Vector{eltype(ltrs)}()
    for l in ltrs
        l ∈ invs && continue
        push!(gens, l)
        push!(invs, inv(A, l))
    end

    return FreeGroup(gens, A)
end

function FreeGroup(n::Integer)
    symbols = Symbol[]
    inverses = Int[]
    sizehint!(symbols, 2n)
    sizehint!(inverses, 2n)
    for i in 1:n
        push!(symbols, Symbol(:f, i), Symbol(:F, i))
        push!(inverses, 2i, 2i - 1)
    end
    return FreeGroup(symbols[1:2:2n], Alphabet(symbols, inverses))
end

Base.show(io::IO, F::FreeGroup) =
    print(io, "free group on $(ngens(F)) generators")

# mandatory methods:
relations(F::FreeGroup) = Pair{eltype(F)}[]
KnuthBendix.ordering(F::FreeGroup) = F.ordering
KnuthBendix.alphabet(F::FreeGroup) = alphabet(KnuthBendix.ordering(F))
rewriting(F::FreeGroup) = alphabet(F)

# GroupsCore interface:
# these are mathematically correct
Base.isfinite(::FreeGroup) = false

GroupsCore.isfiniteorder(g::AbstractFPGroupElement{<:FreeGroup}) =
    isone(g) ? true : false

## FP Groups

struct FPGroup{T,R,S} <: AbstractFPGroup
    gens::Vector{T}
    relations::Vector{Pair{S,S}}
    rws::R
end

relations(G::FPGroup) = G.relations
rewriting(G::FPGroup) = G.rws
KnuthBendix.ordering(G::FPGroup) = KnuthBendix.ordering(rewriting(G))
KnuthBendix.alphabet(G::FPGroup) = alphabet(KnuthBendix.ordering(G))

function FPGroup(
    G::AbstractFPGroup,
    rels::AbstractVector{<:Pair{GEl,GEl}};
    ordering = KnuthBendix.ordering(G),
    kwargs...,
) where {GEl<:FPGroupElement}
    for (lhs, rhs) in rels
        @assert parent(lhs) === parent(rhs) === G
    end
    word_rels = [word(lhs) => word(rhs) for (lhs, rhs) in [relations(G); rels]]
    rws = KnuthBendix.RewritingSystem(word_rels, ordering)

    KnuthBendix.knuthbendix!(rws; kwargs...)

    return FPGroup(G.gens, rels, rws)
end

function Base.show(io::IO, ::MIME"text/plain", G::FPGroup)
    print(io, "Finitely presented group generated by:\n\t{")
    Base.print_array(io, permutedims(gens(G)))
    println(io, " },")
    println(io, "subject to relations:")
    return Base.print_array(io, relations(G))
end

function Base.show(io::IO, G::FPGroup)
    print(io, "⟨")
    Base.print_array(io, permutedims(gens(G)))
    println(io, " | ")
    print(io, "\t ")
    Base.print_array(io, permutedims(relations(G)))
    return print(io, " ⟩")
end

Base.show(io::IO, ::Type{<:FPGroup{T}}) where {T} =
    print(io, FPGroup, "{$T, …}")

## GSymbol aka letter of alphabet

abstract type GSymbol end
Base.literal_pow(::typeof(^), t::GSymbol, ::Val{-1}) = inv(t)

function subscriptify(n::Integer)
    subscript_0 = Int(0x2080) # Char(0x2080) -> subscript 0
    return join([Char(subscript_0 + i) for i in reverse(digits(n))], "")
end


text:
module ADVUMPS

using Zygote
using OMEinsum

export vumps, vumps_env
export Z,magnetisation, energy
export hamiltonian, model_tensor, mag_tensor
export Ising, TFIsing, Heisenberg
export init_ipeps, optimiseipeps

include("hamiltonianmodels.jl")
include("vumps.jl")
include("exampletensors.jl")
include("fixedpoint.jl")
include("autodiff.jl")
include("ipeps.jl")
include("variationalipeps.jl")

end


text:
module ModelClockSystem

    using CFTime
    using Dates

    mutable struct Alarm
        name :: String
        time :: AbstractCFDateTime
        callbacks :: Array{Function, 1}
    end

    mutable struct ModelClock

        time         :: AbstractCFDateTime
        alarms       :: Array{Alarm, 1}

        function ModelClock(time :: AbstractCFDateTime)
            alarms = Array{Alarm}(undef, 0)
            return new(
                time + Dates.Second(0),
                alarms,
            )
        end

    end

    function checkType(time1 :: AbstractCFDateTime, time2 :: AbstractCFDateTime)
        if typeof(time1) != timeof(time2)
            throw(ErrorException("Time type does not match!"))
        end
    end

    function setClock!(
        mc :: ModelClock,
        time :: AbstractCFDateTime,
    )

        checkType(mc.time, time)
        mc.time = time + Second(0)

    end

    function advanceClock!(
        mc :: ModelClock,
        t :: Union{Second, AbstractCFDateTime},
    )

        if typeof(t) <: Second  # inteprete as time interval
            mc.time += t
        elseif typeof(Δt) <: AbstractCFDateTime
            setClock!(mc, t)
        end

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


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

using Markdown
using InteractiveUtils

# ╔═╡ f3a59724-187d-11eb-2f7c-2df2f3339a6e
using Pkg, DrWatson

# ╔═╡ f3a5da02-187d-11eb-24ed-a7c460c0219a
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ f3b3b5aa-187c-11eb-2e75-07bc5598cd16
md"## Turing intro clip 3_quap.jl"

# ╔═╡ f3a66514-187d-11eb-22eb-47f74d2724cd
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
	x̄ = mean(df.weight)
	x = range(minimum(df.weight), maximum(df.weight), length = 100)     # either
end;

# ╔═╡ 2f0ce4de-187e-11eb-296e-8b4e16814350
@model function ppl4_3(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* weights
    for i in eachindex(heights)
        heights[i] ~ Normal(μ[i], σ)
    end
end

# ╔═╡ 33b27276-187e-11eb-26a8-1bb7ad288539
begin
	m4_3t = ppl4_3(df.weight, df.height)
    sampler=NUTS(0.65)
	nsamples=2000
	nchains=4
    chns4_3t = mapreduce(c -> sample(m4_3t, sampler, nsamples), chainscat, 1:nchains)
	Particles(Chains(chns4_3t, [:parameters]))
end

# ╔═╡ 7779262e-187e-11eb-0441-37f564373c1e
q4_3t = quap(m4_3t)

# ╔═╡ 32db0b98-2c08-11eb-145d-a1ea44641e75
rand(q4_3t.distr, 4_000)'

# ╔═╡ 3efadbfe-1882-11eb-0900-b17b00bf8f31

begin
	scatter(df.weight, df.height, leg=false)
	quap4_3t = sample(q4_3t)
	a_map = mean(quap4_3t.a)
	b_map = mean(quap4_3t.b)
	plot!(x, a_map .+ b_map .* x)
end

# ╔═╡ Cell order:
# ╟─f3b3b5aa-187c-11eb-2e75-07bc5598cd16
# ╠═f3a59724-187d-11eb-2f7c-2df2f3339a6e
# ╠═f3a5da02-187d-11eb-24ed-a7c460c0219a
# ╠═f3a66514-187d-11eb-22eb-47f74d2724cd
# ╠═2f0ce4de-187e-11eb-296e-8b4e16814350
# ╠═33b27276-187e-11eb-26a8-1bb7ad288539
# ╠═7779262e-187e-11eb-0441-37f564373c1e
# ╠═32db0b98-2c08-11eb-145d-a1ea44641e75
# ╠═3efadbfe-1882-11eb-0900-b17b00bf8f31


text:
# License is MIT: https://github.com/JuliaString/LaTeX_Entities/LICENSE.md
#
# Portions of this are based on code from julia/base/latex_symbols.jl
#
# Mapping from LaTeX math symbol to the corresponding Unicode codepoint.
# This is used for tab substitution in the REPL.

println("Running LaTeX build in ", pwd())

using LightXML
using StrTables

VER = UInt32(1)

#const dpath = "http://www.w3.org/Math/characters/"
const dpath = "http://www.w3.org/2003/entities/2007xml/"
const fname = "unicode.xml"
#const lpath = "http://mirror.math.ku.edu/tex-archive/macros/latex/contrib/unicode-math/"
const lpath = "https://raw.githubusercontent.com/wspr/unicode-math/master/"
const lname = "unicode-math-table.tex"

const disp = [false]

# Get manual additions to the tables
include("../src/manual_latex.jl")

const datapath = "../data"

const empty_str = ""
const element_types = ("mathlatex", "AMS", "IEEE", "latex")

function get_math_symbols(dpath, fname)
    lname = joinpath(datapath, fname)
    if isfile(fname)
        println("Loaded: ", lname)
        vers = lname
    else
        vers = string(dpath, fname)
        download(vers, lname)
        println("Saved to: ", lname)
    end
    xdoc = parse_file(lname)

    latex_sym  = [Pair{String, String}[] for i = 1:length(element_types)]

    info = Tuple{Int, Int, String, String, String}[]
    count = 0
    # Handle differences in versions of unicode.xml document
    rt = root(xdoc)
    top = find_element(rt, "charlist")
    top == nothing || (rt = top)
    for c in child_nodes(rt)
        if name(c) == "character" && is_elementnode(c)
            ce = XMLElement(c)
            for (ind, el) in enumerate(element_types)
                latex = find_element(ce, el)
                if latex == nothing
                    disp[] && println("##\t", attribute(ce, "id"), "\t", ce)
                    continue
                end
                L = strip(content(latex))
                id = attribute(ce, "id")
                U = string(map(s -> Char(parse_hex(UInt32, s)), split(id[2:end], "-"))...)
                mtch = _contains(L, r"^\\[A-Za-z][A-Za-z0-9]*(\{[A-Za-z0-9]\})?$")
                disp[] &&
                    println("#", count += 1, "\t", mtch%Int, " id: ", id, "\tU: ", U, "\t", L)
                if mtch
                    L = L[2:end] # remove initial \
                    if length(U) == 1 && isascii(U[1])
                        # Don't store LaTeX names for ASCII characters
                        typ = 0
                    else
                        typ = 1
                        push!(latex_sym[ind], String(L) => U)
                    end
                    push!(info, (ind, typ, L, U, empty_str))
                end
            end
        end
    end
    latex_sym, vers, info
end

function add_math_symbols(dpath, fname)
    lname = joinpath(datapath, fname)
    if isfile(fname)
        println("Loaded: ", lname)
        vers = lname
    else
        vers = string(dpath, fname)
        download(vers, lname)
        println("Saved to: ", lname)
    end
    latex_sym = Pair{String, String}[]
    info = Tuple{Int, Int, String, String, String}[]
    open(lname) do f
        for L in eachline(f)
            (isempty(L) || L[1] == '%') && continue
            x = map(s -> rstrip(s, [' ','\t','\n']),
                    split(_replace(L, r"[{}\"]+" => "\t"), "\t"))
            ch = Char(parse_hex(UInt32, x[2]))
            nam = String(x[3][2:end])
            startswith(nam, "math") && (nam = nam[5:end])
            if isascii(ch)
                typ = 0 # ASCII
            elseif Base.is_id_char(ch)
                typ = 1 # identifier
            elseif Base.isoperator(Symbol(ch))
                typ = 2 # operator
            else
                typ = 3
            end
            typ != 0 && push!(latex_sym, nam => string(ch))
            push!(info, (2, typ, nam, string(ch), x[5]))
        end
    end
    latex_sym, vers, info
end

#=

 standard | v7.0    | new   | type
----------|---------|-------|-------------------------------
mscr	  | scr	    | c_    | script/cursive
msans	  | sans    | s_    | sans-serif
Bbb       | bb	    | d_    | blackboard / doublestruck
mfrak     | frak    | f_    | fraktur
mtt	  | tt	    | t_    | mono
mit	  | it	    | i_    | italic
mitsans   | isans   | is_   | italic sans-serif
mitBbb    | bbi	    | id_   | italic blackboard / doublestruct
mbf	  | bf	    | b_    | bold
mbfscr	  | bscr    | bc_   | bold script/cursive
mbfsans   | bsans   | bs_   | bold sans-serif
mbffrak   | bfrak   | bf_   | bold fraktur
mbfit	  | bi	    | bi_   | bold italic
mbfitsans | bisans  | bis_  | bold italic sans-serif
<greek>             | G     | greek
it<greek>           | i_G   | italic greek
bf<greek>           | b_G   | bold greek
bi<greek>	    | bi_G  | bold italic greek
bsans<greek>	    | bs_G  | bold sans-serif greek
bisans<greek>       | bis_G | bold italic sans-serif greek
var<greek>          | V     | greek variant
mitvar<greek>       | i_V   | italic greek variant
mbfvar<greek>       | b_V   | bold greek variant
mbfitvar<greek>	    | bi_V  | bold italic greek variant
mbfsansvar<greek>   | bs_V  | bold sans-serif greek variant
mbfitsansvar<greek> | bis_V | bold italic sans-serif greek variant

i -> imath                     ı
=#
function str_chr(val)
    isempty(val) && return ""
    io = IOBuffer()
    for ch in val
        print(io, hex(ch%UInt32,4), ':')
    end
    String(take!(io))[1:end-1]
end

function str_names(nameset)
    io = IOBuffer()
    allnames = sort(collect(nameset))
    for n in allnames
        print(io, n, " ")
    end
    String(take!(io))
end

function add_name(dic::Dict, val, nam)
    if haskey(dic, val)
        push!(dic[val], nam)
        disp[] && println("\e[s$val\e[u\e[4C$(rpad(str_chr(val),20))", str_names(dic[val]))
    else
        dic[val] = Set((nam,))
        disp[] && println("\e[s$val\e[u\e[4C$(rpad(str_chr(val),20))$nam")
    end
end

function check_name(out::Dict, dic::Dict, val, nam, old)
    oldval = get(dic, nam, "")
    # Check if short name is already in table with same value
    oldval == "" && return (add_name(out, val, nam); true)
    oldval != val && disp[] && println("Conflict: $old => $val, $nam => $oldval")
    false
end

function replace_suffix(out, dic, val, nam, suffix, pref, list)
    for (suf, rep) in list
        suffix == suf && return check_name(out, dic, val, pref * rep, nam)
    end
    false
end

#=
function replace_greek(out, dic, val, nam, off, pref, list)
    for (suf, rep) in list
        if nam[off:end] == suf
            return check_name(out, dic, val, pref * rep, nam) |
                   check_name(out, dic, val, pref[1:end-1] * suf, nam)
        end
    end
    false
end
=#

function replace_all(out, dic, val, nam, suffix, pref)
#    replace_greek(out, dic, val, nam, off, pref * "G_", greek_letters) ||
#    replace_greek(out, dic, val, nam, off, pref * "V_", var_greek) ||
    replace_suffix(out, dic, val, nam, suffix, pref * "_", digits)
end

function shorten_names(names::Dict)
    valtonam = Dict{String,Set{String}}()
    for (nam, val) in names
        # handle combining accents, change from 'accent{X}' to 'X-accent'
        if !startswith(nam, "math") && sizeof(nam) > 3 &&
            nam[end]%UInt8 == '}'%UInt8 && nam[end-2]%UInt8 == '{'%UInt8
            ch = nam[end-1]%UInt8
            if ch - 'A'%UInt8 < 0x1a || ch - 'a'%UInt8 < 0x1a || ch - '0'%UInt8 < 0xa
                # tst = string(nam[end-1], '-', nam[1:end-3])
                # check_name(valtonam, names, val, tst, nam)
                add_name(valtonam, val, nam)
                continue
            end
        end
        # Special handling of "up"/"mup" prefixes
        if startswith(nam, "up")
            # Add it later when processing "mup" prefix if they have the same value
            get(names, "m" * nam, "") == val && continue
        elseif startswith(nam, "mup")
            # If short form in table with same value, continue, otherwise, add short form
            #upval = get(names, nam[2:end], "") # see if "up..." is in the table
            #val == upval && continue # short name is already in table with same value
            oldval = get(names, nam[4:end], "") # see if "..." is in the table
            val == oldval && continue # short name is already in table with same value
            check_name(valtonam, names, val, oldval == "" ? nam[4:end] : nam[2:end], nam)
            continue
        else
            flg = false
            nam in remove_name && continue
            for (oldnam, newnam) in replace_name
                if nam == oldnam
                    flg = check_name(valtonam, names, val, newnam, nam)
                    break
                end
            end
            flg && continue
            if nam[1] in remove_lead_char
                for pref in remove_prefix
                    startswith(nam, pref) || continue
                    flg = true
                    tst = nam[sizeof(pref)+1:end]
                    oldval = get(names, tst, "")
                    oldval == val || (oldval == "" && add_name(valtonam, val, tst))
                    break
                end
            elseif nam[1] in replace_lead_char
                for (pref, rep, repv7) in replace_prefix
                    startswith(nam, pref) || continue
                    suff = nam[sizeof(pref)+1:end]
                    flg = replace_all(valtonam, names, val, nam, suff, rep)
                    #=
                    if rep == "i"
                        flg = replace_all(valtonam, names, val, suff, "i")
                    elseif rep == "t" || rep == "d" || rep == "s" || rep == "c"
                        flg = replace_all(valtonam, names, val, suff, rep)
                    elseif rep == "" || rep[1] != 'b'
                    elseif rep == "b"
                        flg = replace_all(valtonam, names, val, suff, "b")
                    elseif rep == "sb"
                        flg = replace_all(valtonam, names, val, suff, "bs")
                    elseif rep == "ib"
                        flg = replace_all(valtonam, names, val, suff, "bi")
                    elseif rep == "cib"
                        flg = replace_all(valtonam, names, val, suff, "bic")
                    end
                    =#
                    flg || (flg = check_name(valtonam, names, val, rep * "_" * suff, nam))
                    #check_name(valtonam, names, val, repv7 * nam[sizeof(pref)+1:end], nam)
                    #flg = true
                    break
                end
            end
            # Add short forms, if not already handled
            flg && continue
        end
#            replace_suffix(valtonam, names, val, nam, 1, "G_", greek_letters) ||
#                replace_suffix(valtonam, names, val, nam, 1, "V_", var_greek)
        add_name(valtonam, val, nam)
    end
    # Split into two vectors
    syms = Vector{String}()
    vals = Vector{String}()
    for (val, namset) in valtonam
        for nam in namset
            startswith(nam, "math") && length(namset) > 1 && continue
            push!(syms, nam)
            push!(vals, val)
        end
    end
    syms, vals
end

function make_tables()
    sym1, ver1, inf1 = get_math_symbols(dpath, fname)
    sym2, ver2, inf2 = add_math_symbols(lpath, lname)

    latex_sym = [mansym..., sym1[1], sym2, sym1[2:end]...]
    et = (mantyp..., element_types[1], "tex", element_types[2:end]...)

    latex_set = Dict{String,String}()
    diff_set = Dict{String,Set{String}}()

    # Select the first name found, ignore duplicates
    for (ind, sym_set) in enumerate(latex_sym)
        countdup = 0
        countdiff = 0
        for (nam, val) in sym_set
            old = get(latex_set, nam, "")
            if old == ""
                push!(latex_set, nam => val)
            elseif val == old
                countdup += 1
            else
                countdiff += 1
                if haskey(diff_set, nam)
                    push!(diff_set[nam], val)
                else
                    push!(diff_set, nam => Set([old, val]))
                end
            end
        end
        println(countdup, " duplicates, ", countdiff, " overwritten out of ", length(sym_set),
                " found in ", et[ind])
    end
    # Dump out set
    disp[] && println("LaTeX set:\n", latex_set)
    disp[] && println("Differences:\n", diff_set)

    # Now, replace or remove prefixes and suffixes
    symnam, symval = shorten_names(latex_set)

    disp[] && println(length(symval), " distinct entities found\n", symnam)

    # We want to build a table of all the names, sort them, then create a StrTable out of them
    srtnam = sortperm(symnam)
    srtval = symval[srtnam] # Values, sorted the same as srtnam

    # BMP characters
    l16 = Tuple{UInt16, UInt16}[]
    # non-BMP characters (in range 0x10000 - 0x1ffff)
    l32 = Tuple{UInt16, UInt16}[]
    # two characters packed into UInt32, first character in high 16-bits
    l2c = Tuple{UInt32, UInt16}[]

    for i in eachindex(srtnam)
        chrs = srtval[i]
        len = length(chrs)
        len > 2 && error("Too long sequence of characters $chrs")
        ch1 = chrs[1]%UInt32
        if len == 2
            ch2 = chrs[end]%UInt32
            (ch1 > 0x0ffff || ch2 > 0x0ffff) &&
                error("Character $ch1 or $ch2 > 0xffff")
            push!(l2c, (ch1<<16 | ch2, i))
        elseif ch1 > 0x1ffff
            error("Character $ch1 too large")
        elseif ch1 > 0x0ffff
            push!(l32, (ch1-0x10000, i))
        else
            push!(l16, (ch1%UInt16, i))
        end
    end

    # We now have 3 vectors, for single BMP characters, for non-BMP characters, and for 2 BMP chars
    # each has the value and a index into the name table
    # We need to create a vector the same size as the name table, that gives the index
    # of into one of the three tables, in order to go from names to 1 or 2 output characters
    # We also need, for each of the 3 tables, a sorted vector that goes from the indices
    # in each table to the index into the name table (so that we can find multiple names for
    # each character)

    indvec = create_vector(UInt16, length(srtnam))
    vec16, ind16, base32 = sortsplit!(indvec, l16, 0)
    vec32, ind32, base2c = sortsplit!(indvec, l32, base32)
    vec2c, ind2c, basefn = sortsplit!(indvec, l2c, base2c)

    ((VER, string(now()), string(ver1, ",", ver2),
      base32%UInt32, base2c%UInt32, StrTable(symnam[srtnam]), indvec,
      vec16, ind16, vec32, ind32, vec2c, ind2c),
     (ver1, ver2), (inf1, inf2))
end

savfile = joinpath(datapath, "latex.dat")
if isfile(savfile)
    println("Tables already exist")
else
    tup = nothing
    println("Creating tables")
    try
        global tup
        tup = make_tables()
    catch ex
        println(sprint(showerror, ex, catch_backtrace()))
    end
    println("Saving tables to ", savfile)
    StrTables.save(savfile, tup[1])
    println("Done")
end


text:
"""
jFoma.jl

(C) 2018 by Damir Cavar (http://damir.cavar.me/)

An example implementation of Foma access from Julia.

This example code will load a Foma Finite State Tranducer file, and return the analysis strings in
a Vector of Strings.

The assumptions are:

- you have Foma installed and Julia can find the dynamic library in the system-path.

For Foma, see the GitHub page:
https://fomafst.github.io
"""

const LIBNAME = "libfoma"
const MORPHOLOGY = "english.fst"

"""
Loading the network from a file specified in the variable MORPHOLOGY.
"""
function loadNet()
    ccall((:fsm_read_binary_file, LIBNAME), Ptr{Cvoid}, (Cstring,), MORPHOLOGY)
end

"""
Initializes the net and returns the init pointer.
"""
function initNet(net)
    ccall((:apply_init, LIBNAME), Ptr{Cvoid}, (Ptr{Cvoid},), net)
end

"""
Clears the init-pointer in the Foma library.
"""
function clearNet(ah)
    ccall((:apply_clear, LIBNAME), Cvoid, (Ptr{Cvoid},), ah)
end

"""
Frees the memory for the allocated FST.
"""
function destroyNet(net)
    ccall((:fsm_destroy, LIBNAME), Cvoid, (Ptr{Cvoid},), net)
end

"""
Applies UP to the net with the given token string.
"""
function applyUp(ah, token)
    res = Vector{String}()
    r = ccall((:apply_up, LIBNAME), Cstring, (Ptr{Cvoid}, Cstring), ah, token)
    while true
        if r == C_NULL
            break
        end
        push!(res, unsafe_string(r))
        r = ccall((:apply_up, LIBNAME), Cstring, (Ptr{Cvoid}, Cstring), ah, Ptr{UInt8}(C_NULL))
    end
    return res
end

"""
main function.
"""
function main()

    net = loadNet()
    ah = initNet(net)

    if isempty(ARGS)
        # sending the token "calls" through the analyser.
        println("calls:")
        result = applyUp(ah, String("calls"))
        for r in result
            println(r)
        end
    else
        for arg in ARGS
            println("$arg:")
            result = applyUp(ah, arg)
            for r in result
                println(r)
            end
        end
    end

    clearNet(ah)
    destroyNet(net)
end

main()

max_stars_repo_name: dcavar/juliaFoma

max_issues_repo_name: dcavar/juliaFoma

max_forks_repo_name: dcavar/juliaFoma