text:
# this function applies the form of GLS stabilization described in:
# A new finite element formulation for computational fluid dynamics:
# X The compressible Euler and Navier-Stokes Equations
# by Shakib, Hughes, and Johan

# this *should* work for both conservative and entropy variables, I think
function applyGLS2(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh, Tsol, Tres, Tdim}

#  println("----- Entered applyGLS2 -----")
  # extract some constants
  numDofPerNode = mesh.numDofPerNode
  numNodesPerElement = mesh.numNodesPerElement
  w = sbp.w
  # reusable storage

  # flux jacobians: a n x n matrix for each coordinate direction for each node
  A_mats = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim,
                      mesh.numNodesPerElement)
#  dxidx_hat = zeros(Tmsh, Tdim, Tdim, mesh.numNodesPerElement)

  # temporarily hold the transposed q variables for the element
  qtranspose = zeros(Tsol, mesh.numNodesPerElement, mesh.numDofPerNode)
  # temporarily hold the result of D*qtranspose
  qxitranspose = zeros(Tsol, mesh.numNodesPerElement, mesh.numDofPerNode, Tdim)
  # hold un-transposed qxi (temporarily)
  qxi = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement, Tdim)
  # hold non-scaled mapping jacobian
  dxidx = zeros(Tmsh, Tdim, Tdim, numNodesPerElement)
  # hold tau, a n x n matrix for each node
  taus = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                    mesh.numNodesPerElement)

  # hold differentiation matrices D
  D = zeros(Float64, numNodesPerElement, numNodesPerElement, Tdim)

  Dtranspose = zeros(D)  # transpose D because we access it in that order
  # temporary vectors 1 through 4, stored in a matrix so they can be
  # access programaticaly
  tmps = zeros(Tres, mesh.numDofPerNode, 2*Tdim)
  qx = zeros(Tres, numDofPerNode, Tdim)  # accumulation vectors

  # DEBUGGING

  gls_res = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement,
                  mesh.numEl)

  gls_full = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement,
                  mesh.numEl)
  tau_eval_sum = 0.0
  tau_eval_cnt = 0
  res_before = copy(eqn.res)

  # calculate D
  for d=1:Tdim
    smallmatmat!(diagm(1./sbp.w), sview(sbp.Q, :, :, d), sview(D, :, :, d))
    Dtranspose[:, :, d] = D[:, :, d].'
  end

#   for el =1:1  #DEBUGGING
  for el = 1:mesh.numEl
    # get all the quantities for this element
    dxidx_hat = sview(mesh.dxidx, :, :, :, el)
    jac = sview(mesh.jac, :, el)
    aux_vars = sview(eqn.aux_vars, :, :, el)
    getGLSVars(eqn.params, eqn.params_conservative, sview(eqn.q, :, :, el), aux_vars, dxidx_hat, sview(mesh.jac, :, el), D, A_mats, qtranspose, qxitranspose, qxi, dxidx, taus)
#=
    if el == 1
      println("q = ", sview(eqn.q, :, :, el))
      println("dxidx_hat = ", dxidx_hat)
      println("D = ", D)
      println("A_mats = ", A_mats)
      println("qxi = ", qxi)
      println("dxidx = ", dxidx)
    end
=#
    for i=1:numNodesPerElement
      res_i = sview(eqn.res, :, i, el)
      #DEBUGGING
      gls_res_i = sview(gls_res, :, i, el)
      gls_full_i = sview(gls_full, :, i, el)
      for j=1:numNodesPerElement

        # zero  out some things
        fill!(qx, 0.0)
        fill!(tmps, 0.0)
        # trial space part
        # rotate qxi, qeta to qx, qy
        for d1=1:Tdim  # the x-y coordinate direction
          for d2=1:Tdim  # the xi-eta coordinate direction
            for n=1:numDofPerNode
              qx[n, d1] += dxidx[d2, d1, j]*qxi[n, j, d2]
            end
          end
        end

#        println("\n  qx = \n", qx)

        # multiply qx, qy  flux jacobians Ax, Ay
        for d1=1:Tdim
          A_d1 = sview(A_mats, :, :, d1, j)
          q_d1 = sview(qx, :, d1)
          tmp_d1= sview(tmps, :, d1)
          smallmatvec!(A_d1, q_d1, tmp_d1)
        end

#        println("\n  trial space terms = \n", tmps[:, 1:2])

        # now add them together
        tmp1 = sview(tmps, :, 1)
        for d1=2:Tdim
          for n=1:numDofPerNode
            tmp1[n] += tmps[n, d1]
          end
        end

        #DEBUGGING

        for n=1:numDofPerNode
          gls_res_i[n] += tmp1[n]
        end

#        println("\n trial space term = \n", tmp1)
        # now multiply by tau
        tau_j = sview(taus, :, :, j)
#        println("  \ntau = \n", tau_j)
        tmp2 = sview(tmps, :, 2)  # store result of multiplication here
        smallmatvec!(tau_j, tmp1, tmp2)  # this overwrites tmp2

#        println("\n  after multiplication by tau, reduction vector = \n", tmp2)
        # copy tmp2 into tmp1 so the next phase of reductions can proceed in
        # parallel
        # multiply by integration weight (from sbp.w) at the same time
        # the indices should be i not j ???
        for n=1:numDofPerNode
          tmp2[n] *= w[i]/jac[i]
          tmp1[n] = tmp2[n]
        end

#        println("\n  after multiplication by w, reduction vector = \n", tmp2)

        # now do weighting space
         # now multiply by flux jacobians transposed
         for d1=1:Tdim
           A_d1 = sview(A_mats, :, :, d1, j)
           x_d1 = sview(tmps, :, d1)
           b_d1 = sview(tmps, :, d1+Tdim)  # now we use the second half of tmps
           smallmatTvec!(A_d1, x_d1, b_d1)
         end

#         println("\n  after multiplication by flux jacobians, reduction vector = \n", tmps[:, 3:4])

         # rotate the differentiation matrices D into x-y
        for d1=1:Tdim  # the x-y coordinate direction
#          println("d1 = ", d1)
          fac = zero(Tmsh)  # variable to accumulate the factor in
          for d2=1:Tdim  # the xi-eta coordinate direction
#              println("d2 = ", d2)
#              println("dxidx = ", dxidx[d2, d1, i])
#              println("d value = ", Dtranspose[i, j, d2])
              fac += dxidx[d2, d1, i]*Dtranspose[i, j, d2]
          end

#          println("fac = ", fac)

          # now update upper half of tmps with the factor
          for n=1:numDofPerNode
            tmps[n, d1 + Tdim] *= fac
          end
        end

#       println("\n  after multiplication by Dx, Dy, reduction vector = \n", tmps[:, 3:4])

       # now update res
       for d1=(Tdim+1):(2*Tdim)  # loop over the second half of tmps
         @simd for n=1:numDofPerNode
           gls_full_i[n] -= tmps[n, d1]  # DEBUGGING
           res_i[n] -= tmps[n, d1]
         end
       end

      end  # end loop over j

      # DEBUGGING
      Dvals, V = eig(sview(taus, :, :, i))
      tau_eval_sum += real(Dvals[1])
      tau_eval_cnt += 1

    end  # end loop over i

  end  # end loop over elements

  return nothing

end  # end function

@doc """
### EulerEquationMod.getGLSVars

  This function calculates all the quantities needed to calculate the GLS
  stabilization term for an element.

  Although I only worked out the math for 2D, I am pretty sure this will
  work in 3D as well

  Inputs:
    params: ParamType, var_type can be anythinga
    params_c: ParamType, var_type must be :conservative
    q: array of solution variables for the element,
       numDofPerNode x numNodesPerElement
    aux_vars: array of auxiliary variables for the element,
              numAuxVars numNodesPerElement
    dxidx_hat: an array of the mapping jacobian scaled by 1/|J|, ie.
               (dxi/dx)/|J| for the entire element,
               (Tdim x Tdim) x numNodesPerElement
    D : an array of the SBP differentiation matrices in the parametric
        coordinate directions, (numNodesPerElement x numNodesPerElement) x Tdim

  Inputs/Outputs:
    A_mats: array the be populated with flux jacobian for the entire element,
            (numDofPerNode x NumDofPerNode) x Tdim x numNodesPerElement
    qtranspose: array to hold transpose(q) (temporarily)
    qxitranspose: array to hold [Dxi*transpose(q), Deta*transpose(q) ...]
                 numNodesPerElement x numDofPerNode x Tdim  (temporarily)
    qxi: transpose of qxtranspose (reversing the first two dimensions)
    dxidx: dxidx_hat * |J| (so the scaling by 1/|J| is undone)
    tau:  array of tau matrices for each node,
          numDofPerNode x numDofPerNode x numNodesPerElement

"""->
function getGLSVars(params::ParamType{Tdim},
params_c::ParamType{Tdim, :conservative},
q::AbstractArray{Tsol, 2}, aux_vars::AbstractArray{Tsol, 2},
dxidx_hat::AbstractArray{Tmsh, 3},
jac::AbstractArray{Tmsh, 1},
D::AbstractArray{Float64, 3},
A_mats::AbstractArray{Tsol, 4}, # begin outputs
qtranspose::AbstractArray{Tsol, 2},
qxitranspose::AbstractArray{Tsol, 3},
qxi::AbstractArray{Tsol, 3},
dxidx::AbstractArray{Tmsh, 3},
tau::AbstractArray{Tres, 3}) where {Tmsh, Tsol, Tres, Tdim}

  # who cares about performance?
  tau_type = params.tau_type

#  println("----- Entered getGLSVars -----")

  numNodesPerElement = size(q, 2)
  numDofPerNode = size(q, 1)
  # move q into q_transpose
#  for i=1:numDofPerNode
#    for j=1:numNodesPerElement
#      qtranspose[j, i] = q[i, j]
#    end
#  end

  # multiply by Dxi, Deta ...
  for d=1:Tdim
    smallmatmatT!(sview(D, :, :, d), q, sview(qxitranspose, :, :, d))
  end

  # un-transpose while putting into qx
  for d=1:Tdim
    for j=1:numNodesPerElement
      for i=1:numDofPerNode
        qxi[i, j, d] = qxitranspose[j, i, d]
      end
    end
  end

   # get dxidx - not scaled by 1/|J|
  for k = 1:numNodesPerElement
    for j=1:Tdim
      for i=1:Tdim
        dxidx[i, j, k] = dxidx_hat[i,j,k]*jac[k]
      end
    end
  end

  # get conservative variable flux jacobian
  if tau_type == 1
#    println("getting conservative variable flux jacobians")
    q_c = params_c.q_vals
    fill!(q_c, 0.0)
    for k=1:numNodesPerElement
      q_k = sview(q, :, k)
      convertToConservative(params, q_k, q_c)
      A1_k = sview(A_mats, :, :, 1, k)
      calcA1(params_c, q_c, A1_k)
      A2_k = sview(A_mats, :, :, 2, k)
      calcA2(params_c, q_c, A2_k)

      if Tdim == 3  # three cheers for static analysis
        A3_k = sview(A_mats, :, :, 3, k)
        calcA3(params_c, q_c, A3_k)
      end
    end
  else
#    println("getting entropy variable flux jacobians")
    # get entropy variable  flux jacobian
    for k=1:numNodesPerElement
      q_k = sview(q, :, k)
#      println("q_k code = ", q_k)
      A1_k = sview(A_mats, :, :, 1, k)
      calcA1(params, q_k, A1_k)
#      println("A1 code = \n", A1_k)
      A2_k = sview(A_mats, :, :, 2, k)
      calcA2(params, q_k, A2_k)
#      println("A2 code = \n", A2_k)

      if Tdim == 3  # three cheers for static analysis
        A3_k = sview(A_mats, :, :, 3, k)
        calcA3(params, q_k, A3_k)
      end
    end
  end

  # get tau for each node - using conservative variable flux jacobian
  for k=1:numNodesPerElement
    q_k = sview(q, :, k)
    tau_k = sview(tau, :, :, k)
    A_mat_k = sview(A_mats, :, :, :, k)
    dxidx_k = sview(dxidx, :, :, k)

    # dance branch prediction, dance
    if tau_type == 1
     getTau(params, q_k, A_mat_k, dxidx_k, tau_k)
    elseif tau_type == 2
      getTau(params, jac[k], tau_k)
    elseif tau_type == 3
      p = 1
#      println("\n arguments passed to getTau from applyGLS2:")
#      println("q_k = ", q_k)
#      println("A_mat_k = \n", A_mat_k)
#      println("dxidx_k = \n", dxidx_k)
#      println("p_val = ", p)
#      println("tau_k = \n", tau_k)
      getTau(params, q_k, A_mat_k, dxidx_k, p, tau_k)
    else
      println(STDERR, "Warning: unsupported Tau requested for GLS")
    end

  end

  if tau_type == 1
    # get entropy variable  flux jacobian
    for k=1:numNodesPerElement
      q_k = sview(q, :, k)
      A1_k = sview(A_mats, :, :, 1, k)
      calcA1(params, q_k, A1_k)
      A2_k = sview(A_mats, :, :, 2, k)
      calcA2(params, q_k, A2_k)

      if Tdim == 3  # three cheers for static analysis
        A3_k = sview(A_mats, :, :, 3, k)
        calcA3(params, q_k, A3_k)
      end
    end
  end

  return nothing
end  # end function

#=
function getGLSVars(params::ParamType{Tdim},
                    params_c::ParamType{Tdim, :conservative},
                    q::AbstractArray{Tsol, 2}, aux_vars::AbstractArray{Tsol, 2},
                    dxidx_hat::AbstractArray{Tmsh, 3},
                    jac::AbstractArray{Tmsh, 1},
                    D::AbstractArray{Float64, 3},
                    A_mats::AbstractArray{Tsol, 4}, # begin outputs
                    qtranspose::AbstractArray{Tsol, 2},
                    qxitranspose::AbstractArray{Tsol, 3},
                    qxi::AbstractArray{Tsol, 3},
                    dxidx::AbstractArray{Tmsh, 3},
                    tau::AbstractArray{Tres, 3}) where {Tmsh, Tsol, Tres, Tdim}

=#

function applyGLS3(mesh::AbstractMesh{Tmsh}, sbp, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tsol, Tres, Tdim, Tmsh}

  numDofPerNode = mesh.numDofPerNode
  numNodesPerElement = mesh.numNodesPerElement
  w = sbp.w
  # reusable storage

  # flux jacobians: a n x n matrix for each coordinate direction for each node
  A_mats = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim,
                      mesh.numNodesPerElement)
#  dxidx_hat = zeros(Tmsh, Tdim, Tdim, mesh.numNodesPerElement)

  # hold non-scaled mapping jacobian
  dxidx = zeros(Tmsh, Tdim, Tdim, numNodesPerElement)
  # hold tau, a n x n matrix for each node
  taus = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                    mesh.numNodesPerElement)

  # hold differentiation matrices Dxi (parametric directions) and Dx
  # (physical directions)
  Dxi = zeros(Float64, numNodesPerElement, numNodesPerElement, Tdim)
  Dx = zeros(Dxi)

  # temporary vectors 1 through 4, stored in a matrix so they can be
  # access programaticaly
  tmps = zeros(Tres, mesh.numDofPerNode, 2*Tdim)
  trial_term = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)
  middle_term = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)
  complete_term = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)

  # calculate Dxi
  # calculate D
  for d=1:Tdim
    smallmatmat!(diagm(1./sbp.w), sview(sbp.Q, :, :, d), sview(Dxi, :, :, d))
  end

  # DEBUGGING
  gls_res = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement,
                  mesh.numEl)
  middle_terms = zeros(gls_res)
  gls_full = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement,
                  mesh.numEl)
  tau_eval_sum = 0.0
  tau_eval_cnt = 0
  res_before = copy(eqn.res)

#   for el = 1:1  # DEBUGGING
  for el = 1:mesh.numEl

    q_el = sview(eqn.q, :, :, el)
    aux_vars_el = sview(eqn.aux_vars, :, :, el)
    dxidx_hat_el = sview(mesh.dxidx, :, :, :, el)
    jac_el = sview(mesh.jac, :, el)

    fill!(Dx, 0.0)
    fill!(taus, 0.0)
    # get all the variables needed for this element
    getGLSVars3(eqn.params, eqn.params_conservative, q_el, aux_vars_el,
                dxidx_hat_el, jac_el, Dxi, Dx, A_mats, dxidx, taus)

    fill!(trial_term, 0.0)
    fill!(middle_term, 0.0)
    fill!(complete_term, 0.0)

    for i=1:mesh.numNodesPerElement  # row index  (free index)
      gls_res_debug_i = sview(gls_res, :, i, el)

      # variable to accumulate the trial term at node i
      trial_term_i = sview(trial_term, :, i)
      for j=1:mesh.numNodesPerElement  # column index  (summed index)
        fill!(tmps, 0.0)
        q_j = sview(eqn.q, :, j, el)

        ### Trial space term ###

        # multiply q_j by entries of the Dx and store in the lower half of
        # tmps
        for d1 = 1:Tdim
          tmp_d1 = sview(tmps, :, d1)
          d_val = Dx[i, j, d1]
          for n=1:numDofPerNode
            tmp_d1[n] = d_val*q_j[n]
          end
        end

        # now multiply by flux jacobian evaluated at i, storing in upper
        # half of tmps
        for d1 = 1:Tdim
          A_i = sview(A_mats, :, :, d1, i)
          tmp_1 = sview(tmps, :, d1)
          tmp_2 = sview(tmps, :, d1 + Tdim)
          smallmatvec!(A_i, tmp_1, tmp_2)
        end

        # accumulate into trial term vector
        # Perf Note: do parallel reads of tmps vectors
        # is there an AVX instruction for scaling by a constant?
        for n=1:numDofPerNode
          for d1 = (Tdim+1):2*Tdim
            trial_term_i[n] += tmps[n, d1]
            #DEBUGGING
            gls_res_debug_i[n] += tmps[n, d1]
          end
        end

      end  # end sum over j in the trial term
    end  # end loop over i for the trial term

#    println("trial_term = \n", trial_term)

    ### now do middle terms ###

    for i=1:mesh.numNodesPerElement
      # multiply by integration weight, jacobian factor
      trial_term_i = sview(trial_term, :, i)
      middle_term_i = sview(middle_term, :, i)
      middle_terms_debug_i = sview(middle_terms, :, i, el)

      fac = w[i]/jac_el[i]
      for n=1:numDofPerNode
        trial_term_i[n] *= fac
      end

      # multiply by tau
      tau_i = sview(taus, :, :, i)
      smallmatvec!(tau_i, trial_term_i, middle_term_i)

      # DEBUGGING
      for n=1:numDofPerNode
        middle_terms_debug_i[n] += middle_term[n]
      end

    end  # end loop over i for middle terms

#    println("middle_term = ", middle_term)

    ### now do weighting space term ###

    for i=1:mesh.numNodesPerElement  # free index
      complete_term_i = sview(complete_term, :, i)
      res_i = sview(eqn.res, :, i, el)
      gls_full_debug_i = sview(gls_full, :, i, el)

      for j=1:mesh.numNodesPerElement  # summed index
        fill!(tmps, 0.0)
        middle_term_j = sview(middle_term, :, j)

        # multiply by transposed flux jacobian, store in lower half of tmps
        for d1 = 1:Tdim
          A_d1 = sview(A_mats, :, :, d1, j)
          tmp_d1 = sview(tmps, :, d1)
          smallmatTvec!(A_d1, middle_term_j, tmp_d1)
        end

        # multiply by entries of the Dx
        for d1 = 1:Tdim
          d_d1 = Dx[j, i, d1]
          tmp_d1 = sview(tmps, :, d1)
          for n=1:numDofPerNode
            tmp_d1[n] *= d_d1
          end
        end

        # sum into complete_term
        for n=1:numDofPerNode
          for d1=1:Tdim
            tmp_d1 = sview(tmps, :, d1)
            complete_term_i[n] += tmp_d1[n]
          end
        end

      end  # end loop over j

      # update res
      for n=1:numDofPerNode
        res_i[n] -= complete_term_i[n]
        #DEBUGGING
        gls_full_debug_i[n] -= complete_term_i[n]
      end

    end  # end loop over i

#    println("complete_term = \n", complete_term)
#    println("res = \n", eqn.res[:, :, el])

  end  # end loop over elements

#=
  println("element 1 gls_res = \n", gls_res[:, :, 1])
  println("element 1 middle_terms = \n", middle_terms[:, :, 1])
  println("element 1 gls_full = \n", gls_full[:, :, 1])
=#
  #=
  gls_res_vec = zeros(Tres, mesh.numDof)
  middle_vec = zeros(gls_res_vec)
  gls_fullvec = zeros(gls_res_vec)
  array3DTo1D(mesh, sbp, eqn, opts, gls_res_vec
  =#

  return nothing

end

function getGLSVars3(params::ParamType{Tdim},
params_c::ParamType{Tdim, :conservative},
q::AbstractArray{Tsol, 2},
aux_vars::AbstractArray{Tsol, 2},
dxidx_hat::AbstractArray{Tmsh, 3},
jac::AbstractArray{Tmsh, 1}, D::AbstractArray{Float64, 3},
# begin output parameters
Dx::AbstractArray{Tmsh, 3}, A_mats::AbstractArray{Tsol, 4},
dxidx::AbstractArray{Tmsh,3}, tau::AbstractArray{Tres, 3}) where {Tmsh, Tsol, Tres, Tdim}

  numDofPerNode = size(q, 1)
  numNodesPerElement = size(q, 2)

   # get dxidx - not scaled by 1/|J|
  for k = 1:numNodesPerElement
    for j=1:Tdim
      for i=1:Tdim
        dxidx[i, j, k] = dxidx_hat[i,j,k]*jac[k]
      end
    end
  end

  # calculate Dx, Dy, Dz
  for d_phys = 1:Tdim  # loop over x, y, z directions
    D_phys = sview(Dx, :, :, d_phys)
    for d_param = 1:Tdim  # loop for xi, eta, ...
      D_param = sview(D, :, :, d_param)
      for j=1:numNodesPerElement
        for i=1:numNodesPerElement
          fac = dxidx[d_param, d_phys, i]  # get jacobian term
          D_phys[i, j] += fac*D_param[i, j]
        end
      end
    end
  end

  tau_type = params.tau_type
  if tau_type == 1
#    println("getting conservative variable flux jacobians")
    q_c = params_c.q_vals
    fill!(q_c, 0.0)
    for k=1:numNodesPerElement
      q_k = sview(q, :, k)
      convertToConservative(params, q_k, q_c)
      A1_k = sview(A_mats, :, :, 1, k)
      calcA1(params_c, q_c, A1_k)
      A2_k = sview(A_mats, :, :, 2, k)
      calcA2(params_c, q_c, A2_k)

      if Tdim == 3  # three cheers for static analysis
        A3_k = sview(A_mats, :, :, 3, k)
        calcA3(params_c, q_c, A3_k)
      end
    end
  else
#    println("getting entropy variable flux jacobians")
    # get entropy variable  flux jacobian
    for k=1:numNodesPerElement
      q_k = sview(q, :, k)
#      println("q_k code = ", q_k)
      A1_k = sview(A_mats, :, :, 1, k)
      calcA1(params, q_k, A1_k)
#      println("A1 code = \n", A1_k)
      A2_k = sview(A_mats, :, :, 2, k)
      calcA2(params, q_k, A2_k)
#      println("A2 code = \n", A2_k)

      if Tdim == 3  # three cheers for static analysis
        A3_k = sview(A_mats, :, :, 3, k)
        calcA3(params, q_k, A3_k)
      end
    end
  end

  # get tau for each node - using conservative variable flux jacobian
  for k=1:numNodesPerElement
    q_k = sview(q, :, k)
    tau_k = sview(tau, :, :, k)
    A_mat_k = sview(A_mats, :, :, :, k)
    dxidx_k = sview(dxidx, :, :, k)

    # dance branch prediction, dance
    if tau_type == 1
     getTau(params, q_k, A_mat_k, dxidx_k, tau_k)
    elseif tau_type == 2
      getTau(params, jac[k], tau_k)
    elseif tau_type == 3
      p = 1
#      println("\n arguments passed to getTau from applyGLS2:")
#      println("q_k = ", q_k)
#      println("A_mat_k = \n", A_mat_k)
#      println("dxidx_k = \n", dxidx_k)
#      println("p_val = ", p)
#      println("tau_k = \n", tau_k)
      getTau(params, q_k, A_mat_k, dxidx_k, p, tau_k)
    else
      println(STDERR, "Warning: unsupported Tau requested for GLS")
    end

  end

  if tau_type == 1
    # get entropy variable  flux jacobian
    for k=1:numNodesPerElement
      q_k = sview(q, :, k)
      A1_k = sview(A_mats, :, :, 1, k)
      calcA1(params, q_k, A1_k)
      A2_k = sview(A_mats, :, :, 2, k)
      calcA2(params, q_k, A2_k)

      if Tdim == 3  # three cheers for static analysis
        A3_k = sview(A_mats, :, :, 3, k)
        calcA3(params, q_k, A3_k)
      end
    end
  end

  return nothing

end  # end function

function test_GLS(mesh::AbstractMesh{Tmsh}, sbp, eqn::AbstractSolutionData{Tsol, Tres}, opts) where {Tsol, Tres, Tmsh}

  eqn.params.tau_type = 2
  Dxi = diagm(1./sbp.w)*sbp.Q[:, :, 1]
  Deta = diagm(1./sbp.w)*sbp.Q[:, :, 2]

  # create indices
  idx_range = Array{UnitRange{Int64}}(mesh.numNodesPerElement)
  for i=1:mesh.numNodesPerElement
    start_idx = (i-1)*mesh.numDofPerNode + 1
    end_idx = i*mesh.numDofPerNode
    idx_range[i] = copy(start_idx:end_idx)
  end

#    println("idx_range = ", idx_range)

    size_block = mesh.numNodesPerElement*mesh.numDofPerNode
#    println("size_block = ", size_block)

  # testing: only do one element
  for el =1:mesh.numEl
#    println("testing element ", el)
    res_el = sview(eqn.res, :, :, el)

    # constant mapping elements only
    dxidx = zeros(2,2)
    dxidx_hat_el = sview(mesh.dxidx, :, :, 1, el)
    jac_el = sview(mesh.jac, :, el)
    q_el = reshape(copy(eqn.q[:, :, el]), size_block)

    for i=1:2
      for j=1:2
        dxidx[i,j] = dxidx_hat_el[i, j, 1]*jac_el[1]
      end
    end

    # calculate Dx, Dy
    Dx = dxidx[1,1]*Dxi + dxidx[2, 1]*Deta
    Dy = dxidx[1,2]*Dxi + dxidx[2, 2]*Deta

#    println("Dx = \n", Dx)
#    println("Dy = \n", Dy)

    # create block Dx, Dy
    Dx_tilde = zeros(Tmsh, size_block, size_block)
    Dy_tilde = zeros(Dx_tilde)

    for i=1:mesh.numNodesPerElement
      idx_i = idx_range[i]
#      println("idx_i = ", idx_i)
      for j=1:mesh.numNodesPerElement
        idx_j = idx_range[j]
#        println("  idx_j = ", idx_j)
        Dx_tilde[idx_i, idx_j] = Dx[i, j]*eye(mesh.numDofPerNode)
        Dy_tilde[idx_i, idx_j] = Dy[i, j]*eye(mesh.numDofPerNode)
      end
    end

#    println("Dx_tilde = \n", Dx_tilde)
#    println("Dy_tilde = \n", Dy_tilde)

    # create A1 tilde and A2 tilde
    A1_tilde = zeros(Tsol, size_block, size_block)
    A2_tilde = zeros(Tsol, size_block, size_block)

    for i=1:mesh.numNodesPerElement
      idx_i = idx_range[i]
      q_i = q_el[idx_i]

      A1 = sview(A1_tilde, idx_i, idx_i)
      EulerEquationMod.calcA1(eqn.params, q_i, A1)

      A2 = sview(A2_tilde, idx_i, idx_i)
      EulerEquationMod.calcA2(eqn.params, q_i, A2)
    end

#    println("A1 = \n", A1_tilde)
#    println("A2 = \n", A2_tilde)

    # create middle terms, including tau
    middle_tilde = zeros(Tres, size_block, size_block)

    for i=1:mesh.numNodesPerElement

      idx_i = idx_range[i]
      tau = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode)
      EulerEquationMod.getTau(eqn.params, jac_el[i], tau)  # tau number 2

      middle_tilde[idx_i, idx_i] = (sbp.w[i]/jac_el[i])*tau
    end

    # create the operator
    fancy_L  = A1_tilde*Dx_tilde + A2_tilde*Dy_tilde
    gls_operator = fancy_L.'*middle_tilde*fancy_L

    @assert isSymmetric(gls_operator, 1e-12)
    #println("max asymmetry = ", maximum(abs(gls_operator - gls_operator.')))

    gls_test = -gls_operator*q_el

    for i=1:size_block
      res_el[i] += gls_test[i]
    end
  end  # end loop over elements

  return nothing

end  # end function

@doc """
### EulerEquationMod.getTau

  This function computes the tau matrix for a single node, for a steady problem
  as described in Hughes part X.

  Inputs
    params:  a ParamType, var_type can be anything
    A_mat: a matrix containing the flux jacobians of the conservative variables
           at the node
           (numDofPerNode x numDofPerNode) x Tdim
    dxidx: a matrix the mapping jacobian for the node (dxi_i/dx_j), Tdim x Tdim
           This should *not* be scaled by 1/|J|

  Inputs/Outputs:
    tau: the numDofPerNode x numDofPerNode matrix to be populated

"""->
function getTau(params::ParamType{Tdim, var_type, Tsol, Tres, Tmsh},
                q::AbstractVector{Tsol}, A_mat::AbstractArray{Tsol, 3},
                dxidx::AbstractArray{Tmsh, 2}, tau::AbstractArray{Tres, 2}) where {Tdim, var_type, Tsol, Tres, Tmsh}

#  println("----- Entered getTau original-----")

  numDofPerNode = size(A_mat, 1)

  data = params.get_tau_data
  @unpack data AjAk flux_term A0inv tmp_mat

  fill!(AjAk, 0.0)  # unneeded?
  fill!(flux_term, 0.0)

  for k=1:Tdim
    for j=1:Tdim
      Aj = sview(A_mat, :, :, j)
      Ak = sview(A_mat, :, :, k)
      smallmatmat!(Aj, Ak, AjAk)  # Aj*Ak

      # calculate factor of dxidx*dxidx + dxidx*dxidx ...
      jacobian_fac = zero(Tmsh)
      for i=1:Tdim
        jacobian_fac += dxidx[i, j]*dxidx[i, k]
      end

      # accumulate dxidx*dxidx*Aj*Ak for all j, k in tmp2
      for p=1:numDofPerNode
        for r=1:numDofPerNode
          flux_term[p, r] += jacobian_fac*AjAk[p, r]
        end
      end

    end  # end loop over j
  end  # end loop over k

  # add a source term here

  # now take negative square root of the matrix
  # there are more efficient ways of doing this
  D, V = eig(flux_term)
  # check that D contains only positive numbers
  # this should be inside a @debug1
  for i=1:numDofPerNode
    if real(D[i]) < 0.0
      println("warning, D[$i] = ", D[i])
    end
  end

  for i=1:numDofPerNode
    D[i] = D[i]^(-0.5)
  end

  # reconstruct M = (V*D)*inv(V)
  Vinv = inv(V)
  for i=1:numDofPerNode
    entry = D[i]
    for j=1:numDofPerNode
      V[j, i] *= entry
    end
  end

  smallmatmat!(V, Vinv, tmp_mat)

  calcA0Inv(params, q, A0inv)
  smallmatmat!(A0inv, tmp_mat, tau)
#  println("----- Finished getTau -----")
  return nothing
end  # end function

function getTau(params::ParamType, jac::Number, tau::AbstractArray{Tres, 2}) where Tres

  fac = 1.0
  for i=1:size(tau, 1)
    tau[i,i] = fac*1/(jac^(1/2))
  end

end

# not sure if this works with conservative variables
# A_mats must be the entropy variable flux jacobians
function getTau(params::ParamType{Tdim, :entropy},
q::AbstractVector{Tsol}, A_mat::AbstractArray{Tsol, 3},
dxidx::AbstractArray{Tmsh, 2}, p::Integer, tau::AbstractArray{Tres, 2}) where {Tsol, Tres, Tmsh, Tdim}

  data = params.get_tau_data
  @unpack data B_d B_p A_mat_hat A0 tmp_mat2
  fill!(B_p, 0.0)
  calcA0(params, q, A0)
  L = chol(A0)' # is the hermitian transpose right here?

  Linv = full(inv(L))  # convert to full matrices because I think it is more
                       # efficient for small matrices

  # calculate the A_i hats = inv(L)*flux_jacobian_i*inv(L).'
  for d1=1:Tdim
    A_hat_d1 = sview(A_mat_hat, :, :, d1)
    A_d1 = sview(A_mat, :, :, d1)
    smallmatmat!(Linv, A_d1, tmp_mat2)
    smallmatmatT!(tmp_mat2, Linv, A_hat_d1)
  end

#  println("A_mat_hat = \n", A_mat_hat)

  for d1=1:Tdim  # loop over B_d
    fill!(B_d, 0.0)

    for d2=1:Tdim  # summed index dxi_d1/dx_d2 A_d2
      A_hat_d2 = sview(A_mat_hat, :, :, d2)
      dxidx_d2 = dxidx[d1, d2]
      # accumulate into B_d
      for i=1:size(B_d, 1)
        for j=1:size(B_d, 2)
          B_d[i, j] += dxidx_d2*A_hat_d2[i, j]
        end
      end

    end  # end loop d2

 #   println("before make_symmetric, symmetry norm = ", vecnorm(B_d - B_d.'))
    make_symmetric!(B_d)  # make sure it is symmetric so we get real eigenvalues
                          # when using real variables

#    println("B$d1 = \n", B_d)
    D, V = eig(B_d)
#    println("D$d1 = \n", D)
#    println("V$d1 = \n", V)
#    println("p = ", p)
 #   println("before update, B_p = \n", B_p)
    # take absolute value, raise to  power p while accumulating into B_p
    for k=1:length(D)  # for each eigenvalue, do  outer product
      v_k = sview(V, :, k)
      val_k = absvalue(D[k])^p
      for i=1:size(B_p, 1)
        for j=1:size(B_p, 2)
          B_p[i, j] += val_k*v_k[i]*v_k[j]
        end
      end
    end

#    println("after adding b$d1, B_p = \n", B_p)

  end  # end loop d1

#  println("B_p = \n", B_p)
  D2, V2 = eig(B_p)

  # now calculate tau: invert and take the pth root of D_p at the same time
  fac = 1.0
  fill!(tau, 0.0)
  for k=1:length(D2)
    v_k = sview(V2, :, k)
    val_k = D2[k]^(-1/p)
    for i=1:size(B_p, 1)
      for j=1:size(B_p, 2)
        tau[i, j] += fac*val_k*v_k[i]*v_k[j]
      end
    end
  end

#  println("tau = \n", tau)

#  println("----- finished getTau -----")
  return nothing
end


text:
for F in (:+, :-, :abs, :inv, :sqrt, :cbrt)
  @eval begin
    $F(x::XFloat16) = XFloat16($F(Float32(x)))
    $F(x::XFloat32) = XFloat32($F(Float64(x)))
  end
end

for F in (:+, :-, :*, :/, :\, :^)
  @eval begin
    $F(x::XFloat16, y::XFloat16) = XFloat16($F(Float32(x), Float32(y)))
    $F(x::XFloat16, y::Float16) = XFloat16($F(Float32(x), Float32(y)))
    $F(x::Float16, y::XFloat16) = XFloat16($F(Float32(x), Float32(y)))
    $F(x::XFloat16, y::Float32) = XFloat16($F(Float32(x), y))
    $F(x::Float32, y::XFloat16) = XFloat16($F(x, Float32(y)))
    $F(x::XFloat16, y::T) where {T<:Real} = $F(promote(x,y)...)
    $F(x::T, y::XFloat16) where {T<:Real} = $F(promote(x,y)...)

    $F(x::XFloat32, y::XFloat32) = XFloat32($F(Float64(x), Float64(y)))
    $F(x::XFloat32, y::Float32) = XFloat32($F(Float64(x), Float64(y)))
    $F(x::Float32, y::XFloat32) = XFloat32($F(Float64(x), Float64(y)))
    $F(x::XFloat32, y::Float64) = XFloat32($F(Float64(x), y))
    $F(x::Float64, y::XFloat32) = XFloat32($F(x, Float64(y)))
    $F(x::XFloat32, y::T) where {T<:Real} = $F(promote(x,y)...)
    $F(x::T, y::XFloat32) where {T<:Real} = $F(promote(x,y)...)

    $F(x::XFloat16, y::XFloat32) = XFloat32($F(Float64(x), Float64(y)))
    $F(x::XFloat32, y::XFloat16) = XFloat32($F(Float64(x), Float64(y)))
  end
end

ext: jl


text:
# broadcast2arg:
# These functions use broadcasting to handle arrays of different sizes.
# Unless otherwise specified they support:
# (N,N) (N,A) (A,N) (A,A) (A,B)
# where N:Number, A,B arrays of broadcast compatible sizes.

broadcast2arg = [
(:.+, :dy, :dy),                    # extra (A,)
(:.*, :(dy.*x2), :(dy.*x1)),                # extra (A,)
(:.-, :dy, :(-dy)),
#:.% => (:dy,:(dy.*(-trunc(x1./x2)))),  # BUG: WARNING: (:check_grads,(:sum,:.%),:args,([-1.6685861285973334,2.349598738753782],[0.5880954718832765,-0.0010728600840855926]),:exact,([1.0,1.0],[2.0,2190.0]),:numeric,([1.0000000000021103,-9.728600840858691],[1.9999999999997797,-4.863172375468294])), WARNING: (:check_grads,(:sum,:.%),:args,([0.20579984208295538,-0.5521335915808314],[0.14504947039368943,-5.795215813098871e-5]),:exact,([1.0,1.0],[-1.0,-9527.0]),:numeric,([0.9999999999998899,-0.15904316261985962],[-0.9999999999998899,0.5895451080050601]))
(:./, :(dy./x2), :(-dy.*x1./abs2(x2))),
(:.\, :(-dy.*x2./abs2(x1)), :(dy./x1)),
(:.^, :(dxndx(x1,x2,dy)), :(dy.*y.*log(x1))), # domain: x1 >= 0 (unless we use complex args)
#:.<< => :todo,                   # domain: Integers, left bit shift; operators,arraymath,broadcast
#:.>> => :todo,                   # domain: Integers, right bit shift
]

for (f,g1,g2) in broadcast2arg
    @eval @primitive $f(x1,x2),dy,y  unbroadcast(x1,$g1)  unbroadcast(x2,$g2)
    if f==(:.^)
        addtest3(f,(0,Inf))
    else
        addtest3(f,(-Inf,Inf))
    end
end

function dxndx(x1,x2,dy)
    if x2==0
        dy.*0
    elseif x2==1
        dy
    elseif x2==2
        2x1.*dy
    else
        dy.*x2.*x1.^(x2-1)
    end
end

broadcast2cmp = [
:.<,
:.<=,
:.==,
:.>,
:.>=,
]

for f in broadcast2cmp
    @eval begin
        # To avoid conflict at broadcast.jl:414
        $f(x1::AbstractArray,x2::Value)=$f(x1,x2.value)
        $f(x1::Value,x2::AbstractArray)=$f(x1,x2.value)
        @zerograd $f(x1,x2)
    end
end

# Other functions in broadcast.jl:

# eval
# droparg1: Not exported
# longer_tuple: Not exported
# longer_size: Not exported
# broadcast_shape: Not exported
# check_broadcast_shape: Not exported
# gen_broadcast_body_cartesian: Not exported
# gen_broadcast_body_iter: Not exported
# bpack: Not exported
# dumpbitcache: Not exported
# gen_broadcast_body_cartesian_tobitarray: Not exported
# gen_broadcast_body_iter_tobitarray: Not exported
# gen_broadcast_function: Not exported
# gen_broadcast_function_tobitarray: Not exported
# broadcast!
# broadcast
# bitbroadcast
# broadcast!_function
# broadcast_function
# broadcast_getindex
# broadcast_getindex!: Not exported
# broadcast_setindex!
# .*
# .%
# .<<
# .>>
# eltype_plus: Not exported
# .+
# type_minus: Not exported
# .-
# type_div: Not exported
# ./
# .\
# type_rdiv: Not exported
# .//
# type_pow: Not exported
# .^
# $(Expr(:$, :f)): Not a symbol
# $(Expr(:$, :bitf)): Not a symbol
# $(Expr(:$, :cachef)): Not a symbol
# bitcache_pow: Not exported


text:
# Node Traverser
# ==============

struct Traverser{T}
    start::T
end

function Base.eltype(::Type{Traverser{T}}) where T
    return T
end

function Base.IteratorSize(::Type{Traverser{T}}) where T
    return Base.SizeUnknown()
end

function traverse(start::Union{NFANode,DFANode,Node})
    return Traverser(start)
end

function Base.iterate(t::Traverser{T}, state=nothing) where T
    if state == nothing
        state = (visited = Set{T}(), unvisited = Set([t.start]))
    end
    if isempty(state.unvisited)
        return nothing
    end
    s = pop!(state.unvisited)
    push!(state.visited, s)
    for (_, t) in s.edges
        if t ∉ state.visited
            push!(state.unvisited, t)
        end
    end
    return s, state
end


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4
coefs = [0.5;0.25]
dense_stencil = @def_stencil_expression(
            @sum(i,-1,1,
                @sum(j,-1,1,
                    @sum(k,-1,1, c[max(abs(i),abs(j),abs(k))]*D[x+i,y+j,z+k]))))
st_def = CreateStencilDefinition(dense_stencil, coefs)
st_inst1 = NewStencilInstance(st_def, m_step=false)
## Input Data size Definition
radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)

t_steps = 16
## Compile Functions

gpu_out = ApplyFFTstencil(st_inst1, data, t_steps)
CUDA.cuProfilerStart()
NVTX.@range "r1" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,1) end
NVTX.@range "r2" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,2) end
NVTX.@range "r3" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,4) end
NVTX.@range "r4" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,8) end
NVTX.@range "r5" begin
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,16) end
exit()


text:
function electrolyzer_TEA(
    flow,
    product_phase,
    product_mw,
    product_density,
    electrons_product,
    electrons_co2,
    system_potential,
    selectivity,
    hours_year,
    current_density,
    single_pass_conversion,
    electrolyzer_cost,
    electricity_price,
    water_price,
    water_rate,
    crf)

    #_____________________________________________________________________________

    F = 96485 # [C/mol]
    g_kg = 1e3
    sec_hour = 3600
    W_MW = 1e6
    CO2_mw = 44 # [g/mol]
    CO2_density = 1.98 #[kg/m^3]
    l_m3 = 1e3 # [liters / m^3]
    sec_min = 60
    cm2_m2 = 100^2
    min_hour = 60
    h2_mw = 2.015 #[g/mol]
    h2_electrons = 2
    h2_density = 0.08375 #[kg/m^3]
    days_year = hours_year / 24
    hours_day = 24
    watergal_kg = 0.2642
    H2O_electrons = 4 # [mol e- / mol product]
    H2O_mw = 18 # [g/mol]

    # [total current]
    amps = flow * F / hours_year / sec_hour / (product_mw / g_kg) * electrons_product / selectivity

    #[m^2 active area]
    electrolyzer_area = amps / current_density / cm2_m2

    #[kg/day CO2]
    co2_flow_needed = amps / electrons_co2 / F * (CO2_mw / g_kg) * sec_hour * hours_day

    # [kg/hour CO2] inlet flow to electrolyzer
    co2_inlet_flow = co2_flow_needed / hours_day / single_pass_conversion

    # [m^3/hour CO2] outlet flow from electrolyzer
    co2_outlet_flow = (co2_inlet_flow - (co2_flow_needed / hours_day))/ CO2_density#co2_inlet_flow * (1 - single_pass_conversion) / CO2_density

    if product_phase == "Liquid"

        #[m^3/hour]
        gas_product_flow = 0.0
        #[l/min]
        liquid_product_flow = flow / product_density * l_m3 / hours_year / min_hour

    elseif product_phase == "Gas"

        #[m^3/hour]
        gas_product_flow = flow / product_density / hours_year
        #[l/min]
        liquid_product_flow = 0.0

    else
        println("Error in product phase classification!")

    end

    #[gal/day]
    water_flow = amps / H2O_electrons / F * (H2O_mw / g_kg) * watergal_kg * sec_hour * hours_day

    #[m^3/hour]
    h2_outlet_flow = amps * (1 - selectivity) / h2_electrons / F * h2_mw / g_kg / h2_density * sec_hour

    #[m^3/hour]
    gas_outlet_flow = co2_outlet_flow + gas_product_flow + h2_outlet_flow

    #[l/min]
    liquid_outlet_flow = liquid_product_flow

    # [$]
    electrolyzer_capex = electrolyzer_area * electrolyzer_cost

    # [$]
    electrolyzer_BOP = electrolyzer_capex * 0.35

    # [$] https://www.eia.gov/analysis/studies/electricity/batterystorage/pdf/battery_storage.pdf
    battery_capex = 1554 * (amps * system_potential / 1e3)

    #[$]
    PSA_capex = 1323.5 * gas_outlet_flow + 632838

    #[$]
    distil_capex = 2773.8 * liquid_outlet_flow + 1e6

    # [$]
    capex =
    electrolyzer_capex +
    electrolyzer_BOP +
    PSA_capex +
    distil_capex +
    battery_capex

    # [$/year]
    maintenance_cost = capex * 0.1

    # [$/year]
    PSA_opex = 0.25 * gas_outlet_flow * hours_year * electricity_price

    # [$/year]
    distil_opex = 9895 * (gas_outlet_flow / 1000) * days_year

    # [$/year]
    water_opex = water_rate * flow * water_price

    # [$/year]
    opex =
    maintenance_cost +
    PSA_opex +
    distil_opex +
    water_opex

    # [$/year]
    annualized_capex = capex * crf

    # [$/year]
    annual_costs = annualized_capex + opex

    return annual_costs

end


text:
"""
    hasdups(guess::NTuple{N,Char}) where {N}

Returns `true` if there are duplicate characters in `guess`.
"""
function hasdups(guess::NTuple{N,Char}) where {N}
    @inbounds for i in 1:(N - 1)
        gi = guess[i]
        for j in (i + 1):N
            gi == guess[j] && return true
        end
    end
    return false
end

"""
    scorecolumn!(col, guess::NTuple{N,Char}, targets::AbstractVector{NTuple{N,Char}})

Return `col` updated with the scores of `guess` on each of the elements of `targets`

If there are no duplicate characters in `guess` a simple algorithm is used, otherwise
the more complex algorithm that accounts for duplicates is used.
"""
function scorecolumn!(
    col::AbstractVector{<:Integer},
    guess::NTuple{N,Char},
    targets::AbstractVector{NTuple{N,Char}},
) where {N}
    if axes(col) ≠ axes(targets)
        throw(
            DimensionMismatch("axes(col) = $(axes(col)) ≠ $(axes(targets)) = axes(targets)")
        )
    end
    if hasdups(guess)
        onetoN = (1:N...,)
        svec = zeros(Int, N)                 # scores for characters in guess
        unused = trues(N)                    # has a character in targets[i] been used
        @inbounds for i in axes(targets, 1)
            targeti = targets[i]
            fill!(unused, true)
            fill!(svec, 0)
            for j in 1:N                     # first pass for target in same position
                if guess[j] == targeti[j]
                    unused[j] = false
                    svec[j] = 2
                end
            end
            for j in 1:N                     # second pass for match in unused position
                if iszero(svec[j])
                    for k in onetoN[unused]
                        if guess[j] == targeti[k]
                            svec[j] = 1
                            unused[k] = false
                            break
                        end
                    end
                end
            end
            sc = 0                           # similar to undup for evaluating score
            for s in svec
                sc *= 3
                sc += s
            end
            col[i] = sc
        end
    else                                     # simplified alg. for guess w/o duplicates
        @inbounds for i in axes(targets, 1)
            sc = 0
            targeti = targets[i]
            for j in 1:N
                sc *= 3
                gj = guess[j]
                sc += (gj == targeti[j]) ? 2 : Int(gj ∈ targeti)
            end
            col[i] = sc
        end
    end
    return col
end

"""
	scoretype(nchar)

Return the smallest type `T<:Unsigned` for storing the scores from a pool of items of length `nchar`
"""
@inline function scoretype(nchar)
    if nchar ≤ 0 || nchar > 80
        throw(ArgumentError("nchar = $nchar is not in `1:80`"))
    end
    return if nchar ≤ 5
        UInt8
    elseif nchar ≤ 10
        UInt16
    elseif nchar ≤ 20
        UInt32
    elseif nchar ≤ 40
        UInt64
    else
        UInt128
    end
end

"""
	tiles(score::Integer, ntiles::Integer)
    tiles(svec::AbstractVector{<:Integer})

Return a length-`ntiles` `String` tile pattern from the numeric score `score`.
"""
function tiles(sc::Integer, ntiles)
    result = sizehint!(Char[], ntiles)    # initialize to an empty array of Char
    for _ in 1:ntiles                     # _ indicates we won't use the value of the iterator
        sc, r = divrem(sc, 3)
        push!(result, iszero(r) ? '🟫' : (isone(r) ? '🟨' : '🟩'))
    end
    return String(reverse(result))
end


text:
"""
Quasi-stationary vector with iterative methods
"""

"""
qstgs(Q::SparseCSC{Tv,Ti}, xi::Vector{Tv};
      x0::Vector{Tv}=stguess(Q,Tv), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6))

Get a quasi-stationary vector of CTMC.

Parameters:
- Q: CTMC Kernal
- xi: Exit vector
- x0: Initial vector for iteration
- maxiter: The maximum number of iteration. The algorithm stops when the number of iteration becomes maxiter.
- steps: The number of steps to check the convergence
- rtol: the tolerance error. When the relative errors of two successive vectors with steps attains rtol, the algorithm stops.
Return value:
A tuple of
- x: quasi-stationary vector
- gam: minimum eigen value
- conv: A boolean whether the algorithm converges or not
- iter: The number of iterations
- rerror: The relative error when the algorithm stops
"""

function qstgs(Q::SparseMatrixCSC{Tv,Ti}, xi::Vector{Tv}; x0::Vector{Tv}=stguess(Q,Tv),
        maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv,Ti}
    qstgs(SparseCSC(Q), xi, x0=x0, maxiter=maxiter, steps=steps, rtol=rtol)
end

function qstgs(Q::SparseCSC{Tv,Ti}, xi::Vector{Tv}; x0::Vector{Tv}=stguess(Q,Tv),
        maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv,Ti}
    m, n = size(Q)
    @assert m == n
    b = zeros(Tv, n)
    x = copy(x0)
    iter = 0
    conv = false
    rerror::Tv = 0
    gam::Tv = 0
    prevx = similar(x)
    while true
        prevx .= x
        for i in 1:steps
            gam = @dot(x, xi)
            gsstep!(x, Q, b, sigma=-gam)
            x ./= sum(x)
        end
        # rerror = maximum(abs.((x - prevx) ./ x))
        rerror = maximum(abs.(x - prevx)) / maximum(x)
        iter += steps
        if rerror < rtol
            conv = true
            break
        end
        if iter >= maxiter
            break
        end
    end
    return x, gam, conv, iter, rerror
end

"""
qstpower(P::AbstractMatrix{Tv}; x0::Vector{Tv}=stguess(Q,Tv), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6))

Get a quasi-stationary vector of DTMC with power method.

Parameters:
- P: The transition probability matrix for DTMC
- xi: Exit vector
- x0: Initial vector for iteration
- maxiter: The maximum number of iteration. The algorithm stops when the number of iteration becomes maxiter.
- steps: The number of steps to check the convergence
- rtol: the tolerance error. When the relative errors of two successive vectors with steps attains rtol, the algorithm stops.
Return value:
A tuple of
- x: quasi-stationary vector
- gam: maximimum eigen value
- conv: A boolean whether the algorithm converges or not
- iter: The number of iterations
- rerror: The relative error when the algorithm stops
"""

function qstpower(P::AbstractMatrix{Tv}, xi::Vector{Tv};
    x0::Vector{Tv}=stguess(P,Tv), maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv}
    m, n = size(P)
    @assert m == n
    Pdash = P'
    x = copy(x0)
    iter = 0
    conv = false
    rerror::Tv = 0
    prevx = similar(x)
    while true
        prevx .= x
        for i in 1:steps
            x = Pdash * x
            x ./= sum(x)
        end
        # rerror = maximum(abs.((x - prevx) ./ x))
        rerror = maximum(abs.(x - prevx)) / maximum(x)
        iter += steps
        if rerror < rtol
            conv = true
            break
        end
        if iter >= maxiter
            break
        end
    end
    nu = @dot(x, xi)
    return x, nu, conv, iter, rerror
end

# function invC(x::Vector{Float64}, Q::AbstractSparseMatrix{T})
#     v = x
#     for j in 1:Q.n
#         tmpx = 0.0
#         for z in Q.colptr[j]:(Q.colptr[j+1]-1)
#             i = Q.rowval[z]
#             if i == j
#                 tmpx += v[j]
#                 v[j] = tmpx / (-Q.nzval[z])
#                 break
#             else
#                 tmpx -= (-Q.nzval[z]) * v[i]
#             end
#         end
#     end
#     v
# end

# function ex(x::Vector{Float64}, Q::SparseMatrixCSC;
#     maxiter=5000, tol=1.0e-16)::Vector{Float64}
#     b = zeros(Q.n)
#     v = x
#     y = invC(v, Q)
#     while maximum(v) > tol
#         gsstep!(v, Q, b)
#         y += v
#     end
#     y
# end

# function ex2(Q::SparseMatrixCSC, b::Vector{Float64};
#     x0::Vector{Float64}=fill(1/Q.n, Q.n),
#     maxiter=5000, steps=50, rtol=1.0e-6)::Vector{Float64}
#     x = x0
#     iter = 0
#     conv = false
#     rerror = 0.0
#     while true
#         prevx = x
#         for i in 1:steps
#             gsstep!(x, Q, b)
#         end
#         rerror = maximum(abs.((x - prevx) ./ x))
#         iter += steps
#         if rerror < rtol
#             conv = true
#             break
#         end
#         if iter >= maxiter
#             break
#         end
#     end
#     @printf "convergence   : %s\n" conv ? "true" : "false"
#     @printf "iteration     : %d / %d\n" iter maxiter
#     @printf "relative error: %e < %e\n" rerror rtol
#     x
# end

# function test2(Q, x)
#     C = triu(Q, 0)
#     D = tril(Q, -1)
#     y = zeros(size(x))
#     x = (-C)' \ x
#     y += x
#     for i in 1:100
#         x = D' * x
#         x = (-C)' \ x
#         y += x
#     end
#     y
# end

# function test(Q::SparseMatrixCSC;
#     x0::Vector{Float64}=fill(1/Q.n, Q.n),
#     maxiter=5000, steps=50, rtol=1.0e-6)::Vector{Float64}
#     x = x0
#     iter = 0
#     conv = false
#     rerror = 0.0
#     while true
#         prevx = x
#         for i in 1:steps
#             b = x
#             gsstep!(x, Q, b)
#             x /= sum(x)
#         end
#         rerror = maximum(abs.((x - prevx) ./ x))
#         iter += steps
#         if rerror < rtol
#             conv = true
#             break
#         end
#         if iter >= maxiter
#             break
#         end
#     end
#     @printf "convergence   : %s\n" conv ? "true" : "false"
#     @printf "iteration     : %d / %d\n" iter maxiter
#     @printf "relative error: %e < %e\n" rerror rtol
#     x
# end


text:
using SchumakerSpline
from = 0.5
to = 10
x1 = collect(range(from, stop=to, length=40))
y1 = (x1).^2
x2 = [0.5,0.75,0.8,0.93,0.9755,1.0,1.1,1.4,2.0]
y2 = sqrt.(x2)
s1 = Schumaker(x1,y1)
s2 = Schumaker(x2,y2)

crossover_point = get_intersection_points(s1,s2)
abs(evaluate(s1, crossover_point[1]) - evaluate(s2, crossover_point[1]))  < 100*eps()

y1 = (x1 .- 5).^2
x2 = [0.5,0.75,0.8,0.93,0.9755,1.0,1.1,1.4,2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
y2 = 1 .+ 0.5 .* x2
s1 = Schumaker(x1,y1)
s2 = Schumaker(x2,y2)

crossover_points = get_intersection_points(s1,s2)
abs(evaluate(s1, crossover_points[1]) - evaluate(s2, crossover_points[1]))  < 100*eps()
abs(evaluate(s1, crossover_points[2]) - evaluate(s2, crossover_points[2]))  < 100*eps()

y1 = (x1 .- 5).^2
x2 = [0.5,0.75,0.8,0.93,0.9755,1.0,1.1,1.4,2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
y2 = -0.5 .* x2
s1 = Schumaker(x1,y1)
s2 = Schumaker(x2,y2)

crossover_points = get_intersection_points(s1,s2)
length(crossover_points) == 0

# Testing Rootfinder and OptimaFinder
from = 0.0
to = 10.0
x = collect(range(from, stop=to, length=400))
# This should have no roots or optima.

y = log.(x) + sqrt.(x)
spline = Schumaker(x,y)
rootfinder = find_roots(spline)
optimafinder = find_optima(spline)
length(rootfinder.roots) == 1
length(optimafinder.optima) == 0
# But it has a point at which it has a value of four:
fourfinder = find_roots(spline; root_value = 4.0)
abs(evaluate(spline, fourfinder.roots[1]) - 4.0) < 1e-10
# and no points where it is negative four::
negfourfinder = find_roots(spline; root_value = -4.0)
length(negfourfinder.roots) == 0
fourfinder2 = find_roots(spline - 2.0; root_value = 2.0)
abs(fourfinder[:roots][1] - fourfinder2[:roots][1]) < eps()
# And if we strict domain to after 2.5 then we will find one root at 103ish
rootfinder22 = find_roots(spline; interval = (2.5,Inf))
spline(rootfinder22.roots[1]) < 1e-10

# This has a root but no optima:
y = y .-2.0
spline2 = Schumaker(x,y)
rootfinder = find_roots(spline2)
optimafinder = find_optima(spline2)
length(rootfinder.roots) == 1
abs(rootfinder.roots[1] - 1.878) < 0.001
length(optimafinder.optima) == 0

y = (x .- 3).^2 .+ 6 # Should be an optima at x = 3. But no roots.
spline3 = Schumaker(x,y)
rootfinder = find_roots(spline3)
optimafinder = find_optima(spline3)
length(rootfinder.roots) == 0
length(optimafinder.optima) == 1
abs(optimafinder.optima[1] - 3.0) < 1e-2
optimafinder.optima_types[1] == :Minimum

# This is a historical bug - It did not find the root because it happened in the interval right before 4.0
spline = Schumaker{Float64}([0.0, 0.166667, 0.25, 0.277597, 0.5, 0.581563, 0.75, 0.974241, 1.0, 95.1966, 100.0, 116.344, 200.0, 233.333], [1.92913 -9.4624 9.91755; 7.71653 -23.3254 8.39407; 40.1234 -101.466 6.50388; 0.617789 -3.15302 3.73426; 1.84611 -5.68404 3.06358; 0.432882 -2.23154 2.61226; 0.3396 -1.88818 2.24866;
                               25.7353 -51.7173 1.84233; 2.00801e-5 -0.183652 0.527206; 0.00772226 0.288784 -16.594; 0.00225409 0.0127389 -15.0287; 8.60438e-5 -0.0760628 -14.2184; 0.000216017 -0.0446771 -19.9793; 5.40043e-5 -0.0550539 -21.2285])
interval = (1e-14, 4.0)
root_value = 0.0
optima = find_roots(spline; root_value = root_value, interval = interval)
length(optima.roots) == 1
abs(optima.roots[1] - 3.875) < 0.005

spline = Schumaker{Float64}([0.0, 0.166667, 0.25, 0.277597, 0.5, 0.581563, 0.75, 0.913543, 1.0, 1.65249, 2.0, 2.33612, 3.0, 3.5334, 4.0, 4.66667], [-0.645627 2.21521 0.0; -2.58251 2.0 0.351267; -13.4282 1.56958 0.5; -0.206757 0.828427 0.533089; -0.617841 0.736461 0.707107; -0.144873 0.635674 0.763065;
                             -0.155837 0.58687 0.866025; -0.557609 0.535898 0.957836; -0.0193613 0.43948 1.0; -0.068257 0.414214 1.27851; -0.072796 0.366774 1.41421; -0.0186602 0.317837 1.52927; -0.0235395 0.293061 1.73205; -0.030761 0.267949 1.88167; -0.0108734 0.239243 2.0; -0.00271836 0.224745 2.15466])
root_value = 2.0
optima = find_roots(spline; root_value = root_value)
length(optima.roots) == 1
abs(optima.roots[1] - 4.0) < 1e-10

# This covers the case where there is an intercept but it comes after the last value of IntStarts but within the interval.
spline = Schumaker{Float64}([0.0, 0.166667, 0.25, 0.277597, 0.5, 0.581563, 0.75, 0.913543, 1.0, 1.75963, 2.0, 2.66667], [1.92913 -10.7562 11.8412; 7.71653 -28.5007 10.1021; 40.1234 -128.376 7.78064;
                          0.617789 -3.56736 4.2684; 1.84611 -6.92217 3.50557; 0.432882 -2.52187 2.95326; 0.465641 -2.46738 2.54076; 1.66614 -5.5315 2.1497; 0.0496917 -1.00572 1.68391; 0.496301 -1.87488 0.948612;
                          0.0929378 -0.847743 0.526628; 0.0232345 -0.618539 0.00277145])
interval = (1e-14, 4.0)
optima = find_roots(spline; interval = interval)
length(optima.roots) == 1
abs(spline(optima.roots[1])) < 1e-10
# In this case the root is in the last interval which is linear rather than quadratic.
spline = Schumaker{Float64}([-1.0e-10, 0.0, 0.02, 0.03, 0.0460546, 0.1, 0.170253, 0.3, 0.361175, 0.5, 0.598784, 0.75, 0.881918, 1.0, 1.27089, 1.5, 1.8286, 2.0, 2.66667, 4.0], [-0.0 0.0 1.74211; 66.2224 -229.519 34.2346; 264.889 -872.979 29.6707; 182.141 -599.117 20.9674; 16.1324 -56.256 11.3957; 6.4478 -23.5354 8.40789; 1.89036 -8.26419 6.78629; 4.79607 -16.6813 5.74585; 0.931297 -4.40443 4.74333; 1.12989 -4.72817 4.14983; 0.482179 -2.64836 3.69379; 0.577267 -2.76551 3.30434; 0.720484 -2.99668 2.94957; 0.162497 -1.37639 2.60576; 0.227166 -1.44474 2.24484; 0.162381 -1.19387 1.92576; 0.59688 -1.97635 1.55098; 0.0393714 -0.694514 1.22978; 0.00984285 -0.598087 0.784266; -0.0 -0.583443 0.0121747])
interval = (1e-14, 5.0)
optima = find_roots(spline; interval = interval)
length(optima.roots) == 1
abs(spline(optima.roots[1])) < 1e-10

# This spline is NOT continuous. It jumps the root. And hence no root is found.
spline = Schumaker{Float64}([-1.0e-10, 0.0, 0.02, 0.03, 0.0607165, 0.1, 0.207868, 0.3, 0.360524, 0.5, 0.584994, 0.75, 0.859153, 1.0, 1.06934, 1.5, 1.75, 2.0, 2.66667, 4.0],
                            [-0.0 0.0 11.2616; 113.819 -400.816 87.8368; 455.276 -1499.77 79.866; 98.1356 -338.854 64.9138; 60.0001 -210.615 54.598; 10.1658 -47.5069 46.4169; 13.9347 -57.0219 41.4107; 30.028 -103.46 36.2755; 5.6543 -26.5882 30.1237; 9.39851 -35.9132 26.5253; 2.49363 -14.6014 23.5408; 4.42393 -19.0774 21.1993; 2.65699 -13.3898                         19.1697;                     3.03898 -13.6262 17.3365; 0.0787707 -5.70898 16.4063; 0.989202 -7.68514 13.9623; -2.41678 -0.111457 12.1028; -0.132839 -5.68619 11.9239; -0.0332096 -6.00949 8.07405; -0.0 -6.05822 -0.023263])
interval = (1e-14, 5.0)
optima = find_roots(spline; interval = interval)
length(optima.roots) == 0

# There was an issue with this spline not getting 3 roots as two of them are within one interval.
spline = Schumaker{Float64}([-1.0e-10, 0.0, 0.0208333, 0.03125, 0.040147, 0.0625, 0.1023, 0.125, 0.203443, 0.25, 0.393202, 0.5, 0.740062, 1.0, 1.26141, 2.0, 3.0, 4.0, 5.33333, 8.0], [-0.0 0.0 4.79187; 63.3179 -152.462 47.9744; 253.272 -515.484 44.8255; 426.822 -843.087 39.4834; 67.6186 -148.651 32.0162; 22.7273 -60.4593 28.7272; 69.8655 -146.831 26.3569; 12.8026 -37.776 23.0599; 36.3453 -78.2209 20.1754; 8.38371 -25.2826 16.6125; 15.0731 -34.0978 13.1639; 5.82125 -16.024 9.69419; 4.96507 -11.9915 6.18291; 5.60305 -10.2219 3.40135; 0.701903 -1.91131 1.1121; 0.979196 -1.0424 0.0833198; -1.7432 0.749616 0.0201181; -0.283294 -1.67429 -0.973462; -0.0708235 -2.08624 -3.70948; -0.0 -2.22356 -9.77643])
interval = (1e-14, 5.0)
root_value = 0.0
optima = find_roots(spline; root_value = root_value, interval = interval)
length(optima.roots) == 3
abs(spline(optima.roots[1])) < 1e-10
abs(spline(optima.roots[2])) < 1e-10
abs(spline(optima.roots[3])) < 1e-10
# There was also a problem that it is getting roots outside of our interval of interest.
interval = (1e-14, 3.0)
root_value = 0.0
optima = find_roots(spline; root_value = root_value, interval = interval)
length(optima.roots) == 2
# There was also a problem that it is getting roots outside of our interval of interest.
interval = (1e-14, 2.9)
root_value = 0.0
optima = find_roots(spline; root_value = root_value, interval = interval)
length(optima.roots) == 1

x = [1,2,3]
y = [3,5,6]
spline = Schumaker(x,y, ; extrapolation = (Constant,Linear))
interval = (1e-14, Inf)
root_value = 6
optima = find_roots(spline; root_value = 5, interval = interval)
length(optima.roots) == 1


text:
module JellyFish
	using Libdl: find_library

	export jaro_winkler, jaro_distance, hamming_distance,
		levenshtein_distance, damerau_levenshtein_distance, soundex, metaphone,
		match_rating_codex, match_rating_comparison, nysiis
	const _jellyfishlib = find_library(["jellyfish"], [joinpath(@__DIR__, "..", "deps")])

	macro _calls1s2type(func, returntype, s1, s2)
		:(ccall(($(esc(func)), _jellyfishlib), $(esc(returntype)), (Ptr{Cchar}, Ptr{Cchar}), $(esc(s1)), $(esc(s2))))
	end

	function jaro_winkler(s1::String, s2::String)
		ccall((:jaro_winkler, _jellyfishlib), Cdouble, (Ptr{Cchar}, Ptr{Cchar}, Bool), s1, s2, false)
	end

	function jaro_distance(s1::String, s2::String)
		@_calls1s2type(:jaro_distance, Cdouble, s1, s2)
	end

	function hamming_distance(s1::String, s2::String)
		convert(Int, @_calls1s2type(:hamming_distance, Csize_t, s1, s2))
	end

	function levenshtein_distance(s1::String, s2::String)
		@_calls1s2type(:levenshtein_distance, Cint, s1, s2)
	end

	function damerau_levenshtein_distance(s1::String, s2::String)
		@_calls1s2type(:damerau_levenshtein_distance, Cint, s1, s2)
	end

	# for some reason doesn't work:
	function _callstr(func::Symbol, s::String)
		r = ccall((func, _jellyfishlib), Ptr{Cchar}, (Ptr{Cchar},), s)
		unsafe_string(r)
	end

	function soundex(s::String)
		# _callstr(:soundex, s)
		r = ccall((:soundex, _jellyfishlib), Ptr{Cchar}, (Ptr{Cchar},), s)
		unsafe_string(r)
	end

	function metaphone(s::String)
		# _callstr(:metaphone, s)
		r = ccall((:metaphone, _jellyfishlib), Ptr{Cchar}, (Ptr{Cchar},), s)
		unsafe_string(r)
	end

	function match_rating_codex(s::String)
		# _callstr(:match_rating_codex, s)
		r = ccall((:match_rating_codex, _jellyfishlib), Ptr{Cchar}, (Ptr{Cchar},), s)
		unsafe_string(r)
	end

	function match_rating_comparison(s1::String, s2::String)
		@_calls1s2type(:match_rating_comparison, Bool, s1, s2)
	end

	function nysiis(s::String)
		# _callstr(:nysiis, s)
		r = ccall((:nysiis, _jellyfishlib), Ptr{Cchar}, (Ptr{Cchar},), s)
		unsafe_string(r)
	end
end


text:
"A fixed-length bit vector whose data is stored in some manually managed region of memory."
struct ManualBitVector <: AbstractArray{Bool, 1}
    ptr::Manual{UInt64}
    length::Int64
end

struct ManualBit
    ptr::Manual{UInt64}
    mask::UInt64
end

@inline function get_address(pv::ManualBitVector, i::Int)
    (i < 1 || i > pv.length) && throw(BoundsError(pv, i))
    i1, i2 = Base.get_chunks_id(i)
    ManualBit(Manual{UInt64}(pv.ptr.ptr + (i1-1)*sizeof(UInt64)), UInt64(1) << i2)
end

@inline function get_address(pv::Manual{ManualBitVector}, i::Int)
    get_address(unsafe_load(pv), i)
end

@inline function Base.unsafe_load(pb::ManualBit)
    (unsafe_load(pb.ptr) & pb.mask) != 0
end

@inline function Base.unsafe_store!(pb::ManualBit, v::Bool)
    c = unsafe_load(pb.ptr)
    c = ifelse(v, c | pb.mask, c & ~pb.mask)
    unsafe_store!(pb.ptr, c)
end

function unsafe_resize!(pb::ManualBitVector, length::Int64)
    @v pb.length = length
end

@inline function Base.findprevnot(pb::ManualBitVector, start::Int)
    start > 0 || return 0
    start > length(pb) && throw(BoundsError(pb, start))

    # TODO placeholder slow implementation; should adapt optimized
    # BitVector code

    @inbounds while start > 0 && pb[start]
        start -= 1
    end

    start

    # Bc = B.chunks
    #
    # chunk_start = Base._div64(start-1)+1
    # mask = ~Base._msk_end(start)
    #
    # @inbounds begin
    #     if Bc[chunk_start] | mask != _msk64
    #         return (chunk_start-1) << 6 + (64 - leading_ones(Bc[chunk_start] | mask))
    #     end
    #
    #     for i = chunk_start-1:-1:1
    #         if Bc[i] != _msk64
    #             return (i-1) << 6 + (64 - leading_ones(Bc[i]))
    #         end
    #     end
    # end
    # return 0
end

@inline function Base.findnextnot(pb::ManualBitVector, start::Int)
    start > 0 || throw(BoundsError(pb, start))
    start > length(pb) && return 0

    # TODO placeholder slow implementation; should adapt optimized
    # BitVector code

    @inbounds while start <= length(pb) && pb[start]
        start += 1
    end

    start
end

@inline function Base.findnext(pb::ManualBitVector, start::Int)
    # TODO placeholder slow implementation; should adapt optimized
    # BitVector code

    len = length(pb)
    @inbounds while start <= len && !pb[start]
        start += 1
    end

    start
end

# array interface

function Base.size(pv::ManualBitVector)
    (pv.length,)
end

@inline function Base.getindex(pv::ManualBitVector, i::Int)
    unsafe_load(get_address(pv, i))
end

@inline function Base.setindex!(pv::ManualBitVector, v::Bool, i::Int)
    unsafe_store!(get_address(pv, i), v)
end

function Base.IndexStyle(_::Type{ManualBitVector})
    Base.IndexLinear()
end


text:
# Julia program to Find the Number of Digits in a Number

# Extracting the number of digits
function number_of_digits(num)
    count = 0
    while(num != 0)
        count = count + 1
        num = num ÷ 10
    end
    return count
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
count = number_of_digits(abs(num))
println("The number of digits in the given number is $count.")

"""
Time Complexity: O(log(num)), where 'num' is the length of the given number
Space Complexity: O(1)
SAMPLE INPUT AND OUTPUT
SAMPLE 1
Enter the number: -12
The number of digits in the given number is 3.
SAMPLE 2
Enter the number: 43258
The number of digits in the given number is 22.
"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
"""
    mutable struct FiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractMPS

Represents a finite matrix product state

When queried for AL/AR/AC/CL it will check if it is missing.
    If not, return
    If it is, calculate it, store it and return
"""
mutable struct FiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractMPS
    ALs::Vector{Union{Missing,A}}
    ARs::Vector{Union{Missing,A}}
    ACs::Vector{Union{Missing,A}}
    CLs::Vector{Union{Missing,B}}

    function FiniteMPS{A,B}(ALs::Vector{Union{Missing,A}},
                            ARs::Vector{Union{Missing,A}},
                            ACs::Vector{Union{Missing,A}},
                            CLs::Vector{Union{Missing,B}}) where {A<:GenericMPSTensor,B<:MPSBondTensor}
        # Relevant checks:
        # for every site at least one AL/AR/AC should not be missing
        # it should end and start with a trivial space
        # consecutive virtual spaces should 'fit' into eachother
        # checking the site tensor spaces
        # at least one CL/AC

        length(ACs) +1 == length(CLs) || throw(ArgumentError("length mismatch of AC/CL"))

        sum(ismissing.(ACs)) + sum(ismissing.(CLs)) < length(ACs)+length(CLs) || throw(ArgumentError("at least one AC/CL should not be missing"))

        prevD2 = nothing;

        for (i,tup) in enumerate(zip(ALs,ARs,ACs))
            non_missing = filter(!ismissing,tup)
            isempty(non_missing) && throw(ArgumentError("missing site tensor"))

            codom = codomain(non_missing[1])
            dom = domain(non_missing[1])
            for j in 2:length(non_missing)
                (dom == domain(non_missing[j]) && codom == codomain(non_missing[j])) || throw(SectorMismatch("AL/AC/AR should be maps over the same space"))
            end

            D1 = _firstspace(non_missing[1])
            D2 = _lastspace(non_missing[1])

            i == 1 || prevD2 == dual(D1) || throw(SectorMismatch("consecutive space tensors don't fit on the virtual level"))

            prevD2 = D2;
            ismissing(CLs[i]) || domain(CLs[i]) == codomain(CLs[i]) || throw(SectorMismatch("CL isn't a map between identical spaces"))
            ismissing(CLs[i+1]) || domain(CLs[i+1]) == codomain(CLs[i+1]) || throw(SectorMismatch("CL isn't a map between identical spaces"))
            ismissing(CLs[i])  || _lastspace(CLs[i]) == dual(D1) || throw(SectorMismatch("CL doesn't fit"))
            ismissing(CLs[i+1]) || _firstspace(CLs[i+1]) == dual(D2) || throw(SectorMismatch("CL doesn't fit"))

            #i != 1 || D1 == oneunit(D1) || throw(ArgumentError("finite mps should start with a trivial leg"))
            #i != length(ACs) || dual(D2) == oneunit(dual(D2)) || throw(ArgumentError("finite mps should end with a trivial leg"))
        end

        return new{A,B}(ALs,ARs,ACs,CLs);
    end
end

# allow construction with one large tensorkit space
function FiniteMPS(f, elt,P::ProductSpace, args...; kwargs...)
    return FiniteMPS(f, elt, collect(P), args...; kwargs...)
end

# allow construction given only a physical space and length
function FiniteMPS(f,elt, N::Int, V::VectorSpace, args...; kwargs...)
    return FiniteMPS(f, elt,fill(V, N), args...; kwargs...)
end

function FiniteMPS(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S;
                    left::S = oneunit(maxvirtspace),
                    right::S = oneunit(maxvirtspace)) where {S<:ElementarySpace}
    N = length(physspaces)
    virtspaces = Vector{S}(undef, N+1)
    virtspaces[1] = left
    for k = 2:N
        virtspaces[k] = infinum(fuse(virtspaces[k-1], fuse(physspaces[k])), maxvirtspace)
    end
    virtspaces[N+1] = right
    for k = N:-1:2
        virtspaces[k] = infinum(virtspaces[k], fuse(virtspaces[k+1], flip(fuse(physspaces[k]))))
    end
    return FiniteMPS(f, elt,physspaces, virtspaces)
end
function FiniteMPS(f,elt,
                    physspaces::Vector{<:Union{S,CompositeSpace{S}}},
                    virtspaces::Vector{S}) where {S<:ElementarySpace}
    N = length(physspaces)
    length(virtspaces) == N+1 || throw(DimensionMismatch())

    tensors = [TensorMap(f, elt,virtspaces[n] ⊗ physspaces[n], virtspaces[n+1]) for n=1:N]

    return FiniteMPS(tensors)
end

# allow construction with a simple array of tensors
function FiniteMPS(site_tensors::Vector{A}) where {A<:GenericMPSTensor}
    for i in 1:length(site_tensors)-1
        (site_tensors[i],C) = leftorth!(site_tensors[i],alg=QRpos());
        site_tensors[i+1] = _permute_front(C*_permute_tail(site_tensors[i+1]))
    end

    (site_tensors[end],C) = leftorth!(site_tensors[end],alg=QRpos());
    B = typeof(C);

    CLs = Vector{Union{Missing,B}}(missing,length(site_tensors)+1)
    ALs = Vector{Union{Missing,A}}(missing,length(site_tensors))
    ARs = Vector{Union{Missing,A}}(missing,length(site_tensors))
    ACs = Vector{Union{Missing,A}}(missing,length(site_tensors))

    ALs.= site_tensors;
    CLs[end] = C;

    FiniteMPS{A,B}(ALs,ARs,ACs,CLs)
end

Base.copy(psi::FiniteMPS{A,B}) where {A,B} = FiniteMPS{A,B}(copy(psi.ALs), copy(psi.ARs),copy(psi.ACs),copy(psi.CLs));

function Base.getproperty(psi::FiniteMPS,prop::Symbol)
    if prop == :AL
        return ALView(psi)
    elseif prop == :AR
        return ARView(psi)
    elseif prop == :AC
        return ACView(psi)
    elseif prop == :CR
        return CRView(psi)
    else
        return getfield(psi,prop)
    end
end

Base.length(psi::FiniteMPS) = length(psi.ALs)
Base.size(psi::FiniteMPS, i...) = size(psi.ALs, i...)

#conflicted if this is actually true
Base.eltype(st::FiniteMPS{A}) where {A<:GenericMPSTensor} = A
Base.eltype(::Type{FiniteMPS{A}}) where {A<:GenericMPSTensor} = A
bond_type(::Type{FiniteMPS{Mtype,Vtype}}) where {Mtype<:GenericMPSTensor,Vtype<:MPSBondTensor} = Vtype

TensorKit.space(psi::FiniteMPS{<:MPSTensor}, n::Integer) = space(psi.AC[n], 2)
function TensorKit.space(psi::FiniteMPS{<:GenericMPSTensor}, n::Integer)
    t = psi.AC[n]
    S = spacetype(t)
    return ProductSpace{S}(space.(Ref(t), Base.front(Base.tail(TensorKit.allind(t)))))
end

virtualspace(psi::FiniteMPS, n::Integer) =
    n < length(psi) ? _firstspace(psi.AC[n+1]) : dual(_lastspace(psi.AC[n]))

r_RR(state::FiniteMPS{T}) where T = isomorphism(Matrix{eltype(T)},domain(state.AC[end]),domain(state.AC[end]))
l_LL(state::FiniteMPS{T}) where T = isomorphism(Matrix{eltype(T)},space(state.AC[1],1),space(state.AC[1],1))

# Linear algebra methods
#------------------------

#=
    At this moment I don't yet make an effort to convert eltype(finitemps) to make it multiply-able by a
    also don't use lmul! from tensorkit, so this thing copies
=#
function Base.:*(psi::FiniteMPS, a::Number)
    return rmul!(copy(psi),a)
end

function Base.:*(a::Number, psi::FiniteMPS)
    return lmul!(a,copy(psi))
end

function TensorKit.lmul!(a::Number, psi::FiniteMPS)
    psi.ACs .*=a;
    psi.CLs .*=a;
    return psi
end

function TensorKit.rmul!(psi::FiniteMPS,a::Number)
    psi.ACs .*=a;
    psi.CLs .*=a;
    return psi
end

function TensorKit.dot(psi1::FiniteMPS, psi2::FiniteMPS)
    #todo : rewrite this without having to gauge
    length(psi1) == length(psi2) || throw(ArgumentError("MPS with different length"))
    ρr = transfer_right(r_RR(psi2),psi2.AR[2:end],psi1.AR[2:end]);
    return tr(_permute_front(psi1.AC[1])' * _permute_front(psi2.AC[1]) * ρr)
end

#todo : rewrite this without having to gauge
TensorKit.norm(psi::FiniteMPS) = norm(psi.AC[1])
TensorKit.normalize!(psi::FiniteMPS) = rmul!(psi, 1/norm(psi))
TensorKit.normalize(psi::FiniteMPS) = normalize!(copy(psi))

"
    A FiniteMPS/MPO starts and ends with a bond dimension 1 leg
    Its bond dimension also can't grow faster then pspaces (resp. pspaces^2)
    This function calculates the maximal achievable bond dimension
"
function max_Ds(f::FiniteMPS{G}) where G<:GenericMPSTensor{S,N} where {S,N}
    Ds = [1 for v in 1:length(f)+1];
    for i in 1:length(f)
        Ds[i+1] = Ds[i]*prod(map(x->dim(space(f.AC[i],x)),ntuple(x->x+1,Val{N-1}())))
    end

    Ds[end] = 1;
    for i in length(f):-1:1
        Ds[i] = min(Ds[i],Ds[i+1]*prod(map(x->dim(space(f.AC[i],x)),ntuple(x->x+1,Val{N-1}()))))
    end
    Ds
end


text:
export SpecialKey, AdjointKey
export Event, KeyPressedEvent, PasteEvent
export match

@enum SpecialKey begin
    # Function Keys match range 0:11
    F1
    F2
    F3
    F4
    F5
    F6
    F7
    F8
    F9
    F10
    F11
    F12

    # Key Code match range 12:17 (real code: 1:6(8))
    HOME # 1 or 7
    INSERT
    DELETE
    END # 4 or 8
    PAGEUP
    PAGEDOWN

    # Direction Keys
    UP = Int('A')
    DOWN = Int('B')
    RIGHT = Int('C')
    LEFT = Int('D')
    BACKTAB = Int('Z')

    ENTER
    ESC
    TAB
    BACKSPACE
    NULL
end

@enum AdjointKey begin
    SHIFT = 2
    CTRL = 3
    SHIFT_ALT = 4
    ALT = 5
    SHIFT_CTRL = 6
    CTRL_ALT = 7
    SHIFT_CTRL_ALT = 8
    NO_ADJOINT = -1
end

abstract type Event end

KeyboardKey = Union{Char, SpecialKey}

struct KeyPressedEvent <: Event
    key::KeyboardKey
    adjoint_key::AdjointKey
end

KeyPressedEvent(key::KeyboardKey) = KeyPressedEvent(key, NO_ADJOINT)

function Base.show(io::IO, e::KeyPressedEvent)
    adjoint_key = (e.adjoint_key === NO_ADJOINT) ? "" : "$(string.(e.adjoint_key))+"
    Base.print(io, "KeyPressedEvent($(adjoint_key)'$(string(e.key))')")
end

struct PasteEvent <: Event
    content::String
end

Base.show(io::IO, e::PasteEvent) = Base.print(io, "PasteEvent(\"$(e.content)\")")


text:
module KeboolaConnectionComponent

import JSON
include("Config.jl")

"""
`KeboolaConnectionComponent.datadir()` => `String`

Returns the path to the [data directory](https://developers.keboola.com/extend/common-interface/folders/).
"""
function datadir()
    path = get(ENV, "KBC_DATADIR", "/data/")
    if !isdirpath(path)
        path = path * "/"
    end
    path
end

"""
`KeboolaConnectionComponent.config()` => `Config`

Returns the processed configuration.
"""
function config()
    open(joinpath(datadir(), "config.json"), "r") do file
        data = JSON.parse(file)
        Config(
            get(data, "action", "run"),
            get(data, "parameters", Dict()),
            get(data, "storage", Dict()),
            get(data, "image_parameters", Dict()),
            get(data, "authoziation", Dict()),
        )
    end
end

end # module

max_stars_repo_name: keboola/julia-component

max_issues_repo_name: keboola/julia-component

max_forks_repo_name: keboola/julia-component


text:
module ExponentialSmoothing

using Optim
using Statistics
using StatsFuns
using TimeSeries
using RecipesBase
using Colors
using ColorSchemes

include("main.jl")

end # module


text:
#=
# MKL vs. OpenBLAS

With an Intel CPU, [MKL](https://github.com/JuliaLinearAlgebra/MKL.jl) is generally faster
than OpenBLAS. Let's do a quick comparison.
=#

import MLDatasets
import Makie
import CairoMakie
using ValueHistories: MVHistory

# Load MNIST

Float = Float32
train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float}(train_x[:, :, train_y .== 0] .≥ 0.5)
nothing #hide

# Make  sure we are using OpenBLAS first:

using LinearAlgebra
LinearAlgebra.__init__() # use OpenBLAS
BLAS.get_config()

# Number of BLAS threads:

BLAS.get_num_threads()

# Train RBM using OpenBLAS.

import RestrictedBoltzmannMachines as RBMs
rbm = RBMs.BinaryRBM(Float, (28,28), 128)
RBMs.initialize!(rbm, train_x)
history_openblas = MVHistory()
time_0 = time()
RBMs.cd!(
    rbm, train_x; epochs=50, batchsize=128, steps=1,
    callback = function(@nospecialize(args...); @nospecialize(kw...))
        push!(history_openblas, :t, time() - time_0)
    end
)
nothing #hide

# Now load MKL.

using MKL
MKL.__init__() # don't need this on a fresh Julia session
BLAS.get_config()

# Number of BLAS threads:

BLAS.get_num_threads()

# Now let's rerun the RBM training.

RBMs.initialize!(rbm, train_x)
history_mkl = MVHistory()
time_0 = time()
RBMs.cd!(
    rbm, train_x; epochs=50, batchsize=128, steps=1,
    callback = function(@nospecialize(args...); @nospecialize(kw...))
        push!(history_mkl, :t, time() - time_0)
    end
)
nothing #hide

# The epochs should be somewhat faster with MKL.

fig = Makie.Figure(resolution=(600, 400))
ax = Makie.Axis(fig[1,1], xlabel="epoch", ylabel="seconds")
Makie.lines!(ax, get(history_openblas, :t)..., label="OpenBLAS")
Makie.lines!(ax, get(history_mkl, :t)..., label="MKL")
Makie.ylims!(ax, low=0)
Makie.axislegend(ax, position=:rt)
fig


text:
"""
Returns a Dict corresponding to the results of querying `service`, where

## Parameters
* `service`: one of the Google Maps API (e.g. `geocode` or `directions`)

## Keyword Arguments
Often corresponds to parameters that are passed to the REST-ful API.
"""
function requestjson(
        service::String,
        api::String = "https://maps.googleapis.com/maps/api/",
        ext::String = "/json";
        kwargs...
    )
    #result = Requests.get(string(api,service,ext); query = Dict(kwargs))
    result = HTTP.request("GET",string(api,service,ext); query = Dict(kwargs))
    #resultjson = JSON.parse(Requests.readstring(result))
    resultjson = JSON.parse(String(result.body))
    #info("REQUEST:{\"uri\":\"", result.request.value.uri, "\", ",
    #              "\"status\": \"", Requests.statuscode(result), "\"}")
    #Requests.statuscode(result) == 400 && warn("ERROR:",resultjson["error"])
    resultjson
end

function requestimage(service::String; kwargs...)
    result = Requests.get(
        "https://maps.googleapis.com/maps/api/$service"; query = Dict(kwargs)
    )
    info("REQUEST:{\"uri\":\"", result.request.value.uri, "\", ",
                  "\"status\": \"", Requests.statuscode(result), "\"}")
    ImageMagick.readblob(result.data)
end

"""
Decodes a sequence of `lonlats` encoded in the polyline algorithm format.

See https://developers.google.com/maps/documentation/utilities/polylinealgorithm
for details.
"""
function decodepolyline(polyline::Vector{UInt8})
    function f(polyline, index)
        x = shift = 0
        while true
            @inbounds block = polyline[index] - 63 # 63 == Int('?')
            x |= (block & 0x1f) << shift
            index += 1; shift += 5
            block >= 0x20 || break
        end
        index, (x & 1) == 1 ? ~(x >> 1) : (x >> 1)
    end

    index = 1; currlat = currlon = 0
    lonlat = Tuple{Float64,Float64}[]
    polylength = length(polyline)
    while index < polylength
        index, dlat = f(polyline, index)
        index, dlon = f(polyline, index)
        currlat += dlat; currlon += dlon
        push!(lonlat, (currlon / 1e5, currlat / 1e5))
    end
    lonlat
end

decodepolyline(s::String) = decodepolyline(Vector{UInt8}(s))