text:
# differentiated version of sbp_sat_shock_capturing.jl

function calcShockCapturing_diff(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             capture::SBPParabolicSC{Tsol, Tres},
                             shockmesh::ShockedElements,
                             assem::AssembleElementData) where {Tsol, Tres}

  println("\nEntered calcShockCapturing_diff")
  computeGradW(mesh, sbp, eqn, opts, capture, shockmesh,
               capture.entropy_vars, capture.diffusion)

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

  computeFaceTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
                       capture.diffusion, capture.entropy_vars,
                       capture.penalty, assem)

  if shockmesh.isNeumann
    println("computing Neumann boundary condition")
    computeNeumannBoundaryTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
                              capture.diffusion, capture.entropy_vars, assem)
  else
    println("computing Dirichlet boundary condition")
    computeDirichletBoundaryTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
                                      assem)
  end

  #@time computeSharedFaceTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
  #                                 capture.diffusion, capture.penalty)

  return nothing
end

function computeVolumeTerm_diff(mesh, sbp, eqn, opts,
                           capture::SBPParabolicSC{Tsol, Tres},
                           diffusion::AbstractDiffusion,
                           entropy_vars::AbstractVariables,
                           shockmesh::ShockedElements,
                           assembler::AssembleElementData) 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)
  Dx = zeros(mesh.numNodesPerElement, mesh.numNodesPerElement, mesh.dim)
  w_dot = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerElement)
  t1_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode, mesh.dim,
                 mesh.numNodesPerElement, mesh.numNodesPerElement)
  t2_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode, mesh.dim,
                 mesh.numNodesPerElement, mesh.numNodesPerElement)
  gradq_i = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.dim)
  # if dim was the last dimension of t1_dot, we could re-use it, but instead
  # we have to use a 3rd array
  res_jac = eqn.params.calc_volume_integrals_data.res_jac
  op = SummationByParts.Subtract()

  for i=1:shockmesh.numShock
    i_full = shockmesh.elnums_all[i]

    q_i = ro_sview(eqn.q, :, :, i_full)
    w_i = ro_sview(capture.w_el, :, :, i)
    coords_i = ro_sview(mesh.coords, :, :, i_full)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i_full)
    jac_i = ro_sview(mesh.jac, :, i_full)

    # to compute the jacobian, dq/dq = I, so dw/dq = the dw/dq evaluated at
    # each node
    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(eqn.q, :, j, i_full)
      w_dot_j = sview(w_dot, :, :, j)
      getA0inv(entropy_vars, eqn.params, q_j, w_dot_j)
    end

    # apply Dx * w_dot
    calcDx(sbp, dxidx_i, jac_i, Dx)
    applyOperatorJac(Dx, w_dot, t1_dot)

    # need Dx * w for applyDiffusionTensor
    for d=1:mesh.dim
      Dx_d = ro_sview(Dx, :, :, d)
      gradq_d = sview(gradq_i, :, :, d)
      smallmatmatT!(w_i, Dx_d, gradq_d)
    end

    # apply the diffusion tensor to all nodes of the element
    applyDiffusionTensor_diff(diffusion, sbp, eqn.params, q_i, w_i, coords_i,
                              dxidx_i, jac_i, i, gradq_i, t1_dot, t2_dot)

    # apply Qx and sum
    if eqn.params.use_Minv != 1
      @simd for d=1:mesh.dim
        @simd for j=1:mesh.numNodesPerElement
          @simd for i=1:mesh.numNodesPerElement
            Dx[i, j, d] *= sbp.w[i]/jac_i[i]  # Dx is now Qx
          end
        end
      end
    end
    applyOperatorJac(Dx, t2_dot, res_jac)

#    calcQxTransposed(sbp, dxidx_i, Dx)
#    applyOperatorJac(Dx, t2_dot, res_jac, true, op)

    assembleElement(assembler, mesh, i_full, res_jac)
  end

  return nothing
end

function computeFaceTerm_diff(mesh, sbp, eqn, opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements, diffusion::AbstractDiffusion,
                      entropy_vars::AbstractVariables,
                      penalty::AbstractDiffusionPenalty,
                      assem::AssembleElementData) where {Tsol, Tres}

  delta_w = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  delta_w_dotL = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode,
                             mesh.numNodesPerFace, mesh.numNodesPerElement)
  delta_w_dotR = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode,
                             mesh.numNodesPerFace, mesh.numNodesPerElement)

  theta = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  theta_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  theta_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res1L_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res1L_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res1R_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res1R_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res2L_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res2L_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res2R_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res2R_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  data = eqn.params.calc_face_integrals_data
  @unpack data res_jacLL res_jacLR res_jacRL res_jacRR

  #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)
  subtract = SummationByParts.Subtract()
  add = SummationByParts.Add()

  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
    wL = ro_sview(capture.w_el, :, :, iface_red.elementL)
    wR = ro_sview(capture.w_el, :, :, iface_red.elementR)
    qL = ro_sview(eqn.q, :, :, elnumL)
    qR = ro_sview(eqn.q, :, :, elnumR)
    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_diff(eqn.params, capture, diffusion, entropy_vars,
                          sbp, mesh.sbpface,
                          iface_red, wL, wR, qL, qR, gradwL, gradwR,
                          coordsL, coordsR,
                          nrm_face, dxidxL, dxidxR, jacL, jacR,
                          delta_w, delta_w_dotL, delta_w_dotR,
                          theta, theta_dotL, theta_dotR)

    # apply the penalty coefficient matrix
    has_T4 = applyPenalty_diff(penalty, sbp, eqn.params,  mesh.sbpface,
                      diffusion,
                      iface_red, delta_w, delta_w_dotL, delta_w_dotR,
                      theta, theta_dotL, theta_dotR,
                      qL, qR, wL, wR, coordsL, coordsR,  nrm_face, alphas,
                      dxidxL, dxidxR, jacL, jacR,
                      res1L_dotL, res1L_dotR, res1R_dotL, res1R_dotR,
                      t2L, t2R,
                      res2L_dotL, res2L_dotR, res2R_dotL, res2R_dotR)

    # apply Rgk^T, Rgn^T, Dgk^T, Dgn^T
    fill!(res_jacLL, 0); fill!(res_jacRL, 0)
    fill!(res_jacLR, 0); fill!(res_jacRR, 0)
    interiorFaceIntegrate_jac!(mesh.sbpface, iface_red, res1L_dotL, res1R_dotL,
                               res_jacLL, res_jacRL, subtract, subtract, false)
    interiorFaceIntegrate_jac!(mesh.sbpface, iface_red, res1L_dotR, res1R_dotR,
                               res_jacLR, res_jacRR, subtract, subtract, false)

    # apply Dgk^T and Dgn^T
    applyDgkTranspose_diff(capture, sbp, eqn.params, mesh.sbpface, iface_red,
                      diffusion,
                      t2L, res2L_dotL, res2L_dotR,
                      t2R, res2R_dotL, res2R_dotR,
                      qL, qR, wL, wR, coordsL, coordsR,
                      nrm_face, dxidxL, dxidxR, jacL, jacR,
                      res_jacLL, res_jacLR, res_jacRL, res_jacRR,
                      subtract)

    iface_full = replace_interface(iface_red, elnumL, elnumR)
    if !has_T4
      assembleInterfaceVisc(assem, mesh.sbpface, mesh, iface_full, res_jacLL,
                            res_jacLR, res_jacRL, res_jacRR)
    else
      #TODO: the sparsity pattern of the Jacobian may be incorrect for this case
      assembleInterfaceFull(assem, mesh, iface_full, res_jacLL,
                            res_jacLR, res_jacRL, res_jacRR)
    end

  end  # end loop i

  return nothing
end

function computeNeumannBoundaryTerm_diff(mesh::AbstractMesh{Tmsh}, sbp, eqn,
                      opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements,
                      diffusion::AbstractDiffusion,
                      entropy_vars::AbstractVariables,
                      assem::AssembleElementData
                      ) where {Tsol, Tres, Tmsh}

  @assert eqn.params.use_Minv != 1
#  @assert mesh.coord_order == 1  # because of the normal vector

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

  op = SummationByParts.Subtract()

  @unpack capture w_dot Dx t1 t1_dot t2_dot t3_dot t4_dot

  res_jac = eqn.params.calc_face_integrals_data.res_jacLL
  nrm_i = zeros(Tmsh, mesh.dim, mesh.numNodesPerFace)

  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]

    q_i = ro_sview(eqn.q, :, :, elnum_orig)
    w_i = ro_sview(capture.w_el, :, :, bndry_i.element)
    coords_i = ro_sview(mesh.coords, :, :, elnum_orig)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, elnum_orig)
    jac_i = ro_sview(mesh.jac, :, elnum_orig)
    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")
      #nrm_i = ro_sview(mesh.nrm_face, :, :, idx_orig)
      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)
    # compute dw/dq
    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(eqn.q, :, j, elnum_orig)
      w_dot_j = sview(w_dot, :, :, j)
      getA0inv(entropy_vars, eqn.params, q_j, w_dot_j)
    end

    # multiply by Dgk
    calcDx(sbp, dxidx_i, jac_i, Dx)
    applyOperatorJac(Dx, mesh.sbpface, bndry_i.face, w_dot, t1_dot)

    # need t1 = [Dx, Dy] * w
    for d=1:mesh.dim
      Dx_d = ro_sview(Dx, :, :, d)
      t1_d = sview(t1, :, :, d)
      smallmatmatT!(w_i, Dx_d, t1_d)
    end

    applyDiffusionTensor_diff(diffusion, sbp, eqn.params, q_i, w_i, coords_i,
                              dxidx_i, jac_i, bndry_i.element, mesh.sbpface,
                              bndry_i.face, t1, t1_dot, t2_dot)

    # apply N * B * R and do the sum over dimensions
    fill!(t3_dot, 0)
    boundaryFaceInterpolate_jac!(mesh.sbpface, bndry_i.face, t2_dot, t3_dot)

    fill!(t4_dot, 0)  #TODO: tile this
    @simd for q=1:mesh.numNodesPerElement
      @simd for p=1:mesh.numNodesPerFace
        @simd for d=1:mesh.dim
          fac = nrm_i[d, p] #*mesh.sbpface.wface[p]
          @simd for j=1:mesh.numDofPerNode
            @simd for k=1:mesh.numDofPerNode
              t4_dot[k, j, p, q] += fac*t3_dot[k, j, d, p, q]
            end
          end
        end
      end
    end

    # apply R^T
    fill!(res_jac, 0)
    boundaryFaceIntegrate_jac!(mesh.sbpface, bndry_i.face, t4_dot, res_jac, op)
    bndry_orig = replace_boundary(bndry_i, elnum_orig)
    assembleBoundaryFull(assem, mesh, bndry_orig, res_jac)
  end  # end i

  return nothing
end

"""
  Differentiated Dirichlet boundary term
"""
function computeDirichletBoundaryTerm_diff(mesh, sbp, eqn, opts,
                      capture::SBPParabolicSC{Tsol, Tres},
                      shockmesh::ShockedElements,
                      assem::AssembleElementData
                      ) 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_diff(mesh, sbp, eqn, opts, shockmesh, capture,
                     capture.penalty, capture.diffusion,
                     capture.entropy_vars, bc_func, bc_range, assem)
  end

  return nothing
end

function calcBoundaryFlux_diff(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,
                      assem::AssembleElementData,
                      ) where {Tsol, Tres}

  delta_w = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  delta_w_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                            mesh.numNodesPerFace, mesh.numNodesPerElement)
  delta_w_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                             mesh.numNodesPerFace, mesh.numNodesPerElement)
  res1 = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  res1_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                         mesh.numNodesPerFace, mesh.numNodesPerElement)

  res_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                        mesh.numNodesPerElement, mesh.numNodesPerElement)
  res_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                        mesh.numNodesPerElement, mesh.numNodesPerElement)

  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_diff(capture, eqn.params, eqn.params_complex,
                          mesh.sbpface, bndry_i,
                          bc_func, entropy_vars, w_i, q_i, coords_i,
                          nrm_i, delta_w, delta_w_dot)

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

    # apply R^T to T_D * delta_w
    fill!(res_dot, 0)
    boundaryFaceIntegrate_jac!(mesh.sbpface, bndry_i.face, res1_dot, res_dot,
                               op, include_quadrature=false)

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

    for p=1:mesh.numNodesPerFace
      for i=1:mesh.numDofPerNode
        delta_w[i,p] *= -mesh.sbpface.wface[p]
      end
    end

    # the dotR variables aren't needed for boundary terms
    applyDgkTranspose_diff(capture, sbp, eqn.params, mesh.sbpface, bndry_i,
                      diffusion,
                      delta_w, delta_w_dot, delta_w_dotR, q_i, w_i, coords_i,
                      nrm_i, dxidxL, jacL, res_dot, res_dotR, op)

    assembleElement(assem, mesh, elnum_orig, res_dot)
  end

  return nothing
end

"""
  Differentiated version of [`getFaceVariables`](@ref).
"""
function getFaceVariables_diff(params::ParamType,
                          capture::SBPParabolicSC{Tsol, Tres},
                          diffusion::AbstractDiffusion,
                          entropy_vars::AbstractVariables,
                          sbp::AbstractOperator,
                          sbpface::AbstractFace, iface_red::Interface,
                          wL::AbstractMatrix, wR::AbstractMatrix,
                          qL::AbstractMatrix, qR::AbstractMatrix,
                          gradwL::Abstract3DArray, gradwR::Abstract3DArray,
                          coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                          nrm_face::AbstractMatrix,
                          dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                          jacL::AbstractVector, jacR::AbstractVector,
                          delta_w::AbstractMatrix,
                          delta_w_dotL::Abstract4DArray, delta_w_dotR::Abstract4DArray,
                          theta::AbstractArray,
                          theta_dotL::Abstract4DArray, theta_dotR::Abstract4DArray
                         ) where {Tsol, Tres}

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

  getFaceVariables(capture, sbpface, iface_red, wL, wR, gradwL, gradwR,
                   nrm_face, delta_w, theta)

  @unpack capture wL_dot wR_dot Dx t1 t1_dot t2L_dot t2R_dot t3L_dot t3R_dot

  # derivative of delta_w wrt qL and qR
  for i=1:numNodesPerElement
    wL_dot_i = sview(wL_dot, :, :, i)
    wR_dot_i = sview(wR_dot, :, :, i)
    qL_i = ro_sview(qL, :, i)
    qR_i = ro_sview(qR, :, i)
    getA0inv(entropy_vars, params, qL_i, wL_dot_i)
    getA0inv(entropy_vars, params, qR_i, wR_dot_i)
  end

  fill!(delta_w_dotL, 0); fill!(delta_w_dotR, 0)
  interiorFaceInterpolate_jac!(sbpface, iface_red, wL_dot, wR_dot,
                              delta_w_dotL, delta_w_dotR)

  @simd for i=1:length(delta_w_dotR)
    delta_w_dotR[i] = -delta_w_dotR[i]  # wR is negated in wL - wR
  end

  # derivative of theta
  # do volume operations for wL (derivative of wgradL)
  calcDx(sbp, dxidxL, jacL, Dx)
  applyOperatorJac(Dx, sbpface, iface_red.faceL, wL_dot, t1_dot)
  for d=1:dim
    Dx_d = ro_sview(Dx, :, :, d)
    t1_d = sview(t1, :, :, d)
    smallmatmatT!(wL, Dx_d, t1_d)
  end
  applyDiffusionTensor_diff(diffusion, sbp, params, qL, wL, coordsL, dxidxL,
                            jacL, iface_red.elementL, sbpface,
                            iface_red.faceL, t1, t1_dot, t2L_dot)

  # do volume operations for wR (derivative of wgradR)
  calcDx(sbp, dxidxR, jacR, Dx)
  applyOperatorJac(Dx, sbpface, iface_red.faceR, wR_dot, t1_dot)
  for d=1:dim
    Dx_d = ro_sview(Dx, :, :, d)
    t1_d = sview(t1, :, :, d)
    smallmatmatT!(wR, Dx_d, t1_d)
  end
  applyDiffusionTensor_diff(diffusion, sbp, params, qR, wR, coordsR, dxidxR,
                            jacR, iface_red.elementR, sbpface,
                            iface_red.faceR, t1, t1_dot, t2R_dot)

  # interpolate to face and multiply by normal vector
  #TODO: it would be interesting to supply a diagonal matrix to the 5D version
  #      of interiorFaceInterpolate_jac so it can do the reduction and therefore
  #      avoid a second temporary array
  fill!(t3L_dot, 0); fill!(t3R_dot, 0)
  interiorFaceInterpolate_jac!(sbpface, iface_red, t2L_dot, t2R_dot, t3L_dot, t3R_dot)
  fill!(theta_dotL, 0); fill!(theta_dotR, 0)  # TODO: tile this
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerFace
      @simd for d=1:dim
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            theta_dotL[i, j, p, q] += nrm_face[d, p]*t3L_dot[i, j, d, p, q]
            theta_dotR[i, j, p, q] -= nrm_face[d, p]*t3R_dot[i, j, d, p, q]
          end
        end
      end
    end
  end

  return nothing
end

function computeDirichletDelta_diff(capture::SBPParabolicSC{Tsol, Tres},
                               params::ParamType, params_c::ParamType,
                               sbpface::AbstractFace, bndry::Boundary,
                               func::BCType, entropy_vars::AbstractVariables,
                               wL::AbstractMatrix, qL::AbstractMatrix,
                               coords::AbstractMatrix, nrm_face::AbstractMatrix,
                               delta_w::AbstractMatrix,
                               delta_w_dot::Abstract4DArray
                              ) 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)

  # jacobian
  wL_dot = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerElement)
#  wface_dot = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerFace,
#                    numNodesPerElement)
  for i=1:numNodesPerElement
    wL_dot_i = sview(wL_dot, :, :, i)
    qL_i = ro_sview(qL, :, i)
    getA0inv(entropy_vars, params, qL_i, wL_dot_i)
  end
#  println("qL = \n", real(qL))
#  println("wL_dot = \n", real(wL_dot))
#  println("perm = ", sbpface.perm[:, bndry.face])
  fill!(delta_w_dot, 0)
  boundaryFaceInterpolate_jac!(sbpface, bndry.face, wL_dot, delta_w_dot)

  qface = zeros(Tsol, numDofPerNode, numNodesPerFace)
  qgface = zeros(Tres, numDofPerNode, numNodesPerFace)
  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)
    qgface_i = sview(qgface, :, i)  # need this later
    coords_i = sview(coords, :, i)
    nrm_i = sview(nrm_face, :, i)

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

    convertToEntropy(entropy_vars, params, qgface_i, wgface)

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

  # jacobian
  # Compute dqg/qface first, then R^T * dqg/qface to get dqg/dqvolume

  op = SummationByParts.Subtract()
  qg_dot = zeros(Tsol, numDofPerNode, numDofPerNode)
  A0inv = zeros(Tsol, numDofPerNode, numDofPerNode)
  wg_dotface = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerFace)
  q_c = zeros(Complex128, numDofPerNode)
  qg_c = zeros(Complex128, numDofPerNode)
  h = 1e-20
  pert = Complex128(0, h)
  for i=1:numNodesPerFace
    for j=1:numDofPerNode
      q_c[j] = qface[j, i]
    end
    coords_i = sview(coords, :, i)
    nrm_i = sview(nrm_face, :, i)

    # compute jacobian of qgface wrt qface
    for j=1:numDofPerNode
      q_c[j] += pert

      getDirichletState(func, params_c, q_c, aux_vars, coords_i, nrm_i, qg_c)
      for k=1:numDofPerNode
        qg_dot[k, j] = imag(qg_c[k])/h
      end

      q_c[j] -= pert
    end

    # jacobian of convertToEntropy
    getA0inv(entropy_vars, params, sview(qgface, :, i), A0inv)
    wg_dot_i = sview(wg_dotface, :, :, i)
    smallmatmat!(A0inv, qg_dot, wg_dot_i)
  end

  # apply R^T to wg_dotface to get d wg/dqL, and subtract that from delta_w_dot
  # interpolation is linear, so apply the regular operation repatedly
  # It would be better if SBP had a jacobian operation for this
  faceToVolume_jac!(sbpface, bndry.face, wg_dotface, delta_w_dot, op)

  return nothing
end

function applyDgkTranspose_diff(capture::SBPParabolicSC{Tsol, Tres}, sbp,
                       params::ParamType, sbpface, iface::Interface,
                       diffusion::AbstractDiffusion,
                       t2L::AbstractMatrix, t2L_dotL::Abstract4DArray,
                       t2L_dotR::Abstract4DArray,
                       t2R::AbstractMatrix, t2R_dotL::Abstract4DArray,
                       t2R_dotR::Abstract4DArray,
                       qL::AbstractMatrix, qR::AbstractMatrix,
                       wL::AbstractMatrix, wR::AbstractMatrix,
                       coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                       nrm_face::AbstractMatrix,
                       dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                       jacL::AbstractVector, jacR::AbstractVector,
                       resL_dotL::Abstract4DArray, resL_dotR::Abstract4DArray,
                       resR_dotL::Abstract4DArray, resR_dotR::Abstract4DArray,
                       op::SummationByParts.UnaryFunctor=SummationByParts.Add()
                      ) where {Tsol, Tres}

  dim, numNodesPerFace = size(nrm_face)
  numNodesPerElement = size(dxidxL, 3)
  numDofPerNode = size(wL, 1)
  add = SummationByParts.Add()
  subtract = SummationByParts.Subtract()

  @unpack capture temp1L temp1R temp2L temp2R
  fill!(temp2L, 0); fill!(temp2R, 0)
  @unpack capture t3L_dotL t3L_dotR t3R_dotL t3R_dotR Dx
  @unpack capture t4L_dotL t4L_dotR t4R_dotL t4R_dotR
  @unpack capture t5L_dotL t5L_dotR t5R_dotL t5R_dotR

  # 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 N
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerFace
      @simd for d=1:dim
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            # normally t3R would have a - sign because the normal vector
            # is reversed for element nu, but interiorFaceIntegrate -= the
            # second argument, so don't add the minus sign here
            t3L_dotL[i, j, d, p, q] = nrm_face[d, p]*t2L_dotL[i, j, p, q]
            t3L_dotR[i, j, d, p, q] = nrm_face[d, p]*t2L_dotR[i, j, p, q]
            t3R_dotL[i, j, d, p, q] = nrm_face[d, p]*t2R_dotL[i, j, p, q]
            t3R_dotR[i, j, d, p, q] = nrm_face[d, p]*t2R_dotR[i, j, p, q]
          end
        end
      end
    end
  end

  # apply R^T
  fill!(t4L_dotL, 0); fill!(t4L_dotR, 0); fill!(t4R_dotL, 0); fill!(t4R_dotR, 0)
  interiorFaceIntegrate_jac!(sbpface, iface, t3L_dotL, t3R_dotL,
                             t4L_dotL, t4R_dotL, add, subtract, false)
  interiorFaceIntegrate_jac!(sbpface, iface, t3L_dotR, t3R_dotR,
                             t4L_dotR, t4R_dotR, add, subtract, false)

  # multiply by D^T Lambda
  applyDiffusionTensor_diff(diffusion, sbp, params, qL, wL, coordsL, dxidxL,
                            jacL, iface.elementL, temp2L,
                            t4L_dotL, t4L_dotR, t5L_dotL, t5L_dotR)
  # reverse L and R here so the d/dq term will be added to t5R_dotR
  applyDiffusionTensor_diff(diffusion, sbp, params, qR,  wR, coordsR, dxidxR,
                            jacR,iface.elementR, temp2R,
                            t4R_dotR, t4R_dotL, t5R_dotR, t5R_dotL)

  calcDxTransposed(sbp, dxidxL, jacL, Dx)
  applyOperatorJac(Dx, t5L_dotL, resL_dotL, false, op)
  applyOperatorJac(Dx, t5L_dotR, resL_dotR, false, op)

  calcDxTransposed(sbp, dxidxR, jacR, Dx)
  applyOperatorJac(Dx, t5R_dotL, resR_dotL, false, op)
  applyOperatorJac(Dx, t5R_dotR, resR_dotR, false, op)

  return nothing
end

"""
  Jacobian of Dgk * t2L only
"""
function applyDgkTranspose_diff(capture::SBPParabolicSC{Tsol, Tres}, sbp,
                       params::ParamType,
                       sbpface, iface::Union{Interface, Boundary},
                       diffusion::AbstractDiffusion,
                       t2L::AbstractMatrix, t2L_dotL::Abstract4DArray,
                       t2L_dotR::Abstract4DArray,
                       qL::AbstractMatrix,
                       wL::AbstractMatrix,
                       coordsL::AbstractMatrix,
                       nrm_face::AbstractMatrix,
                       dxidxL::Abstract3DArray,
                       jacL::AbstractVector,
                       resL_dotL::Abstract4DArray, resL_dotR::Abstract4DArray,
                       op::SummationByParts.UnaryFunctor=SummationByParts.Add()
                      ) where {Tsol, Tres}

  dim, numNodesPerFace = size(nrm_face)
  numNodesPerElement = size(dxidxL, 3)
  numDofPerNode = size(wL, 1)
  add = SummationByParts.Add()
  subtract = SummationByParts.Subtract()

  @unpack capture temp1L temp2L
  fill!(temp2L, 0);
  @unpack capture t3L_dotL t3L_dotR Dx
  @unpack capture t4L_dotL t4L_dotR
  @unpack capture t5L_dotL t5L_dotR

  # 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(iface), tmp2L, temp1L)
  end

  # multiply by N
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerFace
      @simd for d=1:dim
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            t3L_dotL[i, j, d, p, q] = nrm_face[d, p]*t2L_dotL[i, j, p, q]
            t3L_dotR[i, j, d, p, q] = nrm_face[d, p]*t2L_dotR[i, j, p, q]
          end
        end
      end
    end
  end

  # apply R^T
  fill!(t4L_dotL, 0); fill!(t4L_dotR, 0)
  #TODO: get rid of keyword arguments
  boundaryFaceIntegrate_jac!(sbpface, getFaceL(iface), t3L_dotL, t4L_dotL,
                             include_quadrature=false)
  boundaryFaceIntegrate_jac!(sbpface, getFaceL(iface), t3L_dotR, t4L_dotR,
                             include_quadrature=false)

  # multiply by D^T Lambda
  applyDiffusionTensor_diff(diffusion, sbp, params, qL, wL, coordsL, dxidxL,
                            jacL, getElementL(iface), temp2L,
                            t4L_dotL, t4L_dotR, t5L_dotL, t5L_dotR)

  calcDxTransposed(sbp, dxidxL, jacL, Dx)
  applyOperatorJac(Dx, t5L_dotL, resL_dotL, false, op)
  applyOperatorJac(Dx, t5L_dotR, resL_dotR, false, op)

  return nothing
end


text:
using Pkg
# Julia Kernel
ENV["JUPYTER"]="/opt/conda/bin/jupyter";
Pkg.add("IJulia");
Pkg.build("IJulia");

ext: jl

max_stars_repo_name: colibridigital/jupyter-images

max_issues_repo_name: colibridigital/jupyter-images

max_forks_repo_name: colibridigital/jupyter-images


text:
using Distributed

"""
Use all workers available.
"""
struct ProcParallel <: DistributionStrategy end

function step!(new, old, ::ProcParallel)
    m, n = size(old)
    prepare!(old)
    @sync @distributed for c in 2:TILESIZE:n-1
        for r in 2:TILESIZE:m-1
            rows = r:min(r+TILESIZE-1, m-1)
            cols = c:min(c+TILESIZE-1, n-1)
            I = (rows,cols)
            life_rule!(new, old, I, I)
        end
    end
end


text:
fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
anyerrors = false

println("Running tests:")

for my_test in ["dbi.jl"]
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
push!(LOAD_PATH, joinpath(dirname(@__FILE__), "../.."))
push!(LOAD_PATH, joinpath(dirname(@__FILE__), "..", "deps"))
workspace()
using ArrayFire
using Base.Test
using Base.LinAlg.BLAS

available_backends()
setbackend("cpu")

function getraw(in::AFArray)
  p = device_ptr(in)
  pp = convert(Ptr{Float32}, p)
  pointer_to_array(pp, jl_size(in))
end

x = rand(10, 5)
pointer(x)

a = rand(Float32, 50, 50) |> AFArray
b = rand(Float32, 50, 50) |> AFArray
aa = getraw(a)
bb = getraw(b)

axpy!(0.1f0, aa, bb)
bb

p = device_ptr(a)
pp = convert(Ptr{Float32}, p)
aa = pointer_to_array(pp, (10,5))

p = nothing
finalize(a)

b = rand(Float32, 10, 5)  |> AFArray
a += b
size(a)

a = ArrayFire.af_seq(1.0, 3.0, 1.0)
ArrayFire.af_index_t(a, true, true)
lhs = rand(AFArray, Float32, 10)
rhs = rand(AFArray, Float32, 10)
ArrayFire.assign(lhs, rhs)

ArrayFire.create_indexers()
x = [1,2,3] |> AFArray
x = range(AFArray, Int, 5)
x = rand(Float32,5,4) |> AFArray
i = [3:15] |> AFArray
r = [3.0f0] |> AFArray
ArrayFire.assign(i, x, r)

xs = [rand(Float32,10,5) |> AFArray for i=1:300]
cat(xs, 2)

x = rand(Float32, 10, 5)
xx = AFArray(x)
xxx = AFArray(x)
x[1] = 1.0

x = rand(AFArray{Float32}, 10, 5)

function bench()
  x = rand(AFArray{Float32}, 10, 5, 2, 3)
  for i = 1:10000
    #x = rand(AFArray{Float32}, 100, 10, 10)
    a = af_length(x)
    #unlock(x)
  end
end

gc()
@time bench()

x1 = AFArray(Float32, 1.0, (10, 5))
x1 = rand(AFArray{Float32}, 10, 5)
x2 = rand(AFArray{Float32}, 10, 5)
cat(0, [x1,x2])

tests = ["array"]

for t in tests
    path = joinpath(dirname(@__FILE__), "$t.jl")
    println("$path ...")
    include(path)
end


text:
# let's compare vocoders that are implemented in SynthesisFilters.jl
#
# After you've run the script, you'll get the following synthesized wav files:
# 1. test16k_poledf.wav
# 2. test16k_ltcdf.wav
# 3. test16k_lspdf.wav
# 4. test16k_lmadf.wav
# 5. test16k_mlsadf.wav
# 6. test16k_mglsadf.wav
#
# The original wav file is found in data/test16k.wav.
#
# For more high quality speech waveform syntheis, you can check the `WORLD.jl`
# https://github.com/r9y9/WORLD.jl

using MelGeneralizedCepstrums
using SynthesisFilters
using WAV
using DSP

# Note about excitation
# fs: 16000
# frame period: 5.0 ms
# F0 analysis: esimated  by WORLD.dio and WORLD.stonemask
# Excitation genereration: perioic pulse for voiced segments and gaussian random
# value for un-voiced segments
base_excitation = vec(readdlm(joinpath(dirname(@__FILE__), "data",
                                       "test16k_excitation.txt")))

x, fs = wavread(joinpath(dirname(@__FILE__), "data", "test16k.wav"))
x = vec(x)
framelen = 512
hopsize = 80 # 5.0 ms for fs 16000
noverlap = framelen - hopsize

# Note that mgcep analysis basically assumes power-normalized window
# ŵ(n) = w(n)/|Σₙ w(n)²|¹/²
win = DSP.blackman(framelen) ./ sqrt(sum(abs2, (DSP.blackman(framelen))))

# create windowed signal matrix that each column represents a windowed time slice
as = arraysplit(x, framelen, noverlap)
xw = Array{Float64,2}(framelen, length(as))
for t=1:length(as)
    xw[:,t] = as[t]
end
# col-wise windowing
xw .*= win

function test_poledf_synthesis(; order=25, savepath="test16k_poledf.wav")
    println("testing: poledf_synthesis")
    l = estimate(LinearPredictionCoef(order), xw, use_mgcep=true)
    y = synthesis(base_excitation, l, hopsize)
    wavwrite(y, savepath; Fs=fs)
    println("Dumped to $savepath")
end

function test_ltcdf_synthesis(; order=25, savepath="test16k_ltcdf.wav")
    println("testing: ltcdf_synthesis")
    par = lpc2par(estimate(LinearPredictionCoef(order), xw, use_mgcep=true))
    y = synthesis(base_excitation, par, hopsize)
    wavwrite(y, savepath; Fs=fs)
    println("Dumped to $savepath")
end

function test_lspdf_synthesis(; order=20, savepath="test16k_lspdf.wav")
    println("testing: lspdf_synthesis")
    l = lpc2lsp(estimate(LinearPredictionCoef(order), xw, use_mgcep=true))
    y = synthesis(base_excitation, l, hopsize)
    wavwrite(y, savepath; Fs=fs)
    println("Dumped to $savepath")
end

function test_lmadf_synthesis(; order=25, savepath="test16k_lmadf.wav")
    println("testing: lmadf_synthesis")
    c = estimate(LinearCepstrum(order), xw)
    y = synthesis(base_excitation, c, hopsize)
    wavwrite(y, savepath; Fs=fs)
    println("Dumped to $savepath")
end

function test_mlsadf_synthesis(; order=25, savepath="test16k_mlsadf.wav")
    println("testing: mlsadf_synthesis")
    α = mcepalpha(fs) # automatic α selection
    mc = estimate(MelCepstrum(order, α), xw)
    y = synthesis(base_excitation, mc, hopsize)
    wavwrite(y, savepath; Fs=fs)
    println("Dumped to $savepath")
end

function test_mglsadf_synthesis(; order=25, nstage=6,
                                savepath="test16k_mglsadf.wav")
    println("testing: mglsadf_synthesis")
    α = mcepalpha(fs)
    γ = -1./nstage
    mgc = estimate(MelGeneralizedCepstrum(order, α, γ), xw)
    y = synthesis(base_excitation, mgc, hopsize)
    wavwrite(y, savepath; Fs=fs)
    println("Dumped to $savepath")
end

@time test_poledf_synthesis()
@time test_ltcdf_synthesis()
@time test_lspdf_synthesis()
@time test_lmadf_synthesis()
@time test_mlsadf_synthesis()
@time test_mglsadf_synthesis()


text:
"""
Propagate labels to reach the starting partition. Takes an adjacency matrix
as input and returns a `Modular` object as output.

This approach works better for large graphs (it tends to underperform both the
"random" and "single" approaches for small graphs).

The optimization stops when any label propagation fails to increase the
modularity, measured with `Q`.
"""
function partition_lp(A::Array{Int64, 2})
  nr, nc = size(A)
  min_part = minimum(size(A))
  S = eye(Int64, sum(size(A)))
  M = Modular(A, S)
  old_Q = 0.0
  new_Q = 0.000001
  while new_Q > old_Q
    update_order_top = shuffle([1:nc;])
    update_order_bot = shuffle([1:nr;])
    # Top level
    @inbounds for node_top in update_order_top
       neighbors = [1:nr;][M.A[:,node_top] .> 0] .+ nc
       lab_count = sum(M.S[neighbors,:],1)
       for i = 1:length(lab_count)
          if lab_count[i] == maximum(lab_count)
             M.S[node_top,:] = zeros(Int64, size(M.S)[2])
             M.S[node_top,i] = 1
             break
          end
       end
    end
    # Bottom level
    @inbounds for node_bot in update_order_bot
       neighbors = [1:nc;]'[M.A[node_bot,:] .> 0]
       lab_count = sum(M.S[neighbors,:],1)
       for i = 1:length(lab_count)
          if lab_count[i] == maximum(lab_count)
             M.S[node_bot+nc,:] = zeros(Int64, size(M.S)[2])
             M.S[node_bot+nc,i] = 1
             break
          end
       end
    end
    old_Q = new_Q
    new_Q = Q(M)
  end
  no_empty_modules!(M)
  return M
end


text:
module vector_tile
  include("vector_tile_pb.jl")
end


text:
# Location to monitor pressure at
mon_well_crds = [[-75.0,-75.0]]
# Location of injection and extraction wells
well_crds = [[-50.0,-50.0],[50.0,50.0]]
# Injection rate
Qinj = 0.031688 # [m^3/s] (1 MMT water/yr)
# Time to check pressure
#t = 30.0*86400.0 # 30 days converted to seconds
# Create xy grid for plotting pressures
#xs = ys = -100.0:1.0:100.0
# Training epochs
epochs = 1:100
# Calculate distance between extraction and injection wells
rs = [sqrt((sum((cs-mon_well_crds[1]).^2))) for cs in well_crds]

# Input/output normalization function
lerp(x, lo, hi) = x*(hi-lo)+lo

global nevals = 0
# Set random seed for repeatability
Random.seed!(0)

model = Chain(
    Conv((3, 3), 1=>8, pad=(1,1), relu),
    x -> maxpool(x, (2,2)),

    Conv((3, 3), 8=>16, pad=(1,1), relu),
    x -> maxpool(x, (2,2)),

    Conv((3, 3), 16=>8, pad=(1,1), relu),
    x -> maxpool(x, (2,2)),

    flatten,
    Dense(288, 1),
) |> f64

# Make neural network parameters trackable by Flux
θ = params(model)

# Set ranges for P, T
# Fix pressure target for now. It should be removed, but when I removed it with MultiTheis, the PIML framework didn't work as well. Moral of the story, NN's are squirrelly things that don't always make sense.
PRES = (0.0,0.0) # Pressure target range [m]
TRANS = (-2.0,0.0) # log10 transmissivity [m^2/s]
# Assume that reservoir thickness is 100m (S = S_s*b)
STOR = (-4.0,-1.0) # log10 storativity

# Random T, S, and P generator
sample() = [10^lerp(rand(), TRANS...), 10^lerp(rand(), STOR...), lerp(rand(), PRES...)]

# Wrapper around MultiTheis solve returning pressure head at a point
function head(T::Array{Float64,2},Q1::Float64,Qinj::Float64,rs::Array{Float64,1};train=false)
  Qs = [Q1,Qinj]
  if train
    global nevals = nevals + 1
  end
  num_res = TheisLike.solve_numerical(Qs, T, rs)
end

# Neural network function (T, S, target)->(Q1,Q2)
# 'target' is the target pressure head at the monitoring location
function aim(T)
  T = Flux.unsqueeze(T, 3)
  T = Flux.unsqueeze(T, 3)
  Q1 = model(T)

  # Make sure that Q1 and Q2 are always negative for extraction
  Q1 = Q1[1]
  #Q2 = -softplus(Q2)
  Q1
end

# Physics model vs NN
function mydiff(T, target; train=false)
    #@show myid()
    (head(T, aim(T), Qinj, rs, train=train) - target)
end

function loss(T, target; train=true)
  mydiff(T, target, train=train)^2
end

function cb2()
  # Terminal output
  loss_train = sum(map(x->loss(x,train=false),data_train_batch))
  loss_test = sum(map(x->loss(x,train=false),data_test))
  println(string("callback: train loss: ", loss_train, " test loss: ", loss_test))
end

cb = Flux.throttle(cb2, 10)
opt = ADAM()

loss(x; train=true) = sum(pmap(x->loss(x[1], x[2]; train=train), x) )
mydiff(x; train=true) = sum(pmap(x->mydiff(x[1], x[2]; train=train), x) )

global losses_train = Float64[]
global losses_test = Float64[]
global rmses_train = Float64[]
global rmses_test = Float64[]
global meandiffs_train = Float64[]
global stddiffs_train = Float64[]
global meandiffs_test = Float64[]
global stddiffs_test = Float64[]
global train_time = Float64[]


text:
"""
    StyleTagProxy(io)

This is a transparent proxy that ensures neither `<!--` nor `</style>`
occur in the output stream.

# Examples
```julia-repl
julia> gp = StyleTagProxy(stdout);
julia> print(gp, "valid");
valid
julia> print(gp, "</style>")
ERROR: "Content within a style tag must not contain `</style>`"
```
"""
mutable struct StyleTagProxy{T<:IO} <: IO where {T}
    io::T
    index::Int

    StyleTagProxy(io::T) where T = new{T}(io::T, 0)
end

"""
    StyleTag(data)

This object prints `data` unescaped within a `<style>` tag, wrapped in a
`StyleTagProxy` that guards against invalid style content. Content is
treated as if it had occurred within an attribute value, only that
amperstand escaping is not used.
"""
struct StyleTag
    content
end

Base.print(ep::EscapeProxy, x::StyleTag) =
    print_value(StyleTagProxy(ep.io), x.content)

function scan_for_style(index::Int, octet::UInt8)::Int
    if 1 == index
        return octet == Int('!') ? 12 : octet == Int('/') ? 2 : 0
    elseif 2 == index
        return (octet == Int('S') || octet == Int('s')) ? 3 : 0
    elseif 3 == index
        return (octet == Int('T') || octet == Int('t')) ? 4 : 0
    elseif 4 == index
        return (octet == Int('Y') || octet == Int('y')) ? 5 : 0
    elseif 5 == index
        return (octet == Int('L') || octet == Int('l')) ? 6 : 0
    elseif 6 == index
        return (octet == Int('E') || octet == Int('e')) ? 7 : 0
    elseif 7 == index
        if octet == Int('>')
            throw("Content within a style tag must not contain `</style>`")
        end
    elseif 12 == index
        return octet == Int('-') ? 13 : 0
    elseif 13 == index
        if octet == Int('-')
            throw("Content within a style tag must not contain `<!--`")
        end
    else
        @assert false # unreachable?!
    end
    return 0
end

function Base.write(sp::StyleTagProxy, octet::UInt8)
    if 0 == sp.index
        if octet == Int('<')
            sp.index = 1
        end
        return write(sp.io, octet)
    end
    sp.index = scan_for_style(sp.index, octet)
    return write(sp.io, octet)
end

function Base.unsafe_write(sp::StyleTagProxy, input::Ptr{UInt8}, nbytes::UInt)
    cursor = input
    index = sp.index
    final = input + nbytes
    while cursor < final
        octet = unsafe_load(cursor)
        if octet == Int('<')
            index = 1
        elseif 0 != index
            index = scan_for_style(index, octet)
        end
        cursor += 1
    end
    sp.index = index
    return unsafe_write(sp.io, input, nbytes)
end


text:
struct DecisionBoundaryQss <: ModelBasedQss
    occ::SVDD.OCClassifier
    optimizer::QuerySynthesisOptimizer
    DecisionBoundaryQss(occ; optimizer::QuerySynthesisOptimizer) = new(occ, optimizer)
end

function qs_score_function(qs::DecisionBoundaryQss, data::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Function where T <: Real
    return x -> -abs.(SVDD.predict(qs.occ, x))
end


text:
# Schrodinger equation
SchrodingerEvo(H₀::Operator) = Liouvillian(dims(H₀),-1im*data(H₀))

function SchrodingerEvo(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}})
    dimsmatch(H₀,first.(Hₙ))
    # Unpack time-dependent Hamiltonian terms
    Lₙ,fₙ,pₙ = unpack_operators(Hₙ)
    return Liouvillian(dims(H₀),-1im*data(H₀),-1im.*data.(Lₙ),fₙ,pₙ)
end
SchrodingerEvo(H₀::Operator,Hₙ::Vararg{Tuple}) = SchrodingerEvo(H₀,Hₙ)
SchrodingerEvo(H::Tuple{Operator,Vararg{Tuple}}) = SchrodingerEvo(first(H),tail(H))
SchrodingerEvo(::Any) = throw(ArgumentError("invalid Hamiltonian specification"))

# von Neumann equation
# TODO

# Lindblad form master equation
function LindbladEvo(H₀::Operator, Cₘ::Tuple{Vararg{Operator}})
    dimsmatch(H₀,Cₘ)
    Id = data(qeye(dims(H₀)))
    # Constant Hamiltonian term
    L₀ = -1im*super_vonneumann(H₀,Id)
    # Add constant collapse operator terms
    L₀ += sum_collapse(Cₘ,Id)
    return Liouvillian(dims(H₀),L₀)
end
function LindbladEvo(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, Cₘ::Tuple{Vararg{Operator}})
    dimsmatch(H₀,first.(Hₙ)); dimsmatch(H₀,Cₘ)
    Id = data(qeye(dims(H₀)))
    # Constant Hamiltonian term
    L₀ = -1im*super_vonneumann(H₀,Id)
    # Add constant collapse operator terms
    L₀ += sum_collapse(Cₘ,Id)
    # Unpack time-dependent Hamiltonian terms
    Lₙ,fₙ,pₙ = unpack_operators(Hₙ)
    return Liouvillian(dims(H₀),L₀,-1im.*super_vonneumann.(Lₙ,(Id,)),fₙ,pₙ)
end
LindbladEvo(H₀::Operator,Cₘ::Vararg{Operator}) = LindbladEvo(H₀,Cₘ)
LindbladEvo(H::Tuple{Operator,Vararg{Tuple}},Cₘ::Vararg{Operator}) = LindbladEvo(first(H),tail(H),Cₘ)
LindbladEvo(H::Tuple{Operator,Vararg{Tuple}},Cₘ::Tuple{Vararg{Operator}}) = LindbladEvo(first(H),tail(H),Cₘ)
LindbladEvo(::Any) = throw(ArgumentError("invalid Hamiltonian specification"))

# Schrodinger evo propagator
SchrodingerProp(H₀::Operator, tspan) = SchrodingerProp(H₀,float(tspan[2]-tspan[1]))
function SchrodingerProp(H₀::Operator, Δt::Real)
    # Constant Hamiltonian term
    U = expim(Hermitian(-Array(H₀).*Δt))
    return Propagator(U,Δt,dims(H₀))
end
SchrodingerProp(H₀::Operator, Hₙ::Tuple, tspan, steps::Integer) = SchrodingerProp(H₀,(Hₙ,),tspan,steps)
function SchrodingerProp(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, tspan, steps::Integer)
    dimsmatch(H₀,first.(Hₙ)); F = real(eltype(H₀))
    # Sampling times and spacing dt
    t₁, t₂ = tspan; dt = (t₂-t₁)/steps
    # Unpack time-dependent Hamiltonian terms
    H_tdep = unpack_operators(Hₙ)
    # Multiply sampled propagators together to generate total evolution
    U = Matrix{Complex{F}}(I, size(H₀))
    H = Hermitian(zeros(compute_H_type(H₀,H_tdep),size(H₀)...))
    A = similar(U); B = similar(U); C = similar(H.data); D = similar(U)
    Λ = Vector{F}(undef,size(H,1))
    for i = 1:steps
        step_hamiltonian!(H.data,H₀,H_tdep,(t₁,dt,i)) # calc H for this time step
        expim!(A,H,Λ,C,D) # A = exp(-1im*H*dt)
        mul!(B,A,U) # U = A*U
        U,B = B,U # swappitty swap for the next step
    end
    return Propagator(U,float(t₂-t₁),dims(H₀))
end
SchrodingerProp(H) = throw(ArgumentError("invalid Propagator specification"))

# Lindblad evo propagator
LindbladProp(H₀::Operator, Cₘ::Operator, tspan) = LindbladProp(H₀,(Cₘ,),tspan)
LindbladProp(H₀::Operator, Cₘ::Tuple{Vararg{Operator}}, tspan) = LindbladProp(H₀,Cₘ,tspan[2]-tspan[1])
function LindbladProp(H₀::Operator, Cₘ::Tuple{Vararg{Operator}}, Δt::Real)
    dimsmatch(H₀,Cₘ)
    Id = Diagonal(I, prod(dims(H₀)))
    # Constant Hamiltonian term
    L₀Δt = (-1im*Δt).*(Id⊗Array(H₀) .- transpose(Array(H₀))⊗Id)
    # Add constant collapse operator terms
    L₀Δt .+= Δt.*sum_collapse(dense.(Cₘ),Id)
    # Build constant propagator
    U = LinearAlgebra.exp!(L₀Δt)
    return Propagator(U,Δt,dims(H₀))
end
LindbladProp(H₀::Operator, Hₙ::Tuple, Cₘ::Operator, tspan, steps::Integer) = LindbladProp(H₀,(Hₙ,),(Cₘ,),tspan,steps)
LindbladProp(H₀::Operator, Hₙ::Tuple, Cₘ::Tuple{Vararg{Operator}}, tspan, steps::Integer) = LindbladProp(H₀,(Hₙ,),Cₘ,tspan,steps)
LindbladProp(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, Cₘ::Operator, tspan, steps::Integer) = LindbladProp(H₀,Hₙ,(Cₘ,),tspan,steps)
function LindbladProp(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, Cₘ::Tuple{Vararg{Operator}}, tspan, steps::Integer)
    dimsmatch(H₀,first.(Hₙ)); dimsmatch(H₀,Cₘ); F = real(eltype(H₀));
    Id = Diagonal(I, prod(dims(H₀)))
    # Sampling times and spacing dt
    t₁, t₂ = tspan; dt = (t₂-t₁)/steps
    # Unpack time-dependent Hamiltonian terms
    H_tdep = unpack_operators(Hₙ)
    # Build constant collapse propagator part
    U₀ = LinearAlgebra.exp!(dt.*sum_collapse(dense.(Cₘ),Id))
    # Multiply sampled propagators together to generate total evolution
    U = Matrix{Complex{F}}(I, size(U₀))
    H = Hermitian(zeros(compute_H_type(H₀,H_tdep),size(H₀)...))
    A = Matrix{Complex{F}}(undef,size(H₀)...); B = similar(U); C = similar(H.data); D = similar(A)
    Λ = Vector{F}(undef,size(H,1))
    for i = 1:steps
        step_hamiltonian!(H.data,H₀,H_tdep,(t₁,dt,i))
        expim!(A,H,Λ,C,D) # A = exp(-1im*H*dt)
        invA = LinearAlgebra.inv!(lu(A))
        #
        mul!(B,U₀,U) # TODO: use the Lie product formula here for better results
        mul!(U,transpose(invA⊗Id),B)
        I_kron_mul!(B,A,U)
        U,B = B,U # swappitty swap for the next step
    end
    return Propagator(U,float(t₂-t₁),dims(H₀))
end
LindbladProp(H) = throw(ArgumentError("invalid Propagator specification"))

# Utils
function unpack_operators(t::Tuple{Tuple{Operator,Function,Array},Vararg{Tuple}})
    a, b = first(t), unpack_operators(tail(t))
    return (a[1], b[1]...), (a[2], b[2]...), (a[3], b[3]...)
end
function unpack_operators(t::Tuple{Tuple{Operator,Function},Vararg{Tuple}})
    a, b = first(t), unpack_operators(tail(t))
    return (a[1], b[1]...), (a[2], b[2]...), (Float64[], b[3]...)
end
unpack_operators(t::Tuple{}) = (), (), ()

function super_vonneumann(H,Id)
    # represents the action [H,ρ] as a linear superoperator I⊗H - Hᵀ⊗I
    return Id⊗data(H) - copy(transpose(data(H)))⊗Id
end

function sum_collapse(Cₘ,Id)
    # sum the Lindblad superoperators corresponding to the action of
    # D[Cᵢ](ρ) = CᵢρCᵢ† - 1/2 * (Cᵢ†Cᵢρ + ρCᵢ†Cᵢ) as a linear superoperator
    mapreduce(+,data.(Cₘ)) do Cᵢ
        CdC = Cᵢ'*Cᵢ
        return conj(Cᵢ)⊗Cᵢ - (Id⊗CdC + copy(transpose(CdC))⊗Id)./2
    end
end

function step_hamiltonian!(H,H0,Hn,tspec)
    Hₙ,fₙ,pₙ = Hn
    t₁, dt, i = tspec; tᵢ = t₁ + (i-1)*dt
    H .= -dt.*data(H0) # constant part
    for j = 1:length(Hₙ) # sample function at 3 points
        Hᵢ,fᵢ,pᵢ = data(Hₙ[j]),fₙ[j],pₙ[j]
        H .+= (-dt*(fᵢ(tᵢ,pᵢ) + fᵢ(tᵢ+0.5dt,pᵢ) + fᵢ(tᵢ+dt,pᵢ))/3) .* Hᵢ
    end
end

function compute_H_type(H0,Hn)
    T1 = promote_eltype(H0,Hn[1]...)
    T2 = promote_type(_typeof_H_funs(Hn[2],Hn[3])...)
    return promote_type(T1,T2)
end

@inline function _typeof_H_funs(funs::NTuple{N,Function},ps::NTuple{N,Vector}) where N
    this_T = typeof(first(funs)(0,first(ps)))
    return (this_T, _typeof_H_funs(tail(funs),tail(ps))...)
end

@inline _typeof_H_funs(funs::Tuple{},ps::Tuple{}) = ()


text:
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.NonlinearEquations (ModiaMath/NonlinearEquations/_module.jl)
#

"""
    module KINSOL - Solve nonlinear equation system with Sundials KINSOL
The goal is to solve the same system several times with KINSOL.
"""

module KINSOL

import NLsolve
import ModiaMath
using LinearAlgebra
using NLsolve

mutable struct NonlinearEquationsInfo
    extraInfo                # Model-specific extra information
    name::String             # Name of equation system
    ny::Int                  # Number of equations (length of y-vector)
    getResidues!::Function   # Function of the nonlinear equation system
    y0::Vector{Float64}      # The initial y vector (to print in error messages)
    lastNorm_r::Float64
    lastrScaledNorm_r::Float64

    function NonlinearEquationsInfo(name::String, ny::Int, getResidues!::Function; extraInfo=nothing)
        @assert(ny >= 0)
        new(extraInfo, name, ny, getResidues!, zeros(ny), 1.0, 1.0)
    end
end

function solveNonlinearEquations!(eqInfo::NonlinearEquationsInfo, y::Vector{Float64};
                                  FTOL::Float64=eps(Float64)^(1 / 3),
                                  yScale::Vector{Float64}=ones(length(y)),
                                  rScale::Vector{Float64}=ones(length(y)))
                                  #inputs are the same for backwards compatibility
    eqInfo.y0 .= y
    itnum = Int32(1000.0 * norm(yScale, Inf))
    nsol_f!(r, y) = eqInfo.getResidues!(eqInfo, y, r)
    nsol = nlsolve(nsol_f!, y, xtol = FTOL, ftol=FTOL, iterations=itnum)
    eqInfo.y0 .= nsol.zero
    #println("nsol = $nsol")
    #return nothing
end

end


text:
import Mustache, Highlights
import .WeaveMarkdown
using Compat
using Dates
using Markdown
using REPL.REPLCompletions: latex_symbols

function format(doc::WeaveDoc)
    formatted = AbstractString[]
    docformat = doc.format

    #Complete format dictionaries with defaults
    formatdict = docformat.formatdict
    get!(formatdict, :termstart, formatdict[:codestart])
    get!(formatdict, :termend, formatdict[:codeend])
    get!(formatdict, :out_width, nothing)
    get!(formatdict, :out_height, nothing)
    get!(formatdict, :fig_pos, nothing)
    get!(formatdict, :fig_env, nothing)

    docformat.formatdict[:cwd] = doc.cwd #pass wd to figure formatters
    docformat.formatdict[:theme] = doc.highlight_theme

    #strip header
    if isa(doc.chunks[1], DocChunk)
        if !occursin("pandoc", doc.doctype)
            doc.chunks[1] = strip_header(doc.chunks[1])
        end
    end

    for chunk in copy(doc.chunks)
        result = format_chunk(chunk, formatdict, docformat)
        push!(formatted, result)
    end

    formatted = join(formatted, "\n")
    # Render using a template if needed
    rendered = render_doc(formatted, doc, doc.format)

    return rendered
end

"""
  render_doc(formatted::AbstractString, format)

Render formatted document to a template
"""
function render_doc(formatted, doc::WeaveDoc, format)
  return formatted
end

function highlight(mime::MIME, source::AbstractString, lexer, theme=Highlights.Themes.DefaultTheme)
    return sprint( (io, x) -> Highlights.highlight(io, mime, x, lexer, theme), source)
end

function stylesheet(m::MIME, theme)
  return sprint( (io, x) ->  Highlights.stylesheet(io, m, x), theme)
end

function render_doc(formatted, doc::WeaveDoc, format::JMarkdown2HTML)
  css = stylesheet(MIME("text/html"), doc.highlight_theme)
  path, wsource = splitdir(abspath(doc.source))
  #wversion = string(Pkg.installed("Weave"))
  wversion = ""
  wtime =  string(Date(now()))

  if isempty(doc.css)
    theme_css = read(joinpath(dirname(@__FILE__), "../templates/skeleton_css.css"), String)
  else
    theme_css = read(doc.css, String)
  end

  if isa(doc.template, Mustache.MustacheTokens)
    template = doc.template
  elseif isempty(doc.template)
    template = Mustache.template_from_file(joinpath(dirname(@__FILE__), "../templates/julia_html.tpl"))
  else
    template = Mustache.template_from_file(doc.template)
  end

  return Mustache.render(template; themecss = theme_css,
                          highlightcss = css, body = formatted, header_script = doc.header_script,
                          source = wsource, wtime = wtime, wversion = wversion,
                          [Pair(Symbol(k), v) for (k,v) in doc.header]...)
end

function render_doc(formatted, doc::WeaveDoc, format::JMarkdown2tex)
  highlight = stylesheet(MIME("text/latex"), doc.highlight_theme)
  path, wsource = splitdir(abspath(doc.source))
  #wversion = string(Pkg.installed("Weave"))
  wversion = ""
  wtime =  string(Date(now()))

  if isa(doc.template, Mustache.MustacheTokens)
      template = doc.template
  elseif isempty(doc.template)
    template = Mustache.template_from_file(joinpath(dirname(@__FILE__), "../templates/julia_tex.tpl"))
  else
    template = Mustache.template_from_file(doc.template)
  end

  return Mustache.render(template; body = formatted,
    highlight = highlight,
    [Pair(Symbol(k), v) for (k,v) in doc.header]...)
end

function strip_header(chunk::DocChunk)
  if occursin(r"^---$(?<header>.+)^---$"ms, chunk.content[1].content)
    chunk.content[1].content = lstrip(replace(chunk.content[1].content, r"^---$(?<header>.+)^---$"ms => ""))
  end
  return chunk
end

function format_chunk(chunk::DocChunk, formatdict, docformat)
    return join([format_inline(c) for c in chunk.content], "")
end

function format_inline(inline::InlineText)
    return inline.content
end

function format_inline(inline::InlineCode)
    isempty(inline.rich_output) || return inline.rich_output
    isempty(inline.figures) || return inline.figures[end]
    isempty(inline.output) || return inline.output
end

function ioformat!(io::IOBuffer, out::IOBuffer, fun=WeaveMarkdown.latex)
    text = String(take!(io))
    if !isempty(text)
        m = Markdown.parse(text, flavor=WeaveMarkdown.weavemd)
        write(out, string(fun(m)))
    end
end

function addspace(op, inline)
    inline.ctype == :line && (op = "\n$op\n")
    return op
end

function format_chunk(chunk::DocChunk, formatdict, docformat::JMarkdown2tex)
    out = IOBuffer()
    io = IOBuffer()
    for inline in chunk.content
        if isa(inline, InlineText)
            write(io, inline.content)
        elseif !isempty(inline.rich_output)
            ioformat!(io, out)
            write(out, addspace(inline.rich_output, inline))
        elseif !isempty(inline.figures)
            write(io, inline.figures[end], inline)
        elseif !isempty(inline.output)
            write(io, addspace(inline.output, inline))
        end
    end
    ioformat!(io, out)
    formatdict[:keep_unicode] || return uc2tex(String(take!(out)))
    return String(take!(out))
end

function format_chunk(chunk::DocChunk, formatdict, docformat::JMarkdown2HTML)
    out = IOBuffer()
    io = IOBuffer()
    fun = WeaveMarkdown.html
    for inline in chunk.content
        if isa(inline, InlineText)
            write(io, inline.content)
        elseif !isempty(inline.rich_output)
            ioformat!(io, out, fun)
            write(out, addspace(inline.rich_output, inline))
        elseif !isempty(inline.figures)
            write(io, inline.figures[end])
        elseif !isempty(inline.output)
            write(io, addspace(inline.output, inline))
        end
    end
    ioformat!(io, out, fun)
    return String(take!(out))
end

function format_chunk(chunk::CodeChunk, formatdict, docformat)
    #Fill undefined options with format specific defaults
    chunk.options[:out_width] == nothing &&
        (chunk.options[:out_width] =  formatdict[:out_width])
    chunk.options[:fig_pos] == nothing &&
        (chunk.options[:fig_pos] =  formatdict[:fig_pos])

    #Only use floats if chunk has caption or sets fig_env
    if chunk.options[:fig_cap] != nothing && chunk.options[:fig_env] == nothing
        (chunk.options[:fig_env] =  formatdict[:fig_env])
    end

    if haskey(formatdict, :indent)
        chunk.content = indent(chunk.content, formatdict[:indent])
    end

    chunk.content = format_code(chunk.content, docformat)

    if !chunk.options[:eval]
        if chunk.options[:echo]
            result = "$(formatdict[:codestart])\n$(chunk.content)$(formatdict[:codeend])"
            return result
        else
            r = ""
            return r
        end
    end

    if chunk.options[:term]
        result = format_termchunk(chunk, formatdict, docformat)
    else

    if chunk.options[:echo]
      #Convert to output format and highlight (html, tex...) if needed
        result = "$(formatdict[:codestart])$(chunk.content)$(formatdict[:codeend])\n"
    else
        result = ""
    end

    if (strip(chunk.output)!= "" || strip(chunk.rich_output) != "") && (chunk.options[:results] != "hidden")
        if chunk.options[:results] != "markup" && chunk.options[:results] != "hold"
            strip(chunk.output) ≠ "" && (result *= "$(chunk.output)\n")
            strip(chunk.rich_output) ≠ "" && (result *= "$(chunk.rich_output)\n")
        else
            if chunk.options[:wrap]
                chunk.output = "\n" * wraplines(chunk.output, chunk.options[:line_width])
                chunk.output = format_output(chunk.output, docformat)
            else
                chunk.output = "\n" * rstrip(chunk.output)
                chunk.output = format_output(chunk.output, docformat)
            end

            if haskey(formatdict, :indent)
                chunk.output = indent(chunk.output, formatdict[:indent])
            end
            strip(chunk.output) ≠ "" &&
                (result *= "$(formatdict[:outputstart])$(chunk.output)\n$(formatdict[:outputend])\n")
            strip(chunk.rich_output) ≠ "" && (result *= chunk.rich_output * "\n")
        end
    end

    end

    #Handle figures
    if chunk.options[:fig] && length(chunk.figures) > 0
        if chunk.options[:include]
            result *= formatfigures(chunk, docformat)
        end
    end

    return result
end

function format_output(result::AbstractString, docformat)
  return result
end

function format_output(result::AbstractString, docformat::JMarkdown2HTML)
  return Markdown.htmlesc(result)
end

function format_output(result::AbstractString, docformat::JMarkdown2tex)
  # Highligts has some extra escaping defined, eg of $, ", ...
  result_escaped = sprint( (io, x) -> Highlights.Format.escape(io, MIME("text/latex"), x, charescape=true), result)
  docformat.formatdict[:keep_unicode] || return uc2tex(result_escaped, true)
  return result_escaped
end

function format_code(result::AbstractString, docformat)
  return result
end

function format_code(result::AbstractString, docformat::JMarkdown2tex)
  highlighted = highlight(MIME("text/latex"), strip(result),
        Highlights.Lexers.JuliaLexer, docformat.formatdict[:theme])
  docformat.formatdict[:keep_unicode] || return uc2tex(highlighted)
  return highlighted
  #return "\\begin{minted}[mathescape, fontsize=\\small, xleftmargin=0.5em]{julia}\n$result\n\\end{minted}\n"
end

#Convert unicode to tex, escape listings if needed
function uc2tex(s, escape=false)
    for key in keys(latex_symbols)
        if escape
            s = replace(s, latex_symbols[key] => "(*@\\ensuremath{$(texify(key))}@*)")
        else
            s = replace(s, latex_symbols[key] => "\\ensuremath{$(texify(key))}")
        end
    end
    return s
end

# Make julia symbols (\bf* etc.) valid latex
function texify(s)
    ts = ""
    if occursin(r"^\\bf[A-Z]$", s)
        ts = replace(s, "\\bf" => "\\bm{\\mathrm{") * "}}"
    elseif startswith(s, "\\bfrak")
        ts = replace(s, "\\bfrak" => "\\bm{\\mathfrak{") * "}}"
    elseif startswith(s, "\\bf")
        ts = replace(s, "\\bf" => "\\bm{\\") * "}"
    elseif startswith(s, "\\frak")
        ts = replace(s, "\\frak" => "\\mathfrak{") * "}"
    else
        ts = s
    end
    return ts
end

function format_code(result::AbstractString, docformat::JMarkdown2HTML)
  return highlight(MIME("text/html"), strip(result),
    Highlights.Lexers.JuliaLexer, docformat.formatdict[:theme])
end

function format_code(result::AbstractString, docformat::Pandoc2HTML)
    return highlight(MIME("text/html"), strip(result),
      Highlights.Lexers.JuliaLexer, docformat.formatdict[:theme])
end

function format_termchunk(chunk, formatdict, docformat)
    if chunk.options[:echo] && chunk.options[:results] != "hidden"
        result = "$(formatdict[:termstart])$(chunk.output)\n" * "$(formatdict[:termend])\n"
    else
        result = ""
    end
    return result
end

function format_termchunk(chunk, formatdict, docformat::JMarkdown2HTML)
    if chunk.options[:echo] && chunk.options[:results] != "hidden"
        result = highlight(MIME("text/html"), strip(chunk.output),
          Highlights.Lexers.JuliaConsoleLexer, docformat.formatdict[:theme])
    else
        result = ""
    end
    return result
end

function format_termchunk(chunk, formatdict, docformat::Pandoc2HTML)
    if chunk.options[:echo] && chunk.options[:results] != "hidden"
        result = highlight(MIME("text/html"), strip(chunk.output),
          Highlights.Lexers.JuliaConsoleLexer, docformat.formatdict[:theme])
    else
        result = ""
    end
    return result
end

function format_termchunk(chunk, formatdict, docformat::JMarkdown2tex)
    if chunk.options[:echo] && chunk.options[:results] != "hidden"
        result = highlight(MIME("text/latex"), strip(chunk.output),
                  Highlights.Lexers.JuliaConsoleLexer,
                  docformat.formatdict[:theme])
        #return "\\begin{minted}[mathescape, fontsize=\\small, xleftmargin=0.5em]{julia}\n$result\n\\end{minted}\n"
    else
        result = ""
    end
    return result
end

function indent(text, nindent)
    return join(map(x->
                    string(repeat(" ", nindent), x), split(text, "\n")), "\n")
end

function wraplines(text, line_width=75)
    result = AbstractString[]
    lines = split(text, "\n")
    for line in lines
        if length(line) > line_width
            push!(result, wrapline(line, line_width))
        else
            push!(result, line)
        end
    end

    return strip(join(result, "\n"))
end

function wrapline(text, line_width=75)
    result = ""
    while length(text) > line_width
        result*= first(text, line_width) * "\n"
        text = chop(text, head=line_width, tail=0)
    end
    result *= text
end


text:
module AstroTime

using Reexport

import Dates
import EarthOrientation
import MacroTools

export @timescale

include("TimeScales.jl")
include("Periods.jl")
include("AstroDates.jl")
include("Epochs.jl")

@reexport using .TimeScales
@reexport using .Periods
@reexport using .AstroDates
@reexport using .Epochs

import .Epochs: format

function __init__()
    Dates.CONVERSION_SPECIFIERS['t'] = TimeScale
    Dates.CONVERSION_SPECIFIERS['D'] = Epochs.DayOfYearToken
    Dates.CONVERSION_DEFAULTS[TimeScale] = TimeScales.NotATimeScale()
    Dates.CONVERSION_DEFAULTS[Epochs.DayOfYearToken] = Int64(0)

    Dates.CONVERSION_TRANSLATIONS[Epoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
        TimeScale,
    )

    Dates.CONVERSION_TRANSLATIONS[TAIEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[UTCEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[UT1Epoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TTEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TCGEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TCBEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TDBEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )
end

"""
    @timescale scale [parent[, oneway]]

Define a new time scale and the corresponding `Epoch` type alias.

### Arguments ###

- `scale`: The name of the time scale
- `parent`: The "parent" time scale to which it should be linked (optional)
- `oneway`: If `true`, only the transformation from `parent` to `scale` is
    registered (optional, default: `false`)

# Example

```jldoctest
julia> @timescale GMT UTC

julia> GMT isa TimeScale
true

julia> GMTEpoch
Epoch{GMTScale,T} where T

julia> find_path(TAI, GMT)
3-element Array{TimeScale,1}:
 TAI
 UTC
 GMT
```
"""
macro timescale(scale::Symbol, parent=nothing, oneway=false)
    scale_type = Symbol(scale, "Scale")
    epoch_type = Symbol(scale, "Epoch")
    name = String(scale)

    MacroTools.@esc scale scale_type epoch_type

    if parent === nothing
        reg_expr = MacroTools.@q(register_scale!($scale))
    else
        MacroTools.@esc parent
        reg_expr = MacroTools.@q(link_scales!($parent, $scale, oneway=$oneway))
    end

    return MacroTools.@q begin
        struct $scale_type <: TimeScale end
        const $scale = $scale_type()
        const $epoch_type = Epoch{$scale_type}
        Base.show(io::IO, ::$scale_type) = print(io, $name)
        AstroTime.TimeScales.tryparse(::Val{Symbol($name)}) = $scale

        Dates.CONVERSION_TRANSLATIONS[$epoch_type] = (
            Dates.Year,
            Dates.Month,
            Dates.Day,
            Dates.Hour,
            Dates.Minute,
            Dates.Second,
            Dates.Millisecond,
        )
        Dates.default_format(::Type{$epoch_type}) = Dates.ISODateTimeFormat

        $reg_expr

        nothing
    end
end

function load_eop(files...)
    EarthOrientation.push!(EarthOrientation.EOP_DATA, files...)
end

function load_test_eop()
    finals = joinpath(@__DIR__, "..", "test", "data", "finals.csv")
    finals2000A = joinpath(@__DIR__, "..", "test", "data", "finals2000A.csv")
    load_eop(finals, finals2000A)
end

"""
    AstroTime.update()

Download up-to-date IERS tables from the internet to enable transformations to and from
the [`UT1`](@ref) time scale.
"""
function update()
    EarthOrientation.update()
    nothing
end

end # module


text:
using YaoExtensions, Test, Yao

@testset "gate count, time" begin
    qc = QFTCircuit(3)
    @test qc |> gatecount |> length == 2
    @test qc |> gatecount |> values |> sum == 6
    qc = rand_supremacy2d(4,4,3)
    gc = qc |> gatecount
    @test gc[HGate] == 32
end


text:
module ManualScraper
using CHM
using DataStructures
import JSON
import Base: contains, delete!, show

export CHMSource
export Parser
export @ena_str
export @awg5k_str

export scrape

immutable Parser
    init::AbstractArray
    getset::Regex
    cmd::Regex
    rettype::Regex
    cmd1::AbstractArray
    cmd2::AbstractArray
    arg::Regex
    val::AbstractArray
end

immutable CHMSource{T}
    path::AbstractString
    dirs::AbstractArray
    filter::Regex
    parser::Parser
end
show(io::IO, x::CHMSource) = print(io, "CHM source at: $(x.path)")

# Modules describing how to parse each different .chm file.
include("ENA.jl")
include("AWG5K.jl")

"""
`@awg5k_str(path)`

String macro. Given a path to the Tektronix AWG5000 series' .chm file,
construct a `CHMSource` that contains information on how to parse the .chm file.

Usage: `src = awg5k"/path/to/file.chm"`

When entering a file path using this macro, do not escape spaces in the path.
"""
macro awg5k_str(path)
    CHMSource{AWG5K.SYM}(path, AWG5K.DIRS, AWG5K.FILTER, AWG5K.PARSER)
end

"""
`@ena_str(path)`

String macro. Given a path to the Keysight ENA's .chm file, construct a
`CHMSource` that contains information on how to parse the .chm file.

Usage: `src = ena"/path/to/file.chm"`

When entering a file path using this macro, do not escape spaces in the path.
"""
macro ena_str(path)
    CHMSource{ENA.SYM}(path, ENA.DIRS, ENA.FILTER, ENA.PARSER)
end

"""
`capture_names(r::Regex)`

Returns a ValueIterator of the capture group names for a given Regex.
Test to see if a name is there via: `"abc" in capture_names(r)`

Not really used in the code but I thought it was strange that it is not in Base.
"""
function capture_names(r::Regex)
    values(Base.PCRE.capture_names(r.regex))
end

"""
`scrape(inpath::AbstractString, outpath::AbstractString)`

Function used to scrape information out of the Yokogawa GS200 manual.
The manual should be exported to plain text using Adobe Acrobat as other
choices seem to give poor results.
Eventually this will be restructured to accommodate any plain text file.
"""
function scrape(inpath::AbstractString, outpath::AbstractString)
    f = open(inpath)
    str = readall(f)
    close(f)

    r = r"Function\s+([\w\d\s\(\)/,’'.]+)Syntax\s+(:[A-Za-z:\[\]/?]+)\s+((?:(?:[\[\]<>A-Za-z0-9\s]+\|)+[A-Za-z0-9\[\]]+)?(?:<character string>)?)"
    optional = r"\[:[A-Za-z]+\]"

    entries = eachmatch(r, str)

    docs     = UTF8String[]
    commands = ASCIIString[]
    types = ASCIIString[]
    values   = []

    # Initial parse pass.
    for e in entries
        # Get rid of extraneous newlines.
        push!(docs, replace(e[1], "\n", ""))
        # Strip optional parts of the commands and make everything uppercase.
        c = replace(e[2], optional, "")
        push!(types, cmdtostr(c))
        c = replace(c, r"[a-z]+", "")
        push!(commands, uppercase(c))
        opt = e[3]
        opt = replace(opt, "|", ",")
        opt = replace(opt, r"\s", "")
        opt = replace(opt, "<characterstring>", "s::ASCIIString")
        opt = replace(opt, "1,0,ON,OFF", "v::Bool")
        opt = replace(opt, "<voltage>,<current>", "v::Real")
        opt = replace(opt, "<voltage>", "v::Real")
        opt = replace(opt, "<current>", "v::Real")
        opt = replace(opt, "<time>", "v::Real")
        opt = replace(opt, "<integer>", "v::Integer")
        opt = replace(opt, r"[\[\]]", "")

        # opt = replace(opt, r"[a-z]+", "")
        push!(values, opt)
    end

    symbols  = Array{Dict}(length(values))

    v2 = copy(values)
    for (i,v) in enumerate(v2)
        a = parse(v)
        if isa(a, Symbol)
            # Just a single symbol... unexpected.
            error("Unexpected argument.")
        elseif isa(a, Expr)
            if a.head == :(::)
                # Just a single type
                symbols[i] = Dict()
            elseif a.head == :tuple
                if isa(a.args[1], Expr)
                    # Both values and symbols. Just keep the values.
                    values[i] = string(a.args[1])
                    symbols[i] = Dict()
                    # push!(values, "v::Symbol in symbols")
                    # s = map(string, a.args[2:end])
                    # push!(symbols, Dict(zip(map(upperfirst,s),map(x->replace(x,r"[a-z]+",""),s))))
                    # push!(commands, commands[i])
                    # push!(docs, docs[i])
                    # push!(types, types[i])
                else
                    # Only symbols
                    values[i] = "v::Symbol in symbols"
                    s = map(string, a.args[1:end])
                    symbols[i] = Dict(zip(map(upperfirst,s),map(x->replace(x,r"[a-z]+",""),s)))
                end
            end
        else
            # could be nothing
            symbols[i] = Dict()
        end
    end

    insdict = Dict("module"=>"GS200", "type"=>"InsGS200", "make"=>"Yokogawa",
        "model"=>"GS200", "writeterminator"=>"")
    template(zip(commands,values,symbols,types),
        ["cmd","values","symbols","type"], outpath, insdict=insdict, merge=true)
end

"""
`scrape(src::CHMSource, outpath::AbstractString; merge::Bool=true,
    insdict=emptyinsdict(), breakpoints=[])`

Given a `CHMSource` object named `src` and an `outpath`, scrape information out
of the CHM source into organized arrays, and pass that on to [`template`]({ref})
to assemble a JSON file.

Keyword arguments:

- `merge::Bool=true`: Default to attempting to merge with a previous JSON file
rather than overwriting. See the documentation for [`template`]({ref}) for
further details.
- `insdict=emptyinsdict()`: Default to writing an empty instrument dictionary.
- `breakpoints=[]`: If this array has any command strings (as they appear in
the processed JSON file) then we will pause at that command so the source file
can be examined.

"""
function scrape(src::CHMSource, outpath::AbstractString; merge::Bool=true,
    insdict=emptyinsdict(), breakpoints=[])

    cmds = []
    infixes = []
    types = []
    getargs = []
    setargs = []
    syms = []
    rettypes = []
    docs = UTF8String[]
    writeonly = []
    allfiles = (cmds, infixes, types, getargs, setargs, syms, docs, writeonly)

    flag = true
    f=CHM.open(src.path)
    try
        for directory in src.dirs
            files = CHM.readdir(f, directory)
            deleteat!(files, find(x->(match(src.filter, x) == nothing), files))
            for file in files
                text = CHM.retrieve(f, file)

                # Initial pass of text replacement to correct inconsistencies.
                for x in src.parser.init
                    text = replace(text, x...)
                end

                # Parse for get / set capability if desired
                _get, _set = cangetset(src, text)
                if src.parser.getset != r"" && (length(_get) == 0 || length(_set) == 0)
                    flag && (flag = diagnose(text, "Could not determine read/write capability."))
                    continue
                end

                # Get commands and list of arguments
                _cmds, _types, _args, _infixes = commands(src, text)
                len = length(_cmds)
                if len == 0
                    flag && (flag = diagnose(text, "Could not parse command in file: $file"))
                    continue
                elseif len > 2
                    info("$len matches in file: $file")
                end
                flag2 = false
                for cmd in _cmds
                    if contains(cmd, ["{", "(", "["])
                        if contains(cmd, ["}", ")", "]"])
                            flag && (flag = diagnose(text, "Uncaught braces in $cmd"))
                            flag2 = true
                        else
                            flag && (flag = diagnose(text, "Suspect an incomplete command $cmd"))
                            flag2 = true
                        end
                    end
                end
                flag2 && continue

                for cmd in _cmds
                    if findfirst(breakpoints, cmd) != 0
                        flag && (flag = diagnose(text, "Hit a breakpoint."))
                        continue
                    end
                end

                # Parse arguments and find symbols
                argdict = Dict()
                if src.parser.arg != r""
                    for h in eachmatch(src.parser.arg, text)
                        argdict[h[1]] = h[2]
                    end
                end
                _syms = args!(src, _args, argdict)
                println(argdict)

                # Parse for return type of get methods
                _rettypes = rettype(src, text)
                len = length(_rettypes)
                if len == 0 && reduce(|, _get)
                    flag && (flag = diagnose(text, "Could not parse return types in file: $file"))
                    continue
                end

                # Assemble documentation
                _docs = doc(src, file, text)
                len = length(_docs)
                if len == 0
                    flag && (flag = diagnose(text, "Could not parse for documentation in file: $file"))
                    continue
                end

                # Do a final pass for consistency of array lengths, etc.
                _getargs, _setargs, _writeonly = finishfile!(src,
                    _set, _get, _cmds, _infixes, _types, _args, _syms, _rettypes, _docs)

                # all equal lengths?
                valid = reduce(==, extrema(map(length, (_cmds, _types, _getargs, _setargs, _docs))))
                if !valid
                    warn("Could not process file: $file")
                    warn("cmds: $(length(_cmds)) types: $(length(_types)) get: $(length(_getargs)) set: $(length(_setargs)) docs: $(length(_docs))")
                    println(_cmds, _types, _getargs, _setargs, _docs)
                    continue
                end

                # Keep everything we've learned from this file
                thisfile = (_cmds, _infixes, _types, _getargs, _setargs, _syms, _docs, _writeonly)
                for (x,y) in zip(allfiles, thisfile)
                    for z in y
                        push!(x,z)
                    end
                end
                for x in _cmds
                    info("Processed command: $x")
                end

            end
        end
    finally
        CHM.close(f)
    end

    # Save to JSON template file
    template(zip(cmds,infixes,types,setargs,getargs,docs,writeonly),
            ["cmd","infixes","type","setargs","getargs","doc","writeonly"], syms,
            outpath, insdict=insdict, merge=merge)
end

"""
`cangetset(src::CHMSource, text)`

Parse text for get and set capability of commands in the file.
Called by [`scrape`]({ref}). This is not sufficiently general yet.
"""
function cangetset(src::CHMSource, text)
    p = src.parser

    if p.getset == r""
        return [],[]
    end

    entries = eachmatch(p.getset, text)

    getarray = Bool[]
    setarray = Bool[]
    for e in entries
        get = false
        set = false
        for i in 1:3
            ans = (e[i] == nothing ? "" : e[i])
            get |= (lowercase(ans) == "read" ? true : false)
            set |= (lowercase(ans) == "write" ? true : false)
        end
        push!(getarray, get)
        push!(setarray, set)
    end

    getarray, setarray
end

"""
`commands(src::CHMSource, text)`

Parse text for commands. Called by [`scrape`]({ref}).
"""
function commands(src::CHMSource, text)
    p = src.parser
    # entries = eachmatch(p.cmd, text)

    cmds = ASCIIString[]
    types = ASCIIString[]
    args = []
    syms = []
    infixes = []
    argdict = Dict()

    while (e = match(p.cmd, text)) != nothing
        rawcmd = e[1]

        # We obliterate matches in the text as we go.
        text = replace(text, ASCIIString(e.match), "")

        # First pass of command regex replacement
        for x in p.cmd1
            rawcmd = replace(rawcmd, x...)
        end

        # Make types from first pass results
        push!(types, cmdtostr(rawcmd))

        # Second pass of command regex replacement
        for x in p.cmd2
            rawcmd = replace(rawcmd, x...)
        end

        # Save command we discovered
        push!(cmds, rawcmd)

        # Get infixes from that command (look for lowercase strings)
        _infixes = []
        inf = eachmatch(r"[a-z]+", rawcmd)
        for g in inf
            push!(_infixes, g.match*"::Integer=1") # assume Integers, default 1
        end
        push!(infixes, _infixes)

        # Retain arguments
        vals = e.captures[2:end]
        deleteat!(vals, find(x->x==nothing, vals))
        push!(args, vals)
    end

    cmds, types, args, infixes
end

"""
`args!(src::CHMSource, args)`

`args`: array of arrays, one array entry for each command with multiple arguments.

Called by [`scrape`]({ref}).
"""
function args!(src::CHMSource, args, argdict)
    p = src.parser
    syms = []
    sregex = r"\|"

    # Initial pass of argument replacement
    for i in 1:length(args)
        args[i] = convert(Array{ASCIIString,1}, args[i])
        arr = args[i]
        for j in 1:length(arr)
            for (k,v) in argdict
                arr[j] = replace(arr[j], k, v)
            end
        end
    end
    if p.val != r""
        for i in 1:length(args)
            for x in p.val
                args[i] = map(val->replace(val, x...), args[i])
            end
        end
    end

    for i in 1:length(args)
        s = OrderedDict()

        # Separate entries which are still connected by commas
        args[i] = mapreduce(x->split(x,","), vcat, UTF8String[], args[i])

        arr = args[i]
        for j in 1:length(arr)
            # Make argument names unique
            arr[j] = replace(arr[j], "v::", "v$j::")

            # Replace symbols with appropriate arguments
            if ismatch(sregex, arr[j])
                vs = split(arr[j], "|")
                ks = map(upperfirst, vs)
                vs = map(x->replace(x,r"[a-z]",""), vs)
                arr[j] = "v$j::Symbol in s$j"
                s["s$j"] = OrderedDict(zip(ks,vs))
            end
        end

        push!(syms, s)
    end
    syms
end

function doc(src::CHMSource{ENA.SYM}, file, text)
    d = replace(file, "(", "%28")
    d = lowercase(replace(d, ")", "%29"))
    url = "http://ena.support.keysight.com/e5071c/manuals/webhelp/eng"
    ["- [External documentation]($url$d)"]
end

function doc(src::CHMSource{AWG5K.SYM}, file, text)
    [""]
end

function rettype(src::CHMSource, text)
    p = src.parser
    if p.rettype == r""
        return [nothing]
    end

    entries = eachmatch(p.rettype, text)

    rettypes = []
    for e in entries
        t = e.captures[1:end]
        deleteat!(t, find(x->x==nothing, t))
        push!(rettypes, e)
    end
    rettypes
end

function finishfile!(src::CHMSource{ENA.SYM},
        _set, _get, _cmds, _infixes, _types, _args, _syms, _rettypes, _docs)

    getargs = []
    setargs = []
    writeonly = Bool[]
    set = _set[1]
    get = _get[1]
    !set && !get && warn("Neither read nor write in `finishfile!`.")

    # println(_cmds)
    # println(_types)
    # println(_args)
    # println(_syms)
    # println(_docs)
    # chomp(readline())

    if !set && get# read only
        for i in _args
            push!(getargs, i)
            push!(setargs, [])
            push!(writeonly, false)
        end
    elseif !get && set# write only
        for i in _args
            push!(setargs, i)
            push!(getargs, [])
            push!(writeonly, true)
        end
    else
        # only retain the write commands
        getinds = find(x->x[end]=='?', _cmds)
        setinds = find(x->x[end]!='?', _cmds)
        deleteat!(_cmds, getinds)
        deleteat!(_types, getinds)
        deleteat!(_syms, getinds)
        deleteat!(_infixes, getinds)

        for i in setinds
            length(setargs) < length(_cmds)   && push!(setargs, _args[i])
            length(writeonly) < length(_cmds) && push!(writeonly, false)
        end

        for i in getinds
            length(getargs) < length(_cmds) && push!(getargs, _args[i])
        end
    end

    # Documentation goes on all the commands
    d = _docs[1]
    while length(_docs) < length(_cmds)
        push!(_docs, d)
    end

    getargs, setargs, writeonly

end

function finishfile!(src::CHMSource{AWG5K.SYM},
        _set, _get, _cmds, _infixes, _types, _args, _syms, _rettypes, _docs)

    getargs = []
    setargs = []
    writeonly = Bool[]

    println(_cmds)
    println(_types)
    println(_args)
    println(_syms)
    println(_docs)
    chomp(readline())

    if length(_cmds) == 1
        if _cmds[1][end]=='?'
            (get = true; set = false)
        else
            (get = false; set = true)
        end
    else
        get = true; set = true;
    end

    if !set && get# read only
        for i in _args
            push!(getargs, i)
            push!(setargs, [])
            push!(writeonly, false)
        end
    elseif !get && set# write only
        for i in _args
            push!(setargs, i)
            push!(getargs, [])
            push!(writeonly, true)
        end
    else
        getinds = find(x->x[end]=='?', _cmds)
        setinds = find(x->x[end]!='?', _cmds)
        deleteat!(_cmds, getinds)
        deleteat!(_types, getinds)
        deleteat!(_syms, getinds)
        deleteat!(_infixes, getinds)

        for i in setinds
            length(setargs) < length(_cmds)   && push!(setargs, _args[i])
            length(writeonly) < length(_cmds) && push!(writeonly, false)
        end

        for i in getinds
            length(getargs) < length(_cmds) && push!(getargs, _args[i])
        end
    end

    d = _docs[1]
    while length(_docs) < length(_cmds)
        push!(_docs, d)
    end
    while length(_docs) > length(_cmds)
        d = pop!(_docs)
        warn("Dropping documentation: $d")
    end

    for x in _cmds
        info("Processed command: $x")
    end

    getargs, setargs, writeonly
end

"""
`delete!(inpath::AbstractString, itr)`

Delete entries from the property dictionary at key values given by `itr`.
"""
function delete!(path::AbstractString, itr)
    # splitext(inpath)[2] == ".json" || error("Not a JSON file.")
    r = JSON.parsefile(path, dicttype=OrderedDict)
    d = r["properties"]
    for i in itr
        if haskey(d, i)
            delete(d, i)
        end
    end
    r["properties"] = d
    open(path, "w") do f
        JSON.print(f, r, 4)
    end
end

"""
`allproperties(path::AbstractString, key::AbstractString)`

Return all values for a given `key` from instrument property dictionaries in a
JSON file at `path`.
"""
function allproperties(path::AbstractString, key::AbstractString)
    r = JSON.parsefile(path, dicttype=OrderedDict)
    d = r["properties"]
    cmds = ASCIIString[]
    for (k,v) in d
        push!(cmds, v[key])
    end
    cmds
end

"""
`template(data, labels, ids, outpath::AbstractString;
    insdict=emptyinsdict(),
    merge::Bool=true)`

Writes a JSON file given `data` and `labels`. For example:

```
data = [(":OUTPUT", "v::Bool", ""), (":OTHERCMD", "v::Symbol", "")]
labels = ["cmd", "values", "type"]
template(data, labels, 1:2, "/my/path.json")
```

Optionally `insdict` can be provided for the `"instrument"` field of the resulting
JSON file, or if not provided a template is included to be filled in later.

If `merge` is true, then the "type" and symbol fields of a JSON file at `outpath`
are retained for a given `cmd`. The `cmd` field of a property dictionary is used
to distinguish different instrument properties. The user is prompted to advise
how to proceed if a new command is found that wasn't at `outpath`, and also
if multiple entries are found with the same `cmd` field. It is recommended to
work on a duplicate of the original file.
"""
function template(data, labels, syms, outpath::AbstractString;
    insdict=emptyinsdict(),
    merge::Bool=true)

    keepall = false
    if merge == true
        local r
        try
            r = JSON.parsefile(outpath, dicttype=OrderedDict)
        catch e
            warn("No file found at `outpath` to merge with. Writing a new template.")
            return template(data, labels, outpath, insdict=insdict, merge=false)
        end

        # Keep all of the old instrument key/value pairs.
        if haskey(r, "instrument")
            oldid = r["instrument"]
            for k in keys(oldid)
                insdict[k] = oldid[k]
            end
        end

        # As a starting point go with the old properties array.
        if haskey(r, "properties")
            proparr = r["properties"]
        else
            proparr = []
        end

        # Now process the new properties and merge with the old
        for (i,x) in enumerate(data)
            d = OrderedDict(zip(labels, x))
            for (k,v) in syms[i]
                d[k] = v
            end

            ind = find(y->y["cmd"]==d["cmd"], proparr)
            if length(ind) == 0
                if keepall
                    push!(proparr, d)
                else
                    # Ask if we should add this command.
                    warn("""Command not found in old file: $(d["cmd"])""")
                    println("How do you want to proceed? [a=abort, i=ignore, k=keepall, anything else to keep it]")
                    response = chomp(readline())
                    if response == "a" || response == "abort"
                        error("Aborted!")
                    end
                    if response == "i" || response == "ignore"
                    elseif response == "k" || response == "keepall"
                        keepall = true
                        push!(proparr, d)
                    else
                        push!(proparr, d)
                    end
                end
            else
                if length(ind) > 1
                    warn("""Multiple entries with same command: $(d["cmd"])""")
                    println("How do you want to proceed? [a=abort, anything else to expunge all but first entry]")
                    response = chomp(readline())
                    if response == "a" || response == "abort"
                        error("Aborted!")
                    end
                    deleteat!(proparr, ind[2:end])
                end

                olddict = proparr[ind[1]]

                # Retain old types
                if haskey(olddict, "type")
                    d["type"] = olddict["type"]
                end

                # Retain old symbols
                k1 = collect(keys(olddict))
                symdicts = find(x->ismatch(r"s[0-9]+", x), k1)
                for i in symdicts
                    d[k1[i]] = olddict[k1[i]]
                end

                # Everything else, we take from the current parse pass
                proparr[ind[1]] = d
            end
        end
    else
        # Just start from scratch
        proparr = []
        for (i,x) in enumerate(data)
            d = OrderedDict(zip(labels, x))
            for (k,v) in syms[i]
                d[k] = v
            end
            push!(proparr, d)
        end

        # Now look for duplicate cmds
        for (i,d) in enumerate(proparr)
            ind = find(y->y["cmd"]==d["cmd"], proparr)
            if length(ind) > 1
                warn("""Multiple entries with same command: $(d["cmd"]). Expunging all but first.""")
                deleteat!(proparr, ind[2:end])
            end
        end

    end

    temp = OrderedDict("instrument"=>insdict, "properties"=>proparr)

    # Write the dictionary to JSON
    open(outpath, "w") do f
        JSON.print(f, temp, 4)
    end
end

# # # ### Utility functions ### # # #

"""
`cmdtostr(a::AbstractString)`

Returns a string based on a command, e.g. ":SOURCE:FUNCTION" becomes "SourceFunction".
"""
function cmdtostr(a::AbstractString)
    s = mapreduce(upperfirst, *, "", split(a,":"))
    s = replace(s, r"(\[|\{)[A-Za-z0-9\|\-\[\]]+(\}|\])", "")    # only text
    replace(s, "?", "")
end

"""
`diagnose(text)`

Sets up an interactive prompt for diagnosing issues with files that are not
correctly parsed by the given regular expressions. This is called e.g. whenever
a .chm file is scraped, not directly by the user.
"""
function diagnose(text, message)
    println(text)
    warn(message)
    println("Above are the pre-processed file contents.")
    println("\nHow do you want to proceed? [i=ignoreall, a=abort, anything else to continue]")
    response = chomp(readline())
    if response == "a" || response == "abort"
        error("Aborted!")
    elseif response == "i" || replace(response," ","") == "ignoreall"
        return false
    else
        return true
    end
end

"""
`emptyinsdict()`

Returns a `Dict` intended to be used in an instrument JSON file for the "instrument"
dictionary. It has the expected keys and corresponding empty string values.
"""
emptyinsdict() = OrderedDict("module" => "",
                             "type"   => "",
                             "make"   => "",
                             "model"  => "",
                             "writeterminator" => "")

"""
`upperfirst(a::AbstractString)`

Returns a string based on `a` with the first character only in uppercase.
"""
function upperfirst(a::AbstractString)
    length(a) == 0 && return a
    a = lowercase(a)
    string(uppercase(a[1]))*a[2:end]
end

"""
`contains(haystack::AbstractString, needles::AbstractArray)`

Implementation of `contains` for multiple needles...
"""
function contains(string::AbstractString, strings::AbstractArray)
    mapreduce(x->contains(string,x), |, false, strings)
end

end # module


text:
# Arrange
z = MvNormal([0, 0, 0], eye(3))
examplePrior = ExamplePrior(z)

@testset "Factor Creation" begin
    # Assert
    @test examplePrior != nothing
end

@testset "Testing Sampler and Residual" begin
    # Assert
    @test size(getSample(examplePrior, 100)[1]) == (3,100)
end

@testset "Serialization and Deserialization" begin
    # JSON Pack
    jsonConvert = convert(Dict{String, Any}, examplePrior)
    # JSON Unpack
    backToFactor = convert(ExamplePrior, jsonConvert)
    backToJson = convert(Dict{String, Any}, backToFactor)
    # Assert
    @test JSON.json(backToJson) == JSON.json(jsonConvert)

    # Packing and Uunpacking
    packed = convert(PackedExamplePrior, examplePrior)
    unpacked = convert(ExamplePrior, packed)
    repacked = convert(PackedExamplePrior, unpacked)
    # Assert
    @test JSON.json(packed) === JSON.json(repacked)
end